Identification of Continuous-Discrete Hidden Markov Models with Multiplicative Observation Noise

Abstract

The paper aims to identify hidden Markov model parameters. The unobservable state represents a finite-state Markov jump process. The observations contain Wiener noise with state-dependent intensity. The identified parameters include the transition intensity matrix of the system state, conditional drift and diffusion coefficients in the observations. We propose an iterative identification algorithm based on the fixed-interval smoothing of the Markov state. Using the calculated state estimates, we restore all required system parameters. The paper contains a detailed description of the numerical schemes of state estimation and parameter identification. The comprehensive numerical study confirms the high precision of the proposed identification estimates.

1. Introduction

The parameter identification of the stochastic dynamic systems given the indirect heterogeneous noisy observations is one of the most complicated and demanding problems in system analysis. Usually, this type of estimation appears for either the development and validation of the mathematical models with their subsequent analysis and forecasting or various versions of the stochastic control under incomplete information. Engineers often use the results of identification in mechanical and electronic systems, in navigation, biology and medicine, economics and finance as well as speech and image processing.

There exists both the series of classical monographs and the number of comprehensive surveys, reflecting the state-of-the-art in the stochastic identification theory. The hidden Markov models (HMMs) are the prospective ones [1,2,3,4]. First, they form the mathematical basis for the description of the processes in telecommunications [5,6], finance [7,8,9], biology and medicine [10,11,12,13], tracking and navigation [14,15] and speech and image recognition [16,17,18,19]. Second, the HMM advanced theoretical framework provides a solution to a series of optimal estimation and control problems [20].

The investigation object of this paper represents a subclass of continuous-time HMMs with the states described by the finite-state Markov jump processes (MJPs) given the diffusion observations. Researchers have used various methods and algorithms to identify parameters in these HMMs [21,22,23,24,25,26]; however, the most effective is the maximum likelihood method as realized by the Expectation-Maximization algorithm (EM algorithm) [2,27]. It was successfully used for parameter identification in various stochastic models and particularly for the separation of the Gaussian mixtures [28]. To identify the HMM parameters, the authors of [29,30] introduced the corresponding version of the algorithm and proved the conditions of its convergence to a local maximum of the likelihood function. The essential assumption is the constant intensity of the Wiener observation noises.

The identified parameters include the transition intensity matrix (TIM) of the MJP and the conditional drifts in the observations. Note that the EM-procedure leans on the optimal filtering of the system state and some functionals of their trajectories. The paper aims to develop a new modification of the identification algorithm for the continuous-time HMM parameters with the state-dependent noise intensity in the observations. Such observation noises is also called multiplicative. Traditionally, the estimation problems in such systems are rarely investigated because of the known theoretical difficulties (see [31,32] for the details). The formal theoretical solution to the optimal filtering problem of the MJP state given the diffusion observations with multiplicative noises has been published recently [33,34] along with a complex of the stable numerical algorithms for their realization. The presented work is the natural evolution of that research.

The paper has the following organization. Section 2 contains a formal description of the investigated observation system and the statement of its parameter identification problem. The section also includes a comparative analysis of the various problem formulations and the properties of the EM-estimates, calculated for the adjacent HMMs: the discrete-time systems and the continuous-time ones with the homogeneous observation Wiener noises. Additionally, the section presents some reasons in favor of the proposed version of the identification algorithm. It is based on the high precision fixed interval smoothing estimates of the hidden state.

Section 3 contains a numerically effective “two-filter smoothing” algorithm, which expresses the smoothed estimates as a symmetric function of the forward- and backward-time filtering ones. The calculation of the smoothed estimates represents the E-step of the EM algorithm. M-step itself and the whole recursive EM-procedure is introduced in Section 4. The section contains the conversion of the ${\mathcal{L}}_2$-optimal state estimates to the ${\mathcal{L}}_1$-optimal ones and the formulae for the estimation of the required functionals of the MJP state trajectory. The section also presents pseudocode sketching the structure of the EM algorithm.

The solution to the optimal MJP state filtering given the diffusion observations with the multiplicative noises is a backbone of the identification algorithm. The theoretical estimator, being a variant of the generalized Bayes rule, contains the integrals by an abstract probability measure. They cannot be calculated analytically; hence, Section 5 proposes their numerical approximations and the accuracy characteristics.

Section 6 contains a complex of illustrative numerical examples. Section 6.1 presents the comparison of the filtering and fixed-interval smoothing estimates. The new version of the estimates is also compared with one from [2], initially proposed for the processing of the discrete-time HMMs. Section 6.2 contains diversified numerical analysis of the proposed identification estimates. Section 7 contains our concluding remarks.

2. Problem Formulation

On the probability triplet with filtration $(\Omega ,\mathcal{F},\mathcal{P},{\left\{\mathcal{F}\right\}}_{t\in [0,T]})$, we consider the partially observable HMM

$$\begin{array}{c}{\displaystyle X_t=X_0+{\int }_0^t{\Lambda }^{\top }{X}_{s}ds+M_t^X,}\end{array}$$

(1)

$$\begin{array}{c}{\displaystyle {\mathcal{Y}}_r={\int }_{t_{r-1}}^{t_r}f{X}_{s}ds+{\int }_{t_{r-1}}^{t_r}\sum _{n=1}^N{X}_s^n{g}_n^{\frac{1}{2}}d{W}_s,r\in \{1,\dots ,R\},t_r=r\delta ,T=R\delta ,}\end{array}$$

(2)

where

• $X_t=col(X_t^1,\dots ,X_t^N)\in {\mathbb{S}}^N$ is an unobservable state, which is a finite-state MJP with the state space ${\mathbb{S}}^N\triangleq \{e_1,\dots ,e_N\}$ (${\mathbb{S}}^N$ stands for the set of all unit coordinate vectors of the Euclidean space ${\mathbb{R}}^N$) with the transition intensity matrix (TIM) $\Lambda$ and initial distribution $\pi =col({\pi }^1,\dots ,{\pi }^N)$; the process $M_t^X$ is an  ${\mathcal{F}}_t$-adapted martingale; and
• ${\mathcal{Y}}_r=col({\mathcal{Y}}_r^1,\dots ,{\mathcal{Y}}_r^M)\in {\mathbb{R}}^M$ are the diffusion observations

  $$Y_t={\int }_0^{t}f{X}_{s}ds+{\int }_0^t\sum _{n=1}^N{X}_s^n{g}_n^{\frac{1}{2}}d{W}_s,$$

  (3)

  discretized with a (small) time increment $\delta$; $W_t=col(W_t^1,\dots ,W_t^M)\in {\mathbb{R}}^M$ is an ${\mathcal{F}}_t$-adapted standard Wiener process characterizing the observation noise; and
  f is an $M\times N$-dimensional observation matrix, whose columns define the conditional observation drift, whereas the collection of $M\times M$-dimensional matrices ${\left\{g_n\right\}}_{n=\overline{1,N}}$ represents the conditional observation diffusion coefficient given $X_t=e_n$.

The non-decreasing family of $\sigma$-algebras generated by the discretized observations $\mathcal{Y}$ up to the moment $t_r$ is denoted by ${\mathfrak{Y}}_r\triangleq \{{\mathcal{Y}}_q:q=\overline{1,r}\}$, ${\mathfrak{Y}}_0\triangleq \{\varnothing ,\Omega \}$.

We suppose that the assumptions below hold for the system (1) and (2).

• The probability triplet with filtration $(\Omega ,\mathcal{F},\mathcal{P},{\left\{\mathcal{F}\right\}}_{t\in [0,T]})$ is the Wiener–Poisson one [35].
• All the HMM parameters $\Lambda$, f and ${\left\{g_n\right\}}_{n=\overline{1,N}}$ are unknown nonrandom matrices of appropriate dimensionality.
• The transition intensity matrix $\Lambda$ satisfies the condition

  $$\underset{i,j:i\ne j}{min}{\lambda }_{ij}>0$$

  (4)

  along with the usual identity of the TIMs

  $$\sum _j{\lambda }_{ij}\equiv 0.$$
