Sequential Monte Carlo methods for filtering of unobservable components of multidimensional diffusion Markov processes

Abstract

The problem of filtering of unobservable components x(t) of a multidimensional continuous diffusion Markov process $z\left(t\right)=\left(x\left(t\right),y\left(t\right)\right)$, given the observations of the (multidimensional) process y(t) taken at discrete consecutive times with small time steps, is analytically investigated. On the base of that investigation the new algorithms for simulation of unobservable components, x(t), and the new algorithms of nonlinear filtering with the use of sequential Monte Carlo methods, or particle filters, are developed and suggested. The analytical investigation of observed quadratic variations is also developed. The new closed-form analytical formulae are obtained, which characterize dispersions of deviations of the observed quadratic variations and the accuracy of some estimates for x(t). As an illustrative example, estimation of volatility (for the problems of financial mathematics) is considered. The obtained new algorithms extend the range of applications of sequential Monte Carlo methods, or particle filters, beyond the hidden Markov models and improve their performance.

Keywords:

• nonlinear filtering
• multidimensional diffusion Markov process
• particle filters
• sequential Monte Carlo methods
• simulation
• quadratic variation
• volatility

MSC codes:

• 60G35
• 60J60
• 62M05
• 93E11

Public Interest Statement

The problems of filtering of an unobservable process of interest, x(t), or estimation of the value of some function f(x(t)) (or functional) based on observations of the other random process y(t) related to x(t) arise in many areas: signal detection and filtering in radio physical devices; design of control and intelligent systems; financial mathematics; physical chemistry. In the last two decades, a great deal of works has been devoted to development and investigation of Monte Carlo methods for filtering, or particle filtering, for hidden Markov models where x(t) is a Markov process in itself. In the present paper, the new algorithms for the solution of nonlinear filtering problems (with the use of Monte Carlo calculations) are derived and obtained in explicit and closed analytical form for the general case where (x(t),y(t)) represents a multidimensional continuous diffusion Markov process. The new, effective methods for solution of nonlinear filtering difficult problems are obtained.

1. Introduction

In the last two decades, a great deal of works has been devoted to the development and the investigation of particle filters, or sequential Monte Carlo algorithms, for filtering of an unobservable process $x\left(t\right){=}^{\text{def}}(x_1\left(t\right),\dots ,x_m\left(t\right))$ given observations of the other process $y\left(t\right){=}^{\text{def}}(y_{m+1}\left(t\right),\dots ,y_p\left(t\right))$, taken at discrete times $t_k,t_0<t_1<t_2<\dots <t_k<\dots ,$ with small time steps $t_k-t_{k-1}{=}^{\text{def}}\mathrm{\Delta }{t}_k$ (see, e.g. survey (Crisan & Doucet, 2002), collection (Doucet, Freitas, & Gordon, 2001), works (Carvalho, Del Moral, Monin, & Salut, 1997; Del Moral, 1998; Doucet, Godsill, & Andrieu, 2000), recent complete survey (Doucet & Johansen, 2011), survey (Del Moral & Doucet, 2014), and references wherein). The observations $y(t_k)$ are being obtained consecutively, and an estimate for $x\left(t_k\right)$ should be updated at each time $t_k$. It is assumed that the whole process $z\left(t\right){=}^{\text{def}}\left(x\left(t\right),y\left(t\right)\right){=}^{\text{def}}(z_1\left(t\right),\dots ,z_p\left(t\right))$ is a Markov process.

In order to solve the problem of filtering of x(t_n), given the sequence of observations$y\left(t\right){|}_0^n{=}^{\text{def}}\left(y\left(t_0\right),y\left(t_1\right),\dots ,y\left(t_n\right)\right),$ with the use of Monte Carlo methods, the samples of the random sequences $x^i\left(t\right){|}_0^n{=}^{\text{def}}(x^i\left(t_0\right),x^i\left(t_1\right),\dots ,x^i\left(t_n\right))$ (from the distributions of $x\left(t_k\right)$ when $y\left(t\right){|}_0^k$ and ${x(t)|}_0^{k-1}$ are given) should have been simulated numerically for i = 1, …, N.

In general case, first of all, the problem arises: how to obtain an explicit and compact analytical expression (exactly or approximately) for the conditional probability density $P\left(x\left(t_k\right)|x\left(t\right){{|}_0^{k-1},y\left(t\right)|}_0^k\right)$ in order to simulate samples of $x\left(t_k\right)$ when $y\left(t\right){|}_0^k$ and ${x(t)|}_0^{k-1}$ are given.

Suppose that such sample sequences $x^i{(t)|}_0^k$ are being simulated, for i = 1, …, N. The joint probability density for the given observations $y\left(t\right){|}_0^n$ and the sequence $x^i\left(t\right){|}_0^n$ could be computed:(1) $$P_i\left(n\right){=}^{\text{def}}P\left(x^i\left(t\right){{|}_0^n,y\left(t\right)|}_0^n\right)=P_0(x^i\left(t_0\right),y\left(t_0\right))\prod _{k=1}^{n}P(x^i(t_k),y(t_k)|x^i\left(t_{k-1}\right),y(t_{k-1})),$$(1)

where $P(x\left(t_k\right),y(t_k)|x\left(t_{k-1}\right),y\left(t_{k-1}\right))$ represents the transition probability density of the Markov process $(x\left(t\right),y\left(t\right))$, and $P_0\left(x_0,y_0\right)$ is the probability density for the initial joint distribution of the initial value x(t₀) and the initial observation y(t₀).

In a large number of works devoted to particle filters the various algorithms of resampling were introduced in order to obtain the most a posteriori probable samples. For the problem of searching for the maximum a posteriori probable sample, the following algorithm of resampling can be suggested: introduce the values $W_i(n){=}^{\text{def}}{P}_i(n)/\sum _{j=1}^N{P}_j(n)$, for$i=1,\dots ,N$. Then W_i(n) can be considered as a weight of the sample sequence $x^i\left(t\right){|}_0^n$ (as well as of the sample point xⁱ(t_n)), which characterizes the a posteriori probability (or the importance) of the sample sequence $x^i\left(t\right){|}_0^n$ in comparison with the other samples $x^j\left(t\right){|}_0^n$ given $y\left(t\right){|}_0^n$ . Due to the Markov property of the process$z(t)$, the recursive formulae for $P_i\left(n\right)$ and W_i(n) can be written in general form; we have

(2) $$P_i\left(n+1\right)=P_i\left(n\right)P\left(x^i\left(t_{n+1}\right),y\left(t_{n+1}\right)\text{|}{x}^i\left(t_n\right),y\left(t_n\right)\right).$$(2)

Then, it is possible to calculate recursively the new values P_i(n + 1) and W_i(n + 1), when the new measurement y(t_n+1) is obtained and the new sample point xⁱ(t_n+1) is augmented to the sequence $x^i\left(t\right){|}_0^n$, so that $x^i\left(t\right){|}_0^{n+1}=(x^i\left(t_0\right),\dots ,x^i\left(t_n\right),x^i\left(t_{n+1}\right))$. It is easy to see, that it is not necessary to keep in memory all the sequence $x^i\left(t\right){|}_0^n$, but only the last point $x^i\left(t_n\right).$Then the samples$x^i\left(t\right){|}_0^n$ with small weights W_i(n) could be deleted, but the “more important” samples $x^j\left(t\right){|}_0^n$ should be continued with a few “offsprings”, so that the number of all sample sequences under consideration is equal to N (exactly or approximately). The point $x^{j^{\ast }}(t_n)$ and the sequence $x^{j^{\ast }}\left(t\right){|}_0^n$ that have the maximum weight $W_{j^{\ast }}(n)$ correspond to the point x^j(t_n) and to the sample sequence that is maximally a posteriori probable given the observations $y\left(t\right){|}_0^n$, among all the considered sample points x^j(t_n) and sample sequences $x^j\left(t\right){|}_0^n$. Then the value of $x^{j^{\ast }}(t_n)$ can be taken as the sought estimate of $x\left(t_n\right).$ On the base of the Laws of Large Numbers, it could be proved that this Monte Carlo estimate converges to the maximum a posteriori probability estimate for x(t_n) if N increases (see also Miguez, Crisan, & Djuric, 2013).

In most of works the additional assumption is accepted that the process x(t) is a Markov process in itself, and that the conditional probability density $P(y(t_k)|y\left(t_{k-1}\right),x\left(t_k\right),x\left(t_{k-1}\right))$ (for the observable process $y\left(t\right)$) can be presented in explicit and simple analytical form. Such cases are often referred to as hidden Markov models. Then, in many works, the sample sequences $x^j\left(t\right){|}_0^n$ are being simulated as trajectories of the Markov process $x(t)$ in itself, since such a simulation can be easily done. The joint probability density, corresponding to the constructed sample $x^j\left(t\right){|}_0^n$ and to the observations $y\left(t\right){|}_0^n$, is equal to(3) $$P(x^j\left(t\right){|}_0^n,y\left(t\right){|}_0^n)=P_0\left(x^j\left(t_0\right),y\left(t_0\right)\right)\prod _{k=1}^{n}P(x^j(t_k)|x^j\left(t_{k-1}\right))P\left(y\left(t_k\right)|y\left(t_{k-1}\right),x^j\left(t_k\right),x^j\left(t_{k-1}\right)\right)=P_j\left(n\right).$$(3)

Then the weights W_j(n) (or the other similar weights) could be easily calculated, and the above procedure of sampling and resampling could be realized, in order to obtain the estimate for x(t_n). In some cases it would be better to simulate samples x^j(t_k) with the use of the conditional probability density $P\left(x\left(t_k\right)|x\left(t_{k-1}\right),y\left(t_k\right),y\left(t_{k-1}\right)\right),$ which includes the observed value y(t_k). But this could have required a large amount of computations if we do not have an explicit and compact analytical expression for that conditional probability density. Meantime, if the process x(t) is being simulated just as a Markov process in itself, the sample sequences $x^j\left(t\right){|}_0^n$ represent samples from the a priori distribution for $x\left(t_k\right)$that can be far apart from the a posteriori distribution for $x\left(t_k\right)$ given ${y(t)|}_0^k$. Then resampling with the use of weights and significant increase in the number N are needed in order to find the most a posteriori probable values of $x\left(t_n\right)$ using the algorithm described above. Nevertheless, for hidden Markov models the particle filters with such a simple simulation of the sample sequences $x^j\left(t\right){|}_0^n$, with resampling based on consideration of weights, but with large N, have been implemented and proved to be useful in some applications (Carvalho et al., 1997; Thrun, Fox, Burgard, & Dellaert, 2001).

In general case when x(t) represents only some components of a multidimensional Markov process $(x\left(t\right),y\left(t\right))$ (and x(t) is not a Markov process in itself), it was not shown in the literature how to simulate sample sequences $x^i\left(t\right){|}_0^n$ when the values of the other components, $y\left(t\right){|}_0^n$, are given (i.e. how to simulate sample sequences $x^i\left(t\right){|}_0^n$ without a formidable amount of computations at each time t_k).

Note that the difficulties in obtaining the samples $x^i{(t)|}_0^n$ that correspond to the distributions $P(x(t_k)|x(t){|}_0^{k-1},y\left({t)|}_0^k\right)$ have led to the introduction and the use of some “proposal sampling distributions” (which can be simulated easier) and “auxiliary particle filters” in many works (see, e.g. collection (Doucet, De Freitas, & Gordon, 2001) and survey (Doucet & Johansen, 2011)). But the “optimal proposal sampling distribution” would be equal to the true distribution $P(x(t_k)|x(t){|}_0^{k-1},y\left({t)|}_0^k\right)$ if the latter could have been found in closed form. In the present work the precise, explicit and compact analytical formulae and algorithms for simulation of the sample sequences $x^i\left(t\right){|}_0^n$, when $y\left(t\right){|}_0^n$ is given, and for recursive calculations of P_i(n), weights, and estimates are obtained for the general case of a multidimensional continuous diffusion Markov process $(x\left(t\right),y\left(t\right))$.

For the problem of estimation of a function $f\left(x\left(t_n\right)\right)$ or $\mathit{\phi }(x\left(t\right){|}_0^n)$, the new estimates in explicit and closed analytical form are obtained in the following Section 2.2.

Moreover, in the quotients $P_i(n)/\sum _{j=1}^N{P}_j(n)$ the “scale” of all the P_j(n) is cancelled, and the information that could characterize the smallness of the a posteriori probability of all the generated samples $x^j\left(t\right){|}_0^n$ is lost. In the present work, the new sequential Monte Carlo algorithms are derived that include some tests (developed in Section 2.4) in order to discard the samples of low a posteriori probability before the calculation of all the weights is done. The implementation of the suggested tests guarantees that the samples that remain under consideration belong to the domain where the a posteriori probability density is localized.

In some important cases of hidden Markov models, in the particle filters that were suggested and studied in the literature, where the samples $x^i{(t)|}_0^n$ (with i = 1, …, N) are being simulated as the trajectories of the Markov process x(t) in itself, it is possible that the large part of those N generated samples $x^i{(t)|}_0^n$ would not be localized in the “vicinity” of the “true” realization of the trajectory, $x^{tr}{(t)|}_0^n$, because they are simulated with the use of the a priori transition probability density for x(t), and they could be of low a posteriori probability given ${y(t)|}_0^n$. For example, if the “level of intensity of noise” in observations y(t) is small, the chance to generate randomly the sample ${x(t)|}_0^n$ that belongs to the domain where the a posteriori probability density (given ${y(t)|}_0^n$) is concentrated might be as small as the chance to catch randomly a needle in a haystack. Thus, although theoretically in those cases the Monte Carlo estimates converge as N → ∞, the required number N increases if the “level of intensity of noise” decreases. Meanwhile, a large amount of calculations would be wasted for processing of those samples $x^i{(t)|}_0^n$ of low a posteriori probability. In the present paper, in the following Section 2.6, we shall show (for the problem of nonlinear filtering of a signal) that with the use of the new algorithm (17) of simulation $x^i{(t)|}_0^n$ given ${y(t)|}_0^n$ and with some appropriate change of the mathematical model of description of the process of observations (which could be justified in applications), the new algorithm (17) will generate at once the samples $x^i{(t)|}_0^n$ that are localized in the domain that is a posteriori probable given ${y(t)|}_0^n$, for all i = 1, …, N.

The special case when the diffusion coefficients of the observable process y(t) depend on unobservable components x(t) is considered. In that case “observed quadratic variations” of the process y(t) can be introduced (when observations y(t_k) are taken in discrete time), and they contain a lot of information about x(t). In the present paper, the analytical formulae that characterize the observed quadratic variations are obtained in explicit and closed form. On the base of these results, the algorithms of filtering and estimation with the use of the observed quadratic variations are developed.

Theoretically, the derivative of the process of quadratic variation could be incorporated into the set of observed, known processes. If that derivative process were known, then with the use of the new algorithm of simulation (17), new estimates (26) and algorithms developed in the present paper, the filtering problem would be effectively solved, as shown in Section 2.5. But in practice we cannot assume that this derivative process is straightforwardly observed. However, if additional observations of y(t) are available, namely, y(t₀ + sδt), with δt ≪ Δt, $s=1,2,\dots ,$ then the estimates $\stackrel{\wedge }{u}(t_k)$ for the derivative u(t_k) can be obtained and used for the filtering, as shown in Section 2.5.

The similar systems arise if some additional precisely known observations are available, for example, measurements of some function $H(x\left(t\right),t)$. In that case the solution of the filtering problem, described in Section 2.5, could be implemented.

The obtained new algorithms improve performance of particle filters, or sequential Monte Carlo methods, and extend the range of their applications beyond the hidden Markov models.

