Ivan Gorynin

Senior Member, IEEE Stéphane Derrode

Emmanuel Monfrini

Wojciech Pieczynski

Markov Systems with Applications to Stochastic Volatility

Keywords: Non-linear systems, Stochastic volatility model, Optimal statistical filter, Conditionally Gaussian linear state-space model, Filtering in switching systems, Kalman filter

published or not. The documents may come

I. INTRODUCTION

Let us consider two random sequences X N 1 = (X1, . . . , XN ) and Y N 1 = (Y1, . . . , YN ), taking their values in R m and R q , respectively. X N 1 is hidden, while Y N 1 is observed. In this paper, we focus on the optimal filtering problem, which consists in the sequential search of X N 1 from Y N 1 . More precisely, we present a non-stochastic iterative algorithm which computes the expectation of Xn conditional on Y n 1 . We present a workable approach for filtering in general stationary (or asymptotically stationary) Markov dynamic systems, provided that one can sample a realization of the systems. This approach makes use of a recent switching model, in which fast exact optimal filtering is computationally feasible.

The idea is to approximate a given non-linear non-Gaussian system by a switching Gaussian system. Such ideas are not new, as it is well known that a wide range of probability distributions can be approximated by a Gaussian mixture. In the standard switching models, e.g. in jump Markov linear systems (JMLSs), there is no known fast exact optimal filtering algorithm [START_REF] Andrieu | Efficient particle filtering for jump Markov systems. Application to time-varying autoregressions[END_REF]- [START_REF] Doucet | A tutorial on particle filtering and smoothing: Fifteen years later[END_REF], therefore such approximations are useless in their context. However, there are some recent switching models, in which fast exact optimal filtering is computationally feasible, e.g. the conditionally Markov switching hidden linear model (CMSHLM [START_REF] Pieczynski | Exact filtering in conditionally Markov switching hidden linear models[END_REF]) and the conditionally Gaussian observed Markov switching model (CGOMSM [START_REF] Abbassi | Kalman filtering approximations in triplet Markov Gaussian switching models[END_REF], [START_REF] Abbassi | Optimal filter approximations in conditionally Gaussian pairwise Markov switching models[END_REF]). The novelty of the paper is to use them as an approximation of some given model, then to apply the corresponding filtering algorithms. Besides, let us notice that any standard JMLS can be approximated by a CMSHLM [START_REF] Derrode | Exact fast computation of optimal filter in Gaussian switching linear systems[END_REF], [START_REF] Petetin | A class of fast exact Bayesian filters in dynamical models with jumps[END_REF].

To be more precise, CMSHLMs are Markov triplet

(X N 1 , R N 1 , Y N 1 ) models, where R N 1
is a chain of switches, and where there exists an exact fast filtering algorithm which is as fast as the standard Kalman filter. The main difference between CMSHLMs and JMLSs is that in CMSHLMs (R N 1 , Y N 1 ) is Markovian and (X N 1 , R N 1 ) is not necessarily Markovian, contrary to JMLSs where (X N 1 , R N 1 ) is Markovian and (R N 1 , Y N 1 ) is not necessarily Markovian.

Ivan Gorynin, Emmanuel Monfrini and Wojciech Pieczynski are with SAMOVAR, Telecom Sudparis, CNRS, Université Paris-Saclay, 9 rue Charles Fourier, 91000 Evry, FRANCE. e-mail: ivan.gorynin@telecom-sudparis.eu Stéphane Derrode is with Ecole Centrale de Lyon, LIRIS, CNRS UMR 5205, 69130 Écully, FRANCE. e-mail: stephane.derrode@ec-lyon.fr

The CGOMSM is a sub-model of CMSHLM, in which

(X N 1 , Y N 1 ) is Gaussian conditional on R N 1 .
As the main novelty, deriving from its prior simplified version introduced in [START_REF] Derrode | Fast filter in nonlinear systems with application to stochastic volatility model[END_REF], we put forward a method for approximating any stationary non-linear non-Gaussian Markov model by a CGOMSM, in order to perform fast filtering in this approximation. The other main novelty compared with [START_REF] Derrode | Fast filter in nonlinear systems with application to stochastic volatility model[END_REF] is the design of a new EM (Expectation-Maximization)-based approximation algorithm, as well as several new experiments. In particular, we present consistent results that we obtained for a recent asymmetric stochastic volatility model.

Our method can be used as an alternative to the particle filter (PF) based methods, which are widely used in different areas like finance [START_REF]A survey of sequential Monte Carlo methods for economics and finance[END_REF]- [START_REF] Kim | State-space models with regime switching[END_REF] and tracking [START_REF] Arulampalam | A tutorial on particle filters for online non-linear/non-Gaussian Bayesian tracking[END_REF]- [START_REF] Doucet | A tutorial on particle filtering and smoothing: Fifteen years later[END_REF], [START_REF] Ristic | Beyond the Kalman filter: particle filters for tracking applications[END_REF]. Although PFs are asymptotically optimal, their use may be problematic due to the particle degeneracy phenomenon, or to the need of a large amount of particles when the dimension of the hidden space is high. The method proposed is fundamentally different and is free of these impediments. Besides, our filtering algorithm is exact and is as fast as the standard Kalman filter.

Let us mention some other non-linear non-Gaussian filters which similarly use a mixture model to represent the filtering pdf. Such examples include the Gaussian Sum Filter (GSF) [START_REF] Sorenson | Recursive Bayesian estimation using Gaussian sums[END_REF], [START_REF] Šimandl | Filtering, prediction and smoothing with Gaussian sum representation[END_REF], the Unscented Kalman Filter (UKF) [START_REF] Wan | The unscented Kalman filter for nonlinear estimation[END_REF], the Unscented Gaussian Sum Filter (UGSF) [START_REF] Šimandl | Sigma point Gaussian sum filter design using square root unscented filters[END_REF], and the Gaussian Sum Unscented Kalman Filter (GSUKF) with adaptive scaling parameters [START_REF] Straka | Gaussian sum unscented Kalman filter with adaptive scaling parameters[END_REF]. The main difference between these methods and ours is that they rely on a supplementary approximation to prevent the number of mixands to growth exponentially with time (this number remains constant in our representation). We show through some experiments that the GSF may be disadvantaged due to this additional approximation, whereas our method attains the optimal accuracy. Moreover, contrary to PF and GSUKF, our filter is designed for stationary (or asymptotically stationary) systems only. However, let us notice that the context of stationary systems is relevant for many applications, and different efficient approaches have been recently proposed [START_REF] Feng | Real-time state estimator without noise covariance matrices knowledge-fast minimum norm filtering algorithm[END_REF], [START_REF] Feng | Kalman filter with recursive covariance estimation-sequentially estimating process noise covariance[END_REF].

We provide two series of experiments. The first one is devoted to filtering for a standard stochastic volatility model (SV) [START_REF] Ghysels | Stochastic volatility[END_REF]- [START_REF] Taylor | Modeling Financial Times Series[END_REF]. The second is dedicated to filtering for the asymmetric stochastic volatility model (ASV) [START_REF] Nikolaev | Nonlinear filtering of asymmetric stochastic volatility models and value-at-risk estimation[END_REF]- [START_REF] Takashi | Estimating stochastic volatility models using daily returns and realized volatility simultaneously[END_REF], which extends the SV model by incorporating a leverage effect. We conclude from these experiments that once the approximation is established, our method has the same performance in terms of the mean squared error that the particle filter, while being significantly faster.

The rest of the paper is organized as follows. In the next Section, we describe and justify our approximation setting. In Section three, we recall some crucial properties of the CGOMSM, and then explain our approximation algorithm in Section four. Fifth Section contains experiments and the last Section gives some concluding remarks and prospects for further works.

II. SWITCHING MODEL APPROXIMATION