• The noises in $\mathcal{Y}$ are uniformly nondegenerate [36], i.e., $g_n>0$ for all $n=\overline{1,N}$.
• The collection ${\left\{g_n\right\}}_{n=\overline{1,N}}$ of the observation noise intensity meets the indentifiability condition [33], i.e., all the matrices $\left\{g_n\right\}$ are different.

The HMM identification problem is to find the estimates of $\Lambda$, f and ${\left\{g_n\right\}}_{n=\overline{1,N}}$ given the observations $\{{\mathcal{Y}}_1,\dots ,{\mathcal{Y}}_R\}$.

Remark 1.

Note that condition (4) can be weakened, leaving the MJP X ergodicity only.

The identification problem of the system (1), (2) parameters has its inner history. It represents an approximation of the parameter estimation problem by the continuous-time diffusion observations (3). The authors of [29,30] investigated the continuous-time system (1), (3) with scalar observations and additive noise (i.e., the case $g_n\equiv g$). They developed a procedure of simultaneous MJP state filtering and system parameter identification based on the EM algorithm. The parameter identification requires the estimation of the number of $e_i\to e_j$ transitions

$$N_t^{i,j}={\int }_0^t{X}_{s-}^{i}d{X}_s^j,$$

(5)

and the cumulative occupation time of the value $e_n$

$$O_t^n={\int }_0^t{X}_s^{n}ds.$$

(6)

The definitions given in (5) and (6) refer, respectively, to the transitions and occupation time related to interval $[0,t]$.

The authors of the monograph [2] investigated a discrete-time analog of (1), (3). The available observation is a vector-valued process corrupted by some multiplicative noise. The authors modified the EM algorithm and incorporated the fixed-interval smoothed estimates of $\left\{N_t^{i,j}\right\}$, $\left\{O_t^n\right\}$ and the “level sums”

$$Q_t^{n,q}={\int }_0^t{X}_s^{n}q\left(Y_s\right)d{Y}_s$$

(7)

in it. Note, the mentioned estimates of $\left\{N_t^{i,j}\right\}$, $\left\{O_t^n\right\}$ and $\left\{Q_t^{n,q}\right\}$ are not calculated by the two-filter smoothing algorithm, similar to [37], due to the nature of the underlying processes.

The paper [33] contains a solution to the optimal filtering of the MJP states given the vector observations with multiplicative noises. The authors suggested an observation transform that splits the initial ones into a diffusion with unit intensity, collection of counting processes and a set of random vectors observed at some non-random moments. The theoretical estimate represents the solution to the discrete-continuous stochastic differential system (SDS) with diffusion and counting processes on the right-hand side.

The article declares that identifiability condition 4 is sufficient for the precise restoration of the MJP state given the indirect noisy observations. Unfortunately, the solution to the SDS cannot be realized numerically due to the necessity to perform a double limit passage for the required observation transform. To overcome this obstacle, the authors suggested a collection of numerical filtering algorithms by the descretized observations. They investigated the performance of the proposed approximation and show that the algorithms deliver high precision to the MJP state estimates.

It is difficult to find the optimal estimates of the processes $\left\{N_t^{i,j}\right\}$, $\left\{O_t^n\right\}$ and $\left\{Q_t^{n,q}\right\}$ directly by the discretized observations. Nevertheless, it is possible to construct suboptimal ones based on the high precision two-filter smoothing estimates of the MJP state $X_t$.

Below in the presentation, we use the following notations:

• $\mathbf{1}$ is a row vector of the appropriate dimensionality formed by units.
• I is a unit matrix of the appropriate dimensionality.
• ${\mathbf{I}}_{\mathcal{A}}\left(x\right)$ is a indicator function of the set $\mathcal{A}$.
• ${\parallel \alpha \parallel }_K^2\triangleq {\alpha }^{\top }K\alpha$.
• $N_r^X$ is a random number of the state $X_t$ transitions, occurred on the interval $[t_{r-1},t_r]$.
• ${\mathcal{P}}_{\xi }$ is a contraction of the probability measure $\mathcal{P}$ to the $\sigma$-algebra $\sigma \{X_t,Y_t:0⩽t⩽T\}$.
• ${\tau }_r=col({\tau }_r^1,\dots ,{\tau }_r^N)\triangleq {\int }_{t_{r-1}}^{t_r}{X}_{s}ds$ is a random vector composed of the occupation times of X in each state $e_n$ during the interval $[t_{r-1},t_r]$.
• $\mathcal{D}\triangleq \{u=col(u^1,\dots ,u^N):u^n⩾0,{\sum }_{n=1}^N{u}^n=h\}$ is an $(N-1)$-dimensional simplex in the space ${\mathbb{R}}^M$, $\mathcal{D}$ is a distribution support of the vector ${\tau }_r$.
• $\Pi \triangleq \{\pi =col({\pi }^1,\dots ,{\pi }^N):{\pi }^n⩾0,{\sum }_{n=1}^N{\pi }^n=1\}$ is a “probabilistic simplex” formed by the possible values of $\pi$.
• ${\rho }_r^{k,\ell }\left(du\right)$ is a conditional distribution of the vector $X_{t_r}^{\ell }{\tau }_r$ given $X_{t_{r-1}}=e_k$, i.e.,  for any $\mathcal{G}\in \mathcal{B}\left({\mathbb{R}}^M\right)$ the following equality holds:

  $${\mathsf{E}}_{}\left\{{\mathbf{I}}_{\mathcal{G}}\left({\tau }_r\right)X_{t_r}^{\ell }|X_{t_{r-1}}=e_k\right\}={\int }_{\mathcal{G}}{\rho }_r^{k,\ell }\left(du\right).$$
• ${\rho }_r^{k,\ell ,q}\left(du\right)$ is a conditional distribution of $X_{t_r}^{\ell }{\mathbf{I}}_{\left\{q\right\}}\left(N_r^X\right){\tau }_r$ given $X_{t_{r-1}}=e_k$, i.e.,  for any $\mathcal{G}\in \mathcal{B}\left({\mathbb{R}}^M\right)$ the following equality holds:

  $${\mathsf{E}}_{}\left\{{\mathbf{I}}_{\mathcal{G}}\left({\tau }_r\right){\mathbf{I}}_{\left\{q\right\}}\left(N_r^X\right)X_{t_r}^{\ell }|X_{t_{r-1}}=e_k\right\}={\int }_{\mathcal{G}}{\rho }_r^{k,\ell ,q}\left(du\right).$$
• ${\varrho }_r^{k,\ell }\left(du\right)$ is a conditional distribution of $X_{t_r}^{\ell }{\tau }_{r+1}$ given $X_{t_{r+1}}=e_k$, i.e.,  for any $\mathcal{G}\in \mathcal{B}\left({\mathbb{R}}^M\right)$ the following equality holds:

  $${\mathsf{E}}_{}\left\{{\mathbf{I}}_{\mathcal{G}}\left({\tau }_{r+1}\right)X_{t_r}^{\ell }|X_{t_{r+1}}=e_k\right\}={\int }_{\mathcal{G}}{\varrho }_r^{k,\ell }\left(du\right).$$
• ${\varrho }_r^{k,\ell ,q}\left(du\right)$ is a conditional distribution of $X_{t_r}^{\ell }{\mathbf{I}}_{\left\{q\right\}}\left(N_{r+1}^X\right){\tau }_{r+1}$ given $X_{t_{r+1}}=e_k$, i.e.,  for any $\mathcal{G}\in \mathcal{B}\left({\mathbb{R}}^M\right)$ the following equality holds:

  $${\mathsf{E}}_{}\left\{{\mathbf{I}}_{\mathcal{G}}\left({\tau }_{r+1}\right){\mathbf{I}}_{\left\{q\right\}}\left(N_{r+1}^X\right)X_{t_r}^{\ell }|X_{t_{r+1}}=e_k\right\}={\int }_{\mathcal{G}}{\varrho }_r^{k,\ell ,q}\left(du\right).$$
