miller1997kalman

Kalman �ltering for linear systems with coeÆcients driven by a hidden Markov jump process* Boris M. Miller Institute for Information Transmission Problems Russian Academy of Sciences B. Karetny per. 19, 101447 GSP-4 Moscow, Russia and Wolfgang J. Runggaldier Dipartimento di Matematica Pura e Applicata Universit�a di Padova Via Belzoni 7, I - 35131 - Padova, Italy Abstract. We present an explicit algorithm for the solution of the nonlinear �ltering problem concerning a linear, partially observed di�usion-type model, where the coef- �cients are driven by a hidden Markov jump process, of which one can observe the occurence of a jump. The (recursive) �lter is of a "branching" type with its (�nite) dimension growing for each jump of the hidden Markov process. Key words: linear and nonlinear �ltering, �nite dimensional �lters, hidden Markov models, Markovian switching coeÆcients, Kalman �ltering. 1. INTRODUCTION The �ltering problem for di�usion processes, where the coeÆcients depend on an unobservable Markov jump process, has many applications in various �elds such as telecomunications, image processing, tracking of maneuvering objects, and many others. The traditional problem setting considers a partially observable process (x; y; �) = f(x t ; y t ; � t ; t � 0g, generally described by a linear (in x) Ito equation dx t = a t (� t )x t dt+ b t (� t ) dw (1) t ; x 0 � N (m 0 ;� 0 ); dy t = A t (� t )x t dt+ B t dw (2) t ; y 0 = 0; (1:1) where x t and y t are vectors, and � t is a right continuous Markov process with �nite num- ber of states � t 2 f1; � � � ; ng, having matrix of transition intensities � t = f� i;j t g i;j=1;���;n ; w (i) = fw (i) t ; t � 0g; i = 1; 2 are standard vector Wiener processes, a t (�); b t (�); A t (�); B t are matrix functions of appropriate dimensions, and x 0 ; � t ; w (1) ; w (2) are independent. * Work partially performed during a stay of B.M. Miller in Padova, supported by INTAS Grants 93-2622 and 94-697. 1 It is known that, in general, the �lter equations for (x; �), based on the observations y, namely the equations for Efx t j y t 0 g and EfIf� t = ig j y t 0 g where I(�) denotes the indicator function, are nonlinear and in�nite-dimensional (see [1],[8]) just as in the general nonlinear case ([10]). It is also known that the problem of estimating the state of a Markov process on the basis of continuous observations admits a solution in terms of a �nite-dimensional �lter system only if the observation process y t admits a representation of the form y t = Z t 0 C(� s ) ds+ Z t 0 B s dw s ; (1:2) where C(�) is a function of the current state of the process � t , and w t is a Wiener process independent of � t (see [1],[8]). In [8] the authors consider estimation problems for more general unobserved state processes and this latter monograph covers all known cases admitting a �nite dimensional �lter; still, it does not include models of the form (1.1). On the other hand, in the discrete time case the �ltering problem for models of the type (1.1) can always be reduced to a �nite-dimensional recursive system of equations, which, as shown in [12], has dimension n K , where n denotes the number of possible states of the Markov chain and K the maximal value of the discrete time index. The reason is that in the generic time instant k � 1 the set of possible trajectories of the Markov chain that lead to any of the current n possible states is equal to n k and, consequently, in order to obtain a closed system of �lter equations, it suÆces to have a procedure to compute the conditional probabilitites for each trajectory. As shown in [12] (see also [3]), for given transition probabilities these conditional probabilities can be computed with the use of Bayes formula and as a result one obtains a �nite-dimensional system of �lter equations. These considerations are also at the basis of the present study, where, as in [5], we extend the observation process by adding to it the information about the instants of state transition of � t . This means, we can observe that something has happened, but we do not know exactly what (for a justi�cation of this see [5]). If therefore one can observe the process N t , describing the number of these state transitions up to the generic moment t > 0, then for each given state � t = i; i = 1; � � � ; n; the number of possible trajectories leading to this state is equal to n(n� 1) N t . The problem then consists only in computing the conditional probabilities for each of the possible trajectories which, generally speaking, leads to a smoothing problem. Using the information given by the observation of the process N t , one can construct a system of branching equations for the computation of these probabilities where the branching points coincide with the moments of jump of the process N t . The outline of the paper is as follows. In the next section 2 we formulate the problem and give the basic de�nitions. In the third section we derive a �nite-dimensional system of equations for the estimation of the state of the process � t on the basis of the observations of the process N t together with a continuous process y t of the type as in (1.1) which we may rewrite as y t = Z t 0 g s (!; � s )ds+ w t (1:3) 2 for a suitable adapted process g s (!; � s ) and a Wiener process w t . For simplicity we shall assume y to be scalar and w standard Wiener, the extensions to the vector case implying only technical complications. These equations contain terms of the form Efg t (!; � t ) j F Y t g, where F Y t is the �� algebra generated up to time t by the process-pair Y t = (N t ; y t ), and that can be represented as Efg t (!; � t ) j F Y t g = EfEfg t (!; � t ) j F �;Y t g jF Y t g = X � t 0 g^ t (Y t 0 ; � t 0 )Pf� t 0 j F Y t g; (1:4) where the sum extends over the (�nite) set of admissible trajectories � t 0 that correspond to the observation of the process N t 0 = fN s ; 0 � s � tg, and Pf� t 0 j F Y t g are the conditional probabilities of these trajectories. In the fourth section we derive the system of branching forward equations for the computation of Pf� t 0 j F Y t g leading to a �nite- dimensional (for each t > 0 such that N t < +1) closed system of �lter equations, provided we have an explicit way to compute g^ t (Y t 0 ; � t 0 ). Finally, in the �fth section we obtain the solution of the original problem in the form of a �nite-dimensional (w.P.1) system of �lter equations depending on g^. The original �ltering problem is thus fully solved every time one can get an explicit expression for g^. In the linear case, as in model (1.1), g^ can be computed via a Kalman �lter. 2. PROBLEM SETTING Given a probability space ( ;F ; P ) with a right continuous �ltration F t , consider a nonstationary Markov process � t with �nite state space � t 2 f1; � � � ; ng having transition intensity matrix � t = f� ij t g i;j=1;���;n with deterministic and continuous � ij t . Putting X i t := If� t = ig, where IfAg is the indicator function of the set A, consider the vector X t := fX 1 t ; � � � ; X n t g. Then (see [10; Ch.9]) the process X = fX t ; t � 0g admits the representation X t = X 0 + Z t 0 � � s X s ds+M t ; (2:1) where M = fM t ; t � 0g with M t := fM 1 t ; � � � ;M n t g is a square integrable (P;F t )� martingale with quadratic variation (see [4], [9]) hMi t = � Z t 0 [� � s diagX s + diagX s � s ] ds+ Z t 0 diag (� � s X s ) ds; (2:2) where diagX denotes the matrix with diagonal entries X 1 ; � � � ; X n and � the transpose. Suppose that one can observe the two-component process Y = fY t ; t � 0g, where Y t = � y 1 t = N t y 2 t = R t 0 g s (!; � s )ds+ w t � : (2:3) The process N t in (2.3) denotes the number of state transitions of the Markov process � t that occur up to the moment t, w = fw t ; t � 0g is a standard F t � Wiener process, and g s (!; �) is a given scalar F t � adapted random function, which in particular may correspond to an observation equation as in (1.1) (see also (1.3)). Consider also the 3 sequence of �� algebras F Y t := �fY s ; 0 � s � tg � F t and the stopping times f� k g; k = 0; 1; � � � ; such that � 0 = 0; � k = inffs : N s = kg. Since f� k < tg 2 F Y t , the random variables � k are F Y t � measurable and so with each stopping time � k one can associate the �� algebra F Y � k , generated by the events A : A T f� k � tg 2 F Y t . Remark 2.1. It turns out that F Y � k = F Y � � k . In fact, F Y � k = F Y � � k W �f�N � k g. On the other hand, Pf�N � k = 1g = Pf g = 1 and so F Y � � k di�ers from F Y � k only by events of probability zero. Therefore, after completion of F Y � � k by the null sets, we get the required equality, which also implies F Y t = F Y t� ; 8t: We introduce the following de�nitions : for a random variable h(!), the symbol ^ h t denotes ^ h t = Efh(!) j F Y t g; if h = h t (!); then ^ h t = Efh t (!) j F Y t g. Here we keep the notations from [8]. Notice also that g t (!; � t ) = n X i=1 g t (!; i)If� t = ig = n X i=1 g t (!; i)X i t = hg � t (!); X t i: (2:4) Our purpose is to derive an equation for ^ X t = EfX t j F Y t g: (2:5) Notice that for a model of type (1.1), combining our results with a Kalman �lter (see section 5, in particular (5.4)), from ^ X we can also obtain x^ t = Efx t j F Y t g. 3. DERIVATION OF THE FILTER EQUATIONS To derive the �lter equation in closed form we need the following Theorem 3.1. Let the random processes g t (!; i) be square integrable for any t and i = 1; � � � ; n and continuous in t for �xed ! and i. The process ^ X t = EfX t j F Y t g satis�es then the equation d ^ X t = � � t ^ X t� dt+ 1 t� (dy 1 t + h� � t ; ^ X t� idt) + 2 t (dy 2 t � d z }| { hg � t (!); X t� idt); (3:1) where 1 t =Ifh� � t ; ^ X t� i 6= 0g fdiag(� t )� � � t g ^ X t� h� � t ; ^ X t� i � ^ X t� 2 t = d z }| { diag(g � t (!))X t� � I d z }| { hg � t (!); X t� i ^ X t� ; (3:2) with I standing for the unit matrix of dimension (n � n) ,h�; �i denoting the scalar product, � � s := (� 11 s ; � � � ; � nn s ) ; g � s (!) := (g s (!; 1); � � � ; g s (!; n)): (3:3) and with diag(a) denoting the diagonal matrix with elememts a 1 ; � � � ; a n ; and � meaning transpose. Remark 3.2 : Notice that what Theorem 3.1. really does is generalizing and then combining results from [8; part III] and [2] as well as [6] (concerning 1 ) with 4 results from [10] (concerning 2 ). There is a priori no clear way to do this directly. The results of the theorem can be shown to follow from the theory of �ltration of general semimartingales [11, Th.10.1, Ch.4]. A direct proof is possible but, since it has only technical interest, it is moved to Appendix. The �lter equations (3.1), (3.2) contain terms of the form d z }| { diag(g � t (!))X t� = 0 @ d z }| { g i t (!)X i t� 1 A i=1;���;n = 0 @ d z }| { g i t (!)If� t� = ig 1 A i=1;���;n and d z }| { hg � t (!); X t� i = n X i=1 d z }| { g i t (!)If� t� = ig Analogously to (1.4), recalling also Remark 2.1., we obtain d z }| { g i t (!)If� t� = ig = X � � N t �1 0 g^ t � Y t 0 ; � � N t �1 0 ; � � N t = i � P � � � N t �1 0 j F Y t� ; (3:4) where � � N t �1 0 is any trajectory of the process � t from 0 to � N t �1 compatible with the observations of the process y 1 t = N t , namely having the switching points at the random time instants � k and keeping constant values between them. Generally speaking, the computation of an expression of the form P � � � N t �1 0 j F Y t� leads to a smoothing problem for the Markov process � t . We shall however be able to suggest a recursive procedure allowing the computation of these probabilities by means of a system of forward equations that are analogous to the �lter equations (3.1),(3.2). 4. ALGORITHM FOR THE COMPUTATION OF THE CONDITIONAL PROBABILITY OF A GIVEN TRAJECTORY OF THE HIDDEN MARKOV PROCESS As underlying basis for the algorithm, that will be described below, we give the following auxiliary result Proposition 4.1. Given a time instant t and a sequence i 0 ; � � � ; i N t �1 , let C(!) = Q N t �1 k=0 If� � k = i k g and F Y t := F X;Y � � N t W F Y t for � N t � t < � N t +1 . Then, for � N t � t < � N t +1 , we have that X t C := EfX t C j F Y t g is solution of the following �lter equation X t C = X � N t C + Z t � N t � � s X s Cds+ Z t � N t 1 s� Cd ~ y 1 s + Z t � N t 2 s Cd ~ y 2 s (4:1) where X � N t C =E n X � N t C j F Y � N t o = P n � � N t = i N t j � � N t �1 = i N t �1 ; � � � ; � � 0 = i 0 ; F Y � N t o i N t =1;���;n (4:2) 5 ~ y 1 t = y 1 t � y 1 � N t + Z t � N t h� � s ; X s� Cids; (4:3) ~ y 2 t = y 2 t � y 2 � N t � Z t � N t hg � s (!); X s Cids; (4:4) 1 t = Ifh� � t ; X t� Ci 6= 0g fdiag(� t )� � � t )gX t� C h� � t ; X t� Ci �X t� C; (4:5) 2 t = diag(g � t (!))X t� C � Ihg � t (!); X t� CiX t� C; (4:6) with hg � t (!); X t� Ci = n X i=1 g^ t � Y t 0 ; � � N t �1 0 ; � � N t = i � X i t� C; (4:7) diag(g � t (!))X t� C = � g^ t � Y t 0 ; � � N t �1 0 ; � � N t = i � X i t� C � i=1;���;n ; (4:8) where � � N t �1 0 represents the given sequence i 0 ; � � � ; i N t �1 and g^ t � Y t 0 ; � � N t �1 0 ; � � N t = i � = E n g t � !; � � N t = i � j F Y t� o : (4:9) Proof : The proof of this statement can be obtained in complete analogy to that of Theorem 3.1.(see Appendix), with a suitable adjustment of (2.1) to represent the process X t C on [� N t ; � N t +1 ) and taking also into account that, on one hand, by properties of the indicator function, C 2 = C, on the other that C is F Y � N t � measurable. To determine 1 t ; 2 t one obtains then an equation analogous to (A.9). The Proposition leads to a closed system of �lter equations for X t C on the generic interval [� N t ; � N t +1 ) , that can be solved provided one has the knowledge of the initial condition X � N t C as well as of g^ t . While the problem of determining g^ t will be discussed in the next section 5, the following result shows that, with the help of the system of �lter equations (4.1), one can also obtain the initial condition for the solution of the �lter equations for X t C on the following interval [� N t +1 ; � N t +2 ). Proposition 4.2. Under the same assumptions as for Proposition 4.1, the matrix of conditional transition probabilities of the Markov jump process � t at the moment � N t +1 is D � � � N t +1 ; X � � N t +1 C � = 8 > < > : I; if h� � � N t +1 ; X � � N t +1 Ci = 0 diag(� � N t +1 )�� � � N t +1 h� � � N t +1 ;X � � N t +1 Ci ; if h� � � N t +1 ; X � � N t +1 Ci 6= 0 : (4:10) Proof : The solution of the �lter equations (4.1) at time � N t +1 is X � N t +1 C = 1 � � N t +1 C +X � � N t +1 C = D � � � N t +1 ; X � � N t +1 C � X � � N t +1 C (4:11) 6 which de�nes implicitly D as the matrix of conditional transition probabilities. We show now that it is given by (4.10). To this e�ect notice that the elements of the vector X � � N t +1 C are (taking into account Remark 2.1) X i � � N t +1 C = P n � � N t = i; � � N t �1 = i N t �1 ; � � � ; � � 0 = i 0 j F Y � � N t +1 o = P n � � N t = i; � � N t �1 = i N t �1 ; � � � ; � � 0 = i 0 j F Y � N t +1 o : (4:12) On the other hand, for the components of the vector X j � N t +1 C we have X j � N t +1 C = P n � � N t +1 = j; � � N t �1 = i N t �1 ; � � � ; � � 0 = i 0 j F Y � N t +1 o = n X i=1 P n � � N t +1 = j j � � N t = i; � � N t �1 = i N t �1 ; � � � ; � � 0 = i 0 ; F Y � N t +1 o � P n � � N t = i; � � N t �1 = i N t �1 ; � � � ; � � 0 = i 0 j F Y � N t +1 o ; (4:13) where the term for i = j is zero w.P.1. From here, taking (4.11), (4.12) into account, we obtain the following representation for the elements of the matrix D D ji � � � N t +1 ; X � � N t +1 C � = P n � � N t +1 = j j � � N t = i; � � N t �1 = i N t �1 ; � � � ; � � 0 = i 0 ; F Y � N t +1 o = P n � � N t +1 = j j � � N t = i; � � N t �1 = i N t �1 ; � � � ; � � 0 = i 0 ; F Y � N t +1 o (4:14) since all information about the past history of the process � t up to the moment � � N t +1 is contained in the expression representing the conditioning. This accomplishes the proof of the Proposition. Since on the interval [� 0 ; � 1 ) the initial condition is given by the initial distribu- tion Pf� � 0 = ig i=1;���;n , with the use of Propositions 4.1 and 4.2 we can now derive a (forward) recursive �ltering procedure. We obtain in fact the following Algorithm for the solution of the general �ltering problem - First de�ne the set of random quantities: - on the interval [� 0 ; � 1 ) let C i 0 := If� � 0 = i 0 g; i 0 = 1; � � � ; n, - on the generic interval [� k ; � k+1 ) let C i k ;i k�1 ;���;i 0 := If� � k = i k ; � � � ; � � 0 = i 0 g; i k 6= i k�1 ; i k�1 6= i k�2 ; � � � ; i 0 = 1; � � � ; n. In this way on the k�th interval we obtain n(n� 1) k di�erent random quantities. - The algorithm now starts on the interval [� 0 ; � 1 ), where the elements of the vector X t C coincide with the random variables C i 0 , namely X i 0 t = C i 0 . - The solution of n �lter equations of the type (4.1) leads to the conditional probabilities EfX t C j F Y t g = Pf� � 0 = i 0 j F Y t g i 0 =1;���;n ; for t < � 1 ; 7 - at the time point t = � 1 , in accordance with Proposition 4.2, we obtain n(n�1) initial conditions for the solution of n(n� 1) �lter equations of type (4.1) for the conditional probabilities E � C i 1 i 0 j F Y t � i 0 ; i 1 = 1; � � � ; n i 1 6= i 0 � = E � X i 1 t C i 0 j F Y t � i 0 ; i 1 = 1; � � � ; n i 1 6= i 0 � - The algorithm then proceeds analogously in the obvious way over the various in- tervals [� k ; � k+1 ). In this way, recalling also Remark 2.1., we have obtained a recursive algorithm for the computation of the conditional probabilities Pf� t 0 j F Y t� gthat, by (3.4), are needed for the solution of the general �lter equations (3.1), (3.2). For this recursive algorithm we need at every time point t 2 [� N t ; � N t +1 ) to solve n(n� 1) N t �lter equations of type (4.1), and then to compute the coeÆcients in the equations (3.1), (3.2) by means of (3.4), using the current values of g^ t (Y t 0 ; � t 0 ) as well as of Pf� t 0 j F Y t� g. Remark 4.3 In the case n = 2 the system of �lter equations is