Let us consider two random sequences X N 1 and Y N 1 as previously described. The common model for (X N 1 , Y N 1 ) is the hidden Markov model (HMM), where X N 1 is Markovian, Y1, . . . , YN are independent conditional on X N 1 , and the distribution of each Yn depends only on Xn given X N 1 . We can define the joint pdf p x N 1 , y N 1 via p (x1, y 1 ) and two recursions:

Xn+1 = F(Xn, Un); (1a) Yn = G(Xn, Vn), (1b)
where F, G are appropriate mappings representing the Markov kernel and the observation kernel respectively, and U1, V1, . . . , UN , VN are appropriate independent variables. Equivalently, any HMM is given by p (x1, y 1 ) and conditional densities p (xn+1 |xn ) , p (y n |xn ) for each n in {1, . . . , N -1} (in {2, . . . , N }, respectively). We recall that if, in addition, (X N 1 , Y N 1 ) is Gaussian, then we deal with a Gaussian linear model in which the Kalman filter is statistically optimal.

The HMM can be extended to the "pairwise Markov model" (PMM), where the recursion

(Xn+1, Yn+1) = H(Xn, Yn, Wn), (2)
and thus ( 1) is a particular case of (2) where p (xn+1 |xn,

y n ) = p (xn+1 |xn ) and p y n+1 |xn, y n , xn+1 = p y n+1 |xn+1 .
In our setting, we suppose that X N 1 and Y N 1 are stationary signals, which means that the distributions p xn, y n , xn+1, y n+1 do not depend on n, i.e. for any n, (Xn, Yn, Xn+1, Yn+1) is equal in distribution to (X1, Y1, X2, Y2), what we note by p xn, y n , xn+1, y n+1 = p (x1, y 1 , x2, y 2 ). However, in practice, the algorithm that we put forward applies to asymptotically stationary signals as well, which is detailed in Remark II.1.

The idea proposed in [START_REF] Derrode | Fast filter in nonlinear systems with application to stochastic volatility model[END_REF] is to approximate p x N 1 , y N 1 with a CGOMSM marginal distribution. That is to perform an exact fast filtering afterwards. More precisely, since the model ( 2) is stationary, its distribution derives from p (x1, y 1 , x2, y 2 ), as it provides p (x1, y 1 ) and p xn+1, y n+1 |xn,

y n = p (x2, y 2 |x1, y 1 ) for each n = 1, . . . , N -1. Besides, p (x1, y 1 , x2, y 2 ) can be approximated using a mixture of K 2 components p (x1, y 1 , x2, y 2 ) ≈ 1≤i,j≤K αijpij(x1, y 1 , x2, y 2 ), (4)
where pij(x1, y 1 , x2, y 2 ) are Gaussian distributions which verify some further detailed hypotheses. Then, the scalars αij are seen as a discrete distribution αij = P [R1 = i, R2 = j] of a pair of random variables (R1, R2) taking their values in Ω = {1, . . . , K} and the approximation (4) may be seen as a marginal distribution of

p (x1, y 1 , r1, x2, y 2 , r2) = p (r1, r2) p (x1, y 1 , x2, y 2 |r1, r2 ) .
(5) Then the main idea is to consider the stationary triplet Markov model T

N 1 = (X N 1 , R N 1 , Y N 1 )
, with R N 1 = (R1, . . . , RN ), whose distribution is defined by [START_REF] Pieczynski | Exact filtering in conditionally Markov switching hidden linear models[END_REF] and which would belong to the CGOMSM family. As specified in Remark II.2, such a model approximates [START_REF] Pieczynski | Exact filtering in conditionally Markov switching hidden linear models[END_REF] in that for any n in {1, . . . , N -1} the distribution p xn, y n , xn+1, y n+1 in the CGOMSM is as close to the distribution p xn, y n , xn+1, y n+1 in the PMM as the distribution p (x1, y 1 , x2, y 2 ) in the CGOMSM is close to the distribution p (x1, y 1 , x2, y 2 ) in the PMM.

Let us now specify what are the properties of the Gaussian distributions pij(x1, y 1 , x2, y 2 ) in (4) needed to make the corresponding approximating switching triplet model being a CGOMSM. Let us note Z n = X n , Y n and assume that for each n in {1, . . . , N -1}

p (rn+1 |xn, rn, y n ) = p (rn+1 |rn ) , (6)
which implies the Markovianity of R N 1 . This is equivalent to p (r2 |x1, r1, y 1 ) = p (r2 |r1 ) by the stationarity assumption. Besides, the above equation means that Z1 and R2 are independent conditional on R1. Therefore, we obtain p (x1, y 1 |r1, r2 ) = p (x1, y 1 |r1 ). Since (X1, R1, Y1) and (X2, R2, Y2) are identically distributed by the stationarity assumption for T N 1 , it follows that p (x2, y 2 |r1, r2 ) = p (x2, y 2 |r2 ). Thus, Gaussian distributions p (x1, x2, y 1 , y 2 |r1, r2 ) are given by the variance matrices Γ Z 1 (r1), Γ Z 2 (r2) (note that the mappings Γ Z 1 (.) and Γ Z 2 (.) are equal), and the cross-covariance matrices

Σ Z 1 Z 2 r 2 1 . Let us set Γ r 2 1 = Γ Z 1 (r1) Σ Z 1 Z 2 r 2 1 Σ Z 1 Z 2 r 2 1 Γ Z 2 (r2) , (7)
then we obtain

Z2 = a(r 2 1 )Z1 + b(r 2 1 )W1 + c(r 2 1 ), ( 8
)
where W1 is a standard Gaussian vector, independent of Z1, and

a(r 2 1 ) = Σ Z 1 Z 2 r 2 1 Γ -1 Z 1 (r1) ; (9a) b(r n+1 n )b (r n+1 n ) = Γ Z 2 (r2) -a(r 2 1 )Σ Z 1 Z 2 r 2 1 ; (9b) c(r 2 1 ) = E [Z2 |r2 ] -a(r 2 1 )E [Z1 |r1 ] . (9c)
The matrices

Γ Z 1 (r1), Σ Z 1 Z 2 r 2 1 , a(r 2 1 ), b(r 2 1
) and the vector c(r 2 1 ) may be written in the following block-form:

Γ Z 1 (r1) = Γ X 1 (r1) Σ X 1 Y 1 (r1) Σ Y 1 X 1 (r1) Γ Y 1 (r1) ; Σ Z 1 Z 2 r 2 1 = Σ X 1 X 2 r 2 1 Σ X 1 Y 2 r 2 1 Σ Y 1 X 2 r 2 1 Σ Y 1 Y 2 r 2 1 ; a(r 2 1 ) = a1(r 2 1 ) a2(r 2 1 ) a3(r 2 1 ) a4(r 2 1 ) ; b(r 2 1 ) = b1(r 2 1 ) b2(r 2 1 ) b3(r 2 1 ) b4(r 2 1 )
;

c(r 2 1 ) = c1(r 2 1 ) c2(r 2 1 )
.

Then, by stationarity, we have for n = 1, . . . , N -1

Xn+1 Yn+1 = a1(r n+1 n ) a2(r n+1 n ) a3(r n+1 n ) a4(r n+1 n ) Xn Yn + b1(r n+1 n ) b2(r n+1 n ) b3(r n+1 n ) b4(r n+1 n ) Un+1 Vn+1 + c1(r n+1 n ) c2(r n+1 n ) , (10)
where U1, V1, . . . , UN , VN are independent standard Gaussian vectors.

Definition II.1. The mixture (4) will be said "Conditionally Gaussian Observed Markov Switching Model mixture" (CGOMSM mixture) if it verifies p (r2 |x1, r1, y 1 ) = p (r2 |r1 ) and if each Gaussian pdf pij satisfies one of the two following equivalent properties:

(i) a3(r 2 1 ) = 0 for each r 2 1 ∈ Ω 2 ; (ii) Σ X 1 Y 2 r 2 1 = Σ X 1 Y 1 r 2 1 Γ -1 Y 1 (r1) Σ Y 1 Y 2 r 2 1 for each r 2 1 ∈ Ω 2 .
As explained previously, the above definition also includes the distribution p xn, y n , xn+1, y n+1 of CGOMSM for any n in {1, . . . , N -1}. Our filtering algorithm, which runs within the framework of the CGOMSM, is explained in next Section.

Definition II.2. We call "CGOMSM Approximation Based Filter" (CGOMSM-ABF) the following algorithm: (iii) use these parameters to recover the hidden signal x N 1 from y N 1 . For the step (ii), we propose to use an Expectation-Maximizationbased method, which we detail in Section four. However, there are different alternative approaches that may apply as well.

Remark II.1. In practice, the initial pdf p(x1, y 1 ) may be unspecified. However, if the system has the rapid mixing property, the residual effect of the starting distribution on p(xn, y n ) would be negligible for moderate and large values of n [31]- [START_REF] Douc | Practical drift conditions for subgeometric rates of convergence[END_REF]. Thus, in order to generate a realization z M 1 , it might be possible to begin with some auxiliary pdf p (x1, y 1 ). It is also common in practice to drop the burn-in sample z 1 , . . . , z n .

Remark II.2. The approximation (4) can be obtained with an arbitrary precision. If we denote by p the pdf of the approximation, then p (xn, y n , xn+1, y n+1 ) approximates p xn, y n , xn+1, y n+1 for any n in {1, . . . , N -1} with exactly the same precision. The readers interested by further investigations of the switching model may consult [START_REF] Lanchantin | Unsupervised segmentation of randomly switching data hidden with non-gaussian correlated noise[END_REF] to see that the marginal distribution p x N 1 , y N 1 of a CGOMSM is not necessarily Markovian (it is always not for the reversible switched linear models in [START_REF] Lanchantin | Unsupervised segmentation of randomly switching data hidden with non-gaussian correlated noise[END_REF]). However, in our use of CGOMSM its marginal distribution mimics closely the Markov property.

III. EXACT FAST OPTIMAL FILTERING IN CGOMSM MODELS

We recall in this Section how the fast optimal filter runs in a stationary CGOMSM. Let T

N 1 = (X N 1 , R N 1 , Y N 1 )
E Xn+1 y n+1 1 = r n+1 p rn+1 y n+1 1 E Xn+1 rn+1, y n+1 1 , (11)
and its variance is

Var Xn+1 y n+1 1 = E Xn+1X n+1 y n+1 1 -E Xn+1 y n+1 1 E X n+1 y n+1 1 , with E Xn+1X n+1 y n+1 1 = r n+1 p rn+1 y n+1 1 E Xn+1X n+1 rn+1, y n+1 1 . (12)

A. Some CGOMSM properties

The main property of CGOMSM is p rn+1, y n+1 |xn, rn, y n = p rn+1, y n+1 |rn, y n ,

which is straightforward from (i) in Definition II.1. There are some consequences which may be drawn from the equation above and (6). First, (R N 1 , Y N 1 ) with latent R N 1 is a hidden Markov chain (with correlated noise). Thus, p rn+1, y n+1 |rn, y n = p (rn+1 |rn ) p y n+1 r n+1 n , y n . [START_REF] Eraker | Do stock prices and volatility jump? Reconciling evidence from spot and option prices[END_REF] Second, ( 13) is equivalent to

p xn r n+1 n , y n+1 n = p (xn |rn, y n ) , (15)
which is fundamental for the conception of our fast filter. Finally, since the distribution p xn+1 xn, r n+1 n , y n+1 n is Gaussian, Xn+1 is Gaussian conditional on the pair (Rn, Rn+1) and on a linear combination of Xn, Yn and Yn+1. A similar reasoning holds for p y n+1 y n , r n+1 n and summarizing, we have

Yn+1 = D(r n+1 n )Yn + H(r n+1 n ) + Λ(r n+1 n )Vn+1;
(16a)

Xn+1 = A(r n+1 n )Xn + B(r n+1 n )Yn + C(r n+1 n )Yn+1 + F (r n+1 n ) + Π(r n+1 n )Un+1, (16b)
for some parameters

D(r n+1 n ), H(r n+1 n ), Λ(r n+1 n ), A(r n+1 n ), B(r n+1 n ), C(r n+1 n ), F (r n+1 n ), Π(r n+1 n
) and standard Gaussian vectors U1, V1, . . . , UN , VN . These parameters may be derived as follows.

Let . From [START_REF] Derrode | Fast filter in nonlinear systems with application to stochastic volatility model[END_REF], we find that the conditional mean of Xn+1, Yn+1 is

M X r n+1 M Y r n+1 + a1(r n+1 n ) a2(r n+1 n ) a3(r n+1 n ) a4(r n+1 n ) xn -M X rn y n -M Y rn = M X r n+1 + a1(r n+1 n )(xn -M X rn ) + a2(r n+1 n )(y n -M Y rn ) M Y r n+1 + a3(r n+1 n )(xn -M X rn ) + a4(r n+1 n )(y n -M Y rn ) , (17)
and that the conditional variance matrix of Xn+1, Yn+1

is b(r n+1 n )b (r n+1 n ), written in block-form as b(r n+1 n )b (r n+1 n ) = γ 1 (r n+1 n ) γ 2 (r n+1 n ) γ 2 (r n+1 n ) γ 4 (r n+1 n ) . ( 18
)
Since a3(r 2 1 ) = 0 for each r 2 1 in Ω 2 , equation (16a) holds for

D(r n+1 n ) = a4(r n+1 n ), H(r n+1 n ) = -a4(r n+1 n )M Y rn + M Y r n+1
and for some matrix

Λ(r n+1 n ) such that Λ(r n+1 n )Λ (r n+1 n ) = γ 4 (r n+1 n ). (19)
Likewise, Xn+1 is also normally distributed conditional on Xn,

R n+1 n and Y n+1 n . The conditional variance of Xn+1 is γ 1 (r n+1 n ) - γ 2 (r n+1 n )γ -1 4 (r n+1 n )γ 2 (r n+1 n
), and its conditional mean is

M X r n+1 + a1(r n+1 n )(xn -M X rn ) + a2(r n+1 n )(y n -M Y rn ) + γ 2 (r n+1 n )γ -1 4 (r n+1 n ) y n+1 - (M Y r n+1 + a3(r n+1 n )(xn -M X rn ) + a4(r n+1 n )(y n -M Y rn )
) . Term-by-term identification of (16b) with the equation above gives

C(r n+1 n ) = γ 2 (r n+1 n )γ -1 4 (r n+1 n ) A(r n+1 n ) = a1(r n+1 n ) -C(r n+1 n )a3(r n+1 n ) B(r n+1 n ) = a2(r n+1 n ) -C(r n+1 n )a4(r n+1 n ) F (r n+1 n ) = M X r n+1 -A(r n+1 n )M X rn -B(r n+1 n )M Y rn -C(r n+1 n )M Y r n+1
and Π(r n+1 n

) is a matrix such that

Π(r n+1 n )Π (r n+1 n ) = γ 1 (r n+1 n ) -C(r n+1 n )γ 2 (r n+1 n ). (20)

B. Filtering in the CGOMSM

Let us recall how fast filter runs in a stationary CGOMSM. We sequentially compute p rn+1 y n+1 .

), H(r n+1 n ), Λ(r n+1 n ), A(r n+1 n ), B(r n+1 n ), C(r n+1 n ), F (r n+1 n ), Π(r n+1 n ) as follows. Since (R N 1 , Y N 1 )
(22) Then, E Xn+1 rn+1, y n+1 1 = rn p rn rn+1, y n+1 1 A(r n+1 n )E [Xn |rn, y n 1 ] + L(r n+1 n , y n+1 n ) , (23)
with

L(r n+1 n , y n+1 n ) = B(r n+1 n )y n + C(r n+1 n )y n+1 + F (r n+1 n
). We obtain the equation above by taking the conditional expectation of (16b) and this result is mainly due to [START_REF] Kim | State-space models with regime switching[END_REF] and, consequently, to the fact that