• $\mathcal{N}(y,m,K)\triangleq {\left(2\pi \right)}^{-M/2}{\det }^{-1/2}Kexp\left\{-\frac{1}{2}{\parallel y-m)\parallel }_{K^{-1}}^2\right\}$ is an M-dimensional Gaussian probability density function (pdf) with the expectation m and nondegenerate covariance matrix K.
• $\xi \triangleq \mathrm{vec}(\Lambda ,f,{\left\{g_n\right\}}_{\overline{1,N}})$ is an $\left(N^2+NM+N{M}^2\right)$-dimensional vector of all HMM parameters.
• ${\mathsf{Y}}_t\triangleq \sigma \{Y_s,0⩽s⩽t\}$ is a natural filtration generated by the continuous-time process $Y_s$ available up to the moment t.

3. Auxiliary Problem: Fixed Interval Smoothing Given Discretized Observations

The proposed HMM parameter identification algorithm uses the iterative calculation of the fixed-interval smoothing estimates ${\left\{{\widehat{X}}_r^s\right\}}_{r=\overline{1,R}}$: ${\widehat{X}}_r^s\triangleq {\mathsf{E}}_{}\left\{X_{t_r}|{\mathfrak{Y}}_R\right\}$. For the effective numerical realization of the smoothing procedure, we suggest the two-filtering algorithm [37,38,39].

First, we present an algorithm for the calculation of the forward-time filtering estimate ${\widehat{X}}_t\triangleq {\mathsf{E}}_{}\left\{X_{t_r}|{\mathfrak{Y}}_r\right\}$.

Let us introduce the following positive random numbers and matrices:

$$\begin{array}{cc}{\displaystyle {\theta }_r^{kj}\triangleq {\int }_{\mathcal{D}}\mathcal{N}({\mathcal{Y}}_r,fu,\sum _{p=1}^N{u}^p{g}_p){\rho }_r^{k,j}\left(du\right),}\hfill & {\theta }_r\triangleq {\parallel {\theta }_r^{kj}\parallel }_{k,j=\overline{1,N}}.\hfill \end{array}$$

(8)

The assertion below gives a recursive algorithm to calculate ${\widehat{X}}_t$.

Proposition 1.

The Conditional Expectation (CE) ${\widehat{X}}_r$ can be calculated by the formula

$${\widehat{X}}_r={\left(\mathbf{1}{\theta }_r^{\top }{\widehat{X}}_{r-1}\right)}^{-1}{\theta }_r^{\top }{\widehat{X}}_{r-1},r>0,$$

(9)

with the initial condition

$${\widehat{X}}_0=\pi .$$

(10)

The proof of Proposition 1 is quite similar to one of Lemma 3 in [33].

Second, we introduce the backward-time martingale representation of the MJP $X_t$ along with the backward-time filtering estimate. On the initial probability triplet $(\Omega ,\mathcal{F},\mathcal{P})$ we consider

• The backward-time filtration ${\left\{{\mathcal{F}}_t^r\right\}}_{t\in [0,T]}$: ${\mathcal{F}}_t^r\triangleq \sigma \{X_s,Y_T-Y_s:t⩽s⩽T\}$.
• The backward-time collection of the discretized observations ${\mathfrak{Y}}_r^b\triangleq \sigma \{{\mathcal{Y}}_q:r<q⩽R\}$; ${\mathfrak{Y}}_R^b\triangleq \{\varnothing ,\Omega \}$.

The probability distribution $\pi \left(t\right)={\mathsf{E}}_{}\left\{X_t\right\}$, being the unique solution to the Kolmogorov equation $\pi \left(t\right)=\pi +{\int }_0^t{\Lambda }^{\top }\pi \left(s\right)ds$, has strictly positive components ${\pi }^n\left(t\right)>0$ for any $n=\overline{1,N},t>0$ due to condition 3 above (see [40]). Hence, we can define the backward-time TIM for the initial MJP $X_t$ as follows

$$\Gamma \left(t\right)\triangleq {diag}^{-1}\pi \left(t\right){\Lambda }^{\top }diag\pi \left(t\right)-{diag}^{-1}\pi \left(t\right)diag\left({\Lambda }^{\top }\pi \left(t\right)\right).$$

(11)

The assertion below declares a backward-time representation of the MJP $X_t$.

Proposition 2.

The MJP $X_t$ (1) satisfies the backward time SDS

$$X_t=X_T+{\int }_t^T{\Gamma }^{\top }\left(s\right)X_{s}ds+\left[M_t^X-M_T^X-{\int }_t^T{\left(\Gamma \left(s\right)+\Lambda \right)}^{\top }{X}_{s}ds\right].$$

(12)

The process

$$M_t^b\triangleq M_t^X-M_T^X-{\int }_t^T{\left(\Gamma \left(s\right)+\Lambda \right)}^{\top }{X}_{s}ds$$

represents ${\mathcal{F}}_t^r$-adapted square integrable martingale:

$${\mathsf{E}}_{}\left\{M_t^b|{\mathcal{F}}_s^r\right\}=M_s^b\mathcal{P}-a.s. for any0⩽t⩽s⩽T.$$

The proof of Proposition 2 can be found in [41].

To define the backward-time filtering estimate, we introduce the positive random numbers and matrices:

$$\begin{array}{cc}{\displaystyle {\vartheta }_r^{kj}\triangleq {\int }_{\mathcal{D}}\mathcal{N}({\mathcal{Y}}_{r+1},fu,\sum _{p=1}^N{u}^p{g}_p){\varrho }_r^{k,j}\left(du\right),}\hfill & {\vartheta }_r\triangleq {\parallel {\vartheta }_r^{kj}\parallel }_{k,j=\overline{1,N}}.\hfill \end{array}$$

(13)

One can obtain the backward-time filtering estimates ${\widehat{X}}_r^b\triangleq {\mathsf{E}}_{}\left\{X_{t_r}|{\mathfrak{Y}}_r^b\right\}$, which are calculated by an algorithm similar to the original forward-time one

$${\widehat{X}}_r^b={\left(\mathbf{1}{\vartheta }_r^{\top }{\widehat{X}}_{r+1}^b\right)}^{-1}{\vartheta }_r^{\top }{\widehat{X}}_{r+1}^b,r<R,$$

(14)

with the initial condition

$${\widehat{X}}_R^b=\pi \left(T\right).$$

(15)

Third, in the assertion below, we present a numerically efficient two-filtering algorithm for the calculation of the fixed-interval smoothing MJP state estimate.

Proposition 3.

The fixed-interval smoothing estimate ${\widehat{X}}_r^s=col({\widehat{X}}_r^{1s},\dots ,{\widehat{X}}_r^{Ns})={\mathsf{E}}_{}\left\{X_{t_r}|{\mathfrak{Y}}_R\right\}$, $r=\overline{0,R}$, has the form

$${\widehat{X}}_r^{ns}=\frac{{\tilde{X}}_r^{ns}}{{\sum }_{j=1}^N{\tilde{X}}_r^{js}},n=\overline{1,N},$$

(16)

where the unnormalized estimate ${\tilde{X}}_r^s\triangleq col({\tilde{X}}_r^{1s},\dots ,{\tilde{X}}_r^{Ns})$ is calculated via the filtering estimates ${\widehat{X}}_r$, ${\widehat{X}}_r^b$ and prior MJP distribution $\pi \left(t_r\right)$ as

$${\tilde{X}}_r^{ns}=\frac{{\widehat{X}}_r^n{\widehat{X}}_r^{nb}}{{\pi }^n\left(t_r\right)},n=\overline{1,N}.$$

(17)

The proof of Proposition 3 is similar to the one of Theorem 4 in [39].

4. EM-Like Algorithm of HMM Parameter Identification

Before the description of the parameter identification algorithm, we suggest some hints for it.

Let us consider the continuous-time HMM (1), (3), for which both process $X_t$ and $Y_t$ are observable. Due to the strong law of large numbers for MJP [42] and regenerative processes [43] the following strong limits hold as $T\to \infty$:

$$\begin{array}{c}{\displaystyle \frac{N_T^{ij}}{O_T^i}\to {\lambda }_{ij}\mathrm{for}\mathrm{any}i,j=\overline{1,N}:i\ne j,}\end{array}$$

(18)

$$\begin{array}{c}\frac{1}{O_T^n}{Q}_T^{n,1}\to f{e}_n\mathrm{for}\mathrm{any}n=\overline{1,N}.\end{array}$$

(19)

The identity