The implementation of particle filters, as it is usual for Monte Carlo methods, requires a lot of repeating computations, but it became accessible and useful in some applications, since the speed of computations and the capacity of computers have increased dramatically in the last two decades. Note that the speed can be increased further with the use of parallel computing. The computational cost of implementation of a general scheme of particle filtering on some existing processing units that allow parallel programming was considered e.g. in the paper by Hendeby, Karlsson, and Gustafsson (2010).

2. Derivation of new recursive algorithms for Monte Carlo simulation and filtering of unobservable components of multidimensional diffusion Markov processes

2.1. Simulation of trajectories of unobservable components. Analytical investigation of a multidimensional diffusion Markov process observed at discrete times

Consider the multidimensional diffusion Markov process $z(t)\stackrel{def}{=}(x(t),y(t))\stackrel{def}{=}(z_1\left(t\right),\dots ,z_p\left(t\right))$. The components $(z_1\left(t\right),\dots ,z_m\left(t\right))\stackrel{def}{=}x(t)$ are unobservable, but the other components $(z_{m+1}\left(t\right),\dots ,z_p\left(t\right))\stackrel{def}{=}y(t)$ are available for observations, which are being taken at discrete times t_k.

A diffusion Markov process with continuous trajectories z(t) can be characterized by its drift and diffusion coefficients:(4) $$E\left\{\left(z_j\left(t+\mathrm{\Delta }t\right)-z_j\left(t\right)\right)|z\left(t\right)=z\right\}=A_j\left(z,t\right)\mathrm{\Delta }t+o(\mathrm{\Delta }t),i,j=1,\dots ,p,$$(4) (5) $$E\left\{((z_i\left(t+\mathrm{\Delta }t\right)-z_i\left(t\right))\left(z_j\left(t+\mathrm{\Delta }t\right)-z_j\left(t\right)\right)|z\left(t\right)=z\right\}=b_{ij}\left(z,t\right)\mathrm{\Delta }t+o(\mathrm{\Delta }t),$$(5)

where Δt is a small time step. Denote $\mathrm{\Delta }z\left(t\right)\stackrel{def}{=}z\left(t+\mathrm{\Delta }t\right)-z\left(t\right)\stackrel{def}{=}\left(\mathrm{\Delta }{z}_1\left(t\right),\dots ,\mathrm{\Delta }{z}_p\left(t\right)\right).$

From the assumption that the trajectories of the process z(t) are continuous functions of t it follows (Kolmogorov, 1931/1986, 1933/1986) that(6) $$\underset{\mathrm{\Delta }t\to 0}{lim}\frac{1}{\mathrm{\Delta }t}E\left\{\mathrm{\Delta }{z}_{i_1}\left(t\right)\dots \mathrm{\Delta }{z}_{i_r}\left(t\right)|z\left(t\right)=z\right\}=0,withr>2.$$(6)

The matrix ∥ b_ij ∥  is symmetric and nonnegative definite. Therefore, in general case, that matrix can be represented with the use of its “square root” in the form ∥ b_ij ∥ = ∥ a_ij ∥ ∥ a_ij ∥ ^T, where ∥ a_ij ∥ ^T stands for transpose matrix ∥ a_ij ∥ .

Then the process $z\left(t\right)$ could be constructed as a solution of the system of stochastic differential equations:(7) $$d{z}_i\left(t\right)=A_i\left(z,t\right)dt+\sum _{j=1}^p{a}_{ij}\left(z,t\right)d{w}_j(t),i=1,\cdots ,p,t>t_0,$$(7)

where $w_j\left(t\right)$ are independent Wiener processes. The initial condition at t = t₀ is given as $z\left(t_0\right)=z_0$, where z₀ is a random variable independent of all $w_j\left(t\right),$with given probability density P₀(z₀). The system (7) should be interpreted as the system of stochastic integral equations:(8) $$z_i\left(t\right)=\underset{t_0}{\overset{t}{\int }}{A}_i\left(z\left(s\right),s\right)ds+\sum _{j=1}^p\underset{t_0}{\overset{t}{\int }}{a}_{ij}\left(z\left(s\right),s\right)d{w}_{j(}s),i=1,\dots ,p,$$(8)

with stochastic integrals in the sense of Ito (Doob, 1953). It is assumed that the drift coefficients A_i(z, t) satisfy the Lipschitz condition(9) $$\Vert A_i\left(z^{\left(1\right)},t\right)-A_i(z^{\left(2\right)},t)\Vert \le K\Vert z^{(1)}-z^{(2)}\Vert ,withK=Const.$$(9)

It is also assumed that the diffusion coefficients b_ij(z, t) are continuous and differentiable functions of z and t. The Lipschitz conditions (9) guarantee that trajectories z(t) do not have finite escapes, i.e. $\left|z(t)\right|$ does not tend to infinity when t tends to some finite moments of time.

The system of integral equations (8) with any given continuous trajectory $w\left(t\right)\stackrel{def}{=}(w_1\left(t\right),\dots ,w_p\left(t\right))$ can be solved with the use of successive approximations, which converge and define the continuous trajectory $z\left(t\right)$. Thus, the trajectories of the process z(t) are continuous with the probability 1.

The diffusion Markov process z(t) can be also constructed as the limit solution of the system of finite difference equations:(10) $$z_i\left(t_k\right)-z_i\left(t_{k-1}\right)=A_i\left(z\left(t_{k-1}\right),t_{k-1}\right)\mathrm{\Delta }{t}_k+\sum _{j=1}^p{a}_{ij}(z(t_{k-1}),t_{k-1}){\mathit{\eta }}_j(t_k)\sqrt{\mathrm{\Delta }{t}_k},$$(10)

where the random impacts ${\mathit{\eta }}_j\left(t_k\right)$ are independent random variables (for all $j=1,\dots ,p$, $k=1,2,\dots$), with $E{\mathit{\eta }}_j\left(t_k\right)=0$ and $E{\left({\mathit{\eta }}_j(t_k)\right)}^2=1$; $\mathrm{\Delta }{t}_k\stackrel{def}{=}{t}_k-t_{k-1}$; $i=1,\dots ,p$. In particular, the increments of Wiener processes, $\mathrm{\Delta }{w}_j\left(t_k\right)\stackrel{def}{=}{w}_j\left(t_k\right)-w_j\left(t_{k-1}\right),$ can be used in the finite difference equations (10(10) $$z_i\left(t_k\right)-z_i\left(t_{k-1}\right)=A_i\left(z\left(t_{k-1}\right),t_{k-1}\right)\mathrm{\Delta }{t}_k+\sum _{j=1}^p{a}_{ij}(z(t_{k-1}),t_{k-1}){\mathit{\eta }}_j(t_k)\sqrt{\mathrm{\Delta }{t}_k},$$(10) ) instead of ${\mathit{\eta }}_j(t_k)\sqrt{\mathrm{\Delta }{t}_k}$ :(11) $$z_i\left(t_k\right)-z_i\left(t_{k-1}\right)=A_i\left(z\left(t_{k-1}\right),t_{k-1}\right)\mathrm{\Delta }{t}_k+\sum _{j=1}^p{a}_{ij}(z(t_{k-1}),t_{k-1})\mathrm{\Delta }{w}_j(t_k),$$(11)

Here again it is assumed that the Lipschitz conditions (9) are satisfied, and that the functions A_i(z, t) and b_ij(z, t) are differentiable with respect to z. The construction of diffusion Markov processes by the passage to the limit in the scheme of finite difference equations when time steps Δt_k tend to zero was first introduced and investigated by Academician S.N. Bernstein in 1934, 1938 years. The works (Bernstein, 1934/1964, 1938/1964) were published at first in French. They are republished in Collected works of S.N. Bernstein, vol.4 (in Russian), Publishing House Nauka (Academy of Science of USSR, Moscow, 1964). The results, obtained in these works (Bernstein, 1934/1964, 1938/1964), define analytically all the joint probability distributions for the limit process z(t), for all the random variables $z(t_{i_1})$, $z(t_{i_2})$,…, $z\left(t_{i_r}\right),$ for any r and $t_{i_r}$. By the Theorem of A.N. Kolmogorov on extension (or continuation) of measures in function space (Kolmogorov, 1933/1950), those multidimensional distributions can be continued to the measure on the σ-algebra that corresponds to the process z(t) with continuous time t. The convergence of the linear broken lines, which represent interpolation of the solutions of finite difference stochastic equations (10) or (11), to the continuous trajectories of diffusion Markov processes was also proved in (Kushner, 1971). Thus, the diffusion Markov process $z\left(t\right)$ is well defined by the system of finite difference equations (which could be more general, but similar to (10)) and by the passage to the limit from those Markov processes with discrete time (Bernstein, 1934/1964, 1938/1964). Analytical methods and solutions for some problems of filtering and estimation, based on recursive finite difference stochastic equations, with the passage to the limit (as the time steps tend to zero), similarly to the scheme of S.N. Bernstein, were developed and investigated in (Khazen, 1968, 1971, 1977; Chapter 3; Khazen, 2009)

In the following Section 2.2, the goal is to obtain the estimates of $x\left(t_n\right)$ or $f(x\left(t_n\right))$ with the use of Monte Carlo calculation of some integrals that can be interpreted as mathematical expectations. We should notice the following properties that are important in order to achieve that goal. Denote $P_{\mathrm{\Delta }}\left(z\left(t_{i_1}\right),\dots ,z\left(t_{i_r}\right)\right)$ the joint probability density for the random process z_Δ(t) that represents the solution of the finite difference equations (10) or (11) with small time steps Δt_k = Δ (and with piecewise constant or piecewise linear interpolation on the small time intervals $(t_0+k\mathrm{\Delta },t_0+\left(k+1\right)\mathrm{\Delta })$, $k=0,1,2,\dots$). Denote $P\left(z\left(t_{i_1}\right),\dots ,z\left(t_{i_r}\right)\right)$ the joint probability density for the limit diffusion Markov continuous process z(t) obtained as Δ → 0. It was established in the works by Bernstein (1934/1964, 1938/1964) that $P_{\mathrm{\Delta }}(z\left(t_{i_1}\right),\dots ,z\left(t_{i_r}\right))\to P(z\left(t_{i_1}\right),\dots ,z\left(t_{i_r}\right))$as Δ → 0. Hence, the contribution to the error of the Monte Carlo estimates of the integrals $\int f(x(t_n))P(x(t){|}_0^n,y(t){|}_0^n{)dx(t)|}_0^n$ and $\int P(x(t){|}_0^n,y(t){|}_0^{n}dx(t){|}_0^n$ (which are considered below in (23), (24)), caused by approximation of the limit diffusion process z(t) by its pre-limit finite difference model (11), tends to zero as Δ → 0. For the scheme (11), it was also proved in (Kushner, 1971, Chapter 10) that trajectories z_Δ(t) converge in mean square to the trajectories of the limit diffusion process z(t) as Δ → 0, so that E$\left\{{\left[z\left(t\right)-z_{\mathrm{\Delta }}(t)\right]}^2\right\}\to 0$ . Consequently, $\left|E\left\{f(z\left(t\right))\right\}-E\left\{f(z_{\mathrm{\Delta }}\left(t\right))\right\}\right|\to 0$, if the function f(z) satisfies Lipschitz condition. Meanwhile, in the problems of filtering the best estimate of $f(x\left(t_n\right))$ (by the criterion of minimum mean square error) is $\widehat{f(x(t_n)}$ $=E\left\{f(x\left(t_n\right)){|y(t)|}_0^n\right\}$, and, in general case, the value of the mean square error of filtering, $E\left\{{\left[f\left(x(t_n)\right)-\widehat{f(x(t_n))}\right]}^2\right\}$, remains finite, positive, limited from below, even if the error of calculation of that conditional expectation tends to zero. Thus, the accuracy of the estimate of x(t_n) or $f\left(x\left(t_n\right)\right)$ could not be noticeably improved even if the accuracy of approximation of the limit diffusion process z(t) by its finite difference model (11) increased as Δ decreased. Therefore, it is justified to use the finite difference model (11) for description and simulation of the considered process z(t), in order to obtain the estimates of x(t_n) or $f\left(x\left(t_n\right)\right)$ that provide solutions to the considered filtering problems (with feasible precision) in case if the observations y(t_k) are being taken in discrete consecutive times t_k with small time steps.

` In the present work, we are interested to describe analytically the conditional probability density $P(x(t_k)|x\left(t\right){|}_0^{k-1},y\left(t\right){|}_0^k)$. Consider the local increment Δz of the Markov process $z\left(t\right)$ when the value $z\left(t\right)=z$ is given, with small time step Δt; denote$\mathrm{\Delta }z\stackrel{def}{=}z\left(t+\mathrm{\Delta }t\right)-z\left(t\right)\stackrel{def}{=}(\mathrm{\Delta }x,\mathrm{\Delta }y)$. We have the relations (4)–(6) for the moments of Δz given $z\left(t\right)=z.$ Then the characteristic function for the random value Δz can be presented in the form:

(12) $$F\left(u|z,t\right)\stackrel{def}{=}E\left\{exp\left(i\sum _{k=1}^p{u}_k\mathrm{\Delta }{z}_k\right)|z\left(t\right)=z\right\}\stackrel{def}{=}\int exp\left(i\sum _{k=1}^p{u}_k\mathrm{\Delta }{z}_k\right)P\left(\mathrm{\Delta }z\text{|}z,t\right)d\mathrm{\Delta }z=exp\left\{\mathrm{\Delta }t\sum _{k=1}^p{A}_k\left(z,t\right)\left(i{u}_k\right)+\frac{1}{2}\mathrm{\Delta }t\sum _{k,l=1}^p{b}_{kl}\left(z,t\right)\left(i{u}_k\right)(i{u}_l)+o(\mathrm{\Delta }t)\right\},$$(12)

where $u=\left(u_1,\dots ,u_p\right)$ is a real vector. The last equality follows from the known expression for the characteristic function of a random variable with the use of its moments (see, for example, course of the Theory of Probability (Gnedenko, 1976). The inverse transformation provides the representation(13) $$P\left(\mathrm{\Delta }z|z,t\right)\stackrel{def}{=}{(2\mathit{\pi })}^{-p}\int exp\left(-i\sum _{k=1}^p{u}_k\mathrm{\Delta }{z}_k\right)F\left(u|z,t\right)du={(2\mathit{\pi }\mathrm{\Delta }t)}^{-\frac{p}{2}}{(Det\Vert b_{ij}\Vert )}^{-1/2}exp\left\{-\frac{1}{2\mathrm{\Delta }t}\sum _{i,j=1}^p{r}_{ij}(z,t)(\mathrm{\Delta }{z}_i-A_i\left(z,t\right)\mathrm{\Delta }t)(\mathrm{\Delta }{z}_j-A_j\left(z,t\right)\mathrm{\Delta }t)\right\},$$(13)

where $\Vert r_{ij}\Vert \stackrel{def}{=}\Vert b_{ij}{\Vert }^{-1}$ is the inverse matrix of ∥ b_ij ∥  (with 1$\le i,j\le p).$ The above expressions (12), (13) show that the local increment Δz (in the small, but finite time interval Δt) of the multidimensional diffusion Markov process z(t) with continuous trajectories, when $z\left(t\right)=z$ is given, can be considered as a multidimensional Gaussian random variable.

Using the Theorem on Normal Correlation, we find the following expressions for the first and second moments of the conditional probability distribution of the increments Δx_α (with (α = 1, …, m) provided that the increments Δy_ρ (with ρ = (m + 1), …, p) and the value $z\left(t\right)=z$ are given (Khazen, 2009, Chapter 3, Sections 3.1.2 and 3.3.1, pages 79–81, 101–106):(14) $$E\left\{\mathrm{\Delta }{x}_{\mathit{\alpha }}|\mathrm{\Delta }y,z\right\}=A_{\mathit{\alpha }}\left(z,t\right)\mathrm{\Delta }t+b_{\mathit{\alpha }\mathit{\sigma }}\left(z,t\right)c_{\mathit{\sigma }\mathit{\rho }}(z,t)\left(\mathrm{\Delta }{y}_{\mathit{\rho }}-A_{\mathit{\rho }}\left(z,t\right)\mathrm{\Delta }t\right)+o(\mathrm{\Delta }t),$$(14)

(15) $$E\left\{\mathrm{\Delta }{x}_{\mathit{\alpha }}\mathrm{\Delta }{x}_{\mathit{\beta }}|\mathrm{\Delta }y,z\right\}=\left(b_{\mathit{\alpha }\mathit{\beta }}\left(z,t\right)-b_{a\mathit{\sigma }}\left(z,t\right)c_{\mathit{\sigma }\mathit{\rho }}\left(z,t\right)b_{\mathit{\rho }\mathit{\beta }}(z,t)\right)\mathrm{\Delta }t+o(\mathrm{\Delta }t).$$(15)

Hereinafter the summation is being made over repeated indices, with 1 ≤ α, β ≤ m, $\left(m+1\right)\le \mathit{\sigma },\mathit{\rho }\le p;$ the matrix ∥ c_σρ ∥  represents the inverse or pseudo-inverse (Moore–Penrose) matrix of the matrix of the diffusion coefficients of the observable components y_ρ(t), so that $\Vert c_{\mathit{\sigma }\mathit{\rho }}\Vert {=}^{\text{def}}\Vert b_{\mathit{\sigma }\mathit{\rho }}{\Vert }^{+}.$

For the probability density of the increments Δy, we obtain the following expression:(16) $$P\left(\mathrm{\Delta }y\text{|}z\left(t\right)=z\right)=C(z,t)exp\left\{-\frac{1}{2\mathrm{\Delta }t}{c}_{\mathit{\sigma }\mathit{\rho }}\left(z,t\right)\left(\mathrm{\Delta }{y}_{\mathit{\rho }}-A_{\mathit{\rho }}\left(z,t\right)\mathrm{\Delta }t\right)(\mathrm{\Delta }{y}_{\mathit{\sigma }}-A_{\mathit{\sigma }}\left(z,t\right)\mathrm{\Delta }t)\right\},$$(16)

where C(z, t) is the normalization factor.

Note, that in the theory of filtering of diffusion Markov processes with observations made in continuous time t, i.e. if $y\left(s\right)$ is supposed to be known exactly on the time interval t₀ ≤ s ≤ t, in such case it is assumed that the diffusion coefficients b_σρ, with $\left(m+1\right)\le \mathit{\sigma },\mathit{\rho }\le p$, do not depend on unobservable components $x\left(t\right).$ In the contrary case (as it was pointed out by this author (Khazen, 1977; Chapter 3; Khazen, 2009, Chapter 3), since the diffusion coefficients $b_{\mathit{\sigma }\mathit{\rho }}\left(x(t),y(t),t\right)$ could be (at least, theoretically) precisely restored on the base of a single realization of the observed trajectory y(s) on the small time interval t  − δ ≤ s ≤ t + δ (no matter how small the taken value of δ is), the functions b_σρ(x(t), y(t), t) could have been incorporated into the set of observable functions. In the problem at hand, when filtering of the process x(t) should be obtained on the base of observations $y\left(t_k\right),$taken at discrete times t_k with small, but finite time steps Δt_k, that restriction may be cancelled. We shall consider further (in Sections 2.4, 2.5) how to incorporate some estimates of the values of $b_{\mathit{\sigma }\mathit{\rho }}\left(x(t_k),y(t_k),t_k\right)$ given ${y(t)|}_0^k$ in order to improve filtering of $x\left(t_n\right).$

The matrix ∥ b_αβ − b_ασc_σρb_ρβ ∥  is symmetric and nonnegative definite, and it can be presented in the following form:

∥ b_αβ − b_ασc_σρb_ρβ ∥  ${=}^{\text{def}}$ $\Vert g_{\mathit{\alpha }\mathit{\beta }}\Vert \Vert g_{\mathit{\gamma }\mathit{\delta }}{\Vert }^T,$ with 1 ≤ α, β, γ, δ ≤ m,

Then the samples x^j(t_k), when $x^j{(t)|}_0^{k-1}$ and ${y(t)|}_0^k$ are given, can be simulated as(17) $$x_{\mathit{\alpha }}^j\left(t_k\right)=x_{\mathit{\alpha }}^j\left(t_{k-1}\right)+A_{\mathit{\alpha }}\left(x^j\left(t_{k-1}\right),y\left(t_{k-1}\right)\right)\mathrm{\Delta }{t}_k+g_{\mathit{\alpha }\mathit{\beta }}(x^j\left(t_{k-1}\right),y(t_{k-1})){\mathit{\xi }}_{\mathit{\beta }}^j(k)\sqrt{\mathrm{\Delta }{t}_k}+b_{\mathit{\alpha }\mathit{\sigma }}\left(x^j(t_{k-1)},y\left(t_{k-1}\right)\right)c_{\mathit{\sigma }\mathit{\rho }}(x^j\left(t_{k-1}\right),y(t_{k-1}))\left(y_{\mathit{\rho }}\left(t_k\right)-y_{\mathit{\rho }}(t_{k-1})-A_{\mathit{\rho }}\left(x^j\left(t_{k-1}\right),y\left(t_{k-1}\right)\right)\mathrm{\Delta }{t}_k\right),$$(17)

where ${\mathit{\xi }}_{\mathit{\beta }}^j(k)$ are independent samples from Gaussian standard distribution, with $E{\mathit{\xi }}_{\mathit{\beta }}^j\left(k\right)=0$, $E{({\mathit{\xi }}_{\mathit{\beta }}^j\left(k\right))}^2=1$ (for all β = 1, …, m, $j=1,\dots ,N$, $k=1,2,\dots )$; Δt_k = t_k - t_k-1.

We shall denote shortly $b_{\mathit{\alpha }\mathit{\sigma }}\left(x^j(t_{k-1)},y\left(t_{k-1}\right)\right){=}^{\text{def}}{b}_{\mathit{\alpha }\mathit{\sigma }}^j\left(t_{k-1}\right),c_{\mathit{\sigma }\mathit{\rho }}\left(x^j(t_{k-1)},y\left(t_{k-1}\right)\right){=}^{\text{def}}{c}_{\mathit{\sigma }\mathit{\rho }}^j\left(t_{k-1}\right),$ $A_{\mathit{\sigma }}\left(x^j(t_{k-1)},y\left(t_{k-1}\right)\right){=}^{\text{def}}{A}_{\mathit{\sigma }}^j\left(t_{k-1}\right),$ $c_{\mathit{\sigma }\mathit{\rho }}\left(x\left(t_{k-1}\right),y\left(t_{k-1}\right)\right){=}^{\text{def}}{c}_{\mathit{\sigma }\mathit{\rho }}\left(t_{k-1}\right).$

Note, that the above expressions (14)–(17) show that if all the diffusion coefficients b_αρ vanish (with 1 ≤ α ≤ m, $\left(m+1\right)\le \mathit{\rho }\le p$) and if $A_{\mathit{\alpha }}\left(z,t\right)$, b_αβ(z, t) do not depend on y(t), then ${x(t)|}_1^n$ will be simulated without the use of ${y(t)|}_1^n$. Only the initial measurement $y\left(t_0\right)$ will be taken into account: the initial sample point $x^j\left(t_0\right)$ should be simulated as a random variable with the probability density(18) $$P_0(x(t_0)|y\left(t_0\right)=y_0)=P_0(x\left(t_0\right),y_0)/{\int }^{}{P}_0(x(t_0),y_0)dx(t_0).$$(18)

This particular case do not correspond to all the hidden Markov models: The process x(t) could be a Markov process in itself, but it is still possible, that some coefficients b_αρ are not equal to zero, for example, in case, if there are some “white noises”, which enter simultaneously into equations for x(t) and into “random disturbances or random errors” of the observable process y(t).

In general case of a multidimensional continuous Markov diffusion process $(x\left(t\right),y\left(t\right))$, the influence of the observed data ${y(t)|}_1^n$ is manifested in the expressions (17(17) $$x_{\mathit{\alpha }}^j\left(t_k\right)=x_{\mathit{\alpha }}^j\left(t_{k-1}\right)+A_{\mathit{\alpha }}\left(x^j\left(t_{k-1}\right),y\left(t_{k-1}\right)\right)\mathrm{\Delta }{t}_k+g_{\mathit{\alpha }\mathit{\beta }}(x^j\left(t_{k-1}\right),y(t_{k-1})){\mathit{\xi }}_{\mathit{\beta }}^j(k)\sqrt{\mathrm{\Delta }{t}_k}+b_{\mathit{\alpha }\mathit{\sigma }}\left(x^j(t_{k-1)},y\left(t_{k-1}\right)\right)c_{\mathit{\sigma }\mathit{\rho }}(x^j\left(t_{k-1}\right),y(t_{k-1}))\left(y_{\mathit{\rho }}\left(t_k\right)-y_{\mathit{\rho }}(t_{k-1})-A_{\mathit{\rho }}\left(x^j\left(t_{k-1}\right),y\left(t_{k-1}\right)\right)\mathrm{\Delta }{t}_k\right),$$(17) ), which describe simulation of the random samples $x^j\left(t_k\right)$.

The values P_i(n) and W_i(n) can be easily calculated:

(19) $$P_i\left(n\right)=P_0\left(x^i\left(t_0\right),y\left(t_0\right)\right)\prod _{k=1}^{n}P(\mathrm{\Delta }{x}^i(t_k)\text{|}\mathrm{\Delta }y\left(t_k\right),x^i\left(t_{k-1}\right),y\left(t_{k-1}\right)\text{)}P\left(\mathrm{\Delta }y\left(t_k\right)\text{|}{x}^i\left(t_{k-1}\right),y\left(t_{k-1}\right)\right),$$(19)

where the conditional probability densities for the increments, Δxⁱ(t_k)${=}^{\text{def}}{x}^i\left(t_k\right)-x^i(t_{k-1})$, Δy(t_k) ${=}^{\text{def}}y\left(t_k\right)-y\left(t_{k-1}\right),$ are Gaussian densities with moments determined by (14), (15), and (16). And the above procedure of sampling and resampling, which corresponds to a search for the point x(t_n) that provides maximum to the a posteriori probability density, could be implemented. In the next Sections 2.2–2.4, the results of the further analytical investigation and development of the algorithms are presented.

2.2. Estimating of a function f(x(t_n))

Consider the problem of estimation of some functions $f\left(x\left(t_n\right)\right)$ or $\mathit{\phi }\left(x\left({(t)|}_0^n\right)\right)$ given $y\left(t\right){|}_0^n.$ The general expression for the estimate can be written in the following form:(20) $$\widehat{f(x\left(t_n\right))}=E\left\{f(x\left(t_n\right)){|y(t)|}_0^n\right\}=\frac{{\int }^{}f\left(x\left(t_n\right)\right)P\left(x\left(t\right){{|}_0^n,y\left(t\right)|}_0^n\right)dx\left(t\right){|}_0^n}{{\int }^{}{P(x(t)|}_0^n,y\left(t\right){|}_0^n{)dx(t)|}_0^n}.$$(20)

It is assumed that for the considered functions f(x) or $\mathit{\phi }(x\left(t\right){|}_0^n)$ this conditional expectation exists. In case of a Markov process $(x\left(t_n\right),y(t_n))$, we can write:(21) $$P\left(x\left(t\right){{|}_0^n,y\left(t\right)|}_0^n\right)=P_0\left(x\left(t_0\right)\text{|}y\left(t_0\right)\right)P_0\left(y\left(t_0\right)\right)\prod _{k=1}^{n}P(x(t_k)\text{|}y\left(t_k\right),y\left(t_{k-1}\right),x\left(t_{k-1}\right)\text{)}\times \prod _{s=1}^{n}P(y(t_s)|y\left(t_{s-1}\right),x\left(t_{s-1}\right)).$$(21)

In our case of a multidimensional continuous diffusion Markov process, we have derived the explicit analytical expressions (19), (16). Hence, in this case, we obtain(22) $$\prod _{s=1}^{n}P(y(t_s)\text{|}y\left(t_{s-1}\right),x\left(t_{s-1}\right)\text{)}=C\left(x\left(t\right){{|}_0^{n-1},y\left(t\right)|}_0^{n-1}\right)exp\left\{-\frac{1}{2}\sum _{s=1}^n\frac{1}{\mathrm{\Delta }{t}_s}{c}_{o\mathit{\rho }}(t_{s-1})[\mathrm{\Delta }{y}_{\mathit{\sigma }}\left(t_s\right)\mathrm{\Delta }{y}_{\mathit{\rho }}\left(t_s\right)-2\mathrm{\Delta }{y}_{\mathit{\sigma }}(t_s)A_{\mathit{\rho }}\left(t_{s-1}\right)\mathrm{\Delta }{t}_s+A_{\mathit{\sigma }}(t_{s-1})A_{\mathit{\rho }}(t_{s-1}){(\mathrm{\Delta }{t}_s)}^2]\right\},$$(22)

where $C(x\left(t\right){|}_0^{n-1},y(t){|}_0^{n-1}$) represents the known normalization factor. Denote shortly $C_i\left(n-1\right){=}^{\text{def}}C(x^i\left(t\right){|}_0^{n-1},y(t){|}_0^{n-1})$.

If N samples $x^i{(t)|}_0^n$, with $i=1,\dots ,N$ and N ≫ 1, have been independently taken from the distribution$$P_0\left(x\left(t_0\right)|y\left(t_0\right)\right)\prod _{k=1}^{n}P(x(t_k)|y\left(t_k\right),y\left(t_{k-1}\right),x\left(t_{k-1}\right)),$$

then, in accordance with basics of Monte Carlo methods, the following integrals can be interpreted as mathematical expectations. Then we obtain(23) $$\int f\left(x\left(t_n\right)\right){P(x(t)|}_0^n,y\left(t\right){|}_0^n)dx\left(t\right){|}_0^n\approx P_0(y\left(t_0\right)\frac{1}{N}\sum _{i=1}^{N}f(x^i\left(t_n\right))C_i(n-1)\times exp\left\{-\frac{1}{2}\sum _{s=1}^n\frac{1}{\mathrm{\Delta }{t}_s}{c}_{\mathit{\sigma }\mathit{\rho }}^i(t_{s-1})[\mathrm{\Delta }{y}_{\mathit{\sigma }}\left(t_s\right)\mathrm{\Delta }{y}_{\mathit{\rho }}\left(t_s\right)-2\mathrm{\Delta }{y}_{\mathit{\sigma }}(t_s)A_{\mathit{\rho }}^i\left(t_{s-1}\right)\mathrm{\Delta }{t}_s+A_{\mathit{\sigma }}^i(t_{s-1})A_{\mathit{\rho }}^i(t_{s-1}){(\mathrm{\Delta }{t}_s)}^2]\right\},$$(23)

and similarly(24) $${\int P(x(t)|}_0^n,y\left(t\right){|}_0^n)dx\left(t\right){|}_0^n\approx P_0(y\left(t_0\right)\frac{1}{N}\sum _{i=1}^N{C}_i(n-1)\times exp\left\{-\frac{1}{2}\sum _{s=1}^n\frac{1}{\mathrm{\Delta }{t}_s}{c}_{\mathit{\sigma }\mathit{\rho }}^i(t_{s-1})[\mathrm{\Delta }{y}_{\mathit{\sigma }}\left(t_s\right)\mathrm{\Delta }{y}_{\mathit{\rho }}\left(t_s\right)-2\mathrm{\Delta }{y}_{\mathit{\sigma }}(t_s)A_{\mathit{\rho }}^i\left(t_{s-1}\right)\mathrm{\Delta }{t}_s+A_{\mathit{\sigma }}^i(t_{s-1})A_{\mathit{\rho }}^i(t_{s-1}){(\mathrm{\Delta }{t}_s)}^2]\right\}.$$(24)

Due to the Laws of Large Numbers, the accuracy of those approximate estimates increases as N increases.

Denote shortly(25) $${\mathrm{{\rm Y}}}_i\left(n\right)=C_i(n-1)exp\left\{-\frac{1}{2}\sum _{s=1}^n\frac{1}{\mathrm{\Delta }{t}_s}{c}_{\mathit{\sigma }\mathit{\rho }}^i(t_{s-1})[\mathrm{\Delta }{y}_{\mathit{\sigma }}\left(t_s\right)\mathrm{\Delta }{y}_{\mathit{\rho }}\left(t_s\right)-2\mathrm{\Delta }{y}_{\mathit{\sigma }}(t_s)A_{\mathit{\rho }}^i\left(t_{s-1}\right)\mathrm{\Delta }{t}_s+A_{\mathit{\sigma }}^i(t_{s-1})A_{\mathit{\rho }}^i(t_{s-1}){(\mathrm{\Delta }{t}_s)}^2]\right\}.$$(25)

Then the sought estimate $\widehat{f(x\left(t_n\right))}$ can be written in the following form:(26) $$\widehat{f(x\left(t_n\right))}=\sum _{i=1}^{N}f(x^i(t_n))\stackrel{}{{\tilde{W}}_i}(n),\text{where}\stackrel{}{{\tilde{W}}_i}\left(n\right)=\frac{{\mathrm{{\rm Y}}}_i(n)}{{\sum }_{j=1}^N{\mathrm{{\rm Y}}}_j(n)}.$$(26)

In most of applications, we can assume that there exists the mathematical expectation of the random value $\left|f(x^i\left(t_n\right))\right|{\mathrm{{\rm Y}}}_i(n)$, where ${y(t)|}_0^n$ is fixed and $x^i{(t)|}_0^n$ represent the results of independent random simulations with the use of (17); i.e. $E\{\left|f\left(x^i\left(t_n\right)\right)\right|{\mathrm{{\rm Y}}}_i\left(n\right)\}<\infty$. This assumption is equivalent to the following: $E\left\{\left|f(x\left(t_n\right))\right|{|y(t)|}_0^n\right\}<\infty$. Denote ${\mathit{\zeta }}_i=f(x^i(t_n)){\mathrm{{\rm Y}}}_i\left(n\right).$ The random values ζ_i are independent, and they have one and the same probability distribution, and $E\left\{|{\mathit{\zeta }}_i|\right\}<\infty$. Then the Strong Law of Large Numbers (the Theorem of Kolmogorov) guarantees that the sums (i.e. the arithmetic means $\frac{1}{N}\sum _{i=1}^N{\mathit{\zeta }}_i$) in right-hand side of the formulae (23), (24) converge with probability one (i.e. almost surely) to the integrals in left-hand side if N → ∞.

Note that in case if the diffusion coefficients b_σρ (for $m<\mathit{\sigma },\mathit{\rho }\le p$) of the observable process $y\left(t\right)$ do not depend on unobservable process x(t), then $C_i\left(n-1\right)$ and $c_{\mathit{\sigma }\mathit{\rho }}^i\left(t_{s-1}\right)$ do not depend on i, and cofactors $C_i(n-1)exp\left\{-\frac{1}{2}\sum _{s=1}^n\frac{1}{\mathrm{\Delta }{t}_s}{c}_{\mathit{\sigma }\mathit{\rho }}^i(t_{s-1})\mathrm{\Delta }{y}_{\mathit{\sigma }}\left(t_s\right)\mathrm{\Delta }{y}_{\mathit{\rho }}\left(t_s\right)\right\}$ should be cancelled in the numerator and denominator of the formula (2(2) $$P_i\left(n+1\right)=P_i\left(n\right)P\left(x^i\left(t_{n+1}\right),y\left(t_{n+1}\right)\text{|}{x}^i\left(t_n\right),y\left(t_n\right)\right).$$(2) 6). In that case we can put:(27) $${\mathrm{{\rm Y}}}_i\left(n\right)=exp\left\{-\frac{1}{2}\sum _{s=1}^n{c}_{\mathit{\sigma }\mathit{\rho }}(t_{s-1})[-2A_{\mathit{\rho }}^i\left(t_{s-1}\right)\mathrm{\Delta }{y}_{\mathit{\sigma }}(t_s)+A_{\mathit{\sigma }}^i(t_{s-1})A_{\mathit{\rho }}^i(t_{s-1})\mathrm{\Delta }{t}_s]\right\}.$$(27)

Note that the calculations of the values Y_i(n) (25) or (27) can be implemented recurrently since they are determined by the values of the recursively accumulated sums.

In case of hidden Markov models where the unobservable process $x\left(t\right)$ is a Markov process in itself, and the observable process y(t) is described by a stochastic differential equation (SDE):(28) $$dy\left(t\right)=h\left(x\left(t\right),t\right)dt+\mathit{\sigma }dw\left(t\right),$$(28)

where h(x, t) is a given nonlinear function with respect to x, $\mathit{\sigma }=Const,$ and w(t) is a standard Wiener process, the above estimate $\widehat{f(x\left(t_n\right))}$ takes the form, which is similar to the estimates for $f\left(x\left(t_n\right)\right)$ that were earlier constructed and studied for this particular case in the literature. The novelty of the estimates (23)–(27) and the simulation (17), obtained in the present paper, is that the simulation of unobservable components and estimates $\widehat{f(x\left(t_n\right))}$ are obtained in explicit and closed analytical form for the general case of partially observable multidimensional continuous diffusion Markov processes, and the obtained Monte Carlo estimate (26) converges with probability one to the sought posterior expectation of $f(x\left(t_n\right))$ given ${y(t)|}_0^n$, if N → ∞.

Now we shall demonstrate that the Monte Carlo estimates (23)–(26), (27) for $f\left(x\left(t_n\right)\right)$ (obtained in the present paper) hold true also in the case if the sample sequences $x^i{(t)|}_0^n$ are being generated with the use of some branching sampling procedures.

Consider the following branching resampling. As it was already pointed out in the Introduction, it is purposeful to discard the samples $x^i{(t)|}_0^n$ that have negligibly small weights W_i(n) or $\stackrel{}{{\tilde{W}}_i}(n)$, in order to decrease the amount of all the computations.

Suppose that in the end of each time interval $(T_k,T_{k+1}]$ (with $k=0,1,2,\dots$; T₀ = t₀; T_k+1 = T_k + m_kΔt, and $m_k$and N_k are given numbers) each sample sequence $x^i{(t)|}_0^n$ that still existed at the time T_k+1 is being continued with N_k+1 “offsprings”. At the initial moment of time, T₀, there are N₀ independently taken initial sample points xⁱ(t₀), $i=1,\dots ,N_0$. When t_n = T_k+1 and $t_{n+1}=T_{k+1}+\mathrm{\Delta }t$, sample points $x^{\left(i,j\right)}(t_{n+1})$ are being augmented to the sequence $x^i{(t)|}_0^n,$ with $j=1,\dots ,N_{k+1}$, in order to construct the “offsprings”. The sample points $x^{\left(i,j\right)}(t_{n+1})$ are being independently taken from the distribution $P(x(t_{n+1})|x^i(t){|}_0^n,y\left(t\right){|}_0^{n+1})$. Then those sample sequences (with “offsprings”) will be continued until the next time of branching, T_k+2, except for some sample sequences that are discarded before T_k+2 since their weights become “too small” (for example, smaller than some chosen threshold $W_{cr}$); and so on. For the simplicity of discussion, we can renumber again all the current sample sequences (that still exist at the time t_n) as $x^i{(t)|}_0^n$. Consider at first the case when all the sample sequences are being continued without discarding. Then their number is growing, and at the time T_k+1 it is equal to $N_0\times N_1\times \dots \times N_k$ (with k > 0).

Consider the estimate$$\overline{f(x\left(t_n\right))}=E\left\{f\left(x\left(t_n\right)\right){|y\left(t\right)|}_0^n\right\},$$

under condition that the above branching sampling procedure is being implemented. Denote $(t_{r\left(n\right)},t_{(r+1)\left(n\right)}]$ the time interval (between the moments of branching) that contains t_n, so that $t_{r(n)}=T_{k(n)}<t_n\le T_{k\left(n\right)+1}=t_{\left(r+1\right)(n)}$.

From the Markov property of the random process $(x\left(t\right),y\left(t\right))$ and the factorization (21), (22), and from the “tower property” of conditional expectations:

$$E\left\{f\left(x\left(t_n\right)\right){|y\left(t\right)|}_0^n\right\}=E\left\{E\left\{f(x(t_n)|y\left(t\right){|}_0^n,x(t){|}_0^{r(n)}\right\}\right\}=E\left\{E\left\{E\left\{f(x(t_n)|y\left(t\right){|}_o^n,x(t){|}_0^{r\left(n\right)}\right\}{|x(t)|}_0^{(r-1)\left(n\right)}\right\}\right\}=\dots ,$$

it follows that the formulae (23)–(26), (27) hold true for the Monte Carlo estimate of integral $\overline{f(x\left(t_n\right))}$ in the case if the sample sequences $x^i\left(t\right){|}_0^n$ are generated with the use of the above branching procedure.

We can consider the Monte Carlo estimates for the integrals (23), (24) for each i-th “tree” of the branching sample sequences, which begins at the point $x^i\left(t_0\right)$, with $i=1,\dots ,N_0$. Denote these random values as ξ_i(n), with $i=1,\dots ,N_0$. Then the random values ξ_i(n) are independent of each other, and they have one and the same probability distribution, and $E\left\{\left|{\mathit{\xi }}_i\right|\right\}<\infty$. Then they obey the Strong Law of Large Numbers, so that $\frac{1}{N_0}\sum _{i=1}^{N_0}{\mathit{\xi }}_i$ converges with probability one to the sought integral as N₀ → ∞.

But many of the exponential weights $\stackrel{}{{\tilde{W}}_i}(n)$ are rapidly decreasing with increase in the number n of time steps. Therefore, it is possible to discard some highly a posteriori improbable sample sequences $x^i\left(t\right){|}_0^n$, which do not provide noticeable contribution into the estimates $\widehat{f(x\left(t_n\right))}$. The numbers m_k, N_k, and the threshold W_cr can be adjusted (in practical implementation for a particular class of applications) in order to decrease the amount of all the calculations and make them feasible, and at the same time to keep the current sample points xⁱ(t_n) in the domain where the a posteriori probability density is localized. In practical implementation, it can be also purposeful to begin branching not at the fixed moment of time T_k but at some current moment of time t_n when the number of all existed sample sequences (or sample points) becomes less than some given fraction of N.

It is worth noting that in the following Section 2.4 the tests (47), (48) are obtained that allow detection and rejection of samples xⁱ(t_n) of low a posteriori probability to be done without the use of the weights W_i(n) or $\stackrel{}{{\tilde{W}}_i}(n)$ (for $i=1,\dots ,N)$, i.e. independently for each sample sequence $x^i{(t)|}_0^n$. Thus, in the above algorithm we can use such tests instead of the comparison of the weight $\stackrel{}{{\tilde{W}}_i}(n)$ with the threshold W_cr. Then the algorithm described above can be effectively implemented with the use of parallel computing, since the simulation and continuation (the branching resampling) of each sample sequence $x^i{(t)|}_0^m$ can be done independently of the other samples $x^k{(t)|}_0^m$ with k ≠ i. The values Y_i(n) determined by (25) or (27) will be also calculated recurrently and independently for each sample sequence $x^i{(t)|}_0^n$. Thus, that large amount of calculations could be effectively performed with the use of parallel computing. And only on the other stage all the obtained accumulated values Y_i(n) (with $i=1,\dots ,N$) will be used for the calculation of the weights $\stackrel{}{{\tilde{W}}_i}(n)$, which are needed for the calculation of the sought estimate $\widehat{f(x\left(t_n\right))}$.

The new Monte Carlo algorithm of estimation of $f\left(x\left(t_n\right)\right)$ or $\mathit{\phi }(x\left({t)|}_0^n\right)$ (presented above) is derived straightforwardly from Bayes formula (2(2) $$P_i\left(n+1\right)=P_i\left(n\right)P\left(x^i\left(t_{n+1}\right),y\left(t_{n+1}\right)\text{|}{x}^i\left(t_n\right),y\left(t_n\right)\right).$$(2) 0). The estimates are constructed with the use of the samples xⁱ(t_n) or $x^i{(t)|}_0^n$ that are simulated by (17), and their weights $\widehat{f(x\left(t_n\right))}$ are defined by (23)–(26), (27).

For the particular case of hidden Markov models, another algorithm with random branching resampling for “particle filtering” of unobservable Markov “signal” x(t_n) was developed, that is the Sequential Importance Resampling (SIR) algorithm; it is studied in the works (Crisan & Doucet, 2002; Del Moral, 1998), for the processes with discrete time. In that SIR algorithm with random branching resampling, the a posteriori probability distribution for x(t_n) is approximated by some “cloud of random particles” xⁱ(t_n) with weights w_i(n) = 1/N, $i=1,\dots ,N.$ With the use of our new algorithm of simulation of unobservable components x(t_n) given ${y(t)|}_0^n$ (see Section 2.1, (17)) and the new closed form analytical expressions (14)–(17), (23)–(27), that SIR algorithm with random branching can be generalized for the general case of filtering of unobservable components x(t) of a multidimensional diffusion Markov process $(x\left(t\right),y\left(t\right))$, in the following way.

In the generalized algorithm that we are proposing the sample sequence $x^i{(t)|}_o^m$ should be chosen for continuation at the time of branching T_m = t_r(m) (and kept up to the next time of branching, T_m+1) with probability p_i(m) (which will be determined below) in each of N independent attempts to continue their ensemble, so that the expectation of the random number K_i(m) of its “offsprings” is equal to $E\left\{K_i\left(m\right){|x^1\left(t_{r(m)}\right),\dots ,x^{K\left(m\right)}\left(t_{r(m)}\right),y(t)|}_0^{r(m)}\right\}=p_i\left(m\right)N$, with i = 1, …, K(m), where K(m) denotes the number of all the sample points xⁱ(t_r(m)) that still existed at the time t_r(m). We can suggest that the times of branching T_m+1 = T_m + MΔt, with M ≥ 1, $m=0,1,2,\dots .$ (Note, that M = 1 in the standard SIR algorithms). The probability p_i(m) depends on all the sample points x¹(t_m), …, $x^{K\left(m\right)}(t_m)$ or sample sequences $x^i{(t)|}_0^m$. Such a procedure of random resampling is similar to the standard SIR algorithm, but we have to determine the probability p_i(m) for the general case of multidimensional diffusion Markov processes $(x\left(t\right),y\left(t\right))$.

For the general case, with the use of analytical expressions (16)–(27), we derive the following expression:$$p_i(m)=\frac{Q_i\left(m\right)exp\left\{-R_i(m)\right\}}{{\sum }_{j=1}^{K\left(m\right)}{Q}_j\left(m\right)exp\left\{-R_j(m)\right\}},$$

where

$$R_i(m)=\frac{1}{2}\sum _{s=0}^{M-1}\frac{1}{\mathrm{\Delta }t}{c}_{\mathit{\sigma }\mathit{\rho }}^i(t_{r(m-1)+s})[\mathrm{\Delta }{y}_{\mathit{\sigma }}\left(t_{r(m-1)+s+1}\right)\mathrm{\Delta }{y}_{\mathit{\rho }}\left(t_{r(m-1)+s+1}\right)-2\mathrm{\Delta }{y}_{\mathit{\sigma }}\left(t_{r(m-1)+s+1}\right)A_{\mathit{\rho }}^i\left(t_{r(m-1)+s}\right)\mathrm{\Delta }t+A_{\mathit{\sigma }}^i\left(t_{r(m-1)+s}\right)A_{\mathit{\rho }}^i\left(t_{r(m-1)+s}\right){(\mathrm{\Delta }t)}^2],$$

and Q_i(m) denotes the normalization factor of the probability density $Q_i\left(m\right)exp\left\{-R_i(m)\right\}$. Here the notations are the same as earlier in Sections 2.1, 2.2.

In the case if the diffusion coefficients b_σρ of the observed process y(t) do not depend on the unobservable process x(t), the above expression can be written in the following more concise form:$$p_i(m)=\frac{exp\left\{-V_i(m)\right\}}{{\sum }_{j=1}^{K\left(m\right)}exp\left\{-V_j(m)\right\}},$$

where$$V_i(m)=\frac{1}{2}\sum _{s=0}^{M-1}{c}_{\mathit{\sigma }\mathit{\rho }}(t_{r(m-1)+s})[-2\mathrm{\Delta }{y}_{\mathit{\sigma }}\left(t_{r(m-1)+s+1}\right)A_{\mathit{\rho }}^i\left(t_{r(m-1)+s}\right)+A_{\mathit{\sigma }}^i\left(t_{r(m-1)+s}\right)A_{\mathit{\rho }}^i\left(t_{r(m-1)+s}\right)\mathrm{\Delta }t].$$

In case of hidden Markov models (28), which can be considered as a particular case of a multidimensional diffusion Markov process (x(t), y(t)), from the above expression follows:$$p_i(m)=\frac{exp\left\{-{\mathrm{\Gamma }}_i(m)\right\}}{{\sum }_{j=1}^{K\left(m\right)}exp\left\{-{\mathrm{\Gamma }}_j(m)\right\}},$$

where$${\mathrm{\Gamma }}_i\left(m\right)=\frac{1}{2{\mathit{\sigma }}^2}\sum _{s=0}^{M-1}[-2h(x^i\left(t_{r(m-1)+s}\right))\mathrm{\Delta }y\left(t_{r(m-1)+s+1}\right)+h^2(x^i\left(t_{r(m-1)+s}\right))\mathrm{\Delta }t]$$

For the simplicity of notation, the above formula is written for the case of one-dimensional processes $x\left(t\right)$ and $y\left(t\right)$, which satisfy (28). The last expression for p_i(m) (for the case M = 1) is in agreement with the determination of the probabilities of continuation of the samples xⁱ(t_n) presented in (Crisan & Doucet, 2002; Del Moral, 1998) for the case of hidden Markov models.

With the use of the above procedure of random branching and with our algorithm (17) of simulation of $x\left(t_n\right)$ given ${y(t)|}_0^n$, we can generalize the SIR algorithm of particle filtering with random branching, which was developed for hidden Markov models, to the general case of multidimensional diffusion Markov processes $(x\left(t\right),y\left(t\right))$.

Thus, in the present paper a few different possible versions of algorithms of particle filters are provided for filtering of unobservable components x(t_n) given $y\left(t\right){|}_0^n$, which are justified theoretically. The new algorithm developed above, with branching sampling and recurrent calculation of the values (25) or (27), appears to be preferable for implementation with the use of parallel computing. The further practical implementation and comparison for various applications could be achieved in the future works.

2.3. Analytical investigation of observed quadratic variations

Consider the case when some diffusion coefficients b_σρ of the observable process y(t) depend on the unobservable process x(t). We assume throughout this paper that the diffusion coefficients $b_{ij}\left(z,t\right)$ are continuous and differentiable functions with respect to z. For simplicity of notation, consider at first the case when y(t) is one-dimensional. For example, consider the following model that plays an important role in financial mathematics (Hull, 2000):(29) $$dx=A\left(x,t\right)dt+{\mathit{\sigma }}_1\left(x,t\right)d{w}_1\left(t\right),dy=\left(\mathit{\mu }-\frac{x^2}{2}\right)dt+{\mathit{\sigma }}_{2}xd{w}_2\left(t\right).$$(29)

Here $w_1\left(t\right),w_2(t)$ are independent Wiener processes, μ and σ₂ are known constants, and the functions A(x, t) and σ₁(x, t) are continuous and differentiable with respect to x. In the famous Black–Scholes models (Hull, 2000), the observable process y(t) in (29) represents the natural logarithm of a stock price, S(t), so that $y(t){=}^{\text{def}}logS(t)$, while the process $\left|x(t)\right|$ corresponds to volatility, which is to be estimated given the path y(s), t₀ ≤ s ≤ t. If the process y(s) were available for observations continuously over the time interval t₀ ≤ s ≤ t, the value $\left|x(t)\right|$ could have been restored precisely, at least theoretically.

It is well known that a diffusion Markov process z(t) is not differentiable at any time t. It can be proved that its quadratic variations(30) $$\sum _{k=1}^M[z_i\left(t_k\right)-z_i\left(t_{k-1}\right)]\left[z_j\left(t_k\right)-z_j\left(t_{k-1}\right)\right]\underset{max\left\{\mathrm{\Delta }{t}_k\right\}\to 0}{\to }\underset{t-\mathit{\delta }}{\overset{t+\mathit{\delta }}{\int }}{b}_{ij}\left(z\left(s\right),s\right)ds,$$(30)

if time steps Δt_k tend to zero, while M → ∞; here $\left(t-\mathit{\delta }\right)=t_0<t_1<\dots <t_k<\dots <t_M=(t+\mathit{\delta })$, with given δ > 0, and the value of δ can be arbitrarily small. Note that the relation (30(30) $$\sum _{k=1}^M[z_i\left(t_k\right)-z_i\left(t_{k-1}\right)]\left[z_j\left(t_k\right)-z_j\left(t_{k-1}\right)\right]\underset{max\left\{\mathrm{\Delta }{t}_k\right\}\to 0}{\to }\underset{t-\mathit{\delta }}{\overset{t+\mathit{\delta }}{\int }}{b}_{ij}\left(z\left(s\right),s\right)ds,$$(30) ) will be also obtained as a consequence of the new analytical estimates that are developed below.

For the system (29), $b_{22}\left(x(t),y(t),t\right)={({\mathit{\sigma }}_2)}^2{x}^2(t)$, and from (30) it follows that the value of x²(t) could be restored with any desirable accuracy. Then $\left|x(t)\right|=+\sqrt{x^2(t)}$, and the filtering problem would be solved exactly.

But the observations $y\left(t_k\right)$ are being taken at discrete times t_k with small, but finite time steps $\mathrm{\Delta }{t}_k$, so that it is impossible to observe precisely the values of $b_{\mathit{\sigma }\mathit{\rho }}(x\left(t_k\right),y\left(t_k),t_k\right).$ In general case, we can calculate the “observed quadratic variations”,(31) $$V_{\mathit{\sigma }\mathit{\rho }}^{\left(obs\right)}(t_k){=}^{\text{def}}\sum _{s=0}^M[y_{\mathit{\sigma }}\left(t_{k-s}\right)-y_{\mathit{\sigma }}\left(t_{k-s-1}\right)]\left[y_{\mathit{\rho }}\left(t_{k-s}\right)-y_{\mathit{\rho }}\left(t_{k-s-1}\right)\right],$$(31)

and consider it as an estimate for $b_{\mathit{\sigma }\mathit{\rho }}(x\left(t_k\right),y\left(t_k),t_k\right).$ The similar estimates for volatilities were introduced and considered in (Hull, 2000, chapter 15). It is more convenient to use the recursive averaging instead of the moving averaging (31). Then(32) $$V_{\mathit{\sigma }\mathit{\rho }}^{\left(obs\right)}\left(t_k\right){=}^{\text{def}}{e}^{-\mathit{\lambda }\mathrm{\Delta }{t}_k}{V}_{\mathit{\sigma }\mathit{\rho }}^{\left(obs\right)}\left(t_{k-1}\right)+\left[y_{\mathit{\sigma }}\left(t_k\right)-y_{\mathit{\sigma }}\left(t_{k-1}\right)\right]\left[y_{\mathit{\rho }}\left(t_k\right)-y_{\mathit{\rho }}\left(t_{k-1}\right)\right],$$(32)

with $\mathit{\lambda }>0.$ In the limit, if all Δt_k decreased and tended to zero, we would have obtained

(33) $$V_{\mathit{\sigma }\mathit{\rho }}^{\left(obs\right)}\left(t\right)\underset{max\left\{\mathrm{\Delta }{t}_k\right\}\to 0}{\to }\underset{t_0}{\overset{t}{\int }}{e}^{-\mathit{\lambda }\left(t-s\right)}{b}_{\mathit{\sigma }\mathit{\rho }}(x\left(s\right),y\left(s\right),s)ds.$$(33)

We can incorporate the observed quadratic variations $V_{\mathit{\sigma }\mathit{\rho }}^{\left(obs\right)}(t_k)$ into the set of observed data. We shall show further, in Section 2.4, how to use that “observed quadratic variations” in order to reject at once the samples $x^i\left(t\right){|}_0^k$ that are highly improbable given $y\left(t\right){|}_0^k$.

The recursive formula (32(32) $$V_{\mathit{\sigma }\mathit{\rho }}^{\left(obs\right)}\left(t_k\right){=}^{\text{def}}{e}^{-\mathit{\lambda }\mathrm{\Delta }{t}_k}{V}_{\mathit{\sigma }\mathit{\rho }}^{\left(obs\right)}\left(t_{k-1}\right)+\left[y_{\mathit{\sigma }}\left(t_k\right)-y_{\mathit{\sigma }}\left(t_{k-1}\right)\right]\left[y_{\mathit{\rho }}\left(t_k\right)-y_{\mathit{\rho }}\left(t_{k-1}\right)\right],$$(32) ) implies that(34) $$V_{\mathit{\sigma }\mathit{\rho }}^{\left(obs\right)}\left(t_k\right)=\sum _{s=1}^k{e}^{-\mathit{\lambda }(t_k-t_s)}\left[y_{\mathit{\sigma }}\left(t_s\right)-y_{\mathit{\sigma }}\left(t_{s-1}\right)\right][y_{\mathit{\rho }}\left(t_s\right)-y_{\mathit{\rho }}\left(t_{s-1}\right)].$$(34)

We may assume that the considered realization of the process z(t), $\left(x\left(t_0\right),y\left(t_0\right)\right),\left(x\left(t_1\right),y\left(t_1\right)\right),\dots ,\left(x\left(t_s\right),y\left(t_s\right)\right),\dots$, is being obtained consecutively in accordance with finite difference stochastic equations (11(11) $$z_i\left(t_k\right)-z_i\left(t_{k-1}\right)=A_i\left(z\left(t_{k-1}\right),t_{k-1}\right)\mathrm{\Delta }{t}_k+\sum _{j=1}^p{a}_{ij}(z(t_{k-1}),t_{k-1})\mathrm{\Delta }{w}_j(t_k),$$(11) ). The value of observed quadratic variation $V_{\mathit{\sigma }\mathit{\rho }}^{\left(obs\right)}(t_k)$ is being accumulated in parallel, along with that realization of z(t), so that the random increment $\mathrm{\Delta }{y}_{\mathit{\rho }}\left(t_s\right){=}^{\text{def}}{y}_{\mathit{\rho }}\left(t_s\right)-y_{\mathit{\rho }}(t_{s-1})$ will be produced after $(x\left(t_{s-1}\right),y\left(t_{s-1}\right))$ is realized. We are interested to describe the properties of the observed quadratic variations $V_{\mathit{\sigma }\mathit{\rho }}^{\left(obs\right)}\left(t_k\right)$ under condition that the realization of unobservable components, ${x(t)|}_0^k$, is fixed, although it is unknown, and that the next measurement y(t_s) will arrive after the previous measurement $y\left(t_{s-1}\right)$ is already given.

The conditional expectation$$E\left\{\mathrm{\Delta }{y}_{\mathit{\sigma }}\left(t_s\right)\mathrm{\Delta }{y}_{\mathit{\rho }}\left(t_s\right)|x\left(t_{s-1}\right),y(t_{s-1})\right\}=b_{\mathit{\sigma }\mathit{\rho }}\left(x\left(t_{s-1}\right),y\left(t_{s-1}\right),t_{s-1}\right)\mathrm{\Delta }{t}_s.$$

In general case, as it was demonstrated above, the increment Δy = (Δy_m+1, …, Δy_p) in small time step Δt, given $z\left(t\right)=z$, can be described as the Gaussian multidimensional random variable with the probability density (16). The increments Δy(t_k) are independent from the past history before t_k-1 if z(t_k-1) is given. The following properties of Gaussian distributions will be useful in the sequel: If ξ is a Gaussian random variable with Eξ = 0, then Eξ⁴ = 3(Eξ²)². If (ξ₁, ξ₂, ξ₃, ξ₄) is a four-dimensional Gaussian random variable with $E\left\{{\mathit{\xi }}_i\right\}=0$ (for i = 1, …, 4), then$$E\left\{{\mathit{\xi }}_1{\mathit{\xi }}_2{\mathit{\xi }}_3{\mathit{\xi }}_4\right\}=E\left\{{\mathit{\xi }}_1{\mathit{\xi }}_2\right\}E\left\{{\mathit{\xi }}_3{\mathit{\xi }}_4\right\}+E\left\{{\mathit{\xi }}_1{\mathit{\xi }}_3\right\}E\left\{{\mathit{\xi }}_2{\mathit{\xi }}_4\right\}+E\left\{{\mathit{\xi }}_1{\mathit{\xi }}_4\right\}E\left\{{\mathit{\xi }}_2{\mathit{\xi }}_3\right\}.$$

In general case of multidimensional process y(t), we obtain$$E\left\{{(\mathrm{\Delta }{y}_{\mathit{\sigma }}(t_k)\mathrm{\Delta }{y}_{\mathit{\rho }}(t_k)-b_{\mathit{\sigma }\mathit{\rho }}(t_k)\mathrm{\Delta }{t}_k)}^2|x\left(t_{k-1}\right),y(t_{k-1})\right\}=[b_{\mathit{\sigma }\mathit{\sigma }}\left(x\left(t_{k-1}\right),y\left(t_{k-1}\right),t_{k-1}\right)\mathrm{\Delta }{t}_k][b_{\mathit{\rho }\mathit{\rho }}\left(x\left(t_{k-1}\right),y\left(t_{k-1}\right),t_{k-1}\right)\mathrm{\Delta }{t}_k]+{[b_{\mathit{\sigma }\mathit{\rho }}\left(x\left(t_{k-1}\right),y\left(t_{k-1}\right),t_{k-1}\right)\mathrm{\Delta }{t}_k]}^2.$$

Denote shortly $b_{o\mathit{\rho }}(x\left(t_k\right),y\left(t_k\right),t_k){=}^{\text{def}}{b}_{\mathit{\sigma }\mathit{\rho }}(t_k)$. Introduce the following value:(35) $$\sum _{s=1}^k{e}^{-\mathit{\lambda }\left(t_k-t_s\right)}{b}_{\mathit{\sigma }\mathit{\rho }}(t_{s-1})(t_s-t_{s-1}){=}^{\text{def}}\stackrel{\sim }{b_{\mathit{\sigma }\mathit{\rho }}}(t_k).$$(35)

Formula (35) can be written in recursive form, similar to (32):(36) $$\stackrel{\sim }{b_{\mathit{\sigma }\mathit{\rho }}}(t_k)=e^{-\mathit{\lambda }(t_k-t_{k-1})}\stackrel{\sim }{b_{\mathit{\sigma }\mathit{\rho }}}(t_{k-1})+b_{\mathit{\sigma }\mathit{\rho }}\left(t_{k-1}\right)(t_k-t_{k-1}).$$(36)

Consider deviation

(37) $${\mathit{\psi }}_{\mathit{\sigma }\mathit{\rho }}\left(t_k\right){=}^{\text{def}}{V}_{\mathit{\sigma }\mathit{\rho }}^{\left(obs\right)}\left(t_k\right)-\stackrel{\sim }{b_{\mathit{\sigma }\mathit{\rho }}}(t_k).$$(37)

For simplicity of discussion suppose that all the time steps $\mathrm{\Delta }{t}_s=\mathrm{\Delta }t.$, we obtain: $E{\mathit{\psi }}_{\mathit{\sigma }\mathit{\rho }}\left(t_k\right)=o(\mathrm{\Delta }t)$. That notation means that this small value is proportional to a value which decreases faster than Δt if Δt decreases.

The increments Δy_ρ(t_s) and Δy_ρ(t_q) with t_q < t_s are independent when $(x\left(t_{s-1}\right),y\left(t_{s-1}\right))$ is given. We obtain that the dispersion of the deviation ψ_σρ(t_k) (37) (under condition that ${x(t)|}_0^k$ is fixed) is equal to(38) $$Var\left\{{\mathit{\psi }}_{\mathit{\sigma }\mathit{\rho }}\left(t_k\right)\right\}{=}^{\text{def}}E\left\{{\left({\mathit{\psi }}_{\mathit{\sigma }\mathit{\rho }}\left(t_k\right)\right)}^2{|x(t)|}_0^k\right\}=\sum _{s=1}^k{e}^{-2\mathit{\lambda }(t_k-t_s)}E\left\{\left(b_{\mathit{\sigma }\mathit{\sigma }}\left(t_{s-i}\right)b_{\mathit{\rho }\mathit{\rho }}\left(t_{s-1}\right)+{\left[b_{\mathit{\sigma }\mathit{\rho }}\left(t_{s-1}\right)\right]}^2\right){|x(t)|}_0^{s-1}\right\}{(t_s-t_{s-1})}^2.$$(38)

In case if the diffusion coefficients of the observable process, $b_{\mathit{\sigma }\mathit{\rho }}\left(z,t\right),$are the functions of x and t only, so that $b_{\mathit{\sigma }\mathit{\rho }}\left(z\left(t\right),t\right)=b_{\mathit{\sigma }\mathit{\rho }}\left(x\left(t\right),t\right)$ (with $\left(m+1\right)\le \mathit{\sigma },\mathit{\rho }\le p)$we obtain(39) $$Var\left\{{\mathit{\psi }}_{\mathit{\sigma }\mathit{\rho }}\left(t_k\right)\right\}=\sum _{s=1}^k{e}^{-2\mathit{\lambda }\left(t_k-t_s\right)}\left(b_{\mathit{\sigma }\mathit{\sigma }}\left(t_{s-1}\right)b_{\mathit{\rho }\mathit{\rho }}\left(t_{s-1}\right)+{\left[b_{\mathit{\sigma }\mathit{\rho }}\left(t_{s-1}\right)\right]}^2\right){(t_s-t_{s-1})}^2{=}^{\text{def}}\stackrel{}{{\tilde{d}}_{\mathit{\sigma }\mathit{\rho }}}(x(t){|}_0^k,t_k){=}^{\text{def}}\stackrel{}{{\tilde{d}}_{\mathit{\sigma }\mathit{\rho }}}(t_k).$$(39)

The value (39) also can be calculated recursively, similarly to (36). Note, that the value $\stackrel{}{{\tilde{d}}_{\mathit{\sigma }\mathit{\rho }}}(t_k)$ is a value proportional to Δt. If σ = ρ, the above expression takes more concise form:

(40) $${\tilde{d}}_{\mathit{\rho }\mathit{\rho }}(t_k)=2\sum _{s=1}^k{e}^{-2\mathit{\lambda }(t_k-t_s)}{[b_{\mathit{\rho }\mathit{\rho }}\left(t_{s-1}\right)]}^2{\left(t_s-t_{s-1}\right)}^2.$$(40)

In general case, if $b_{\mathit{\sigma }\mathit{\rho }}(z\left(t\right),t)$ depend on y(t), we can find the estimates $E\left\{b_{\mathit{\sigma }\mathit{\rho }}\left(t_s\right){|x(t)|}_0^{s-1}\right\}<K_{\mathit{\sigma }\mathit{\rho }}<\infty$, where $K_{\mathit{\sigma }\mathit{\rho }}$may be determined either as some functions of ${x(t)|}_0^{s-1}$ or as some constants $K_{\mathit{\sigma }\mathit{\rho }}\le K_1$. Then for the variance (38), we obtain the following estimates:$$Var\left\{{\mathit{\psi }}_{\mathit{\sigma }\mathit{\rho }}(t_k)\right\}<2\sum _{s=1}^k{e}^{-2\mathit{\lambda }(t_k-t_s)}(K_{\mathit{\sigma }\mathit{\rho }}{(x\left({t)|}_o^{s-1}\right))}^2{(\mathrm{\Delta }{t}_s)}^2,$$(41) $$Var\left\{{\mathit{\psi }}_{\mathit{\sigma }\mathit{\rho }}(t_k)\right\}<{(K_1)}^2(\frac{1}{\mathit{\lambda }})\mathrm{\Delta }t.$$(41)

The process z(t) is considered on the finite time interval t₀ ≤ s ≤ T, so that the values of $K_{\mathit{\sigma }\mathit{\rho }}(x(t){|}_0^{s-1})$ are limited. Thus, the variance of the deviation (37) is a value proportional to Δt, and if Δt → 0 the observed quadratic variation tends to the limit (33). We shall demonstrate below, that the probability distribution for the deviation ψ_ρρ or ${\mathit{\psi }}_{\mathit{\sigma }\mathit{\rho }}$tends to the Gaussian one as Δt → 0. It follows that the deviations converge to zero if Δt → 0.

Besides, we obtained new analytical formulae (38)–(41), which describe the dispersion of the deviation (37) of the observed quadratic variation when observations are taken in discrete times, with small time steps.

For the system (29), the general analytical formulae (38)–(40), obtained in the present paper, provide characterization of the accuracy of the estimate of volatility $\left|x(t_k)\right|$, which is constructed with the use of the observed quadratic variation. If this accuracy satisfies requirements, the problem of estimation is solved. In that case the value of the given time step Δt can be considered as “small enough”.

In general case, if the value of x(t_k) cannot be uniquely recovered only on the base of observed quadratic variations, the further filtering may be needed.

Finally, we are going to demonstrate that the probability density of the random value ψ_ρρ(t_k) can be approximately described as a Gaussian one. The deviation ψ_ρρ(t_k) contains the sum of squares of the increments $\mathrm{\Delta }{y}_{\mathit{\rho }}(t_s)$, which are independent of the past history before t_s-1 given y(t_s-1), x(t_s-1). The increments $\mathrm{\Delta }{y}_{\mathit{\rho }}(t_s)$ (given y(t_s-1), x(t_s-1)) can be described as the Gaussian random variables with probability density (16). The sum of squares of Gaussian random variables is not a Gaussian random variable. However, the probability density for the random deviation ψ_ρρ(t_k) can be approximately described as the Gaussian probability density. Consider the case when $b_{\mathit{\sigma }\mathit{\rho }}\left(z\left(t\right),t\right)$ (with $\left(m+1\right)\le \mathit{\rho },\mathit{\sigma }\le p)$are the functions of x(t) and t only. Then the variance of deviation ψ_ρρ(t_k) is equal to $\stackrel{}{{\tilde{d}}_{\mathit{\rho }\mathit{\rho }}}(t_k)$ (40), which represents a value proportional to Δt (here the time step Δt can be small, but Δt is finite).

Note, that for a Gaussian random variable ξ with Eξ = 0 the following relations hold true: $E\left\{{\mathit{\xi }}^{2r}\right\}=Q\left(r\right){\left(E\left\{{\mathit{\xi }}^2\right\}\right)}^r$, where Q(r) is a constant, and for the odd moments $E\left\{{\mathit{\xi }}^{2r+1}\right\}=0$, with $r=1,2,3,\dots$. Those relations can be useful in estimating of the moments of $\mathrm{\Delta }{y}_{\mathit{\rho }}\left(t_s\right)$ when (y(t_s-1),x(t_s-1)) is given. Then the consideration of the higher moments of the deviation proves that $E{\left({\mathit{\psi }}_{\mathit{\rho }\mathit{\rho }}(t_k)\right)}^q=o(\mathrm{\Delta }t)$, if q > 2. Then the characteristic function of the random deviation ψ_ρρ(t_k) takes the form similar to (12), which implies the Gaussian expression for its inverse transformation, similar to (13). Hence, the probability density of that deviation can be approximately described as a Gaussian probability density. The smaller Δt the higher the precision of that approximation, although the deviation itself tends to zero if Δt → 0. But in practice the time step should not be chosen too small since in our mathematical model y(t) is considered as a diffusion random process, which is not differentiable at any time t. That is not the case in practice where we have a smooth trajectory $\tilde{y}(t)$, and the quadratic variation of $\tilde{y}(t)$ is equal to zero. Diffusion approximation can be accepted when the time steps Δt_s ≥ τ_cor, where τ_cor represents the characteristic time span of decay of correlation, so that correlation between the random values $\mathrm{\Delta }\stackrel{}{\tilde{y}}\left(t_s\right)$ and $\mathrm{\Delta }\stackrel{}{\tilde{y}}(t_{s-1})$ (given $x\left(t_{s-1}\right),\stackrel{}{\tilde{y}}\left(t_{s-1}\right)$) would be negligibly small (Khazen, 1968, 1971, 1977; Chapters 2, 3; Khazen, 2009).

2.4. Detection and rejection of highly a posteriori improbable samples

In the process of resampling, the a posteriori improbable samples $x{(t)}_0^n$ are being rejected when their weights become small. Besides, it is possible to discard some highly a posteriori improbable samples before all the weights $W_i\left(n\right)$ or $\stackrel{}{{\tilde{W}}_i}(n)$ are calculated. It can be done with the use of the following tests that are based on analysis of the obtained analytical expressions (19(19) $$P_i\left(n\right)=P_0\left(x^i\left(t_0\right),y\left(t_0\right)\right)\prod _{k=1}^{n}P(\mathrm{\Delta }{x}^i(t_k)\text{|}\mathrm{\Delta }y\left(t_k\right),x^i\left(t_{k-1}\right),y\left(t_{k-1}\right)\text{)}P\left(\mathrm{\Delta }y\left(t_k\right)\text{|}{x}^i\left(t_{k-1}\right),y\left(t_{k-1}\right)\right),$$(19) ), (21), (22(22) $$\prod _{s=1}^{n}P(y(t_s)\text{|}y\left(t_{s-1}\right),x\left(t_{s-1}\right)\text{)}=C\left(x\left(t\right){{|}_0^{n-1},y\left(t\right)|}_0^{n-1}\right)exp\left\{-\frac{1}{2}\sum _{s=1}^n\frac{1}{\mathrm{\Delta }{t}_s}{c}_{o\mathit{\rho }}(t_{s-1})[\mathrm{\Delta }{y}_{\mathit{\sigma }}\left(t_s\right)\mathrm{\Delta }{y}_{\mathit{\rho }}\left(t_s\right)-2\mathrm{\Delta }{y}_{\mathit{\sigma }}(t_s)A_{\mathit{\rho }}\left(t_{s-1}\right)\mathrm{\Delta }{t}_s+A_{\mathit{\sigma }}(t_{s-1})A_{\mathit{\rho }}(t_{s-1}){(\mathrm{\Delta }{t}_s)}^2]\right\},$$(22) ).

The measurements ${y(t)|}_0^k$ are being obtained under condition that there is an unobservable realization of the process $x(t),x^{tr}{(t)|}_0^n$, which is unknown. (Here the superscript “tr” stands for “true”). The first and the second moments of the random variables Δy_ρ(t_s) given $x^{tr}{(t)|}_0^k$ are equal to(42) $$E\left\{\mathrm{\Delta }{y}_{\mathit{\rho }}\left(t_s\right)|x^{tr}\left(t_{s-1}\right),y(t_{s-1})\right\}=A_{\mathit{\rho }}(x^{tr}\left(t_{s-1}\right),y\left(t_{s-1}\right),t_{s-1})\mathrm{\Delta }{t}_s{=}^{\text{def}}{A}_{\mathit{\rho }}^{tr}(t_{s-1})\mathrm{\Delta }{t}_s,$$(42) $$E\left\{\mathrm{\Delta }{y}_{\mathit{\sigma }}\left(t_s\right)\mathrm{\Delta }{y}_{\mathit{\rho }}(t_s)|x^{tr}\left(t_{s-1}\right),y(t_{s-1})\right\}=b_{\mathit{\sigma }\mathit{\rho }}\left(x^{tr}\left(t_{s-1}\right),y\left(t_{s-1}\right),t_{s-1}\right)\mathrm{\Delta }{t}_s+o\left(\mathrm{\Delta }{t}_s\right){=}^{\text{def}}{b}_{\mathit{\sigma }\mathit{\rho }}^{tr}\left(t_{s-1}\right)\mathrm{\Delta }{t}_s+o(\mathrm{\Delta }{t}_s)$$

For simplicity of discussion we can assume that all the time steps are equal to Δt.

Consider at first the case when the diffusion coefficients b_σρ of the observable process y(t) do not depend on x(t); then $b_{\mathit{\sigma }\mathit{\rho }}^i=b_{\mathit{\sigma }\mathit{\rho }}^{tr}=b_{\mathit{\sigma }\mathit{\rho }}$ .

The value P_i(k) (19) contains the following cofactor:

(43) $$Cexp\left\{-\frac{1}{2}\sum _{s=1}^k\frac{1}{\mathrm{\Delta }{t}_s}{c}_{o\mathit{\rho }}(t_{s-1})[\mathrm{\Delta }{y}_{\mathit{\sigma }}\left(t_s\right)\mathrm{\Delta }{y}_{\mathit{\rho }}\left(t_s\right)-2\mathrm{\Delta }{y}_{\mathit{\sigma }}(t_s)A_{\mathit{\rho }}^i\left(t_{s-1}\right)\mathrm{\Delta }{t}_s+A_{\mathit{\sigma }}^i(t_{s-1})A_{\mathit{\rho }}^i(t_{s-1}){(\mathrm{\Delta }{t}_s)}^2]\right\},$$(43)

where the normalization factor C does not depend on $x^i{(t)|}_0^k$ in that case.

Consider the hypothesis H_i that the sample $x^i{(t)|}_0^k$ is situated in the vicinity of the true realization $x^{tr}{(t)|}_0^k$, and the “closeness” means that $|A_{\mathit{\rho }}^i\left(t_{s-1}\right)-A_{\mathit{\rho }}^{tr}\left(t_{s-1}\right)|\le D_{\mathit{\rho }}(x^i\left(t_{s-1}\right),y\left(t_{s-1}\right),t_{s-1}){=}^{\text{def}}{D}_{\mathit{\rho }}^i(t_{s-1})$. Here the function D_ρ(x, y, t) can be determined on the basis of some preliminary analysis of the considered dynamical system, for example, we can put $D_{\mathit{\rho }}\left(x,y,t\right)=\left|\frac{\mathit{\partial }{A}_{\mathit{\rho }}(x,y,t)}{\mathit{\partial }{x}_{\mathit{\alpha }}}\right|R_{\mathit{\alpha }}$, where R_α are constants, 1 ≤ α ≤ m. In order to distinguish between two hypotheses, H_i and its negation $\stackrel{}{{\overline{H}}_i}$, consider the following value:(44) $$Q^i\left(t_k\right){=}^{\text{def}}\sum _{s=1}^k{c}_{\mathit{\sigma }\mathit{\rho }}\left(t_{s-1}\right)A_{\mathit{\rho }}^i\left(t_{s-1}\right)\mathrm{\Delta }{y}_{\mathit{\sigma }}(t_s).$$(44)

Note, that Qⁱ(t_k) can be recursively calculated. Denote

(45) $$G^i\left(t_k\right){=}^{\text{def}}\sum _{s=1}^k{c}_{\mathit{\sigma }\mathit{\rho }}\left(t_{s-1}\right)A_{\mathit{\sigma }}^i\left(t_{s-1}\right)A_{\mathit{\rho }}^i(t_{s-1})\mathrm{\Delta }{t}_s.$$(45)

The larger the difference $\left|Q^i(t_k)-G^i(t_k)\right|$ the farther the sample $x^i{(t)|}_0^k$ can be from $x^{tr}{(t)|}_0^k$ Denote $\stackrel{}{\mathrm{\Delta }{\tilde{y}}_{\mathit{\sigma }}}\left(t_s\right)=\mathrm{\Delta }{y}_{\mathit{\sigma }}\left(t_s\right)-A_{\mathit{\sigma }}^{tr}(t_{s-1})\mathrm{\Delta }{t}_s$, $U^i\left(t_k\right){=}^{\text{def}}\sum _{s=1}^k{c}_{\mathit{\sigma }\mathit{\rho }}\left(t_{s-1}\right)A_{\mathit{\rho }}^i\left(t_{s-1}\right)\stackrel{}{\mathrm{\Delta }{\tilde{y}}_{\mathit{\sigma }}}(t_s)$.

Consider the variance of the random value Uⁱ(t_k) when the realization $x^{tr}{(t)|}_0^k$ is given. We obtain$$Var\{U^i\left(t_k\right)\}=\sum _{s=1}^{k}E\{c_{\mathit{\sigma }\mathit{\rho }}\left(t_{s-1}\right)b_{\mathit{\sigma }\mathit{\eta }}\left(t_{s-1}\right)c_{\mathit{\eta }\mathit{\tau }}\left(t_{s-1}\right)A_{\mathit{\rho }}^i\left(t_{s-1}\right)A_{\mathit{\tau }}^i(t_{s-1}){|x^{tr}(t)|}_0^s\text{}}\mathrm{\Delta }{t}_s.$$

If the matrix ∥ b_σρ ∥  is not degenerated, then b_σηc_ητ = δ_στ, where δ_στ is Kronecker symbol, δ_στ = 1 if σ = τ and δ_στ = 0 if σ ≠ τ. In that case, the above expression for $Var\{U^i(t_k)\}$can be written in more concise form:(46) $$Var\{U^i\left(t_k\right)\}=\sum _{s=1}^{k}E\{c_{\mathit{\sigma }\mathit{\rho }}\left(t_{s-1}\right)A_{\mathit{\rho }}^i\left(t_{s-1}\right)A_{\mathit{\sigma }}^i\left(t_{s-1}\right){|x^{tr}(t)|}_0^s\text{}}\mathrm{\Delta }{t}_s$$(46)

In order to reject the samples ${x(t)|}_0^k$, which are highly improbable given ${y(t)|}_0^k$, we can use approximately the following test. If the considered dynamic system is stable, the following value can be used as an approximate estimate for the above expression (46), if the hypothesis H_i is true and the measurements ${y(t)|}_0^k$ are obtained:$$\tilde{V}a{r}^{}\{U^i\left(t_k\right)\}=\sum _{s=1}^k{c}_{\mathit{\sigma }\mathit{\rho }}\left(t_{s-1}\right)A_{\mathit{\rho }}^i\left(t_{s-1}\right)A_{\mathit{\sigma }}^i\left(t_{s-1}\right)\mathrm{\Delta }{t}_s.$$

The probability of the following event:(47) $$\left|Q^i\left(t_k\right)-G^i(t_k)\right|>\text{[}\sum _{s=1}^k{c}_{\mathit{\sigma }\mathit{\rho }}\left(t_{s-1}\right)A_{\mathit{\rho }}^i\left(t_{s-1}\right)D_{\mathit{\sigma }}^i\left(t_{s-1}\right)\mathrm{\Delta }{t}_s+3\sqrt{\stackrel{}{\tilde{V}ar}\left\{U^i(t_k)\right\}}\text{}}$$(47)

is small under condition that H_i is true. Hence, if the inequality (47) is satisfied for the considered sample $x^i{(t)|}_0^k$, the hypothesis $H_i$should be rejected, and that sample $x^i{(t)|}_0^k$ should be discarded. If all the considered sample sequences are discarded, the new sample points $x\left(t_k\right)$ should be generated as initial points, i.e. as samples from the initial distribution $P_0(x(t_k)|y\left(t_k\right))$ (18). For the samples $x^i{(t)|}_0^k$ that remain under consideration the values P_i(k) and W_i(k) or $\stackrel{}{{\tilde{W}}_i}(k)$ will be recursively computed. The samples that were not discarded can be continued with a few “offsprings”. Due to the test (47) the samples $x^i{(t)|}_0^k$ that are highly a posteriori improbable (when the sequence of observations $y\left(t\right){|}_0^k$ is given) will be rejected at once, before the weights W_j(k) or $\stackrel{}{{\tilde{W}}_i}(k)$ are computed.

Consider the case when some diffusion coefficients depend on x(t). In some cases, as it was shown above, the unobservable value of $x\left(t\right)$ can be restored with the use of the observed quadratic variations. In other cases, if x(t) cannot be uniquely recovered, it is purposeful to use observed quadratic variations $V_{\mathit{\rho }\mathit{\rho }}^{\left(obs\right)}(t_k)$ in order to reject the samples $x^i{(t)|}_0^k$ that are highly a posteriori improbable. That test is similar to the above test (47). Assume that if hypothesis H_i is true, the sample $x^i{(t)|}_0^k$ is “close” to the true realization $x^{tr}{(t)|}_0^k$, and the “closeness” means that$$\left|b_{\mathit{\rho }\mathit{\rho }}\left(x^i\left(t_s\right),y\left(t_s\right),t_s\right)-b_{\mathit{\rho }\mathit{\rho }}(x^{tr}\left(t_s\right),y\left(t_s\right),t_s)\right|\le B_{\mathit{\rho }}(x^i\left(t_s\right),y\left(t_s\right),t_s){=}^{\text{def}}{B}_{\mathit{\rho }}^i(t_s).$$

Here the function B_ρ(x, y, t) can be determined as $B_{\mathit{\rho }}\left(x,y,t\right)=\left|\frac{\mathit{\partial }{b}_{\mathit{\rho }\mathit{\rho }}(x,y,t)}{\mathit{\partial }{x}_{\mathit{\alpha }}}\right|R_{\mathit{\alpha }}$, where R_α are constants. Denote shortly (similarly to (35), (36), (38), (39)):$${\tilde{b}}_{\mathit{\sigma }\mathit{\rho }}^i(t_k){=}^{\text{def}}\stackrel{}{{\tilde{b}}_{\mathit{\sigma }\mathit{\rho }}}(x^i(t){|}_0^k,y(t){|}_0^k,t_k){\tilde{d}}_{\mathit{\sigma }\mathit{\rho }}^i(t_k){=}^{\text{def}}\stackrel{}{{\tilde{d}}_{\mathit{\sigma }\mathit{\rho }}}(x^i(t){|}_0^k,y(t){|}_0^k,t_k),$$$${\tilde{B}}_{\mathit{\sigma }\mathit{\rho }}^i(t_k){=}^{\text{def}}\stackrel{}{{\tilde{B}}_{\mathit{\rho }}}\left(x^i\left(t\right){{|}_0^k,y(t)|}_0^k,t_k\right){=}^{\text{def}}\sum _{s=1}^k{e}^{-\mathit{\lambda }\left(t_k-t_s\right)}{B}_{\mathit{\rho }}\left(x^i\left(t_s\right),y\left(t_s\right),t_s\right)(t_s-t_{s-1})$$

$$\stackrel{}{{\tilde{B}}_{\mathit{\rho }}^i}(t_k)=e^{-\mathit{\lambda }\left(t_k-t_{k-1}\right)}\stackrel{}{{\tilde{B}}_{\mathit{\rho }}^i}(t_{k-1})+B_{\mathit{\rho }}^i(t_k)(t_k-t_{k-1})$$

Consider the case when $b_{\mathit{\sigma }\mathit{\rho }}\left(z\left(t\right),t\right)=b_{\mathit{\sigma }\mathit{\rho }}(x\left(t\right),t)$ (with $\left(m+1\right)\le \mathit{\sigma },\mathit{\rho }\le p$). Then the probability of the following event(48) $$\left|V_{\mathit{\rho }\mathit{\rho }}^{\left(obs\right)}\left(t_k\right)-{\tilde{b}}_{\mathit{\rho }\mathit{\rho }}^i(t_k)\right|>\left[\stackrel{}{{\tilde{B}}_{\mathit{\rho }}^i}\left(t_k\right)+3\sqrt{\stackrel{}{{\tilde{d}}_{\mathit{\rho }\mathit{\rho }}^i}(t_k)}\right]$$(48)

is small if the hypothesis H_i is true. Hence, if the inequality (48) is satisfied at least for one number ρ, with (m + 1) ≤ ρ ≤ p, then the hypothesis H_i should be rejected, and the sample $x^i{(t)|}_0^k$ should be discarded at once, before the calculation of all the weights W_j(k) or $\stackrel{}{{\tilde{W}}_j}(k)$ is done.

In case if $b_{\mathit{\sigma }\mathit{\rho }}\left(z\left(t\right),t\right)=b_{\mathit{\sigma }\mathit{\rho }}(x\left(t\right),y\left(t\right),t)$, with $E\{{\left(b_{\mathit{\sigma }\mathit{\rho }}\right)}^2|x^{tr}(t){|}_0^k\}\le K_1^2=Const$, the following test can be used instead of (48):(49) $$\left|V_{\mathit{\rho }\mathit{\rho }}^{\left(obs\right)}\left(t_k\right)-{\tilde{b}}_{\mathit{\rho }\mathit{\rho }}^i(t_k)\right|>\left[\stackrel{}{{\tilde{B}}_{\mathit{\rho }}^i}\left(t_k\right)+3K_1\sqrt{\frac{1}{\mathit{\lambda }}\mathrm{\Delta }t}\right]$$(49)

2.5. Incorporation of the derivative of the process of quadratic variation into the set of observed data. Systems with some additional precise observations

In case of continuous time of observation, the process of quadratic variation can be included into the set of observed processes, as it was demonstrated in Section 2.3.

For simplicity of notation, consider the system of one-dimensional processes x(t), y(t):

(50) $$dx=m\left(x,t\right)dt+{\mathit{\sigma }}_1\left(x,t\right)d{w}_1(t),$$(50)

(51) $$dy=h\left(x,t\right)dt+{\mathit{\sigma }}_2\left(x,t\right)d{w}_2(t),$$(51)

where w₁(t), w₂(t) are independent Wiener processes, the functions m(x, t), h(x, t) satisfy Lipschitz condition, and the functions $h\left(x,t\right),{\mathit{\sigma }}_1(x,t)$, σ₂(x, t) are twice differentiable with respect to x. The process

(52) $$V(t)=\underset{t_0}{\overset{t}{\int }}{e}^{-\mathit{\alpha }\left(t-s\right)}{\mathit{\sigma }}_2^2\left(x\left(s\right),s\right)ds$$(52)

can be considered theoretically as an observed process, given the trajectory $y\left(s\right)$for t₀ ≤ s ≤ t. Here $\mathit{\alpha }=Const\ge 0.$The process V(t) obeys the ordinary differential equation$$\frac{dV}{dt}=-\mathit{\alpha }V+{\mathit{\sigma }}_2^2(x\left(t\right),t).$$

For continuous time of observation and multidimensional diffusion processes $(x\left(t\right),y\left(t\right))$, in some cases the possible values of x(t) are precisely determined by the equalities similar to (52) (and, consequently, x(t) can be restored precisely, at least theoretically), and in other cases the trajectory x(s) must be localized in some “layer”, where those equalities are satisfied, given the observed path of the process of quadratic variation, V(s), for t₀ ≤ s ≤ t.

Theoretically, it is possible to consider the derivative $\frac{dV(t)}{dt}{=}^{\text{def}}u(t)$ as an observed process, and incorporate it into the set of observed components: $(y\left(t\right),u\left(t\right))$. For the system (50), (51), the process u(t) satisfies the following stochastic differential equation (in the sense of Ito):(53) $$du=-\mathit{\alpha }udt+2{\mathit{\sigma }}_2\left(x,t\right)\frac{\mathit{\partial }{\mathit{\sigma }}_2}{\mathit{\partial }x}\left[m\left(x,t\right)dt+{\mathit{\sigma }}_1\left(x,t\right)d{w}_1\left(t\right)\right]+\left[{\mathit{\sigma }}_2\frac{{\mathit{\partial }}^2{\mathit{\sigma }}_2}{\mathit{\partial }{x}^2}+{\left(\frac{\mathit{\partial }{\mathit{\sigma }}_2}{\mathit{\partial }x}\right)}^2\right]{\mathit{\sigma }}_1^2\left(x,t\right)dt.$$(53)

The above SDE (53) could be augmented to the system of SDE (50)–(51) for x(t), y(t).

But in practice we cannot assume that the process u(t) (or, equivalently, the process ${\mathit{\sigma }}_2^2(x\left(t\right),t)$) is straightforwardly observed and known precisely. Practically, the process V(t_k) is not precisely given. Even in case if the deviations of “observed quadratic variation” V^(obs)(t_k) from V(t_k) are small, in obtaining the estimates for derivative, $\frac{dV(t)}{dt}{|}_{t=t_k}{=}^{\text{def}}u(t_k)$, the sought values could be inevitably corrupted by errors.

However, if additional observations of y(t) are available, namely, y(t₀ + sδt), with δt ≪ Δt, $s=1,2,\dots ,$ then the estimates $\stackrel{}{\widehat{u}}(t_k)$ for the values of u(t_k) can be improved. Suppose that there is a special digital device that takes measurements of $y(t_s^{\prime })$ at times $t_s^{\prime }=t_0+s\mathit{\delta }t$, with $t_{k-1}<t_s^{\prime }\le t_k$, and use them in order to calculate the estimate $\stackrel{}{\widehat{u}}\left(t_k\right)$. For the moment, consider the case if we had $\stackrel{}{\widehat{u}}\left(t_k\right)=u(t_k)$. Then the problem of filtering (based on observations y(t_k), u(t_k) taken at discrete times t_k, with small time steps Δt_k) could be solved effectively with the use of the new algorithm of simulation (17) and the estimates (26) and algorithms suggested in the present paper. For the system (50), (51), and (53), we obtain in (17): b₁₁ - b_1σc_σρb_ρ1 ≡ 0. That means that the sample paths $x^i{(t)|}_1^n$ would be deterministically continued given the observations $u\left(t\right){|}_1^n$ and the initial values xⁱ(t₁), with i = 1, …, N. It is purposeful to choose the probability distribution of the initial random samples xⁱ(t₁) in such way that $E\left\{{\mathit{\sigma }}_2^2(x^i\left(t_1\right),t_1)\right\}=u(t_1)$, with some random scattering. Meanwhile, the weights $\stackrel{}{{\tilde{W}}_i}(n)$ are being defined by (25), (26) with the use of observations ${y(t)|}_1$. Thus, the a posteriori distribution is being determined not only by $u\left(t\right){|}_1^n$ but also by $y\left(t\right){|}_1^n$. Then, as it was shown in Section 2.2, the estimate (26) converges with probability one (due to the Strong Law of Large Number, as N → ∞) to the sought a posteriori expectation that provides solution to the filtering problem. That algorithm could be implemented with the use of parallel computing.

The random errors of the obtained estimates $\widehat{u}(t_k)$ can be taken into account in the following way. The random deviations $\widehat{u}(t_k)-u(t){=}^{\text{def}}\mathit{\epsilon }(t_k)$ can be approximately described as independent of each other (for different k) Gaussian random values. Then $\widehat{u}(t_k)$ can be considered as measurements (taken at discrete times ) of the process $\widehat{u}(t_k)$ that satisfies the following SDE:(54) $$d\stackrel{}{\widehat{u}}=-\mathit{\alpha }\stackrel{}{\widehat{u}}dt+2{\mathit{\sigma }}_2\left(x,t\right)\frac{\mathit{\partial }{\mathit{\sigma }}_2}{\mathit{\partial }x}\left[m\left(x,t\right)dt+{\mathit{\sigma }}_1\left(x,t\right)d{w}_1\left(t\right)\right]+\left[{\mathit{\sigma }}_2\frac{{\mathit{\partial }}^2{\mathit{\sigma }}_2}{\mathit{\partial }{x}^2}+{\left(\frac{\mathit{\partial }{\mathit{\sigma }}_2}{\mathit{\partial }x}\right)}^2\right]{\mathit{\sigma }}_1^2\left(x,t\right)dt+{\mathit{\sigma }}_{3}d{w}_3(t),$$(54)

where $w_3(t)$ is an independent Wiener process, and ${\mathit{\sigma }}_3^2(t_k)$ is proportional to the variance of $\mathit{\epsilon }\left(t_k\right)$, i.e. $E({\mathit{\epsilon }}^2(t_k))$. Thus, σ₃ can be determined as some function of x, t. Then the solution of the filtering problem for the system (50), (51), (54) can be obtained in the similar way as above.

Note that the systems that are similar to (50), (51), (53) arise if the measurements of some given function $H(x\left(t\right),t)$ are known without errors. Here the function H(x, t) is supposed to be twice differentiable with respect to x. Then the solution of the filtering problem (given observations of y(t) and $H(x\left(t\right),t)$ at discrete times t_k = t₀ + kΔt, with small time step Δt) can be obtained with the use of the algorithm of simulation (17) and the estimates (26), similarly to the solution described above for the system (50), (51), (53).

2.6. The problem of nonlinear filtering of a signal. Change of the mathematical model that describes the process of observations

For simplicity of notations, consider one-dimensional processes x(t), y(t). In most of works devoted to nonlinear filtering problems the following hidden Markov model is accepted for description of a signal process x(t) and an observed process y(t):(55) $$dx=m\left(x,t\right)dt+{\mathit{\sigma }}_1\left(x,t\right)d{w}_1\left(t\right),$$(55) (56) $$dy=h\left(x(t),t\right)dt+{\mathit{\sigma }}_{2}d{w}_2\left(t\right),$$(56)

where w₁(t), w₂(t) are independent Wiener processes, σ₂ ≡ Const, the functions m(x, t), $h\left(x,t\right),$ σ₁(x, t) satisfy Lipschitz condition, and the functions h(x, t), σ₁(x, t) are twice differentiable with respect to x.

In many applications, it would be possible to describe an observed process Y(t) as follows:(57) $$Y\left(t\right)=h\left(x\left(t\right),t\right)+u(t),$$(57)

where u(t) represents the errors of measurements. The “noise process” u(t) can be considered as a Gaussian Markov process, which can be described by the following stochastic differential equation:(58) $$du=-\mathit{\beta }udt+{\mathit{\sigma }}_{3}d{w}_3(t),$$(58)

where β and σ₃ are constant, and w₃(t) is an independent Wiener process. Then $E\left\{u\right\}=0$, the variance $E\left\{u^2\right\}={\mathit{\sigma }}_3^2/2\mathit{\beta }$, and the correlation function $E\left\{u\left(t\right)u(t+\mathit{\tau })\right\}=\left({\mathit{\sigma }}_3^2/2\mathit{\beta }\right)e^{-\mathit{\beta }\mathit{\tau }}$.

Denote $z_1\left(t\right){=}^{\text{def}}x(t)$, $z_2(t){=}^{\text{def}}u(t)$, $z_3(t){=}^{\text{def}}Y(t)$.

Then the model (55), (57), (58) of the signal process x(t) and the observed process Y(t) can be described as follows:(59) $$d{z}_1=m\left(z_1,t\right)dt+{\mathit{\sigma }}_1\left(z_1,t\right)d{w}_1\left(t\right),d{z}_2=-\mathit{\beta }{z}_{2}dt+{\mathit{\sigma }}_{3}d{w}_3(t),d{z}_3=\frac{\mathit{\partial }h(z_1,t)}{\mathit{\partial }t}dt+\frac{\mathit{\partial }h}{\mathit{\partial }{z}_1}\left[m\left(z_1,t\right)dt+{\mathit{\sigma }}_1\left(z_1,t\right)d{w}_1\right]+\frac{1}{2}\frac{{\mathit{\partial }}^{2}h}{\mathit{\partial }{z}_1^2}{\mathit{\sigma }}_1^2\left(z_1,t\right)dt-\mathit{\beta }{z}_{2}dt+{\mathit{\sigma }}_{3}d{w}_3.$$(59)

Here z₃(t) represents the observed process,  and z₁(t), z₂(t) are unobservable components of the process $z\left(t\right){=}^{\text{def}}(z_1\left(t\right),z_2\left(t\right),z_3\left(t\right))$. It is easy to see that the diffusion coefficients $b_{11}={\mathit{\sigma }}_1^2(z_1,t)$, $b_{13}=\frac{\mathit{\partial }h}{\mathit{\partial }{z}_1}{\mathit{\sigma }}_1^2(z_1,t)$, b₁₂ = 0, $b_{22}={\mathit{\sigma }}_3^2$, $b_{23}={\mathit{\sigma }}_3^2$, $b_{33}={\left(\frac{\mathit{\partial }h}{\mathit{\partial }{z}_1}{\mathit{\sigma }}_1\left(z_1,t\right)\right)}^2+{\mathit{\sigma }}_3^2$. Consider the matrix ∥ r_αβ ∥ = ∥ b_αβ - b_ασc_σρb_ρβ ∥ , which was introduced in (15): This is the variance–covariance matrix for the conditional probability density of the increments $\mathrm{\Delta }{z}_1=z_1\left(t+\mathrm{\Delta }t\right)-z_1(t)$, $\mathrm{\Delta }{z}_2=z_2\left(t+\mathrm{\Delta }t\right)-z_2(t)$, given z(t) and the increment $\mathrm{\Delta }{z}_3=z_3\left(t+\mathrm{\Delta }t\right)-z_3(t)$. For the system (59) we find: $r_{11}={\mathit{\sigma }}_1^2\left[1-\frac{1}{1+{\mathit{\sigma }}_3^2/{\left(\frac{\mathit{\partial }h}{\mathit{\partial }{z}_1}{\mathit{\sigma }}_1\right)}^2}\right]$,$$r_{12}=-{\mathit{\sigma }}_3^2{\mathit{\sigma }}_1^2\frac{\mathit{\partial }h}{\mathit{\partial }{z}_1}{\left[{\left(\frac{\mathit{\partial }h}{\mathit{\partial }{z}_1}{\mathit{\sigma }}_1\right)}^2+{\mathit{\sigma }}_3^2\right]}^{-1},r_{22}={\mathit{\sigma }}_3^2\left[1-{\mathit{\sigma }}_3^2{\left[{\left(\frac{\mathit{\partial }h}{\mathit{\partial }{z}_1}{\mathit{\sigma }}_1\right)}^2+{\mathit{\sigma }}_3^2\right]}^{-1}\right],andDet\Vert r_{\mathit{\alpha }\mathit{\beta }}\Vert \equiv 0.$$

The case when σ₃ is small corresponds to the small errors of measurements. All the entries of the matrix ∥ r_αβ ∥  (shown above) is small if σ₃ is small. That means that most of the samples $(z_1\left(t\right),z_2\left(t\right)){|}_0^n$ (given $z_3{(t)|}_0^n$) simulated with the use of the algorithm (17) continue to be localized in the domain where the a posteriori probability density is concentrated. Besides, the test (48) can be implemented (unless the function $h\left(x,t\right)$ is linear with respect to $x$and σ₁ ≡ Const) in order to discard a posteriori improbable samples with the use of “observed quadratic variation” (for the observations made in discrete time) of the observed process z₃(t).

Thus, although the dimensionality of the unobservable process increased for the model of description (59) in comparison with the model (55), (56), the problem of nonlinear filtering (with observations made in discrete time) will be solved effectively with the use of the new Monte Carlo estimates (23)–(27) and with the use of the new algorithm of simulation (17), derived in the present paper.

Note that the model (59) is a hidden Markov model where unobservable process $(z_1\left(t\right),z_2(t))$ is a Markov process in itself. But if the samples $(z_1\left(t\right),z_2\left(t\right)){|}_0^n$ were simulated as trajectories of the Markov process $(z_1\left(t\right),z_2\left(t\right))$ in itself (as it is suggested in particle filters, developed and studied in the literature), then, in case of the low “level of intensity of noise” (i.e. if σ₃ is small) the large part of generated samples would be localized in the domain of low a posteriori probability (and that part is growing when ${\mathit{\sigma }}_3$ decreases). Then, in order to obtain some feasible accuracy of filtering, the number N should be increased significantly if the value of σ₃ decreased. Meanwhile, the large amount of calculations would be wasted in vain for processing of the samples of low a posteriori probability. Thus, this example (the system (59)) also proves the advantage of the new algorithms, developed in the present paper, in comparison with the known algorithms for hidden Markov models.

3. Conclusion

In the Conclusion, the new results obtained and presented by this author are shortly summarized and underscored.

(1)	It was analytically proved by this author at first in the book (Khazen, 2009) that the increment of a multidimensional continuous diffusion Markov process z(t), $\mathrm{\Delta }z=z\left(t+\mathrm{\Delta }t\right)-z(t)$, in small time interval Δt, given z(t), obeys asymptotically to the multidimensional Gaussian probability distribution (with the first and second moments that are determined by its drift and diffusion coefficients; see expression (13)). As a corollary, the analytical expressions (14)–(16) are obtained, which describe Gaussian conditional probability density for the increment Δx of the unobservable components given observation of Δy and given the value z(t); here Δz = (Δx, Δy). The Gaussian probability density for Δy is also obtained (see (16)). In the present paper the new, precise algorithm (in closed analytical form) is obtained for simulation of samples x(tk) when ${x(t)|}_0^{k-1}$ and ${y(t)|}_0^k$ are given (see (17)). It is important that the influence of the current “next” observation $y\left(t_k\right)$ is taken into account when samples $x\left(t_k\right)$ are being simulated.
(2)	The new Monte Carlo estimates of functions $f(x\left(t_n\right))$ given ${y(t)|}_0^n$ are obtained and presented in explicit, precise, closed analytical form, under the condition that the samples ${x(t)|}_0^n$ are being simulated with the use of the new proposed algorithm (17) (see (26), (25), (27)). The convergence of the Monte Carlo estimates (as the number N of samples increases, N → ∞) is guaranteed by the Strong Law of Large Number (the Theorem of Kolmogorov). For the first time, the estimates, (26), (25), (27), are obtained for the general case of multidimensional diffusion Markov process $(x\left(t\right),y\left(t\right))$. In the particular case when x(t) is a Markov process in itself and the process $(x\left(t\right),y\left(t\right))$ represents the simple Hidden Markov Model described by (28), the estimate (26) takes the form that is in agreement with the estimate that was obtained earlier in this particular case.
(3)	For the first time, some tests are developed in order to discard at once the sample sequences ${x(t)|}_0^n$ that are highly improbable given the observation ${y(t)|}_0^n$ (see Section 2.4). This is important in order to prevent the “degeneration” of the set of considered samples. These tests are being performed independently for each sample ${x(t)|}_0^n$, and that is important also for the implementation with the use of parallel computing.
(4)	Some branching sampling procedures are developed. The standard Sequential Importance Resampling (SIR) procedure is generalized for the general case of partially observed multidimensional diffusion Markov process, when the new algorithm of simulating of samples (17) should be used. Some new versions of branching resampling procedures are proposed, which can be easier implemented with the use of parallel computing.
(5)	Analytical investigation of observed quadratic variations is developed (in case when diffusion coefficients of observed components, y(t), depend on unobservable components x(t)) (see Section 2.5). For the first time, the analytical formulae (in explicit and closed form) that determine the dispersions of deviations of observed quadratic variations are obtained (see (38), (39)), and it is also proved that the deviations obey asymptotically (when Δt → 0) to the Gaussian distribution.
(6)	The important particular case of nonlinear filtering of a signal is considered in Section 2.6. Significant advantage of the new algorithms and estimates obtained in the present paper in comparison with the particle filters suggested earlier in the literature for Hidden Markov Models is demonstrated.
(7)	In Section 2.5 it is demonstrated that the new filtering algorithms, developed in the present paper, provide opportunity to incorporate the observed quadratic variations into the set of observed data. They also provide opportunity to use additional precise observations in case if such observations are available. These results are new, and they also confirm the advantage of the proposed new solution of the nonlinear filtering problem.

The obtained new algorithms and estimates extend the range of applications of particle filtering beyond Hidden Markov Models and improve performance.

Additional information

Funding

Funding. The author received no direct funding for this research.

Notes on contributors

Ellida M. Khazen

Ellida M. Khazen graduated from Department of Mathematics of Lomonosov Moscow State University (Moscow, Russia, USSR) with Honor Diploma in 1959 (Master Degree). She received her PhD Degree in Mathematics from Department of Mathematics of Lomonosov Moscow State University in 1962. She received her Doctor Degree in Applied Mathematics in 1994. She was working as Senior Scientist Researcher Mathematician at Moscow Scientific Research Institute of Device Automation in 1962–1996, and also as a Visiting Lecturer at Department of Mathematics of Lomonosov Moscow State University and as a Visiting Lecturer at Moscow Institute of Radio Engineering, Electronics and Automation (MIREA). She has published two scientific monographs (in Russian and in English), chapters in books (in Russian), and more than 50 scientific research papers in the areas of the theory of turbulence onset, the theory of random processes, filtering and signal detection, optimal statistical decisions and statistical sequential analysis, informational estimation of risks, and the theory of optimal control.