E Xn r n+1 n , y n+1 1 = E [Xn |rn, y n 1 ] . (24)
A similar approach applies to the computation of E Xn+1X n+1 rn+1, y n+1 1 . It may be shown that

E Xn+1X n+1 r n+1 n , y n+1 1 = A(r n+1 n )E [XnX n |rn, y n 1 ] A (r n+1 n ) + A(r n+1 n )E [Xn |rn, y n 1 ] L (r n+1 n , y n+1 n ) + L(r n+1 n , y n+1 n )E [X n |rn, y n 1 ] A (r n+1 n ) + L(r n+1 n , y n+1 n )L (r n+1 n , y n+1 n ) + Π(r n+1 n )Π (r n+1 n ) (25) Then E Xn+1X n+1 rn+1, y n+1 1 = rn p rn rn+1, y n+1 1 E Xn+1X n+1 r n+1 n , y n+1 1 . (26)
Finally, our fast optimal filter is: Algorithm III.1.

IV. APPROXIMATING A NON-LINEAR NON-GAUSSIAN SYSTEM WITH THE CGOMSM

We recall that according to our CGOMSM-ABF principle, we contemplate estimation of the CGOMSM parameters through a signal (x M 1 , y M 1 ) which we sample by using (1), [START_REF] Arulampalam | A tutorial on particle filters for online non-linear/non-Gaussian Bayesian tracking[END_REF] or any other timeinvariant system cf. Remark IV.1.

To this end, we use a common algorithm known as "Expectation-Maximization" [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] (EM). Let us then parametrize the CGOMSM by

Θ = µ i , Γ i , p ij , A ij , B ij , C ij , D ij , F ij , H ij , Π ij , Λ ij (i, j) ∈ Ω 2
, where:

• for each i in Ω, µ i and Γ i define the pdf p (x1, y 1 |r1 ); [START_REF] Ristic | Beyond the Kalman filter: particle filters for tracking applications[END_REF].

• for each i, j in Ω, p ij = P (r1 = i, r2 = j) and A ij , B ij , C ij , D ij , F ij , H ij , Π ij , Λ ij are described in
To make the text easier to read, we will denote the artificial signal by

(x M 1 , y M 1 ) instead of (x M 1 , y M 1 )
. Thus, y M 1 stands for the presampled data and should not be confused with y N 1 , which is the noisy input signal.

The crux point is that Xn+1 and Yn+1 are Gaussian conditional on the pair (Rn, Rn+1) and on a linear combination of Xn, Yn so an appropriate linear regression per each pair (Rn, Rn+1) enables to determines the coefficients of the combination and to solve the Mstep of the EM algorithm. We detail below the EM algorithm suited to our context (cf. section II).

Algorithm IV.1. EM estimation (a) Make an initial guess

Θ 0 = µ 0 i , Γ 0 i , p 0 ij , A 0 ij , B 0 ij , C 0 ij , D 0 ij , F 0 ij , H 0 ij , Π 0 ij , Λ 0 ij (i, j)
∈ Ω 2 as follows:

(i) Apply the K-means clustering method to x M 1 . We will denote by κn(i) the function which assigns 1 if xn is within the i th cluster, and 0 otherwise. We also note δn(i, j) = κn(i)κn+1(j); (ii) For each i in Ω, µ 0 i and Γ 0 i are given by

µ 0 i = M n=1 znκn(i) M n=1 κn(i) ; (27a)
Γ 0 i = M n=1 zn -µ 0 i zn -µ 0 i κn(i) M n=1 κn(i) , (27b)
where z n = [x n , y n ], and for each (i, j) in Ω 2 , p 0 ij is given by

p 0 ij = 1 M -1 M -1 n=1 δn(i, j). ( 28
) (iii) Compute intermediate matrices E 0 ij , S 0 ij , χ 0 ij , Φ 0 ij , G 0 ij , P 0
ij , ξ 0 ij and T 0 ij as follows:

E 0 ij = 1 p 0 ij M -1 n=1 zn y n+1 δn(i, j); (29a) S 0 ij = 1 p 0 ij M -1 n=1
znz n zny n+1 y n+1 z n y n+1 y n+1 δn(i, j); (29b)

χ 0 ij = 1 p 0 ij M -1 n=1
xn+1z n xn+1y n+1 δn(i, j); (29c) For our experiments, we decided to reduce the number of terms as follows: when ξn becomes greater than r, we keep the r mixands in ξn j=1 αnjN (xn; xnj, Pnj) which have the greatest weight coefficients αnj, and we discard the remaining. Therefore, we impose the constraint that ξn = r. We found out empirically that GSF behaves asymptotically for r ≥ 3, but does not attain the optimal MSE. We note that since the model ( 36) is linear, there is no reason for considering the extensions of the GSF for non-linear systems, such as Gaussian Sum Unscented Kalman Filter [START_REF] Šimandl | Sigma point Gaussian sum filter design using square root unscented filters[END_REF] or estimates within the ASV model, whose parameters are given in [START_REF]SV mixture models with application to S&P 500 index returns[END_REF] and reported in Table IV. Our result is shown in Figure 3. We find that the volatility estimates produced by our algorithm are consistent with the log-return process: as we can see in Figure 3, the intervals where the fluctuation of log-returns are low (e.g. between 1991 and 1995) match the intervals where the log-volatility is low, and vice versa. Moreover, we calculated the mean squared distance between our volatility estimates and those of the particle filter, and we find that this distance is negligible compared to the variance of the log-volatility process. Furthermore, when the number K of classes in our algorithm increases, this distance decreases as shown in Table V.

VI. CONCLUSION

We proposed a method to estimate the latent variables in non-linear and non-Gaussian systems, called CGOMSM-ABF for 'conditionally Gaussian observed Markov switching model -Approximation Based Filter'. It makes use of a single global approximation of the system. Our method is very general and has several advantages over existing techniques. Its performance has been examined on synthetic samples of the standard and asymmetric stochastic volatility, as well as on real data (SPX index returns). Indeed, we found that our method attains the asymptotic performances of the particle filter, what could not be obtained with the Gaussian Sum Filter [START_REF] Sorenson | Recursive Bayesian estimation using Gaussian sums[END_REF], [START_REF] Šimandl | Filtering, prediction and smoothing with Gaussian sum representation[END_REF].

The filtering procedure, which is the object of the paper, is applicable in general stationary (or asymptotically stationary) Markov dynamic systems, provided that one can sample its realizations. It is as fast as the standard Kalman filter, however, one has to adjust the filter to a particular model by using, for instance, the EM-based algorithm that we have supplied. Our future prospects include filtering for various stochastic volatility models, including those with high dimension of the state space, i.e. multifactor volatility models. Another important topic would be the unsupervised estimation of parameters of the initial PMM model in eq. ( 2) [START_REF] Cappé | Online EM algorithm for hidden Markov models[END_REF], [START_REF] Ozkan | Recursive maximum likelihood identification of jump Markov nonlinear system[END_REF].

Ivan Gorynin received the M.Sc. degree in statistics from Telecom SudParis, France, in 2014. He is currently a PhD student within the CITI Department, at Telecom SudParis. His research interests include term structure and stochastic volatility modelling, financial Econometrics, machine learning, Big Data, non-linear time series analysis, mathematical statistics and hidden Markov models.

(i) generate a realization z M 1 =

 1 (z 1 , . . . , z M ) from the prior model (2); (ii) infer the CGOMSM parameters based on the observation z M 1 assuming that the chain of switches R N 1 is missing;