$$\frac{1}{O_T^n}\left[Q_T^{n,1}{\left(Q_T^{n,1}\right)}^{\top }-{\int }_0^T{Q}_s^{n,1}d{\left(Q_s^{n,1}\right)}^{\top }-{\int }_0^{T}d{Q}_s^{n,1}{\left(Q_s^{n,1}\right)}^{\top }\right]=g_n$$

(20)

is also true $\mathcal{P}$-a.s. for any $n=\overline{1,N}$ on the set $\{\omega \in \omega :O_T^n>0\}$. In the above formulae,

• The processes $N_t^{ij}$ and $O_t^n$ are determined by formulae (5) and (6).
• The “level sum” $Q_t^{n,1}$ is a particular case of (7):

  $$Q_t^{n,1}={\int }_0^t{X}_s^{n}d{Y}_s.$$

  (21)

One can treat identification formulae (18) and (19) as some application of the method of moments [44]. Furthermore, [33], if  $X_t$ is the unobservable state and identifiability condition 5 above holds, then it is possible to restore $X_t$ precisely given the diffusion observations $Y_t$:

$${\overline{X}}_t\triangleq {\mathsf{E}}_{}\left\{X_t|{\mathsf{Y}}_{t+}\right\}=X_t\mathcal{P}-a.s.\mathrm{for}\mathrm{any}0⩽t<T.$$

(22)

Thus, in “the ideal situation” the CE ${\overline{X}}_t$ allows identifying the HMM parameters $\xi$. However, two facts block the direct implementation of this idea. The estimate ${\mathsf{E}}_{}\left\{X_t|{\mathsf{Y}}_{t+}\right\}$ represents a solution to the SDS with both the diffusion and counting processes on the right-hand side. The SDS uses some transformation of the original observations (3), which is a result of some double limit passage. There is no effective numerical procedure for its realization. Further, the CE ${\overline{X}}_t={\overline{X}}_{\xi }(\omega ,t)$ depends on the exact values of $\Lambda$, f and $\left\{g_n\right\}$, which are unknown.

We perform several steps to overcome this obstacle and develop a suboptimal identification algorithm. First, under rather mild conditions [45], there exists a continuous dependency of the solution to the SDS on the parameters, i.e., if the sequence of parameters ${\left\{{\xi }_j\right\}}_{j\in \mathbb{N}}$ converges to $\xi$ as $j\to \infty$, then the sequence of the corresponding SDS solutions ${\overline{X}}_{{\xi }_j}\stackrel{p}{\to }{\overline{X}}_{\xi }$ as $j\to \infty$. This means the performance of the estimate ${\overline{X}}_{{\xi }_j}$ is close to the CE ${\overline{X}}_{\xi }$ when the used parameters ${\xi }_j$ are close enough to the exact values $\xi$.

Second, in the case of the known parameters $\xi$ and the sequence of the nested partitions ${\left\{t_r\right\}}_{r=\overline{1,R}}$, due to Levy theorem [40], ${\widehat{X}}_r\to {\overline{X}}_{t_r}$$\mathcal{P}-$ a.s. as $R\to \infty$. This means that, under the identifiability condition 5 above, the optimal filtering estimate ${\widehat{X}}_r$ calculated by the discretized observations $\mathcal{Y}$ is close to the exact MJP state $X_{t_r}$, and the optimal smoothed estimate ${\widehat{X}}_r^s$ is even more so.

Third, both the trajectories of the state $X_{t_r}$ and CE ${\overline{X}}_t$ given the diffusion observations $Y_t$ (3) are piecewise constant ones with values in ${\mathbb{S}}^N$. By contrast, the CE ${\widehat{X}}_r$ calculated by the discretized observations $\mathcal{Y}$ (2) belongs to the “probabilistic simplex” $\Pi$. One can neutralize this issue rather easily, replacing the estimate ${\widehat{X}}_r$ with the maximum a posteriori probability (MAP) estimate ${\mathcal{X}}_r$:

$${\widehat{\mathcal{X}}}_r={argmax}_n\left(e_n^{\top }{\widehat{X}}_r\right).$$

(23)

The CE ${\widehat{X}}_r$ is the ${\mathcal{L}}_2$-optimal estimate, whereas ${\widehat{\mathcal{X}}}_r$ is the ${\mathcal{L}}_1$-optimal one.

Fourth, Formulae (18)–(20) with the processes $N^{i\ell }$, $O^i$ and $Q^{i,1}$ substituted by their CEs, can be the base for the iterative identification procedure. Indeed, on the first step of the loop, the smoothed estimate ${\widehat{\mathcal{X}}}_{{\xi }_j}$ is calculated given the previous values of the HMM parameters ${\xi }_j$. Based on ${\widehat{\mathcal{X}}}_{{\xi }_j}$, one can compute the required auxiliary estimates ${\widehat{N}}^{ik}$, ${\widehat{O}}^n$ and ${\widehat{Q}}^{n,1}$:

$${\widehat{N}}_R^{ik}=\sum _{r=1}^R{\widehat{\mathcal{X}}}_{r-1}^{is}({\widehat{\mathcal{X}}}_r^{ks}-{\widehat{\mathcal{X}}}_{r-1}^{ks}),$$

(24)

$${\widehat{O}}_R^n=\delta \sum _{r=1}^R{\widehat{\mathcal{X}}}_{r-1}^{ns},$$

(25)

$${\widehat{Q}}_R^{n,1}=\sum _{r=1}^R{\widehat{\mathcal{X}}}_{r-1}^{ns}{\mathcal{Y}}_r.$$

(26)

Using them, in the second step, one can recalculate the HMM parameters ${\xi }_{j+1}$ by the modifications of (18), (19) and (20):

$${\widehat{\lambda }}_{ik}=\frac{{\widehat{N}}_R^{i,k}}{{\widehat{O}}_R^i},$$

(27)

$$\widehat{f}{e}_n=\frac{1}{{\widehat{O}}_R^n}{\widehat{Q}}_R^{n,1},$$

(28)

$${\widehat{g}}_n=\frac{1}{{\widehat{O}}_R^n}\sum _{r=1}^R{\widehat{\mathcal{X}}}_r^{ns}({\mathcal{Y}}_r-\delta \widehat{f}{e}_n){({\mathcal{Y}}_r-\delta \widehat{f}{e}_n)}^{\top }.$$

(29)

The authors of [2,29,30] presented the versions of EM algorithm for the continuous-time HMM (1), (3) parameters in the case of additive observation noise ($g_n\equiv g$). It is shown that the formulae (18) and (19) establish the maximization step. Moreover, in [29,30], it is shown that, under the condition ${\mathcal{P}}_{{\xi }^{\prime }}\sim {\mathcal{P}}_{{\xi }^{\prime \prime }}$ for all admissible parameter values $\xi \prime$ and $\xi ″$, the recursive identification procedure converges to a local maximum of the conditional likelihood function.

Below (see Algorithm 1) is a sketch of the EM like algorithm of the discrete-continuous HMM (1), (2) parameter identification. The input parameters include the initial values of parameters $\stackrel{\circ }{\xi }$, maximal number of iterations $N_{iter}$ and minimal relative parameter variation per iteration ${\epsilon }_0$.

Algorithm 1: The algorithm for HMM parameter identification.

