{'fraction_non_alphanumeric': 0.03523346590072175, 'fraction_numerical': 0.01116512004713507, 'mean_word_length': 4.336923197861803, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 40, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Kalman filter demystified: from intuition to probabilistic graphical model to real case in financial markets

In this paper, we revisit the Kalman filter theory. After giving the intuition on a simplified financial markets example, we revisit the maths underlying it. We then show that Kalman filter can be presented in a very different fashion using graphical models. This enables us to establish the connection between Kalman filter and Hidden Markov Models. We then look at their application in financial markets and provide various intuitions in terms of their applicability for complex systems such as financial markets. Although this paper has been written more like a self contained work connecting Kalman filter to Hidden Markov Models and hence revisiting well known and establish results, it contains new results and brings additional contributions to the field. First, leveraging on the link between Kalman filter and HMM, it gives new algorithms for inference for extended Kalman filters. Second, it presents an alternative to the traditional estimation of parameters using EM algorithm thanks to the usage of CMA-ES optimization. Third, it examines the application of Kalman filter and its Hidden Markov models version to financial markets, providing various dynamics assumptions and tests. We conclude by connecting Kalman filter approach to trend following technical analysis system and showing their superior performances for trend following detection.


Introduction
One of the most challenging question in finance is to be able from past observation to make some meaningful forecast. Obviously, as no one is God, no one would ultimately be able to predict with one hundred percent accuracy where for instance the price of a stock will end. However, between a pure blind guess and a perfect accuracy forecast, there is room for improvement and reasoning. If in addition, we are able somehow to model the dynamics of the stock price and factor in some noise due to unpredictable human behavior, we can leverage this model information to make an informed estimate, in a sense an educated guess that cannot miss the ballpark figure. It will not be 100 percent accurate but it is much better than a pure random appraisal. Indeed, this scientific question of using a model and filtering noise has been extensively examined in various fields: control theory leading to Kalman filter, Markov processes leading to hidden Markov models and lately machine learning using Bayesian probabilistic graphical models. In this work, we revisit these three fields to give a didactic presentation of this three approaches and emphasizing the profound connection between Kalman filter and Hidden Markov Models (HMM) thanks to Bayesian Probabilistic Graphical models. We show in particular that the derivation of Kalman filter equations is much easier once we use the more general framework of Bayesian Probabilistic Graphical models. In particular, we show that Kalman filter equations are just a rewriting of the sum product algorithm (also referred to as the Viterbi algorithm for HMM). We then provide various dynamics assumptions for financial markets and comment their overall performance. We compare these trading system with simpler moving average trend detection trading systems and show that they provide better performances.

Intuition
In a nutshell, a Kalman filter is a method for predicting the future state of a system based on previous ones. It was discovered in the early 1960's when Kalman introduced the method as a different approach to statistical prediction and filtering (see Kalman (1960) and Kalman and Bucy (1961)). The idea is to estimate the state of a noisy system. Historically, it was build to monitor the position and the speed of an orbiting vehicle. However, this can be applied to non physical system like economic system. In this particular settings, the state will be described by estimated equations rather than exact physical laws. We will emphasize this point in the section 5 dealing with practical application in finance.
To build an intuition, let us walk through a simple example -if you are given the data with green dots (that represent, by the way, the predicted Kalman filtered predicted price of the S&P500 index on October 17, 2018), it seems reasonable to predict that the orange dot should follow, by simply extrapolating the trend from previous samples and inferring some periodicity of the signal. However, how assured would you be anticipating the dark red point on the right (that is just 10 minutes later)? Moreover, how certain would you be about predicting the orange point, if you were given the black series (that represent the real prices) instead of the green one?
From this simple example, we can learn three important rules: • Predicting far ahead in the future is less reliable than near ahead.
• The reliability of your data (the noise) influences the reliability of your forecast. • It's not good enough to give a prediction -you also want to provide a confidence interval.

Fig 1: A simple example from financial markets to build an intuition of Kalman filter
From this simple presentation, we can conceptualize the approach and explain what a Kalman filter is. First, we need a state. A state is a representation of all the parameters needed to describe the current system and perform its prediction. In our previous example, we can define the state by two numbers: the price of the S&P 500 at time t, denoted by p(t) and the slope of the price at time t, s(t). In general, the state is a vector, commonly denoted x. Of course, you can include many more parameters if you wish to accommodate more complex systems.
Second, we need a model. The model describes how the system behaves. It provides the underlying equations that rules our system. It may be an ideal representation or a simplified version of our system. It may be justified by some physical laws (in the case for instance of a Kalman filter for a GPS system), or by some empirical analysis (in the case of finance). In the standard Kalman filter, the model is always a linear function of the state. In extended Kalman filter, it is a non linear function of the state. In our previous example, our model could be: • the price at time t, p(t), is obtained as the previous price p(t − 1) plus the slope at time t − 1: s(t − 1) p(t) = p(t − 1) + s(t − 1) (2.1) • the slope is evolving over time with a periodic sinusoidal function ψ(t) = a sin(bt + c) s(t) = s(t − 1) + ψ(t) (2.2) Usually, one expresses this model in matrix form to make it simple. We have In the above equation, following standard practice, we have denoted by x t = (p(t), s(t)) T , the state vector that combines the price p(t) and the slope s(t). We have denoted by u t = (0, ψ(t)) T , the control vector. The idea underlying the control vector is that it can control somehow the dynamic of the state vector. The previous equation (2.3) is called the state equation and writes Of course, our model is too simple, else we wouldn't need a Kalman Filter! To make it more realistic, we add an additional term -the noise process, w t that represents anything else that we do not measure in our model. Although we don't know the actual value of the noise, we assume we can estimate how "noisy" the noise is and can make assumptions about the noise distribution.
In standard Kalman modeling, we keep things simple and assumed this noise to be normally distributed. To measure how "noisy" the noise is, we refer to the variance of the normal distribution underpinning the noise w t−1 . The state equation is modified into: x t = Fx t−1 + B t · u t + w t−1 (2.5) Third, we need to use measurement to improve our model. When new data arrived, we would like to change the value of our model parameters to reflect our improved understanding of the state dynamic. We would proceed in two steps: make a prediction and then a correction based on what has happened in time t. There is a subtleties here. What we measure does not have to be exactly the states. It has to be related but not the same. For instance, in a Kalman filter system in a GPS, the state could be the GPS's acceleration, speed and position, while what we measure may be the car's position and the wheel velocity. In our financial example, we may only be able to observe the price but not the slope. Our measure would therefore be limited to the price in this particular setting. Hence, the measurement would be represented as follows: In general, the measurement should be a vector, denoted by z, as we may have more than one number for our measurement. Measurements should also be noisy to reflect that we do not measure perfectly. Hence, the equation specify measurement would write: Where v t is the measurement noise, and H is a matrix with rows equal to the number of measurement variables and columns equal to the number of state variables. Fourth, we need to work on prediction as we have set up the scene in terms of modeling. This is the key part in Kalman filter. To build intuition, let us assume, for the time being, noise is equal to zero. This implies our model is perfect. How would you predict the state at time t, x t knowing the previous state (at time t − 1)? This is dammed simple. Just use the state equation and compute our prediction of the state, denoted byx t as: Hold on! In the mean time, we do some measurements reflected by our measurement equation: These are measurements so what we measure might be slightly different. Intuitively, the measurement error computed as z t −ẑ t would be slightly different from 0. We would call this difference y = z t −ẑ t , the innovation. It represents the bias of our measurement estimation. Obviously, if everything were perfect, the innovation should be zero, as we do not have any bias! To incorporate the innovation in our model, we would add it to our state. But we would not add blindly in the state equation. We should multiply this innovation by a matrix factor that reflects somehow the correlation between our measurement bias and our state bias. Hence, we would compute a correction to our state equation as follows:x The matrix K t is called the Kalman gain. We would see in section 3 how to determine it but for the sake of intuition of Kalman filter, we skip this details. What really matters is to build an intuition. We can easily understand the following rules of thumb: 1. The more noisy our measurement is, the less precise it is. Hence, the innovation that represents the bias in our measurement may not be real innovation but rather an artifact of the measurement noise. Noise or uncertainty is directly captured by variance. Thus the larger the measurement variance noise, the lower the Kalman gain should be. 2. The more noisy our process state is, the more important the innovation should be taken into account. Hence, the larger the process state variance, the larger the Kalman gain should be.
Summarizing these two intuitions, we would expect the Kalman gain to be: Fifth, noise itself should be modeled as we have no clue about real noise. Noise in our model is represented by variance, or more precisely by the covariance of our state. Traditionally, we denote by P t the covariance matrix of the state: Using the state equation, and using the fact that for any matrix A, Cov(A.X) = A Cov(X)A , we can derive P t from its previous state: At this stage, we can make the model even more realistic. Previous equation (2.13) assumes that our process model is perfect. But keep in mind that we temporarily assumed no noise. However, real world is more complex and we should now use our real state equation (2.5). In particular, if we assume that the state noise w t is independent from the state x t and if the state noise w t is assumed to be distributed as a normal distribution with covariance matrix Q t , the prediction for our state covariance matrix becomes: (2.14) Likewise, we could compute the covariance matrix for the measurement and model the evolution of noise that we will denote by S t . The last source of noise in our system is the measurement. Following the same logic we obtain a covariance matrix forẑ t and denoting by v t the normally distributed independent noise for our measurements, with covariance matrix given by R t , we get Let us come back to our Kalman gain computation. Intuitively, we found that it should be larger for larger process noise and lower for large measurement noise, leading to an equation of the type K t ∼ Process Noise Measurement Noise . As process noise is measured P t and measurement noise by S t , we should get something like K t = P t S −1 t . We shall see in the next section that the real equation for the Kalman gain is closely related to our intuition and given by

The maths
In this section, we will prove rigorously all the equations provided in section 2. To make notations more robust and make the difference for a time dependent vector or matrix A between its value at time t, knowing information up to time t − 1, and its value at time t, knowing information up to time t, we will denote these two different value A t|t−1 and A t|t .

Kalman filter modeling assumption
The Kalman filter model assumes the true state at time t is obtained from the state at the previous time t − 1 by the following state equation The noise process w t is assumed to follow a multi dimensional normal distribution with zero mean and covariance matrix given by Q t : w t ∼ N (0, Q t ).
At time t, we make a measurement (or an observation) z t of the true state x t according to our measurement equation: Like for the state equation, we assume that the observation noise v t follows a multi dimensional normal distribution with zero mean and covariance matrix given by R t : v t ∼ N (0, R t ). In addition, the initial state, and noise vectors at each step x 0 , w 1 , . . . , w t , v 1 , . . . , v t are assumed to be all mutually independent.

Properties
It is immediate to derive the prediction phase as the estimated value of the statex t|t−1 is simply the expectation of the state equation (2.5). Similarly, it is trivial to derive the estimated value for the error covariance. This provides the prediction phase that is summarized beloŵ The correction phase that consists in incorporating the measurements to correct our prediction imsart-generic ver. 2014/10/16 file: article_KalmanFilterDemystified.tex date: December 14, 2018 is slightly more tricky. It consists in the following equations: Proposition 1. Under the assumptions stated in subsection 3.1, the Kalman gain that minimizes the expected squared error defined as the square of the Euclidean (or L 2 ) norm of the error vector, representing the error between the true state and our estimated state : x t −x t|t is given by It is referred to as the optimal Kalman gain. For any Kalman gain (and not necessarily the optimal one), the estimate covariance updates as follows: For the optimal Kalman filter, this reduces to the usual Kalman filter Updated estimate covariance as follows: Proof. The derivation of all these formulae consists in three steps. First, we derive the posteriori estimate of the covariance matrix. Then we compute the Kalman gain as the minimum square error estimator. Third, we show that in the particular case of the optimal gain the equation for the updated estimate, covariance reduces to the formula provided in equation (3.9). All these steps are provided in appendix A to make the reading of this article smooth.

Kalman as Graphical model
Historically, Hidden Markov Model (HMM) and Kalman filter were developed in distinct and unconnected research communities. Hence, their close relationship has not always been widely emphasized and appreciated. Part of the explanation lies also to the fact that the general framework for unifying these two approaches, namely graphical models came much later than HMM and Kalman filter. Without Bayesian graphical framework, the two algorithms underlying the inference calculation look rather different and unrelated. However, their difference is simply a consequence of the differences between discrete and continuous hidden variables and more specifically between multinomial and normal distribution. These details, important as they may be in practice, should not obscure us from the fundamental similarity between these two models. As we shall see in this section, the inference procedure for the state space model (SSM) shall prove us shortly that HMM and Kalman filter's model are cousin and share the same underlying graphical model structure, namely a hiddent state space variable and an observable variable. The interest of using Bayesian Probabilistic Graphical model is multiple. First, it emphasizes the general graphical model architecture underpinning both HMM and Kalman filter. Second, it provides modern computational tools used commonly in machine learning to do the inference calculation. It shows how to generalize Kalman filter in the case of non Gaussian assumptions. It is interesting to realize that graphical models have been the marriage between probability theory and graph theory. They provide a natural tool for dealing with two problems that occur throughout applied mathematics and engineering uncertainty and complexity and in particular they are playing an increasingly important role in the design and analysis of machine learning algorithms.
From a literature point of view, Hidden Markov models were discussed as early as Rabiner and Juang (1986), and expanded on in Rabiner (1989). The first temporal extension of probabilistic graphical models is due to Dean and Kanazawa (1989), who also coined the term dynamic Bayesian network. Much work has been done on defining various representations that are based on hidden Markov models or on dynamic Bayesian networks; these include generalizations of the basic framework, or special cases that allow more tractable inference. Examples include mixed memory Markov models (see Saul and Jordan (1999)); variable-duration HMMs (Rabiner (1989)) and their extension segment models (Ostendorf et al. (1996)); factorial HMMs (Ghahramani and Jordan (1994)); and hierarchical HMMs (Fine et al. (1998) and Bui et al. (2002)). Smyth et al. (1997) is a review paper that was influential in providing a clear exposition of the connections between HMMs and DBNs. Murphy and Paskin (2001) show how hierarchical HMMs can be reduced to DBNs, a connection that provided a much faster inference algorithm than previously proposed for this representation. Murphy (2002) (and lately the book Murphy (2013)) provides an excellent tutorial on the topics of dynamic Bayesian networks and related representations as well as the non published book of Jordan (2016)

State space model
The state space model as emphasized for instance in Murphy (2013) is described as follows: • there is a continuous chain of states denoted by (x t ) t=1,...,n that are non observable and influenced by past states only through the last realization. In other words, x t is a Markov process, meaning P(x t | x 1 , x 2 , . . . , x t−1 ) = P(x t | x t−1 ) Using graphical model, this can also be stated as given a state at one point in time, the states in the future are conditionally independent of those in the past. • for each state, we can observe a space variable denoted by z t that depends on the non observable space x t Compared to the original presentation of the Kalman filter model, this is quite different. We now assume that there is an hidden variable (our state) and we can only measure a space variable. Since at each step, the space variable only depends on the non observable, there is only one arrow or edge between two latent variables horizontal nodes. This model is represented as a graphical model in figure 2.
Obviously, the graphical model provided in figure 2 can host both HMM and Kalman filter model. To make our model tractable, we will need to make additional assumptions. We will emphasize the ones that are common to HMM and Kalman filter models and the ones that differ. We impose the following conditions: .   Fig 2: State Space model as a Bayesian Probabilistic Graphical model. Each vertical slice represents a time step. Nodes in white represent unobservable or latent variables called the states and denoted by x t while nodes in gray observable ones and are called the spaces and denoted by z t . Each arrow indicates that there is a relationship between the arrow originating node and the arrow targeting node. Dots indicate that there is many time steps. The central dot line is to emphasize the fundamental difference between latent and observed variables • The relationship between the state and the space is linear (this is our measurement equation in the Kalman filter and this is common to Kalman filter and HMM models): where the noise term v t is assumed to follow a multi dimensional normal distribution with zero mean and covariance matrix given by R t . • We will make the simplest choice of dependency between the state at time t − 1 and t and assume that this is linear (this is our state equation and this is common to Kalman filter and HMM models) where the noise term w t is assumed to follow a multi dimensional normal distribution with zero mean and covariance matrix given by Q t and where B t u t is an additional trend term (that represents our control in Kalman filter). This control term is not common in HMM model but can be added without any difficulty. This results in slightly extended formula that we will signal. These formula are slight improvement of the common one found in the literature. Although, from a theoretical point of view, this control term may seem a futility, it is very important in practice and makes a big difference in numerical applications. • We will assume as in the Kalman filter section that the initial state, and noise vectors at each step x 0 , w 1 , . . . , w t , v 1 , . . . , v t are all mutually independent (this is common to HMM and Kalman filter models). • Last but not least, we assume that the distribution of the state variable x t follows a multi dimensional normal distribution. This is Kalman filter specific. For HMM, the state is assumed to follow a multinomial distribution.
The above assumptions restrict our initial state space model to a narrower class called the Linear-Gaussian SSM (LG-SSM). This model has been extensively studied and more can be found in Durbin and Koopman (2012) for instance.
Before embarking into the inference problem for the SSM, it is interesting to examine the unconditional distribution of the states x t . Using equation (4.2) The unconditional mean of x t is computed recursively with the implicit assumption that empty product equals 1: In the specific case of a null control term (B t = 0), the latter equation simplifies into The unconditional co-variance is computed as follows Using our assumptions on independence as well as the state equation (4.2), we can compute it easily as: This last equation remains unchanged in case of a non zero control term as the control term is deterministic. This last equation provides a dynamic equation for the unconditional variance and is referred to as the Lyapunov equation. It is also easy to checked that the unconditional covariance between neighboring states x t and x t+1 is given by F t+1 P t F T t+1 .

Inference
The inference problem consists in calculating the posterior probability of the states given an output sequence. This calculation can be done both forward and backward. By forward, we mean that the inference information (called the evidence) at time t consists of the partial sequence of outputs up to time t. The backward problem is similar except that the evidence consists of the partial sequence of outputs after time t. Using standard graphical model terminology, we distinguish between filtering and smoothing problem. In filtering, the problem is to calculate an estimate of the state x t based on a partial output sequence z 0 , . . . , z t . That is, we want to calculate P(x t | z 0 , . . . , z t ). This is often referred to as the alpha recursion in HMM models (see for instance Rabiner and Juang (1986) and Rabiner (1989)).
Using standard Kalman filter notations, we shall denote bŷ Under this settings, it is fairly easy to derive the following property that provides the conditional posterior distribution in the forward recursion. And to recover traditional results of Kalman filter, we shall decompose our time propagation into two steps: • measurement update: This can be represented nicely in terms of graphical models by the figure 3 below.
Proposition 2. Conditioned on past outputs z 0 , . . . , z t , the variables x t+1 and z t+1 have a joint Gaussian distribution with mean and covariance matrix given by: Proof. This is trivial as the considered state space model is the Linear Gauss State Space Model (LGSSM).
As simple as it may seem, the previous proposition makes our graphical model a full powerhouse as it provides the building block to start the inference. Indeed, knowing the conditional posterior distribution at step (a) of figure 3 gives us the first bullet to conclude for the step (b). Moreover, using results from simpler graphical models like the factor analysis graphical model, we can immediately conclude that the second step is given as follows Proposition 3. Conditioned on past outputs z 0 , . . . , z t+1 , we have the following relationship Proof. This is a direct consequence of the spatial model and can be found for instance in Murphy (2013) or Jordan (2016). We provide a self contained proof in appendix B Summarizing all these results leads to the seminal recursions of the Kalman filter as provided in the section 3 and given by the following proposition.
Proposition 4. Kalman filter consists in the following recursive equations: Proof. This is a trivial consequence of proposition 3. Equation (4.11) (resp. (4.12)) is the same as the one provided in (3.3) (resp. (3.4)) but using graphical models induction.
Remark 4.1. If we introduce the following intermediate variables already provided in section 3: The equations (4.13) and (4.14) transform into equationŝ which are equations (3.8) and (3.9). This proves that the derivation using graphical models and control theory are mathematically equivalent! Remark 4.2. There is nothing new to fancy at this stage except that we have shown with graphical models the Kalman filter recursion. And we can check this is way faster, easier and more intuitive.
It is worth noticing that we have also multiple ways to write the gain matrix. Using the Sherman-MorrisonWoodbury formula, we can also derive various forms for the gain matrix as follows: These forms may be useful whenever the reduced form (which is the last equation) is numerically unstable. Equation (4.23) is useful as it relates K t+1 to P t+1|t+1 . It is interesting to notice the two form of the Kalman gain imsart-generic ver. 2014/10/16 file: article_KalmanFilterDemystified.tex date: December 14, 2018 Remark 4.3. The Kalman filter appeals some remarks. We can first note that Kalman filtering equations can be interpreted differently. Combining equations (3.8) and (3.3), we retrieve an error correcting algorithm as follows: or equivalently, regrouping the termx t−1|t−1 This shows us two things: • Kalman filter can be seen as the discrete version of an Ornstein Uhlenbeck process • Kalman filter can be seen as Auto Regressive process in discrete times Since, recursive equation do appear similarly in Recursive Least Square (RLS) estimates, we also see here a connection with RLS. It is striking that these connections are not often made, mostly because Kalman filter was originally a control problem and not a statistician one.
As nice as the Kalman filter equation may look like, they have one major problem. It is the numerical stability of the filter. If the process noise covariance Q t is small, round-off error would lead to obtain numerically a negative number for small positive eigenvalues of this matrix. As the scheme will propagate round off errors, the state covariance matrix P t will progressively become only positive semi-definite (and hence indefinite) while it is theoretically a true positive definite matrix fully invertible.
In order to keep the positive definite property, we can slightly modify the recursive equation to find recursive equations that preserve the positive definite feature of the two covariance matrices. Since any positive definite matrix S d can be represented and also reconstructed by its upper triangular square root matrix R with S d = R t R T , with the strong property that representing in this form will guarantee that the resulting matrix will never get a negative eigen value, it is worth using a square root scheme with square root representation. Alternatively, we can also represent our covariance matrix using the so called unit diagonal (U-D) decomposition form, with S d = UDU T where U' is a unit triangular matrix (with unit diagonal), and D' is a diagonal matrix. This form avoids in particular many of the square root operations required by the matrix square root representation. Moreover, when comparing the two approaches, the U-D form has the elegant property to need same amount of storage, and somewhat less computation. This explains while the U-D factorization is often preferred Thornton (1976). A slight variation is the LDL decomposition of the innovation covariance matrix. The LDL decomposition relies on decomposing the covariance matrix S d with two matrices: a lower unit triangular (unitriangular) matrix L, and D a diagonal matrix. The algorithm starts with the LU decomposition as implemented in the Linear Algebra PACKage (LAPACK) and further it factors into the LDL form. Any singular covariance matrix is pivoted so that the first diagonal partition is always non-singular and well-conditioned (for further details see Bar-Shalom et al. (2002)).

Connection to information filter
It is worth showing the close connection with particle and information filter. This is a direct consequence of the fact that a multivariate Gaussian belongs to the exponential family and as such admits canonical parameters. Hence, we can rewrite the filter in terms of the later instead of the initial usage of moment parameters. This is an interesting rewriting of the Kalman filter as it makes it numerically more stable. The canonical parameters of a multi variate Gaussian distribution, denoted by Λ and η, are obtained from the moment parameters Σ and µ as follows: Λ = Σ −1 and η = Σ −1 µ. We are interested in deriving the canonical parameters of x t first, at the prediction phase, conditioned on z 1 , . . . , z t−1 .
In the Kalman filter settings, the covariance Σ is denoted by P while the moment is given by x with the relationship η = P −1 x. Hence, we will write the precision matrix as Λ with the relationship with the covariance matrix given by Λ = P −1 . We shall write our new variables to be consistent with previous development as Λ t|t−1 and η t|t−1 . At the correction or measurement phase, we are interested in the same parameters but now conditioned on z 1 , . . . , z t . We shall write them Λ t|t and η t|t . We can easily derive the following recursive scheme that is given by the proposition below: Proposition 5. The filter equations are given by the following recursive scheme: Proof. The derivation is easy and given in section B.3 Remark 4.4. The recursive equation for the information filter are very general. They include the control term B t u t that is often neglected in literature introduced as early as 1979 in Anderson and Moore (1979). It is worth noticing that we can simplify computation by pre-computing a term M t as follows: These equations are more efficient than the ones provided in proposition 5. As for the initialization of this recursion, we define the initial value as followsη 1|0 =η 1 and Λ 1|0 = Λ 1 . It is interesting to note that the Kalman filter and the information filter are mathematically equivalent. They both share the same assumptions. However, they do not use the same parameters. Kalman filter (KF) uses moment parameters while particle or information filter (IF) relies on canonical parameters, which makes the later numerically more stable in case of poor conditioning of the covariance matrix. This is easy to understand as a small eigen value of the covariance translates into a large eigen value of the precision matrix as the precision matrix is the inverse of the covariance matrix. Reciprocally, a poor conditioning of the information filter should convert to a more stable scheme for the Kalman filter as the he condition number of a matrix is the reciprocal of the condition number of its inverse. Likewise, for initial condition as they are inverse, a small initial state in KF should translate to a large state in IF.

Smoothing
Another task we can do on dynamic Bayesian network is smoothing. It consists in obtaining estimates of the state at time t based on information from t on-wards (we are using future information in a sense). Like for HMM, the computation of this state estimate requires combining forward and backward recursion, either by starting by a backward-filtered estimates and then a forwardfiltered estimates (an 'alpha-beta algorithm"), or by doing an algorithm that iterates directly on the filtered-and-smoothed estimates (an "alpha-gamma algorithm"). Both kinds of algorithm are available in the literature on state-space models (see for instance Koller and Friedman (2009) and Cappe et al. (2010) for more details), but the latter approach appears to dominate (as opposed to the HMM literature, where the former approach dominates) . The "alpha-gamma" approach is referred to as the "Rauch-Tung-Striebel (RTS) smoothing algorithm" (although it was developed using very different tool, namely control theory as early as 1965: see Rauch et al. (1965)) while the other approach is just the "alpha-beta" recursion.

Rauch-Tung-Striebel (RTS) smoother
The Rauch-Tung-Striebel (RTS) smoother developed in Rauch et al. (1965) relies precisely on the idea of of doing first a backward estimate and then a forward filter. Smoothing should not be confused with smoothing in times series analysis that is more or less a convolution. Smoothing for a Bayesian network means infering the distribution of a node conditioned on future information. It is the reciprocal of filtering that has also two meanings. Filtering for time series means doing a convolution to filter some noise. But filtering for Bayesian network means infering the distribution of a node conditioned on past information.
We can work out the RTS smoother and find the recursive equations provided by the following proposition Proposition 6. The RTS smoothing algorithm works as follows: with an initial condition given bŷ Proof. The proof consists in writing rigorously the various equations and is given in section B.4

Alternative to RTS filter
It is worth noting that we can develop an alternative approach to the RTS filter that relies on the alpha beta approach without any observation. This approach is quite standard for HMM models (see for instance Rabiner and Juang (1986) or Russell and Norvig (2009)). Following Kitagawa (1987), this approach has been called in the literature the two-filter algorithm. It consists in combining the forward conditional probability P(x t | z 1 , . . . , z t ) with the backward conditional probability P(x t | z t+1 , . . . , z T ). The intuition is illustrated by the figure 4. The graphical model provided by figure 4 can be easily characterized. The joint probability distribution of (x t , x t+1 ) is a multi variate normal whose Lyapunov equation (equation of the covariance matrix) is given by Hence the covariance matrix of (x t , x t+1 ) is given by The underlying idea in this approach is to invert the arrows in the graphical model leading to a new graphical model given by . . . Mathematically, this imply to change the relationship and now solve for P t in terms of P t+1 using equation (4.38). We get: where we have assumed that F t+1 is invertible. As a matter of fact, if it is not the case, we can rewrite everything with an approaching matrix to F t+1 (in the sense of the Frobenius norm) that is invertible, whose distance is less than ε, derive all the relationship below and then take the limit as tends to zero. We will therefore assume in the following that F t+1 is invertible as this is not a strict condition for our derivation. Hence we can rewrite the covariance matrix now as follows: (4.41) If we define the new matrix F t+1 as: we can simplify the covariance matrix as follows: (4.43) As this looks like a forward covariance matrix. Recall that the forward dynamics is defined by: The following proposition gives the inverse dynamics: Proposition 7. The inverse dynamics is given by: and with w t+1 independent of the past information x t+1 , . . . , x T and with the covariance matrix of w t+1 given by and with the forward Lyapunov equation given by: Proof. The proof is straightforward and given in B.5.
The reasoning followed here gives us a new way to generate an information filter but now backward. We should note that the conversion between the state space and the observable variable is unchanged and given by We can now apply the information filter developed previously to get the new filtering equations using canonical instead of moment parameterization. We get the following new information filter : Proposition 8. The modified BrysonFrazier as presented in Bierman (2006) is given by: Proof. The proof consists in combining all previous results and is given in B.6.
Remark 4.5. If we want to convert to to the moment representation, we can use the inverse transform given byx t|t+1 = S −1 t|t+1η t|t+1 and P t|t+1 = S −1 t|t+1

Inferring final posterior distribution
This problem consists in inferring the final posterior distribution P(x t | z 1 , . . . , z T ). We can easily derivex t|T = P t|T (P −1 t|tx t|t + P −1 t|t+1x t|t+1 ) (4.55) and 4.6. Parameter estimation 4.6.1. Intuition Since the connection between HMM and Kalman filter, it has been established that the best way to estimate the initial parameters of a Kalman filter and any of its extension is to rely on EM estimation. The expectationmaximization (EM) algorithm is an iterative method to find the maximum likelihood or maximum a posteriori (MAP) estimates of the parameters of a statistical model. This approach relies on the fact that the model depends on non observable (also called latent) variables. The EM algorithm alternates between performing an expectation (E) step, which provides the expectation of the conditional log-likelihood evaluated using the current estimate for the parameters, and a maximization (M) step, which computes parameters maximizing the expected log-likelihood found on the E step. These parameter-estimates are then used to determine the distribution of the latent variables in the next E step.
The EM algorithm was explained and given its name in a classic 1977 paper by Dempster et al. (1977). The intuition is to find a way to solve the maximum likelihood solution for a model with latent variables. Since variables are hidden, solving the maximum likelihood solution typically requires taking the derivatives of the likelihood function with respect to all the unknown values and then solving the resulting equations. Because of latent variables, this function is unknown and the problem can not be solved. Instead, we compute the expected maximum likelihood with respect to the latent variables. We can then solve this explicit maximum likelihood function with respect to the observable variables.

EM Algorithm
We assume that our statistical model has a set X of observed data and a set of unobserved latent Z that depend on unknown parameters θ that we want to determine thanks to maximum likelihood. We obviously know the likelihood function L(θ; X, Z) = p(X, Z|θ) given by The later is often intractable as the number of values for Z is vary large making the computation of the expectation (the integral) extremely difficult. In order to avoid the exact computation, we can work as follows: First compute the expectation explicitly (called the E step) by taking the expected value of the log likelihood function of θ with respect to latent variables Z, denoted by Q(θ|θ (t) : We then maximize the resulted function with respect to the parameters θ: The EM method works both for discrete or continuous random latent variables. It takes advantage of algoirthm like the Viterbi algorithm for hidden Markov models or recursive equation fo Kalman filter. The various steps are the following: Algorithm 1 EM algorithm: The convergence of this algorithm is given by the following proposition.
Proposition 9. The EM algorithm converges monotonically to a local minimum for the marginal log likelihood.
Proof. There are at least two ways of proving this results that are provided in D.1 and D.2.

Application to Kalman filter
The EM algorithm was initially developed for mixture models in particular Gaussian mixtures but also other natural laws from the exponential family such as Poisson, binomial, multinomial and exponential distributions as early as in Hartley (958). It was only once the link between latent variable and Kalman filter models was made that it became obvious that this could also be applied to Kalman and extended Kalman filter (see Cappe et al. (2010) or Einicke et al. (2010)). The EM method works as folllows: Algorithm 2 EM algorithm for Kalman filter init: Set the Kalman filter parameters to some random values. while Not Converged do Use Kalman filter recursive equations to obtained updated state estimates (E-step) Use the filtered or smoothed state estimates within maximum-likelihood calculations to obtain updated parameter estimates. For instance for the standard Kalman filter, the distributions are assumed to be normal. The corresponding maximum likehood estimate for the updated measurement noise variance estimate leads to a variance given byσ wherex k are the output estimates calculated by the Kalman filter for the measurements z k . Likewise, the the updated process noise variance estimate is calculated aŝ wherex k andx k+1 are scalar state estimates calculated by a filter or a smoother. The updated model coefficient estimate is obtained viaF end while

CMA-ES estimation
Another radically difference approach is to minimize some cost function depending on the Kalman filter parameters. As opposed to the maximum likelihood approach that tries to find the best suitable distribution that fits the data, this approach can somehow factor in some noise and directly target a cost function that is our final result. Because our model is an approximation of the reality, this noise introduction may leads to a better overall cost function but a worse distribution in terms of fit to the data. Let us first introduce the CMA-ES algorithm. Its name stands for covariance matrix adaptation evolution strategy. As it points out, it is an evolution strategy optimization method, meaning that it is a derivative free method that can accommodate non convex optimization problem. The terminology covariance matrix alludes to the fact that the exploration of new points is based on a multinomial distribution whose covariance matrix is progressively determined at each iteration. Hence the covariance matrix adapts in a sense to the sampling space, contracts in dimension that are useless and expands in dimension where natural gradient is steep. This algorithm has led to a large number of papers and articles and we refer to Varelas et al. (2018) CMA-ES relies on two main principles in the exploration of admissible solution for our optimization problem. First, it relies on a multi variate normal distribution as this is the maximum entropy distribution given the first two moments. The mean of the multi variate distribution is updated at each step in order to maximize the likelihood of finding a successful candidate. The second moment, the covariance matrix of the distribution is also updated at each step to increase the likelihood of successful search steps. These updates can be interpreted as a natural gradient descent. Intuitively, the CMA ES algorithm conducts an iterated principal components analysis of successful search steps while retaining all principal axes.
Second, we retain two paths of the successive distribution mean, called search or evolution paths. The underlying idea is keep significant information about the correlation between consecutive steps. If consecutive steps are taken in a similar direction, the evolution paths become long. The evolution paths are exploited in two ways. We use the first path is to compute the covariance matrix to increase variance in favorable directions and hence increase convergence speed. The second path is used to control step size and to make consecutive movements of the distribution mean orthogonal in expectation. The goal of this step-size control is to prevent premature convergence yet obtaining fast convergence to a local optimum.
In order to make it practical, we assume that we have a general cost function that depends on our Bayesian graphical model denoted by Φ(θ) where θ are the parameters of our Kalman filter. Our cost function is for instance the Sharpe ratio corresponding to a generic trend detection strategy whose signal is generated by our Bayesian graphical model that is underneath a Kalman filter. This approach is more developed in a companion paper Benhamou (2018) but we will give here the general idea. Instead of computing the parameter of our Bayesian graphical model using the EM approach, we would like to find the parameters θ max that maximize our cost function Φ(θ). Because our cost function is to enter a long trade with a predetermined target level and a given stop loss whenever our Bayesian graphical model anticipates a price risen and similarly to enter a short trade whenever our prediction based on Bayesian graphical model is a downside movement, our trading strategy is not convex neither smooth. It is a full binary function and generates spike whenever there is a trade. Moreover, our final criterium is to use the Sharpe ratio of the resulting trading strategy to compare the efficiency of our parameters. This is way too complicated for traditional optimization method, and we need to rely on Black box optimization techniques like CMA-ES. Before presenting results, we will also discuss in the following section computational issues and the choice of the State space model dynamics.
We will describe below the most commonly used CMA -ES algorithm referred to as the µ/µ w , λ version. It is the version in which at each iteration step, we take a weighted combination of the µ best out of λ new candidate solutions to update the distribution parameters. The algorithm has three main loops: first, it samples new solutions, second, it reorders the sampled solutions based on their fitness and third it updates the state variables. A pseudo code of the algorithm is provided just below: In this algorithm, for an n dimensional optimization program, the five state variables at iteration step k are:

Algorithm 3 CMA ES algorithm
1. m k ∈ R n , the distribution mean and current best solution, 2. σ k > 0, the variance step-size, 3. C k , a symmetric positive-definite n × n covariance matrix initialized to I n , 4. p σ ∈ R n , the isotropic evolution path, initially set to null, 5. p c ∈ R n , the anisotropic evolution path, initially set to null.
It is worth noting that the order of the five update is important. The iteration starts with sampling λ > 1 candidate solutions x i ∈ R n from a multivariate normal distribution N (m k , σ 2 k C k ), that is x i ∼ m k + σ k N (0, C k ) for i = 1, ..., λ.
We then evaluate candidate solutions x i for the objective function f : R n → R of our optimization.
We then sort candidate solution according to their objective function value: It is worth noticing that we do not even need to know the value. Only the ranking is important for this algorithm.
The new mean value is updated as a weighted mean of our new candidate solutions where the positive (recombination) weights w 1 ≥ w 2 ≥ · · · ≥ w µ > 0 sum to one. Typically, µ ≤ λ/2 and the weights are chosen such that µ w := 1/ µ i=1 w 2 i ≈ λ/4. The step-size σ k is updated using cumulative step-size adaptation (CSA), sometimes also denoted as path length control. The evolution path p σ is updated first as it is used in the update of the step-size σ k : where • c −1 σ ≈ n/3 is the backward time horizon for the evolution path p σ and larger than one (c σ 1 is reminiscent of an exponential decay constant as −1 is the variance effective selection mass and 1 ≤ µ w ≤ µ by definition of is the unique symmetric square root of the inverse of C k , and • d σ is the damping parameter usually close to one. For d σ = ∞ or c σ = 0 the step-size remains unchanged. The step-size σ k is increased if and only if p σ is larger than √ n (1 − 1/(4 n) + 1/(21 n 2 )) and decreased if it is smaller (see Hansen (2006) for more details).
Finally, the covariance matrix is updated, where again the respective evolution path is updated first.
where T denotes the transpose and c −1 c ≈ n/4 is the backward time horizon for the evolution path p c and larger than one, α ≈ 1.5 and the indicator function 1 [0,α √ n] ( p σ ) evaluates to one if and only if p σ ∈ [0, α √ n] or, in other words, p σ ≤ α √ n, which is usually the case, makes partly up for the small variance loss in case the indicator is zero, c 1 ≈ 2/n 2 is the learning rate for the rank-one update of the covariance matrix and c µ ≈ µ w /n 2 is the learning rate for the rank-µ update of the covariance matrix and must not exceed 1 − c 1 .
The covariance matrix update tends to increase the likelihood function for p c and for (x i:λ − m k )/σ k to be sampled from N (0, C k+1 ). This completes the iteration step. The stop condition is quite standard and makes sure we stop if the best solution does not move or if we reached the maximum iterations number.
Remark 4.6. The number of candidate samples per iteration, λ, is an important and influential parameter. It highly depends on the objective function and can be tricky to be determined. Smaller values, like 10, tends to do more local search. Larger values, like 10 times the dimension n makes the search more global. A potential solution to the determination of the number of candidate samples per iteration is to repeatedly restart the algorithm with increasing λ by a factor of two at each restart (see Auger and Hansen (2005)) Besides of setting λ (or possibly µ instead, if for example λ is predetermined by the number of available processors), the above introduced parameters are not specific to the given objective function. This makes this algorithm quite powerful as the end user has relatively few parameters to set to do an efficient optimization search.

Computational issues
Although the computational cost of Kalman filter is very small and its complexity dammed simple, there are room for improvements. The two dominant costs in the Kalman filter are the matrix inversion to compute the Kalman gain matrix, K t , which takes O(|x t | 3 ) time where |x t | stands for the dimension of x t ; and the matrix-matrix multiplication to compute the estimated covariance matrix which takes O(|z t | 2 ) time.
In some applications (e.g., robotic mapping), we have |z t |x t |, so the latter cost dominates. However, in such cases, we can sometimes use sparse approximations (see for instance Thrun et al. (2005)). In the opposite cases where |x t | |z t , we can precompute the Kalman filter gain K t , since, it does not depend on the actual observations (z 1 , . . . , z t ). This is a an unusual property that is specific to linear Gaussian systems.
The iterative equations for updating the estimate covariance in equation (3.9) or in (A.6) are Riccati equations as they are first order ordinary differential recursive equations quadratic in Kalman gain. Standard theory state that for time invariant systems, they converge to a fixed point. Hence, it may be wised to use this steady state solution instead of using a time-specific gain matrix in our filtering. In practice however, this does not provide better results. More sophisticated implementations of the Kalman filter should be used, for reasons of numerical stability. One approach is the information filter also called particle filter, which recursively updates the canonical parameters of the underlying Gaussian distribution, namely the precision matrix defined as the inverse of the covariance matrix Λ t|t = P t|t and the first order canonical parameter given by η t|t = Λ t|t x t|t . These parameters instead of the traditional moment parameters makes the inference stronger.
Another approach is the square root filter, which works with the Cholesky decomposition or the U DU factorization of the estimated error covariance: P t|t . This is much more numerically stable than directly updating the estimated error covariance P t|t . Further details can be found at the reference site http://www.cs.unc.edu/ welch/kalman/ and the refernce book for Kalman filter Simon (2006).

In practice
Financial markets are complex dynamics models. Notoriously, one does not know the exact dynamics as opposed to a noise filtering problem for a GPS. We are doomed to make assumptions about the dynamics of financial systems. This can seriously degrade the filter performance as we are left with some un-modeled dynamics. As opposed to noise that the filter can accommodate (it has been designed for it), non modeled dynamics depends on the input. They can bring the estimation algorithm to instability and divergence. It is a tricky and more experimental issue to distinguish between measurement noise and non modeled dynamics.

Numerical experiments
In order to test our results, we look at the following trend following algorithm based on our Kalman filter algorithm where we enter a long trade if the prediction of our kalman filter is above the close of the previous day and a short trade if the the prediction of our kalman filter is below the close of the previous day. For each comparison, we add an offset µ to avoid triggering false alarm signals. We set for each trade a pre-determined profit and stop loss target in ticks. These parameters are optimized in order to provide the best sharpe ratio over the train period together with the Kalman filter parameters. The pseudo code of our algorithm is listed below Our resulting algorithm depends on the following parameters p 1 , . . . , p n the Kalman filter algorithm, the profit target, the stop loss and the signal offset µ. We could estimate the Kalman filter parameters with the EM procedure, then optimize the profit target, the stop loss and the signal offset µ. However, if by any chance the dynamics of the Kalman filter is incorrectly specified, the noise generated by this wrong specification will only be factored in the three parameters: the profit target, the stop loss and the signal offset µ. We prefer to do a combined optimization of all the parameters. We use daily data of the S&P 500 index futures (whose CQG code is EP ) from 01Jan2017 to 01Jan2018. We train our model on the first 6 months and test it on the next six months. Deliberately, our algorithm is unsophisticated to keep thing simple and concentrate on the parameter estimation method. The overall idea is for a given set of parameter to compute the resulting sharpe ratio over the train period and find the optimal parameters combination. For a model like the fourth one in the table 1, the optimization encompasses 18 parameters: p 1 , . . . , p 15 , the profit target, the stop loss and the signal offset µ, making it non trivial. We use the CMA-ES algorithm to find the optimal solution. In our optimization, we add some penalty condition to force non meaningful Kalman filter parameters to be zero, namely, we add a L1 penalty on this parameters.
Results are given below     We compare our algorithm with a traditional moving average crossover algorithm to test the efficiency of Kalman filter for trend detection. The moving average cross over algorithm generates a buy signal when the fast moving average crosses over the long moving average and a sell signal when the former crosses below the latter. A d period moving average is defined as the arithmetic average of the daily close over a d period, denoted by SMA(d). Our algorithm is given by the following pseudo code We can now compare moving average cross over versus Kalman filter algorith. The table 7 compares the two algorithms. We can see that on the train period, the two algorithms have similar performances : 5, 260 vs 5, 086. However on the test period, moving average performs very badly with a net profit of 935 versus 4, 266 for the bayesian graphical model (the kalman filter) algorithm.

Conclusion
In this paper, we have revisited the Kalman filter theory and showed the strong connection between Kalman filter and Hidden Markov Models. After discussing various recursive scheme for the underlying Bayesian graphical model, we provided various insight on the parameter estimation.
The traditional method is to use the expectation maximization (EM) method. This method aims to find the best fit distribution according to our modeling assumptions. This may work well when the state space model is ruled by physical law. However, in finance, there is no such a thing as the physical law of markets. One makes assumptions about the price dynamics that may be completely wrong. Hence, we present a new method based on CMA-ES optimization that does not decouple the Bayesian graphical model parameters from the other parameters of the financial strategy. We show on a numerical example on the S&P 500 index futures over one year of data (from Jan 1st 2017 to Jan 1st 2018) that the combined estimation of the Bayesian graphical model works well and that there is no much performance deterioration between train and test. This suggests little over fitting. We also show that the traditional moving average cross over method does not perform that well compared to our method and show strong over fitting bias as the net performance over the test data set is dramatically lower than the one on the train, suggesting over fitting.

Appendix A: Kalman filter proof
The proof is widely known and presented in various books (see for instance Gelb (1974) or Brown and Hwang (1997)). We present it here for the sake of completeness. It consists in three steps: • Derive the posteriori estimate covariance matrix • Compute the Kalman gain • Simplify the posteriori error covariance formula A.1.
Step 1: A posteriori estimate covariance matrix Proof. The error covariance P t|t is defined as the covariance between x t and x t+1 : P t|t = Cov x t −x t|t . We can substitute in this equation the definition ofx t|t given by equation (3.8) to get: Using the definition ofỹ t from equation (3.5), we get: Using the measurement equation for z t (equation (3.2)), this transforms into: By collecting the error vectors, this results in: Since the measurement error v k is not correlated with the other terms, this simplifies into: By property of covariance for matrix: Cov(Ax t ) = A Cov(x t )A T and using the definition for P t|t−1 and R t = Cov (v t ), we get the final result: This last equation is referred to as the Joseph form of the covariance update equation. It is true for any value of Kalman gain K t no matter if it is optimal or not. In the special case of optimal Kalman gain, this further reduces to equation (3.9) which we will prove in A.3

A.2. Step 2: Kalman gain
Proof. The goal of the Kalman filter is to find the minimum mean-square error estimator. It is obtained by minimizing the error of the a posteriori state given by x t −x t|t . This error is stochastic, hence we are left with minimizing the expected value of the square of the L 2 norm of this vector given by E x t −x t|t 2 . Since the error of the a posteriori state follows a normal distribution with zero mean and covariance given by P t|t , the optimization program is equivalent to minimizing imsart-generic ver. 2014/10/16 file: article_KalmanFilterDemystified.tex date: December 14, 2018 the matrix trace of the a posteriori estimate covariance matrix P t|t .
Expanding out and collecting the terms in equation (A.6), we find: This minimization program is very simple. It is quadratic in the Kalman gain, K t . The optimum is obtained when the derivative with respect to K t is null. The equation for the critical point writes as: ∂ tr(P t|t ) We can solving for K t to find the optimal Kalman gain as follows: A.3.
Step 3: Simplification of the posteriori error covariance formula Proof. We can work on our formula further to simplify equation (A.6) as follows. A trick is to remark that in equation (A.9), multiplying on the right both sides by S t K T t , we find that: Injecting this equality in equation (A.7) and noticing that the last two terms cancel out, we get: This last formula appeals some remarks. It makes Kalman filter dammed simple. It is computationally cheap and much cheaper than the equation (A.7). Thus it is nearly always used in practice. However, this equation is only correct for the optimal gain. If by back luck, arithmetic precision is unusually low due to numerical instability, or if the Kalman filter gain is non-optimal (intentionally or not), equation A.11 does not hold any more. Using it would lead to inaccurate results and can be troublesome. In this particular case, the complete formula given by equation (A.7) must be used.
where W is also distributed as a a multi dimensional normal distribution W ∼ N (0, Ψ) independent of X. This is obviously a simplified version of our Linear Gaussian State Space Model and can provide building blocks for deriving the Kalman filter. Trivially, we can compute the unconditional distribution of Z as a normal distribution whose mean is E[Z] = E[µ + ΛX + W ] = µ. Likewise, it is immediate to compute the unconditional covariance matrix of Z as The covariance between X and Z is also easy to compute and given by Hence, collecting all our results, we have shown that the joint distribution of X and Z is a Gaussian given by We shall be interested in computing the conditional distribution of X given Z. Using lemma provided in B.2, we obtained that X | Z is a multi dimensional Gaussian whose mean is: Using the Woodbury matrix inversion formula (see formula (158) in Petersen and Pedersen (2012)), this is also The two formula are equivalent and one should use the one that is easier to compute depending on the dimension of X and Z. Our partitioned Gaussian lemma B.2 gives us also the conditional variance of X: imsart-generic ver. 2014/10/16 file: article_KalmanFilterDemystified.tex date: December 14, 2018 where in the last equation, we have used again the Woodbury matrix inversion formula (Woodbury variant formula (157) in Petersen and Pedersen (2012)).

B.4. Proof of RTS recursive equations
Proof. We start by writing down the distribution of x t and x t+1 conditioned on z 1 , . . . , z t . Aŝ x t+1|t = F t+1 x t + B t+1 u t+1 and using the fact that the control term B t+1 u t+1 is deterministic, we have Thus the vector distribution of x t|t x t+1|t has its mean given by x t|t x t+1|t and its covariance matrix given by P t|t P t|t F T t+1 F t+1 P t|t P t+1|t Let us now work on the backward recursion. We are interested in computing the distribution of x t , conditioned on x t+1 and z 1 , . . . , z t . Using the Factor analysis model, it is easy to see that E[x t | x t+1 , z 1 , . . . , z t ] =x t|t + P t|t F T t+1 P −1 t+1|t (x t+1 −x t+1|t ) (B.33) Similarly, we have Var [x t | x t+1 , z 1 , . . . , z t ] = P t|t − P t|t F T t+1 P −1 t+1|t F t+1 P t|t (B.34) But using the Markov property that states that E[x t | x t+1 , z 1 , . . . , z T ] = E[x t | x t+1 , z 1 , . . . , z t ] (B.35) and Var[x t | x t+1 , z 1 , . . . , z t ] = Var(x t | x t+1 , z 1 , . . . , z T ] (B.36) We have E[x t | x t+1 , z 1 , . . . , z T ] =x t|t + P t|t F T t+1 P −1 t+1|t (x t+1 −x t+1|t ) (B.37) Similarly, we have Var[x t | x t+1 , z 1 , . . . , z T ] = P t|t − P t|t F T t+1 P −1 t+1|t F t+1 P −1 t|t (B.38) imsart-generic ver. 2014/10/16 file: article_KalmanFilterDemystified.tex date: December 14, 2018 Using conditional properties states in appendix section C.2, we can work the recursion as follows: x t|T E[x t | z 1 , . . . , z T ] (B.39) The latter equation is the basic update equation in the RTS smoothing algorithm. This equation provides an estimate of x t based on the filtered estimatex t|t corrected by the convolution of L t with the error term x t+1|T −x t+1|t that represents the difference between the smoothed estimate of x T and the filtered estimatex t+1|t . The matrix L t can be interpreted as a gain matrix that depends only on forward information. and can be computed in the forward pass.
As for the conditional variance, we havê We can summarize the RTS smoothing algorithm as follows: x t|T =x t|t + L t (x t+1|T −x t+1|t ) (B.53) P t|T = P t|t + L t (P t+1|T − P t+1|t )L T t (B.54) with an initial condition given byx

B.5. Proof of the inverse dynamics
Proof. We have which provides the first result 4.45.
Using the forward dynamics (equation (4.44)) that relates x t+1 and w t+1 , we get that which proves 4.48.
The independence between w t+1 and the past information x t+1 , . . . , x T is obvious for x t+2 , . . . , x T and inferred from the fact that the two processes are Gaussian with zero correlation: E[ w t+1 x t+1 ] = −F −1 t+1 (Q t+1 − Q t+1 ) = 0. Last but not least, using the independence, the forward Lyapunov equation is trivially derived, which concludes the proof.

B.6. Proof of modified BrysonFrazier equations
Proof. Applying the information filter given in proposition 5 with the forward dynamics equation 4.45 gives: where in our series of equation, the third equality follows from the fact that the conditional probability given by Pr [ W = w, Y = y | Y = y] is Pr[W = w | Y = y] for y = y and 0 otherwise and the sixth equality comes from Bayes' rule.
The second equality is trivial and is a consequence of the first and the law of total variance.
Appendix D: EM Algorithm convergence proofs

D.1. First Proof
Proof. Expectation-maximization improves the expected complete log likelihood: Q(θ|θ (t) ) rather than the complete log likehood: log p(X|θ). We have for any Z with non-zero probability p(Z|X, θ) that the log likehood conditional to the parameters θ can be split into two parts: log p(X|θ) = log p(X, Z|θ) − log p(Z|X, θ) .
which states that the marginal likelihood at each step is non-decreasing, and hence concludes the proof.

D.2. Second Proof
Proof. Another proof relies on the fact that the EM algorithm can be viewed as a coordinate ascent method, which is proved to converges monotonically to a local minimum of the function we maximize (see for instance Hastie et al. (2009) or Neal and Hinton (1999)). This works as follows.
Let us consider the function: where q is an arbitrary probability distribution over latent data Z and H(q) is the Entropy of the distribution q. This can also be written as: F (q, θ) = −Div KL q p Z|X (·|x; θ) + log L(θ; x), where p Z|X (·|x; θ) is the conditional distribution of the latent data given the observed data X and Div KL is the KullbackLeibler divergence.
In the EM algorithm, the expectation step is to choose q to maximize F : The maximization step is to choose the parameter θ to maximize F : θ (t+1) = arg max θ F (q (t) , θ) which concludes the proof by showing that the EM algorithm is a coordinate ascent method.