1 ,E

 1 Xn+1 rn+1, y n+1 1 and E Xn+1X n+1 rn+1, y n+1 1 using p (rn |y n 1 ), E [Xn |rn, y n 1 ], E [XnX n |rn, y n 1 ] and parameters D(r n+1 n

(i) Use p (rn |y n 1 ) and y n+1 to obtain p rn+1 y n+1 1 and p rn rn+1, y n+1 1 by

 111 [START_REF] Straka | Gaussian sum unscented Kalman filter with adaptive scaling parameters[END_REF],[START_REF] Feng | Real-time state estimator without noise covariance matrices knowledge-fast minimum norm filtering algorithm[END_REF]; (ii) Use E [Xn |rn, y n 1 ], p rn rn+1, y n+1 1 and y n+1 to obtain E Xn+1 rn+1, y n+1 1 by (23); (iii) Use E [XnX n |rn, y n 1 ], E [Xn |rn, y n 1 ] and y n+1 to compute E Xn+1X n+1 r n+1 n , y n+1 1 by (25). Then use E Xn+1X n+1 r n+1 n , y n+1 1 and p rn rn+1, y n+1 1 to obtain E Xn+1X n+1 rn+1, y n+1 1 by (26).

Fig. 1 .

 1 Fig. 1. Up: simulated log-volatility trajectory with a SV model (red, plain), simulated log-returns (black, dotted). Down: log-volatility estimates computed using K = 2 classes (blue, dotted), and K = 5 classes (green, dashed).

Fig. 3 .

 3 Fig. 3. Trajectories of the SPX index log-returns (down) and log-volatility estimates (up). The x-axis represents the dates for both trajectories, the y-axis labelling on the left concerns the log-volatility values, and the y-axis labelling on the right is related to the log-return values.

Stephane

  Derrode served as an Associate Profesin Signal and Image processing first at École Nationale Supérieure de Physique de Marseille in 2001 and then at École Centrale de Marseille in 2005, France. He joined École Centrale de Lyon, France, in 2014 where he is currently Professor in Computer Sciences within the Mathematics and IT Department. His research interests include nonlinear, non-Gaussian state-space systems for time series analysis, with applications to UAV, object tracking in video and dynamic biometry. Emmanuel Monfrini is an Associate Professor at T à c l à c com SudParis, France, since 2008. He received his PhD thesis in mathematical statistics from Paris VI University in 2002. His research interest is centered around statistical modeling, with a particular interest for Markovian models and their application in several fields including signal and image processing. Wojciech Pieczynski received the Doctorat d' État degree in mathematical statistics from the Université Pierre et Marie Curie, Paris, France, in 1986. He is currently Professor at the Telecom SudParis and Head of the Communications, Images, and Information Processing (CITI) department. His research interests include mathematical statistics, stochastic processes, and statistical image processing.

  It follows from (10) that Xn+1, Yn+1 is normally distributed conditional on Xn, Yn and R n+1 n

us set M X rn = E [Xn |rn ] and M Y rn = E [Yn |rn ].

  is a pairwise Markov model,

				p rn+1, y n+1 |rn, y n p (rn |y n 1 )
	p rn+1 y n+1 1	=	rn	1 ) p r n+1 , y n+1 |rn, y n p (rn |y n	;
			rn,r n+1
				(21)
	p rn rn+1, y n+1 1	=	p rn+1, y n+1 |rn, y n p (rn |y n 1 ) p rn+1, y n+1 |r
			r n

n , y n p (r n |y n 1 )

TABLE IV THE

 IV PARAMETERS OF THE ASV MODEL FOR THE STOCK MARKET DATA.

	Param.	µr	µ	φ	σ	ρ	β
	Value	7 • 10 -5	-9.54	0.98	0.17	-0.43	1.00

TABLE V MEAN

 V SQUARE DISTANCES BETWEEN OUR VOLATILITY ESTIMATES AND THOSE FROM THE PF, WITH DIFFERENT NUMBER OF CLASSES.

xn+1δn(i, j); (29d)

y n δn(i, j); (29e)

y n y n δn(i, j); (29f)

y n+1 y n δn(i, j); (29g)

y n+1 δn(i, j). (29h)

and Λ 0 ij are given by

y n+1 y n+1 δn(i, j)-

xn+1x n+1 δn(i, j)-

(b) Find the new set of parameters Θ q+1 as follows:

(i) For each i in Ω, compute the posterior probabilities φ q n (i) = p rn = i x M 1 , y M 1 , Θ q , and for each i, j in Ω compute ψ q n (i, j) = p rn = i, rn+1 = j x M 1 , y M 1 , Θ q (computational details are provided below); (ii) For each i in Ω, compute µ q+1 i and Γ q+1 i by substitution φ q n (i) for κn(i) in ( 27); (iii) For each i, j in Ω, p q+1 ij is given by

Then compute intermediate matrices

Details on the derivation of Algorithm IV.1 will be made available upon request to the corresponding author.

Let us recall the formulas for φ q n (i) and ψ q n (i, j). Let us set tn = (xn, rn, y n ), αn(rn) = p (rn, z n 1 ) and βn(rn) = p z N n+1 |tn . Then, the forward-backward algorithm computes recursively αn(rn) and βn(rn) as follows:

Thus,

with

Finally, the CGOMSM-ABF implementation, which we use in experiments, runs as follows: Here we use one to approximate a non-linear non-Gaussian system. Another possible situation is that the system is linear Gaussian, and there are random switches. Although quite different, these contexts can be dealt with using the same model, as the CGOMSM used in this paper. Of course, there also exists non-Gaussian switching systems [START_REF] Hamilton | Autoregressive conditional heteroskedasticity and changes in regime[END_REF], [START_REF] Shumway | Dynamic linear models with switching[END_REF].

V. EXPERIMENTS ON STOCHASTIC VOLATILITY MODELS

A stochastic volatility model describes the variance of a stochastic process. Examples of stochastic volatility models include the Heston model, the GARCH model and many more [START_REF] Kim | State-space models with regime switching[END_REF], [START_REF] Bollerslev | Generalized autoregressive conditional heteroskedasticity[END_REF], [START_REF] Engle | Autoregressive conditional heteroskedasticity with estimates of the variance of U.K. inflation[END_REF]. In this Section, we provide two series of Monte-Carlo experiments and an example of estimating the volatility of a series of stock returns. In the first series, we focus on a standard Stochastic Volatility (SV) model [START_REF] Ghysels | Stochastic volatility[END_REF]- [START_REF] Taylor | Modeling Financial Times Series[END_REF], that is usually presented as follows:

(33a)

where U1, V1, . . . , UN , VN are independent standard Gaussian vectors and µ, φ, β, σ are fixed. Since this model is of the form (1), it is an HMM. The second series of experiments is devoted to a recent extension, known as the asymmetric SV (ASV) model [START_REF] Nikolaev | Nonlinear filtering of asymmetric stochastic volatility models and value-at-risk estimation[END_REF]- [START_REF] Takashi | Estimating stochastic volatility models using daily returns and realized volatility simultaneously[END_REF], which may be presented as follows:

+ λUn+1 ; (34b)

Since both SV and ASV are generative models, they fit well to our CGOMSM-ABF approach. If x(y n 1 ) is a filtered signal obtained from y N 1 , and x N 1 is the "ground-truth", then the mean squared error (MSE) of the filter can be estimated by

which is a useful performance criterion for comparing filtering effectiveness. All the results presented below are means of 100 equivalent independent experiments, each of them being computed using N = 1000 simulated data points.

A. SV model

Here we compare the performance of our method with respect to the Particle Filtering (PF) and the Gaussian Sum Filtering (GSF), using data simulated by different stationary SV models. We set µ = 0.5, β = 0.5, and consider four different cases for φ and σ such that φ 2 + σ 2 = 1 (that is to ensure that the common variance of the variables Xn is unitary). The results are reported in Table I.

Here are the details of each filtering method used in the experiments

• For our method, we try out different values of K (cf. Remark V.1) and we infer the CGOMSM from some independently generated sample (x M 1 , y M 1 ) of size M = 20000, performing 100 EM iterations. See Figure 1 for an example of trajectories.

• The PF implementation uses 1500 particles. We found out empirically that PF behaves asymptotically for this particle number.

• In order to use the GSF, let us linearise the SV model by taking the log of both sides of eq. (33c), to get

where

are standard Gaussians. Then, the common pdf of V n is p (v) = exp v 2 N exp v 2 ; 0, 1 , which we approximate by a Gaussian mixture of r components using the EM algorithm: p (v n ) ≈ r m=1 γnN (v n ; v m , Rm). We found that when r ≥ 5, the approximation is accurate enough to achieve a (µ = 0.5, β = 0.5, AND σ 2 = 0.75), FOR φ = 0.5. (µ = 0.5, β = 0.5, AND σ 2 = 0.75), FOR φ = 0.8. Gaussian Sum Unscented Kalman Filter with adaptive scaling parameters [START_REF] Straka | Gaussian sum unscented Kalman filter with adaptive scaling parameters[END_REF]. Contrary to our filtering method, which makes use of a single global approximation, the GSUKF relies on multiple approximations:

• an approximation of the noise terms with a Gaussian mixture;

• some reduction technique to keep the number of mixands of the filtering pdf within reasonable bounds. Additionally, when the model is non-linear, the GSUKF uses the Unscented Transform (UT) for computing the approximate means and covariances. And the UT relies, in turn, on its scaling parameters. Our experiments show that computing a single global approximation may be advantageous and helps to avoid the residual cumulative effects. However, unlike our method, the GSUKF may be used for non-stationary systems.

B. Asymmetric SV model

Regarding the ASV model, we compare the performance between our method and PF only, since the GSF and GSUKF ignore the impact of the volatility asymmetry coefficient ρ and therefore are not suitable for this model (recalling that the ASV is not an HMM). The experimental configuration is identical to the previous one. For the sake of consistency with the asymmetric volatility phenomenon [START_REF] Black | Studies of stock price volatility changes[END_REF], ρ should be assumed negative.

We set µ = 0.5, β = 0.5, σ 2 = 0.75, and consider five different cases for ρ and λ such that ρ 2 + λ 2 = 1 (to ensure that the common variance of the variables Xn is unitary). The results are reported in Table II for φ = 0.5 and in Table III for φ = 0.8. Figure 2 shows an ASV trajectory, and its restoration with the CGOMSM-ABF algorithm for K = 2 and K = 5 classes.

According to our results, the approach we propose is efficient for both SV and ASV models, and attains the same asymptotic performances as the PF. Regarding the processing time, we find that:

• after having it adjusted to the SV model, our filter is nearly five times faster than the PF; • the model approximation algorithm may take time depending on the system complexity.

Remark V.1. At the moment, we have no computational technique to select the minimum number of classes allowing to obtain asymptotic performances. We note only the trade-off between the computational cost and the variance of the resulting estimates. Indeed, with a greater number of classes the former increases, while the latter decreases.

In practice, five classes would be enough for most situations.

C. Recovering volatility from stock market data

Let us remind that if Pn-1 denotes the stock price at the beginning of the previous trading day and if Pn denotes the stock price at the beginning of the current trading day, then :

is the current daily return on the stock investment;

• un = log (1 + Rn) = log Pn P n-1 is the continuously compounded daily return. It is also often called the log-return. To see why un is called the continuously compounded return, take the exponential of both sides to get exp (un) = Pn P n-1 . Rearranging, we get Pn = Pn-1 exp (un) so that un is the continuously compounded growth rate in prices between the beginning of the previous and the current trading days. This has to be contrasted with Rn, which is the simple growth rate in prices Pn-1 and Pn without any compounding.

Following [START_REF] Eraker | The impact of jumps in volatility and returns[END_REF]- [START_REF]SV mixture models with application to S&P 500 index returns[END_REF] to examine the performance of our method on the stock market data, we compute the log-returns un over daily SPX index data from Jun. [START_REF] Feng | Kalman filter with recursive covariance estimation-sequentially estimating process noise covariance[END_REF]1980 to Aug. 30, 2002 (N = 5604), then we calculate yn = u n -µr, where µr is given in [START_REF]SV mixture models with application to S&P 500 index returns[END_REF] and u n denotes pre-processed log-return with an ARMA(2,1) model [START_REF] Durham | Monte Carlo methods for estimating, smoothing, and filtering one-and two-factor stochastic volatility models[END_REF]- [START_REF] Andersen | An empirical investigation of continuous-time equity return models[END_REF]. Next, we use our algorithm to compute the filtered volatility


Zied Bouyahia
email: bouyahiazied@gmail.com

Stéphane Derrode
email: stephane.derrode@ec-lyon.fr

Wojciech Pieczynski
email: wojciech.pieczynski@telecom-sudparis.eu

An exact smoother in a fuzzy jump Markov switching model

In this paper, we address the exact smoothing problem of Conditionally Gaussian Observed Markov Switching Model (CGOMSM). The proposed approach tackles the discontinuity feature in switching regime models by incorporating fuzzy switches instead of hard jumps. Fuzzy switched based approach is more adapted to real-world application in which regime continuity is an intrinsic property. We define, with respect to a measure, the local density of the switches and we show how fast smoothing equations should be adapted to cope with the fuzzy model. Finally, we show through several experiments the interest of the fuzzy switches model.

Introduction

Let X N 1 = {X 1 , . . . , X N }, Y N 1 = {Y 1 , . . . , Y N } and R N 1 = {R 1 , . . . , R N } be three random sequences taking values in R m , R q and Ω = 1, . . . , K. Let X N 1 be a hidden process and Y N 1 be an observed process. We consider a switching regime model represented by the sequence of switches R N 1 . We address the smoothing problem consisting in an recursive search of the unobserved process X N 1 and the switches sequence R N 1 knowing the observed sequence Y N 1 . A fast Bayesian processing can be carried out by assuming that the distribution of (X N 1 , Y N 1 ) is within the framework of hidden Gaussian Markov Model. The non-linearity can be modeled by a switching regime system. Then, the idea is to approximate a non-linear non-Gaussian system by a regime switching Gaussian system. Some recent switching models have been proposed with efficient fast exact filtering schemes [START_REF] Gorynin | Fast filtering in switching approximations of non-linear Markov switching systems with application to stochastic volatility[END_REF] and satisfactiry computation time. These switching Gaussian models include conditionally Markov switching hidden linear models (CMSHLM) [START_REF] Pieczynski | Exact filtering in conditionally Markov switching hidden linear models[END_REF] and the conditionally Gaussian observed Markov switching model (CGOMSM) which is a sub-model of CMSHLM defined as follows: p(r n+1 |x n , r n , y n ) = p(r n+1 |r n ) and 1 depicts the dependence graph of the CGOMSM model. We assume that R N 1 takes its values in a discrete finite set of K switches Ω = {1, ..., K}. This hard jumps model has been widely used in several contexts dealing with switching regime Markov systems. Its success comes from its ability to represent non-linear dynamic patterns which is an inherent property in several applications (analysis of economic and finance time series [START_REF] Koko | Application of Markov-Switching Model to Stock Returns Analysis[END_REF], sustainable energy [START_REF] Yang | Support Vector Machine Enhanced Markov Model for Short Term Wind Power Forecast[END_REF], robotics [START_REF] Baltzakis | A Hybrid Framework for Mobile Robot Localization: Formulation Using Switching State-Space Models[END_REF] [7] [START_REF] Le Corff | Online Expectation Maximization algorithm to solve the SLAM problem[END_REF], etc.). However, this model does not take into account the intrinsic imprecision of the switches in real-world applications. In fact, hard jumps induce discontinuity in the dynamic behavior of the studied system. This transitory imprecision can be modeling fuzzy modeling which consists in allowing each switch to take its value as a mixture of many components simultaneously. Fuzzy modelling has been already incorporated in several applications dealing with Markov models [START_REF] Carincotte | Unsupervised image segmentation based on a new fuzzy hidden Markov chain model[END_REF][START_REF] Salzenstein | Parameter estimation in hidden fuzzy Markov random fields and image segmentation[END_REF][START_REF] Carincotte | Unsupervised Change Detection on SAR Images using Fuzzy Hidden Markov Chains[END_REF]. In this paper, we present a new method to approximate non-linear Markov systems using a new variant of CGOMSM using fuzzy switches (hereafter called CGOMFSM). The remaining of this paper is organized as follows. In the second section, we detail the formulation of fuzzy switching model. The third section describes the adaptation of CGOMSM algorithms for parameters estimation and for posterior marginal probabilities computation the fuzzy counterpart. The fourth section presents experimental results, and the last one draws conclusions and future work.

-T N 1 = (X N 1 , R N 1 , Y N 1 ) is Markov chain; X1 Y1 R1 X2 Y2 R2 XN YN RN XN-1 YN-1 RN-1
X n+1 Y n+1 = A xx n+1 (R n+1 n ) A xy n+1 (R n+1 n ) 0 A yy n+1 (R n+1 n ) X n Y n + B xx n+1 (R n+1 n ) B xy n+1 (R n+1 n ) B yx n+1 (R n+1 n ) B yy n+1 (R n+1 n ) U n+1 V n+1 + M X (R n+1 ) -A xx n+1 (R n+1 n )M X (R n ) -A xy n+1 (R n+1 n )M Y (R n ) M Y (R n+1 ) -A yy n+1 (R n+1 n )M Y (R n ) with U N 1 and V N 1 a Gaussian unit-variance white noise vector, M X (R n ) and M Y (R n ) are respective means of X N 1 and Y N 1 in each state (independently from n). The CGOMSM is then defined by matrices A(R n+1 n ), B(R n+1 n ), transition matrix denoted as t such that t(i, j) = p(r n+1 = j|r n = i) and vector M = [M (R n ); M (R n+1 )] such that M (R n ) = [M X (R n ); M Y (R n )]. Figure

Fuzzy switching model with K hard classes

In the fuzzy switches system we assume that each jump r n in the random process

R N 1 is a mixture of the K classes. Let ε i n be the contribution of the i th hard component ω i such that ε i n ∈ [0, 1].
Then each switch r n is given by the vector ε n :

r n = ε n = (ε 1 n , ε 2 n , ..., ε K n ), (1)
The normalization condition yields:

K k=1 ε k = 1 (2)
In the remainder, for the sake of simplicity, we refer to each hard component by its index in the discrete set i.e. class i corresponds to the hard component ω i .

Measure of local density

The distribution of each switch is defined by a density h r with respect to a Lebesgue measure υ. This measure includes K discrete components represented by K Dirac functions {δ k, 1≤k≤K } as well as the continuous components representing the combination of the hard switches.

Let S be a subset of 0 < l ≤ K hard components in Ω. If a switch r n is a mixture of the hard components in S ⊂ Ω then r n belongs to the hyperplane of dimension l corresponding to

l i=1 ε i n = 1 such that ε i n > 0, ∀1 ≤ i ≤ l.
The corresponding measure υ is given by :

υ = K k=1 δ k + S⊆Ω,S =∅ µ S (3)

Example

In this example we consider the particular case K = 3 to illustrate the local density measure. The three dimensional space corresponding to the three hard components is depicted in figure 2.1. The hyperplane corresponding to the normalization condition (eq. ( 2)) coincides with the shaded triangle ∆. The edges of ∆ represent the mixture of two hard classes while the vertices of ∆ correspond to the three hard components. When K = 3, eq. ( 3) becomes:

υ = δ ω1 + δ ω2 + δ ω3 + µ ω1ω2 + µ ω1ω3 + µ ω2ω3 + µ ω1ω2ω3

Joint densities

Let H be the hyperplane of dimension K -1 and f (r n , r n+1 ) be the joint density of the pair (r n , r n+1 ) ∈ H × H. The joint density f is defined as follows:

f (r n , r n+1 ) = α ij if both switches are hard βφ(r n , r n+1 ) + θ otherwise . ( 4
) We choose φ(r n , r n+1 ) = 1 -rn+1-rn √ 2 r
, r ∈ R with r n+1 -r n is the distance between two consecutive switches given by the quadratic norm. Then:

φ(r n , r n+1 ) =       1 -     K i=1 (ε i n -ε i n+1 ) 2 2     1 2       r , r ∈ R (5)
The normalization condition yields:

K i=1 α ii + H H θ + βφ(r n , r n+1 )(ν, µ)dνdµ = 1 (6)
Example:

When the number of hard switches equals 3, the expression of normalization condition gives:

3 i=1 α ii + (β a + β b + β c ) 1 0 1 0 (1 -|u -v|) r dudv + 2(β a + β b + β c ) 1 0 (1 -|ε|) r + (1 -|1 -ε|) r dε = 1 (7)

Parameters interpolation

The model matrices of the fuzzy switching model can be calculated by linear interpolation using the following formula:

A(ε 1 , ε 2 ) = 1≤i,j≤K ε i 1 ε j 2 A(i, j), B(ε 1 , ε 2 ) = 1≤i,j≤K ε i 1 ε j 2 B(i, j), M (ε 1 , ε 2 ) = 1≤i,j≤K ε i 1 ε j 2 M (i, j),
where A(i, j) and B(i, j) are the model matrices corresponding to the hard components i and j. M (i, j) is means model vector for hard switches i and j.

The implementation of the fuzzy switching model can be performed by adequate quantification of the interval [0, 1] into F discrete fuzzy levels. The larger F is, the more accurate the representation of data would be. However, choosing a large number of fuzzy levels will lead to high computation time. For example, when the number of crisp components equals three, setting F = 3 yields 15 switches and setting F = 4 gives 21 switches.

Fuzzy switching model with two hard components

In this remaining of the paper, we consider the case of two hard switches Ω = {0, 1}. To model fuzzy switches, we consider that each random variable R n in R N

1 takes its values in the continuous interval [0, 1], instead of the set {0, 1}. Let us denote the pair (ε 0 n , ε 1 n ) ∈ [0, 1], in which ε i n represents the contribution of the hard component i to the switch r n . Without loss of generality, let , where notation R n+1 n represents the couple (R n , R n+1 ). p R n+1 n is defined with respect to the measure product ν ⊗ ν, under normalisation condition

ε n = ε 1 n = 1 -ε 0 n . Then we have R n = ε n : -ε n = 0 if the switch is the hard component 0. -ε n ∈]0, 1[ if the switch is fuzzy. -ε n = 1 if
[0,1] 2 p (u, v) d(ν ⊗ ν)(u, v) = 1. ( 8
)
This distribution can be represented by a density h1 which involves three types of components:

-Hard components when both R n and R n+1 are in {0, 1}; -Composite components when either R n and R n+1 is in {0, 1} and the other is in ]0, 1[; -Fuzzy components when both switches are in fuzzy area ]0, 1[ 2 .

Let assume that the density h is uniform in each of its components. It is so possible to define the density h by a set of 9 parameters explicitly expressing

h(r n , r n+1 ) = p R n+1 n = r n+1 n = π ij . But since we have π 0F = h(0, u) = h(u, 0), (9)
π 1F = h(1, u) = h(u, 1), (10)
with u ∈]0, 1[, only 7 parameters are required to entirely define h : π 00 , π 01 , π 10 , π 11 , π 0F , π 1F and π F (see Fig. 3). With these notations, the normalisation eq. ( 8) writes

π 00 + π 01 + π 10 + π 11 + 2π 0F + 2π 1F + π F = π 00 + π 01 + π 10 + π 11 + 2 ]0,1[ p(0, u)du + 2 ]0,1[ p(1, u)du + ]0,1[ 2 p(u, v)dudv = 1 (11)
Let now assume that the shape of h is defined according to the simple following shape

h(ε 1 , ε 2 ) = aΦ(ε 1 , ε 2 ) + b, (12)
whatever with r ∈ R. Hence, this model for h relies 3 parameters: a, b, and r. Analytic computation of the joint prior densities can be worked out by quantifying the interval [0, 1] into F equal-length sub-intervals i F , i+1 F as described in Figure 3. Using this scheme, the normalization condition in eq. ( 11) yields:

ε 1 , ε 2 ∈ [0, 1], with (a, b) ∈ R 2 and Φ(ε 1 , ε 2 ) = (1 -|ε 1 -ε 2 |) r , (13)
π 00 + π 01 + π 10 + π 11 + 1 2F F -1 i=0 (a(1 -ε i ) r + b) + 1 2F F -1 i=0 (aε r i + b) + 1 2F 2 F -1 i=0 F -1 j=0 (a (1 -|ε i -ε j |) r + b) = 1. ( 14
)
Each sub-interval can be represented by its medium value 2i+1 2F . So, in this discrete approximate scheme, the joint a priori density can be defined by a (2 + F ) × (2 + F ) matrix. Under the assumptions of fuzzy switches, we can define the matrices of the incorporated model using a bi-linear function as follows:

A(ε 1 , ε 2 ) = [(1 -ε 1 )A(0, 0) + ε 1 A(1, 0)](1 -ε 2 ) +[(1 -ε 1 )A(0, 1) + ε 1 A(1, 1)]ε 2 (15)
B(ε 1 , ε 2 ) = [(1 -ε 1 )B(0, 0) + ε 1 B(1, 0)](1 -ε 2 ) +[(1 -ε 1 )B(0, 1) + ε 1 B(1, 1)]ε 2 (16)
The means vectors of the fuzzy model are calculated using the following equations: Hence, the CGOMFSM is entirely defined by the parameters of the corresponding deterministic hard switching model, the number of fuzzy levels F , and parameters a, b and r.

M (ε i ) = [(1 -ε i )M (0) + ε i M (1)] (17)
The parameter r specifies the homogeneity of the switching model. The larger r is, the larger the probability of having two similar consecutive switches is. Figure 5 represents and example of simulation of (X,Y,R) using the set of parameters of a fuzzy switching model defined in table 1. Simulations were performed using the following transition matrix:

t =          
0.99 0.01 0 0 0 0 0 0 0.99 0.01 0 0 0 0 0 0 0.99 0.01 0 0 0 0 0 0 0.99 0.01 0 0 0 0 0 0 0.99 0.01 0 0 0 0 0 0 0.99 0.01 0 0 0 0 0 0 1.00

         
This simulation shows the imprecision between hard switches 0 and 1 as presented in the trajectory of Y N 1 . The choice of the transition matrix allow a progressive regime switching from parameters set corresponding to hard switch 0 to hard switch 1. 1 ] and the model parameters using the procedure detailed in [START_REF] Gorynin | Fast filtering in switching approximations of non-linear Markov switching systems with application to stochastic volatility[END_REF]. The main difference between CGOMSM and CGOMFSM is that in the case of fuzzy switches we involve continuous integration, requiring to be quantified with respect to the number of discrete fuzzy levels F .