1: $\epsilon \leftarrow {\epsilon }_0$            ▹ Required accuracy index: initialization 2: $N_{iter}\leftarrow N_0$        ▹ Maximal number of iterations: initialization    3: for $n_1=1$ to N do              ▹ Parameter initialization: begin 4:       for $m_1=1$ to M do 5:           ${\widehat{f}}_{cur}^{m_1,n_1}\leftarrow {\stackrel{\circ }{f}}^{m_1,n_1}$ 6:           for $m_1=2$ to M do 7:              ${\widehat{g}}_{cur}^{m_1,m_2}\leftarrow {\stackrel{\circ }{g}}^{m_1,m_2}$ 8:           end for 9:       end for 10:     for $n_2=1$ to N do If $(n_1<>n_2)$ 11:         ${\widehat{\lambda }}_{cur}^{n_1,n_2}\leftarrow {\stackrel{\circ }{\lambda }}^{n_1,n_2}$ 12:     end for 13:     ${\widehat{\lambda }}_{cur}^{n_1,n_2}\leftarrow -{\sum }_{m:m\ne n_1}{\widehat{\lambda }}_{cur}^{n_1,m}$ 14: end for                  ▹ Parameter initialization: end       15: $N_{cur}\leftarrow 0$    16: repeat                    ▹ Identification loop: begin    17:     ${\widehat{X}}_0\leftarrow \pi$ 18:     ${\widehat{X}}_0^b\leftarrow \pi \left(T\right)$ 19:     $\widehat{\xi }\leftarrow {\widehat{\xi }}_{cur}$ 20:     $N_{cur}\leftarrow N_{cur}+1$    21:     for $r=1$ to R do            ▹ Forward and backward time filtering: begin 22:         ${\widehat{X}}_r\leftarrow \widehat{X}({\widehat{X}}_{r-1},{\mathcal{Y}}_r,\widehat{\xi })$ formulae (8), (9) 23:         ${\widehat{X}}_{R-r}^b\leftarrow {\widehat{X}}^b({\widehat{X}}_{R-r+1}^b,{\mathcal{Y}}_{R-r+1},\widehat{\xi })$ formulae (13), (14) 24:     end for               ▹ Forward and backward time filtering: end    25:     for $r=0$ to R do        ▹ Smoothing and CE to MAP conversion: begin 26:         ${\widehat{X}}_r^s\leftarrow {\widehat{X}}^s({\widehat{X}}_r,{\widehat{X}}_r^b)$ formulae (16), (17) 27:         ${\widehat{\mathcal{X}}}_r^s\leftarrow {\widehat{\mathcal{X}}}^s\left({\widehat{X}}_r^s\right)$ formula (23) 28:     end for             ▹ Smoothing and CE to MAP conversion: end    29:     for $n_1=1$ to N do          ▹ HMM Parameter recalculation: begin 30:         for $m_1=1$ to M do 31:            ${\widehat{f}}_{cur}^{m_1,n_1}\leftarrow {\widehat{f}}_{cur}^{m_1,n_1}({\widehat{\mathcal{X}}}^s,\mathcal{Y})$ formulae (24)–(26), (28) 32:            for $m_2=1$ to M do 33:                ${\widehat{g}}_{cur}^{m_1,m_2}\leftarrow {\widehat{g}}_{cur}^{m_1,m_2}({\widehat{\mathcal{X}}}^s,\mathcal{Y})$ formulae (24)–(26), (29) 34:            end for 35:         end for 36:         for $n_2=1$ to N do If($n_1<>n_2$) 37:            ${\widehat{\lambda }}_{cur}^{n_1,n_2}\leftarrow {\widehat{\lambda }}_{cur}^{n_1,n_2}({\widehat{\mathcal{X}}}^s,\mathcal{Y})$ formulae (24)–(27) 38:         end for 39:         ${\widehat{\lambda }}_{cur}^{n_1,n_2}\leftarrow -{\sum }_{m:m\ne n_1}{\widehat{\lambda }}_{cur}^{n_1,m}$ 40:     end for                ▹ HMM Parameter recalculation: end    41: until$\left(\left(\frac{\parallel {\widehat{\xi }}_{cur}-\widehat{\xi }\parallel }{\parallel {\widehat{\xi }}_{cur}\parallel }⩽{\epsilon }_0\right)\mathrm{OR}(N_{cur}>=N_{iter})\right)$                ▹ Identification loop: end   return ${\widehat{\xi }}_{cur}$

5. Numerical Schemes of Filtering Estimate Calculation

Neither the forward-time conditional distribution ${\rho }^{k,\ell }$ nor the backward-time one ${\varrho }^{k,\ell }$ are absolutely continuous with respect to the Lebesgue measure; hence, the practical calculation of the integrals (8) and (13) is nontrivial. Due to the law of the total probability,

$${\theta }_r^{kj}\triangleq {\int }_{\mathcal{D}}\mathcal{N}({\mathcal{Y}}_r,fu,\sum _{p=1}^N{u}^p{g}_p){\rho }_r^{k,j}\left(du\right)=\sum _{q=0}^{\infty }{\int }_{\mathcal{D}}\mathcal{N}({\mathcal{Y}}_r,fu,\sum _{p=1}^N{u}^p{g}_p){\rho }_r^{k,j,q}\left(du\right).$$

The results of [33] hint that, for small enough time increments $\delta$, the sum in the last formula can be limited by only two summands ($q=0$ and 1). In this case, we replace the estimate ${\widehat{X}}_r$ with its analytic approximation of order 1:

$${\widehat{X}}_r={\left(\mathbf{1}{\varkappa }_r^{\top }{\widehat{X}}_{r-1}\right)}^{-1}{\varkappa }_r^{\top }{\widehat{X}}_{r-1},r>0,{\widehat{X}}_0=\pi ,$$

(30)

where

$$\begin{array}{cc}{\displaystyle {\varkappa }_r^{kj}\triangleq \sum _{q=0}^1{\int }_{\mathcal{D}}\mathcal{N}({\mathcal{Y}}_r,fu,\sum _{p=1}^N{u}^p{g}_p){\rho }_r^{k,j,q}\left(du\right),}\hfill & {\varkappa }_r\triangleq {\parallel {\varkappa }_r^{kj}\parallel }_{k,j=\overline{1,N}}.\hfill \end{array}$$

(31)

From the probabilistic point of view, this means that we take into account at most one state transition occurring at the discretization interval $[t_{r-1},t_r]$. If $\overline{\lambda }\triangleq {max}_{1⩽n⩽N}\left|{\lambda }_{nn}\right|$ and $\delta$ is so small that $\overline{\lambda }\delta <1$, then (see [33])

$$C_1\triangleq e^{-\overline{\lambda }\delta }\frac{2}{{\left(\overline{\lambda }\delta \right)}^2}{\displaystyle \sum _{k=2}^{\infty }}\frac{{\left(\overline{\lambda }\delta \right)}^k}{k!}<1,$$

(32)

and the filtering approximation error is characterized by the inequality

$${\mathsf{E}}_{}\left\{{∥{\widehat{X}}_r-{\widehat{X}}_r∥}_1\right\}⩽2-2{\left(1-C_1\frac{{\left(\overline{\lambda }\delta \right)}^2}{2}\right)}^r\sim {\overline{\lambda }}^2{t}_r\delta as\delta \to 0.$$

The explicit form of ${\varkappa }_r^{kj}$ is

$${\varkappa }_r^{kj}={\delta }_{kj}{e}^{{\lambda }_{jj}\delta }\mathcal{N}({\mathcal{Y}}_r,\delta f^j,\delta g_j)+(1-{\delta }_{kj}){\lambda }_{kj}{e}^{{\lambda }_{jj}\delta }{\int }_0^{\delta }{V}^{kj}({\mathcal{Y}}_r,u)du,$$

(33)

where ${\delta }_{kj}$ is the Kronecker symbol, and

$$V^{kj}(y,u)\triangleq e^{({\lambda }_{kk}-{\lambda }_{jj})u}\mathcal{N}(y,u{f}^k+(\delta -u)f^j,u{g}_k+(\delta -u)g_j).$$

(34)

The integral in (33) cannot be calculated analytically; thus, as the analytic approximation ${\widehat{X}}_r$, we have to calculate it numerically, appending some error. We suggest using the composite method of the midpoint rectangles with the substep ${\delta }^{1+\alpha }$ for some $\alpha >0$:

$${\varkappa }_r^{kj}\approx {\psi }_r^{kj}\triangleq {\delta }_{kj}{e}^{{\lambda }_{jj}\delta }\mathcal{N}({\mathcal{Y}}_r,\delta f^j,\delta g_j)+(1-{\delta }_{kj}){\lambda }_{kj}{\delta }^{1+\alpha }\sum _{i=1}^{\left[{\delta }^{-\alpha }\right]}{V}^{kj}\left({\mathcal{Y}}_r,{\delta }^{1+\alpha }\left(i-\frac{1}{2}\right)\right).$$

(35)

Denoting ${\psi }_r\triangleq {\parallel {\psi }_r^{kj}\parallel }_{k,j=\overline{1,N}}$, we present the numerical approximation of order $min(\alpha ,1)$ based on the composite midpoint rectangle scheme:

$${\widehat{\mathbf{X}}}_r={\left(\mathbf{1}{\psi }_r^{\top }{\widehat{\mathbf{X}}}_{r-1}\right)}^{-1}{\psi }_r^{\top }{\widehat{\mathbf{X}}}_{r-1},r>0,{\widehat{\mathbf{X}}}_0=\pi .$$

(36)

The total approximation error of ${\widehat{\mathbf{X}}}_r$ is characterized by the inequality [34]

$${\mathsf{E}}_{}\left\{{∥{\widehat{X}}_r-{\widehat{\mathbf{X}}}_r∥}_1\right\}⩽C_2{t}_r(\delta +C_3{\delta }^{\alpha })$$

with some positive constants $C_2$ and $C_3$. Thus, the order of the numerical approximation is at most 1.

Above in this section, we discuss the numerical approximation of the forward-time filtering estimate. A similar framework is applicable for the numerical realization of backward-time filtering. As in the forward-time, we assume at most one MJP state transition occurred on the time increment. We also have to consider the nonhomogeneous nature of the MJP $X_t$ viewed in the reversed time. The analogue of formula (33) has the form

$${\zeta }_r^{kj}={\delta }_{kj}exp\left({\int }_{t_r}^{t_{r+1}}{\gamma }_{kk}\left(u\right)du\right)\mathcal{N}({\mathcal{Y}}_{r+1},\delta f^k,\delta g_k)+(1-{\delta }_{kj}){\int }_0^{\delta }{U}_r^{kj}({\mathcal{Y}}_{r+1},u)du,$$

(37)

where

$$\begin{array}{c}{U}_r^{kj}(y,u)\triangleq {\gamma }_{kj}\left(u\right)exp\left({\int }_{t_r}^u{\gamma }_{kk}\left(s\right)ds+{\int }_u^{t_{r+1}}{\gamma }_{jj}\left(s\right)ds\right)\times \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill \\ \hfill \times \mathcal{N}(y,u{f}^k+(\delta -u)f^j,u{g}_k+(\delta -u)g_j).\end{array}$$

(38)

Both the first summand in (37) and the function $U_r^{kj}$ (38) contain integrals; therefore, we replace them with the following midpoint approximations:

$$exp\left({\int }_{t_r}^{t_{r+1}}{\gamma }_{kk}\left(u\right)du\right)\approx exp\left(\delta {\gamma }_{kk}\left(\frac{t_r+t_{r+1}}{2}\right)\right),$$

$$\begin{array}{c}{\tilde{U}}_r^{kj}(y,u)\triangleq {\gamma }_{kj}\left(u\right)exp\left((u-t_r){\gamma }_{kk}\left(\frac{u+t_r}{2}\right)+(t_{r+1}-u){\gamma }_{jj}\left(\frac{t_{r+1}+u}{2}\right)\right)\times \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill \\ \hfill \times \mathcal{N}(y,u{f}^k+(\delta -u)f^j,u{g}_k+(\delta -u)g_j).\end{array}$$

(39)

Finally, the backward-time analogue of ${\psi }^{kj}$, which is valid for the numerical realization, takes the form

$${\zeta }_r^{kj}\approx {\varphi }_r^{kj}={\delta }_{kj}exp\left(\delta {\gamma }_{kk}\left(\frac{t_r+t_{r+1}}{2}\right)\right)+(1-{\delta }_{kj}){\delta }^{1+\alpha }\sum _{i=1}^{\left[{\delta }^{-\alpha }\right]}{\tilde{U}}^{kj}\left({\mathcal{Y}}_{r+1},{\delta }^{1+\alpha }\left(i-\frac{1}{2}\right)\right).$$

(40)

Denoting ${\varphi }_r\triangleq {\parallel {\varphi }_r^{kj}\parallel }_{k,j=\overline{1,N}}$, we present the numerical approximation of order $min(\alpha ,1)$ based on the composite midpoint rectangle scheme for backward-time filter:

$${\widehat{\mathbf{X}}}_r^b={\left(\mathbf{1}{\varphi }_r^{\top }{\widehat{\mathbf{X}}}_{r+1}^b\right)}^{-1}{\varphi }_r^{\top }{\widehat{\mathbf{X}}}_{r+1}^b,r<R,{\widehat{\mathbf{X}}}_R^b=\pi \left(T\right).$$

(41)

6. Numerical Experiments

To illustrate the properties of the proposed estimation algorithms, we performed a complex of numerical experiments. We fixed the values for some basic parameters of HMM (1), (2). Varying some of them, we fulfilled the calculation, and the obtained results help us to highlight various performance characteristics of the suggested estimates.

We chose an HMM with the following parameters: $N=3$, $M=1$,

$$\Lambda =\left[\begin{array}{ccc}-18.0& 12.0& 6.0\\ 9.0& -18.0& 9.0\\ 6.0& 12.0& -18.0\end{array}\right],\pi =\left[\begin{array}{c}0.3\\ 0.4\\ 0.3\end{array}\right],f=[0.0,0.0,0.0],\begin{array}{c}{g}_1=0.1,\\ g_2=0.2,\\ g_3=0.3.\end{array}$$

(42)

6.1. Filtering vs. Smoothing: Estimation Performance Analysis

The HMM state filtering numerical algorithm is a backbone of subsequent HMM parameter identification; thus, the high precision of the proposed filtering estimates is a vital requirement.

The optimal filtering algorithm of the discrete-time HMM states given the observations with state-dependent noises is presented in ([2], chap. 3). Its adaptation to the time-discretized observations is a specific case of (30), with the matrix $\theta$ calculated other than by (35).

We compare the proposed filtering and smoothing estimates, calculated by the composite midpoint rectangle scheme and the ones from [2]. Both numerical schemes are stable. The difference can be exposed clearly, when the time-discretization step is sufficiently large. The numerical experiment was conducted for the HMM with the basic parameter values (42) and the following auxiliary parameter values: the time interval $[0,T]=[0,1.0]$, the time-discretization step $\delta =0.1$, the MJP state simulation step ${\delta }_{sim}=0.001$ and the step of the numerical integration in (35) ${\delta }_{int}=0.01$. The estimation error variance was calculated by the Monte-Carlo method over the sample of the size $N_{MC}={10}^8$.

Figure 1 contains the following plots:

• ${\widehat{D}}_t$, ${\widehat{D}}_t^b$, ${\widehat{D}}_t^s$—the sample variances of the forward-time filtering estimate, backward-time and fixed-interval smoothing one, calculated by the algorithms proposed in the paper.
• ${\tilde{D}}_t$, ${\tilde{D}}_t^b$, ${\tilde{D}}_t^s$—the sample variances of the forward-time filtering estimate, backward-time and fixed-interval smoothing one, calculated by the algorithms proposed in [2].

Figure 1. Sample error variances of the proposed estimates and the ones from [2].

Figure 1. Sample error variances of the proposed estimates and the ones from [2].

Figure 2. Available discretized observations, true MJP state and its filtering and smoothing estimates.

Figure 2. Available discretized observations, true MJP state and its filtering and smoothing estimates.

The chosen initial distribution of the MJP $X_t$ coincides with the stationary one, hence, the vector of component variances is $\left[0.21,0.24,0.21\right]$. The combination of the high transition intensity and the large discretization step (${max}_j\left(\right|{\lambda }_{jj}\left|\right)\delta \approx 1$), along with the absence of the observation drift (f is a zero vector) and the slightly different noise intensities $\left\{g_n\right\}$ establish poor conditions for the state estimation; hence, one can expect for the variances of estimate errors to be close to the ones of the MJP itself. In this numerical experiment, both the filtering and smoothing algorithms are better than the adaptation of the discrete-time algorithms [2]. Moreover, looking at the sample variances of the second MJP component, one can see that the computational errors become so high that the smoothing demonstrates worse performance than the filtering does.