Fig. 1 .N 1 and R N 1 .

 111 Fig. 1. Directed graph representing dependencies between random sequences X N 1 Y N 1 and R N 1 . Circles represent continuous process and diamond represents discrete process.

Fig. 2 .

 2 Fig. 2. Example of measure density for a fuzzy model with 3 discrete classes.

Fig. 3 .

 3 Fig. 3. Density of p R n+1 n

Fig. 4 .

 4 Fig. 4. Subdivision of the interval [0, 1] into F = 5 equal-length fuzzy sub-intervals.

RFig. 5 .

 5 Fig. 5. Trajectories of simulated CGOMFSM (X N 1 , Y N 1 , R N 1 ).

Table 1 .4E X n+1 y N 1 =p r n+1 = ν y N 1 E 1 and E X n+1 r n+1 , y n+1 1 . 1 from p r n y N 1 and

 1111111 Example of fuzzy switching model with 5 fuzzy switches. Fast smoothing in a CGOMFSMLet us denote byT N 1 the triplet (X N 1 , R N 1 , Y N 1 ). The smoothing problem consists in computing : X n+1 r n+1 = ν, y n+1 1 dν, (19) from p r n+1 y N The optimal smoother computes recursively p r n+1 y N 1 and E X n+1 r n+1 , y n+1 E [X n |r n , y n

  the switch is the hard component 1.

	So this model able to represent signals with both discrete (hard) and continuous
	(fuzzy) components. Let ν define the measure associated to random variable R n
	(ν is defined by 2 diracs and a Lebesgue measure).
	Let us now precisely define the joint a priori density p R n+1 n

Density h does not depend on index n since the model we consider in this paper is assumed stationary.

Since E X n r n+1 n , y n+1

and from (19), we can derive the following recursive equation:

with F n+1 (r

and

βn+1(υ)p(rn+1, y n+1 |rn, y n )dυ (25)

Using forward-backward probabilities, we can compute the smoothed and the filtered probabilities as follows::

Posterior marginal probabilities are calculated using the normalized Baum-Welch algorithm. The algorithm computes recursively the forward and backward probabilities. In the case of fuzzy switches, these probabilities are defined as follows:

with t n (θ) = (x n , y n , r n = θ). Then:

Experiments

In this section, we present two series of experiments to illustrate the smoother, in the case of scalar data (m = q = 1). In the first series we assess the performance of the fuzzy model with synthetic fuzzy signals; in the second series we apply our algorithm to smooth simulated Stochastic Volatility (SV) data. In both experiments, parameters estimation is carried out using EM algorithm using training samples denoted by (x T 1 , y T 1 ) of size T . Then we repeatedly generate, according to the considered model, synthetic sequences of size S denoted by (x S 1 , y S 1 ). Smoothing algorithm is then performed using estimated parameters to generate xS 1 from the observed sequence y S 1 . The criterion used to assess the efficiency of smoothing algorithms is the mean squared error (MSE) defined as follows:

Smoothing synthetic fuzzy signals

We generate a fuzzy signal with 5 discrete fuzzy switches. Then we estimate the model considering different values of F ranging from 1 to 5. For each set of data, we consider 3 cases for the value of r ∈ {2, 5, 20}. For each value of r, we perform 10 independent experiments with generated signals of size S = 1000. For each experiment, we perform 100 EM iterations on independently generated samples (X T 1 , Y T 1 ) of size T = 20000. Figure 6 depicts an example of simulated data and the optimal (but approximated) smoothing output, and Table 2 reports the MSE results for different cases of fuzzy models and different values of r. The chief finding of this series of experiments is that when the number of fuzzy levels increases, the smoothed signal is closer to the "ground-truth" hidden signal: e.g. in case 3, for r = 20, the MSE decreases from 0.741 for F = 0 to 0.577 for F = 5.

Experiments on stochastic volatility models

Stochastic volatility (SV) models are widely used to highlight the variance of stochastic processes [START_REF] Ghysels | Stochastic volatility[END_REF]. Several variants of SV models have been studied (Henston, CEV, GARCH, Chen, etc.). In this paper, we consider standard SV models defined as follows:

CGOMFSM trajectory, together with the automatically restored states in magenta (dashed). where U i , V i are independent standard Gaussian vectors. The SV models is defined by the set of parameters σ, µ, and α. The main conclusion is that when the number of discrete fuzzy states increase, the model approaches the results of the optimal (but time consuming) particle smoother.

Conclusion

In this paper, we presented a novel approach to approximate non-linear Markov system using Conditionally Gaussian Observed Markov Fuzzy Switching Model (CGOMFSM). The chief novelty of this work is the introduction of fuzzy jumps instead of classical crisp states. This model still allows exact (up to required quantification) and fast smoothing equations. The fuzzy jumps allows transient modification of parameters, which is more appropriate for real-world applications. Future work includes the evaluation of the model for real data.