To compare the performance of the filtering estimates and the smoothing one, the next numerical experiment is fulfilled with the following time steps: $\delta =0.2\times {10}^{-4}$, ${\delta }_{sim}=0.2\times {10}^{-6}$ and ${\delta }_{int}=0.2\times {10}^{-5}$.

Figure 2 presents specific trajectories of the available observations, components of the accurate MJP state, their forward- and backward-time filtering estimates and the fixed-interval smoothing ones. One can see the superiority of the smoothing estimate performance.

Condition 5 is valid for the HMM with parameters (42); hence, we can expect convergence of both the filtering and smoothing estimates to the accurate state value as the time discretization step vanishes. Figure 3 presents sample error variances for all types of estimates (forward-time filtering ${\widehat{\mathbf{D}}}_t$, backward-time-filtering ${\widehat{\mathbf{D}}}_t^b$ and fixed-interval smoothing ones ${\widehat{\mathbf{D}}}_t^s$) calculated with various time discretization steps $\delta ={10}^{-2},{10}^{-3},{10}^{-4}$. One can see the expected convergence.

The results of the first stage of the numerical experiment allow us to make the following conclusions.

• Both the proposed numerical approximation based on the composite midpoint rectangle scheme (35) and the adaptation of the algorithm from [2] demonstrate the stable character: both the estimates have non-negative components and satisfy the normalization property. Their performances look similar; however, (35) demonstrates more stability in the case of the large discretization step (${max}_j\left(\right|{\lambda }_{jj}\left|\right)\delta \approx 1$).
• The drift f is absent in the observations. However, all conditional noise intensities are different; therefore, identifiability condition 5 is valid for this HMM. The calculation results confirm the convergence of the proposed estimates to the real MJP state as $\delta \to 0$.
• The quality of the fixed-interval smoothing estimates is significantly better than the one of the filtering estimates: the error variances are distinguished about twice.

6.2. Parameter Identification

We performed the complex of the identification numerical experiments for the HMM with parameters (42), where the observation drift vector f is replaced by $f=\left[-10.0,5.0,20.0\right]$. This change allows us to illustrate the identification performance of all HMM parameters, including the drift f. The observations are available on the time interval $[0.0,100.0]$. We simulate the “continuous-time” state with the small step ${\delta }_{sim}={10}^{-7}$, and use the step ${\delta }_{int}={10}^{-6}$ in the numerical integration scheme (40).

To analyze how the time discretization of observations affects the parameter identification performance, we choose the following time step values: $\delta ={10}^{-2},2\times {10}^{-4},{10}^{-4}$and ${10}^{-5}$. The number of iteration is bounded: $N_{iter}=30$ with the minimal relative parameter variation per iteration ${\epsilon }_0={10}^{-5}$.

Table 1 contains the identification results of the TIM $\Lambda$, Table 2 stands for the results of the drift f identification, and Table 3 stands for the ones of the diffusion g. All the tables include

• The exact parameter values ($\Lambda$, f and g).
• The starting points for the iterative identification procedure (${\widehat{\Lambda }}_0$, ${\widehat{f}}_0$ and ${\widehat{g}}_0$).
• The obtained estimates ($\widehat{\Lambda }$, $\widehat{f}$ and $\widehat{g}$).

Table 1. Identification results of the TIM $\Lambda$.

$\mathit{\delta }$	$\Lambda$			${\tilde{\Lambda }}_{\mathit{X}}$			${\widehat{\Lambda }}_0$			$\widehat{\Lambda }$			Relative Errors (%)
	−18.0	12.0	6.0	−18.6	11.6	7.0	−1.0	0.5	0.5	−14.6	7.4	7.2	18.9	38.3	20.0
10−2	9.0	−18.0	9.0	8.7	−17.2	8.5	0.5	−1.0	0.5	10.0	−16.0	6.0	11.1	11.1	33.3
	6.0	12.0	−18.0	6.0	11.2	−17.2	0.5	0.5	−1.0	9.7	5.3	−15.0	61.7	55.8	16.7
	−18.0	12.0	6.0	−18.9	11.8	7.1	−1.0	0.5	0.5	−16.6	10.1	6.5	7.8	15.8	8.3
2 × 10−4	9.0	−18.0	9.0	8.9	−17.5	8.6	0.5	−1.0	0.5	7.6	−14.6	7.0	15.6	18.9	22.2
	6.0	12.0	−18.0	6.0	11.5	−17.5	0.5	0.5	−1.0	5.5	9.8	−15.3	8.3	18.3	15.0
	−18.0	12.0	6.0	−18.9	11.8	7.1	−1.0	0.5	0.5	−17.4	11.1	6.3	3.3	7.5	5.0
10−4	9.0	−18.0	9.0	8.9	−17.5	8.6	0.5	−1.0	0.5	8.2	−15.6	7.4	8.9	13.3	17.8
	6.0	12.0	−18.0	6.0	11.5	−17.5	0.5	0.5	−1.0	5.6	9.8	−15.4	6.7	18.3	14.4
	−18.0	12.0	6.0	−18.9	11.8	7.1	−1.0	0.5	0.5	−18.4	11.4	7.0	2.2	5.0	16.7
10−5	9.0	−18.0	9.0	8.9	−17.5	8.6	0.5	−1.0	0.5	8.5	−16.9	8.4	5.6	6.1	6.7
	6.0	12.0	−18.0	6.0	11.5	−17.5	0.5	0.5	−1.0	6.1	11.0	−17.1	1.7	8.3	5.0

Table 1. Identification results of the TIM $\Lambda$.

$\mathit{\delta }$	$\Lambda$			${\tilde{\Lambda }}_{\mathit{X}}$			${\widehat{\Lambda }}_0$			$\widehat{\Lambda }$			Relative Errors (%)
	−18.0	12.0	6.0	−18.6	11.6	7.0	−1.0	0.5	0.5	−14.6	7.4	7.2	18.9	38.3	20.0
10−2	9.0	−18.0	9.0	8.7	−17.2	8.5	0.5	−1.0	0.5	10.0	−16.0	6.0	11.1	11.1	33.3
	6.0	12.0	−18.0	6.0	11.2	−17.2	0.5	0.5	−1.0	9.7	5.3	−15.0	61.7	55.8	16.7
	−18.0	12.0	6.0	−18.9	11.8	7.1	−1.0	0.5	0.5	−16.6	10.1	6.5	7.8	15.8	8.3
2 × 10−4	9.0	−18.0	9.0	8.9	−17.5	8.6	0.5	−1.0	0.5	7.6	−14.6	7.0	15.6	18.9	22.2
	6.0	12.0	−18.0	6.0	11.5	−17.5	0.5	0.5	−1.0	5.5	9.8	−15.3	8.3	18.3	15.0
	−18.0	12.0	6.0	−18.9	11.8	7.1	−1.0	0.5	0.5	−17.4	11.1	6.3	3.3	7.5	5.0
10−4	9.0	−18.0	9.0	8.9	−17.5	8.6	0.5	−1.0	0.5	8.2	−15.6	7.4	8.9	13.3	17.8
	6.0	12.0	−18.0	6.0	11.5	−17.5	0.5	0.5	−1.0	5.6	9.8	−15.4	6.7	18.3	14.4
	−18.0	12.0	6.0	−18.9	11.8	7.1	−1.0	0.5	0.5	−18.4	11.4	7.0	2.2	5.0	16.7
10−5	9.0	−18.0	9.0	8.9	−17.5	8.6	0.5	−1.0	0.5	8.5	−16.9	8.4	5.6	6.1	6.7
	6.0	12.0	−18.0	6.0	11.5	−17.5	0.5	0.5	−1.0	6.1	11.0	−17.1	1.7	8.3	5.0

Table 2. Identification results of the observation drift elements f.

$\mathit{\delta }$	f			${\widehat{\mathit{f}}}_0$			$\widehat{\mathit{f}}$			Relative Errors (%)
${10}^{-2}$	−10.0	5.0	20.0	−1.0	0.0	1.0	−9.9	5.1	20.1	1.0	2.0	0.5
$2\times {10}^{-4}$	−10.0	5.0	20.0	−1.0	0.0	1.0	−10.1	5.0	20.3	1.0	0.0	1.5
${10}^{-4}$	−10.0	5.0	20.0	−1.0	0.0	1.0	−10.1	5.0	20.2	1.0	0.0	1.0
${10}^{-5}$	−10.0	5.0	20.0	−1.0	0.0	1.0	−10.0	4.9	20.0	0.0	2.0	0.0

Table 2. Identification results of the observation drift elements f.

$\mathit{\delta }$	f			${\widehat{\mathit{f}}}_0$			$\widehat{\mathit{f}}$			Relative Errors (%)
${10}^{-2}$	−10.0	5.0	20.0	−1.0	0.0	1.0	−9.9	5.1	20.1	1.0	2.0	0.5
$2\times {10}^{-4}$	−10.0	5.0	20.0	−1.0	0.0	1.0	−10.1	5.0	20.3	1.0	0.0	1.5
${10}^{-4}$	−10.0	5.0	20.0	−1.0	0.0	1.0	−10.1	5.0	20.2	1.0	0.0	1.0
${10}^{-5}$	−10.0	5.0	20.0	−1.0	0.0	1.0	−10.0	4.9	20.0	0.0	2.0	0.0

Table 3. Identification results of the observation diffusion elements g.

$\mathit{\delta }$	g			${\widehat{\mathit{g}}}_0$			$\widehat{\mathit{g}}$			Relative Errors (%)
${10}^{-2}$	0.1	0.2	0.3	0.05	0.15	0.4	0.1	0.2	0.3	0.0	0.0	0.0
2 × 10−4	0.1	0.2	0.3	0.05	0.15	0.4	0.1	0.2	0.3	0.0	0.0	0.0
${10}^{-4}$	0.1	0.2	0.3	0.05	0.15	0.4	0.1	0.2	0.3	0.0	0.0	0.0
${10}^{-5}$	0.1	0.2	0.3	0.05	0.15	0.4	0.1	0.2	0.3	0.0	0.0	0.0

Table 3. Identification results of the observation diffusion elements g.

$\mathit{\delta }$	g			${\widehat{\mathit{g}}}_0$			$\widehat{\mathit{g}}$			Relative Errors (%)
${10}^{-2}$	0.1	0.2	0.3	0.05	0.15	0.4	0.1	0.2	0.3	0.0	0.0	0.0
2 × 10−4	0.1	0.2	0.3	0.05	0.15	0.4	0.1	0.2	0.3	0.0	0.0	0.0
${10}^{-4}$	0.1	0.2	0.3	0.05	0.15	0.4	0.1	0.2	0.3	0.0	0.0	0.0
${10}^{-5}$	0.1	0.2	0.3	0.05	0.15	0.4	0.1	0.2	0.3	0.0	0.0	0.0

Additionally, we compare the TIM estimate $\widehat{\Lambda }$ with an estimate ${\tilde{\Lambda }}_X$, calculated by (18) given the MJP state $X_t$. Clearly, $X_t$ is unobservable directly; hence, ${\tilde{\Lambda }}_X$ can be treated as some “ideal estimate”.

Below, we present the evolution of the HMM parameter estimates depending on the iteration number. The results are obtained for the discretization step $\delta ={10}^{-4}$. Figure 4 contains the estimates of the off-diagonal TIM $\Lambda$ coefficients in comparison with their exact values, Figure 5 contains the analogous information concerning the drift coefficients f, and Figure 6 presents the identification results of the diffusion coefficients g.

The identification quality can be illustrated in another way. One can build the histogram using the time-discretized observations as a sample. If the time step $\delta$ is small enough, then the theoretical pdf can be approximated by a mixture of the Gaussians ${\sum }_{n=1}^N{\pi }_{\infty }^n\mathcal{N}(y,\delta f{e}_n,\delta g_n)$, where the vector ${\pi }_{\infty }=\mathrm{col}({\pi }_{\infty }^1,\dots ,{\pi }_{\infty }^N)$ is the stationary distribution of the MJP $X_t$ with the TIM $\Lambda$.

In this case, the approximation precision is $o\left(\delta \right)$. We cannot construct the “theoretical pdf” because the true values of HMM parameters are unknown. Nevertheless, we can build its approximation ${\sum }_{n=1}^N{\widehat{\pi }}_{\infty }^n\mathcal{N}(y,\delta \widehat{f}{e}_n,\delta {\widehat{g}}_n)$, where $\widehat{f}$ and $\left\{{\widehat{g}}_n\right\}$ are identified values of the drift and diffusion and ${\widehat{\pi }}_{\infty }=\mathrm{col}({\widehat{\pi }}_{\infty }^1,\dots ,{\widehat{\pi }}_{\infty }^N)$ is the stationary distribution of a MJP with the TIM $\widehat{\Lambda }$. Figure 7 contains all the pdf estimates mentioned above, and the calculation is performed for the time increment $\delta ={10}^{-2}$.

The results of the second stage of the numerical experiment allow us to make the following conclusions.

• All estimates demonstrate high precision. Even for the relatively large time increment $\delta ={10}^{-2}$ and the corresponding small observation amount $\frac{T}{\delta }={10}^4$, the obtained estimates are precise.
• The iterative identification procedure converges quickly. In all numerical experiments, the iteration number does not exceed 15, and the identification procedure is terminated by the decline of the minimal relative parameter variation per iteration less than the threshold ${\epsilon }_0$.
• The precision of the parameter estimates depends on the time increment $\delta$. The most impact of the $\delta$ declining to the identification performance can be seen in the estimates of the drift f and diffusion g parameters.
• The estimation of the transition intensity matrix $\Lambda$ was found to be a challenging problem. Its identification precision demonstrates low dependency on the time increment $\delta$ value. Under a fixed observation interval $[0,t]$, the estimate $\Lambda$ cannot be precise more than the “ideal estimate” ${\tilde{\Lambda }}_X$ calculated by the precise knowledge of the state trajectory.

7. Conclusions

The paper presented a complete technological circuit to identify continuous-time HMM parameters. We included the identification problem formulation, which is proper from the mathematical point of view. The article contains a theoretical foundation of the solution and a numerical algorithm of the HMM parameter identification. We tested the proposed identification procedure over the extended complex of the numerical experiments. The numerical study confirmed the applicability of the current version of the identification algorithm and the high precision of the calculated estimates.

Nevertheless, we consider some directions for the further refinement of the algorithm and its theoretical basis. First, the authors of [29,30] have proven the convergence of the identification procedure to a local maximum of the conditional likelihood function. Let us remember that they considered the stochastic differential HMM with a scalar observable process and additive observation noises (the case $g_n\equiv g$). The keystone of the proof is the absolute continuity of the probability measures connected using a Girsanov transform. The continuity does not hold in the case of the continuous-time HMM with the observations corrupted by the multiplicative noise. Nevertheless, it is highly likely that this property will be provable for the time-discretized ones. Potentially, this would prove the convergence of the identification procedure and control the value of the conditional likelihood function during the iterations.

Second, the proposed algorithm has a significant drawback: it does not admit direct parallelization. This is an inherent feature of all identification algorithms based on fixed-interval smoothing, including the classic Baum–Welch algorithm [21,22]. Nevertheless, the proposed two-filter version of the fixed-interval smoothing admits double parallelization of the forward- and backward-time filtering procedures in each iteration step. Furthermore, the identification process could converge to a local maximum of the conditional likelihood function, depending on the starting point. The identification procedure admits the parallelization in the following sense: several exemplars of the identification procedure are performed on the independent processors given the same observation bulk and different starting points. Then, we can choose one of the calculated results, using either the likelihood function or some additional information concerning the estimated parameters.

Third, the current version of the algorithm can be modified to incorporate various constraints and describing the admissible set of the HMM parameters. We can illustrate this idea with the above-stated problem. It contains inequality (4), which postulates strong positivity for all off-diagonal elements of the TIM $\Lambda$. This constraint defines an open admissible set of $\Lambda$ values. Theoretically, this leads to the absence of the maximum likelihood TIM estimate for certain samples. Practically, this appears as the convergence of the estimates to a matrix with some zero off-diagonal elements. To avoid this situation, one should replace the initial open set of admissible parameter values with some closed subsets, using additional barriers. They can be either stationary or variable during the iterative process to preserve the calculated estimates from leaving the admissible set.