kalman filtering in R(R语言实现卡尔曼滤波）

【实例简介】
卡尔曼滤波（平方根算法、序贯算法、线性卡尔曼滤波、扩展卡尔曼滤波、区间卡尔曼滤波、系统误差和观测误差有关的卡尔曼滤波）在R语言中的实现，详细讲述了原理和其中工具包的使用，英文原版。
Journal of statistical softwaro computational effort. One possible advantage is that they provide a natural way to specify complete uncertainty about the initial value of a component of the state: we can set the corresponding diagonal term in the information matrix to zero. With a covariance filter, we have to set the corresponding variance in Pol-1 to a"large" number, or else use exact diffuse initialization an option described belo Direct transcription of the equations making up the Kalman filter into computer code is straightforward. I was soon oliced, Chough, thal the resulting programs sullered Irom numerical instability: see for instance Bucy and Joseph(1968). In particular, buildup of floating point crrors in cquation( 8) may eventually yicld non symmetric or non positive delinile Pt matrices. An alternalive lo equation(8)(more expensive Iron the computational point of view)is Pt=(I-KtztPtt-1I(I-KtZt+Kthtkt ith Kt= Plt-1Z1 Ft but cvcn this(" Joscph stabilized form,, scc Bucy and J Joseph (1968), p. 175) Imlay prove insullicient to prevent roundoll' error degeneracy in the liller. A detailed reference describing the pitfalls associated with the numerical implementation of the Kalman filter is Bicrman(1977); scc also Grewal and Andrcws(2001), Chap 6 and Anderson and Moore(1979), Chap. 6. Square root algorithms, that we discuss next, go a long way in improving the numerical stability of the filter 2.2. Square root algorithms Consider a matrix St such that Pt=StSt: square root covariance filter(SRCF) algorithms propagate St instead of Pt with two main benelits(cf. Anderson and Moore(1979),$6.5) (i)Re-constitution of Pt from St will always yield a symmetric non negative matrix, and (ii) The numerical condition of St will in general be much better than that of Pt. If instead of the covariance maTrix we choose to actor the information malrix P 1 we flave a square root information filter algorithm(SRIF) It is casy to producc a rcplaccmcnt for the timc update cquation(6). Consider an orthogonal matrix G such that (Q where M is an upper triangular matrix. Matrix G can be constructed in a variety of ways including repeated application of Householder or Givens transforms. Multiplication of the transpose of (10) by itself produces, taking into account the orthogonality of G 2 2TR MM= TESt-1Sf+RtQt'(Qt) TtPt-1lt+Rt Qt (12) comparison with(6)shows that M is a possiblc choice for Stlt- 1. Likewise, it can be shown (Simon 2006, Section 6.3. 1) that the orthogonal matrix G* such that (Ht+ZtPtt-1Zi)Kt T t-1 t produces in the left-hand side of (13 a block M"that can be taken as St, and thus performs the measurement update. Kalman Filtering in R Although the matrices performing the block triangularization described in equations(10) and (13) can be obtained quite efficiently(Lawson and Hanson 1974; Gentle 2007), clearly the computational effort is greater than that required by the time and measurement updates in equations(6)and(8); square root filtering does have a cost equations above G and G can be chosen so that M and M are Cholesky factors of the corresponding covariance matrices. This needs not be so, and other factorizations are possible. In particular, the factors in the singular value decomposition of Pt-1 can be propagated: see for instance Zhang and Li(1996) and Appendix B of Petris et al. (2009). Also we may notc that timc and mcasurcmcnt updates can be merged in a singlc triangularization (see Anderson and Moore 1979, p. 148, and Vanbegin and Verhaegen 1989) 2.3. Sequential processing In the particular casc where Ht is diagonal, the componcnts of yt ypt arc uncorrelated. We may pretend that we observe one yit at a time and perform p univariate measurement updates similar to(7)-(8), followed by a time update(5)-(6)-see Durbin and Koopman(2001, Scction 6.4)or Andcrson and Moorc(1979) for details The advantage of sequential processing is that Ft becomes 1 x 1 and the inversion of a p p matrix in equation(7)-( 8)is avoided. Clearly, the situation where we stand to gain most from this strategy is when the dimension of the observation vector yt is large Sequential processing can be combined with square root covariance and information filters although in the latter case the computational advantages seem unclear(Anderson and Moore 1979, P. 142; see also Bierman 1977) Asidc from the casc where At is diagonal. sequential proccssing may also bc used when it is block diagonal- in which case we can perform a sequence of reduced dinension updates or when it can be reduced to diagonal form by a linear transformation. In the last case assuming full rank, let Ht 1/2 be a square root of h-1 ultiplying (2\ rough by H,/2 we get yt =di t stat+ et (14 with ElE** =I. If matrix Ht is time-invariant, the same linear transformation wi decorrelate the measurements at all times 2. 4. Smoothing and the simulation smoother The algorithms presented produce predicted atlt-1 or filtered at values of the state vector at. Somctimcs it is of interest to estimate atin for o<ts N, i.C., E[atlyo, .. yN], thc value of the stale vector given all past and future observations. It Turns out that this can be done running once the Kalman filter and then a recursion backwards in time(Durbin and Koopman 2001, Scction 4.3, Harvcy 1989, Scction 3.6) In some cases, and notably for the Bayesian analysis of the state space model, it is of interest to generate random samples of state and disturbance vectors, conditional on the observations y0,..., JN. Friihwirth-Schnatter(1994) and Carter and Kohn(1994) provided algorit,hms to that purpose, improved by de Jong(1995). Durbin and Koopman(2001, Section 1.7)draws on work of the last author, the algorithm they present is fairly involved. Recently. Durbin and Koopman(2002) have provided a. much simpler algorithm for simulation smoothing of the Gaussian state space model; see also Strickland et al.(2009) Journal of statistical Softwarc 2.5. Exact diffuse initial conditions As mentioned above, the Kalman filter iteration needs starting values a-1 and P-1. When nothing is known about the initial state, a customary practice has been to set P-1 with large clements along the main diagonal, to reflect our uncertainty. This practicc has bccn criticized on the ground that "While it can be useful for approximate exploratory work, it, is not recommended for general use, since it can lead to large rounding errors (Durbin and Koopman 2001, p. 101). Grewal and Andrews(2001, Section 6.3.2) show how this may come about when elements of P-1 are set to values" too large relative to the mea surement variances. An alternative is to use an information filter with the initial information matrix(or some diagonal elements of it )set to zero. Durbin and Koopman(2001, Section 5.3) advocate a diffcrcnt approach: scc also Koopman and Durbin(2003)and Koopman(1997) The last reference discusses several alternatives and their respective merits 2.6. Maximum likelihood estimation An inportant special case of the stale space nodel is that in which sOme or all of the matrices Tt, Rt, Zt, Qt and Ht in equations(1)-(2)are time invariant and depend on a vector of parameters, 0 that we seek to estimate Assuming co N(ao. Po)with both ao and Po known, the log likelihood is given by C(8)- logp(yo, .. yN0 N+121g2-2>(gE+r'e) where et= 3t- Ztat.(Durbin and Koopman 2001, Section 7. 2 the result goes back to Schweppe 1965). Except for Ft, all other quantities in(15) are computed when running the Kalman filter (cf. equations (5)-( 8)); therefore, the likelihood is easily computed. (If we resort to sequential processing Section 2.3, the analog of equation(15) does not require computation of determinants nor matrix inversions Maximum likelihood estimation can therefore be implemented easily: it suffices to write a routine computing(15) as a function of the parameters 0 and use R functions such as optim or nlminb to perform the optimization. An alternative is to use the e algorithm which quite naturally adapts to this likelihood maximization problem, but that approach has generally been found slower than the use of quasi-Newton methods, see for instance Shumway and Stoffer(2006),p. 345 3. Kalman filtering in R There are several packages available from the Comprehensive R Archive Network ( CRAN offering general Kalman filter capabilities, plus a number of functions scattered in other packages which cater to special models or problems. We describe five of those packages in chronological order of first appearance on CRAN Kalman Filtering in R 3.1. Package dse Package dse( Gilbert 2011) is the R offering including Kalman filtering that has been longest in existence. Versions for R in the cran archives date back at least to vear 2000(initially its content was split in packages dsel and dse2 Before that, the software existed since at least 1996, at which time it ran on S-PLUS and R. a separate version forR has been in existence since 19 9 8. The version used here is 2009 10-2 dse is a largc package, with fairly cxtcnsivc functionality to dcal with multivariatc timc scrics The author writes While the software does Inally slandard time-series things, it is really intended for doing some non-standard things. Package dse] is designed for working with multivariate time series and for studying estimation techniques and forecasting mode」 dse implements three ($3)classes of objects: TSdata, TSmodel and TSestModel, which stand for data. model and estimated model. Methods take. for instance. a TSmodel and TSdata object to return a TSestModel. There is a wide range of models available, covering just about any variety of the(V)(AR)(MA)(X) family. As regards state space models, two representations are supporled, the 11on-innovalions form Fat-1 I Gut t Qnt (16) Hαt+Ret (17) and the innovations form F +Gut+key (18 t Hat+ Re (19) For the relative advantages of each form. and a much wider discussion on the equivalence of different representations, see Gilbert (1993 The first representalion is similar lo(1)-(2), with Gut taking the place of Ct. There is a neat separation in TSdata objects between inputs ut and outputs yt. Either representation is specified with time-invariant matrices, F,G,etc. There is no provision for time-varying matrices, nor lor missing dala All system matrices can be made of any mixture of constants and parameters to estimate. a vcry uscful default is that all entries arc considered paramcters to estimate, with the exception of those set at 0.0 or 1.0; completely general specification of fixed constants and parameters to estimate in the system matrices can be obtained by specifying a list of logical matrices with the samc sizes and namcs as the systcm matrices, and a logical TRue in place of the constants. Consider. for instance the measurements of the annual How of the nile river at Ashwan, 1871-1970(in datasets). We can set a local level model, I nt yi (2 as follows: For a description of S and $4 classes the reader may turn respectively to Chambers and Hastie(1992) and Chambers(2008) Journal of statistical softwaro R> m1 dse < dse:: SS(F matrix(1, 1, 1),Q= matrix(40, 1, 1) H= matrix(1, 1, 1), R= matrix(130, 1, 1) z0= matrix(0,1,1),P0= matrix(10^5,1,1)) Matrices Q and R in(16)-(17)are 1 x 1 and have been set to values different from 0.0 and implying they are to be estimated. Should we wish to fix the value of Q to 40 and estimate only r, thie specification would become R> m1b dse <-dse:: SS(F matrix(1, 1, 1),Q= matrix(40, 1, 1) H= matrix(1, 1, 1),R= matrix(130, 1, 1) constants list(Q= matrix(TRUE, 1, 1), P0 matrix(TRUE, 1, 1)) z0= matrix(0,1,1),P0= matrix(105,1,1)) (The scope resolution operator: will in general be unnecessary; it is used here because both dse and spir include a function named SS. There is a genera l purpose function l to estimate models, which dispatches to specialized functions according to the type of model. It takes two arguments, a TSmodel and a TSdata(or objects that can be coerced to those). with state-space models, l calls 1. SS, which by default returns a TSestModel(but can be made to rcturn other things, like a computed likelihood, which makes the output from 1.Ss a suitable input lo an optimization routine dse is a nixed language package, wilh soine routines coded in Fortran. In particular, function l SS invokes alternatively R code or Fortran compiled code, according to whether argument compiled is FALSE or TRUE. Both the R and Fortran implementations arc covariance filters following closely the algorithm described in equations(5)-( 8), with(5) substituted in(7)and (6)in( 8), and additional steps of the form PtIt-1-(PtIt-1+Ptt_1)/2 and Pt-(Pt+Pt to prevent rounding crror from destroying the symmetry of Ptlt-1 and Pt. Thc Fortran implementalion uses LAPACK(Anderson el al. 1999)routines Routines for different estimation methods are provided in the package. A convenient function is estMaxLik, which performs maximum likelihood for Gaussian models. The following code estimates the univariate local level model in equations(16)(17)(defined in m1 dse above) with the well known nile data R> data("Nile", package ="datasets") R>m1b dse est <-estMaxLik(m1b dse, ISdata(output Nile)) (15), but rathcr a largc samplc approximation(scc for instance Harvey 1989, Scction jx o However, it is important to realize that the likelihood which is optimized is not that given For Chis reason, package dse is dropped Iron some comparisons in Section 4.2, as the resulls differ from those obtained with other packages. The Kalman filtering routine l SS, however still compared in terms of speed with its counterparts in other packages As is the case with other packages reviewed in the following, dse goes well beyond what has been described. It is geared towards the use of alternative representations and evaluation of different estimation procedures (package EvalEst, Gilbert 2010, provides a complement of Functions for the last purpose). Aside from that, it offers useful utility functions to check stability observability compute the McMillan degree, solve Riccati equations, and more Kalman Filtering in R 3.2. Package sapir The first version of spir (Dethlelsen el al. 2009)on Cra dales back lo 2005; he version we comment upon is version 0.2. 8 of 2009. It does not contain documentation other than thc manual pages, but is woll described in Dcthlofscn and Lundbyc-Christenscn(2006). In addition, the package is used in Cow pertwail and Metcalfe(2009) sPir is a pure R package, i.e., contains no compiled code. It uses an object-oriented approach A Gaussian state space model as in (1)-(2)(but with no deterministic inputs ct, dt nor matrix Rt) is represented by an object of (s3) class SS. There is a constructor function of the same name that yields an SS object when invoked with the required information. The covariance matrices Qt, Ht(Wmat and Vmat in the notation of the package) can be constant matrices or functions, while measurement and transition matrices Zt and Tt(Fmat and gmat in the notation of the package) have to be functions of time, a parameter veclor and possibly a set of covariates. For instance, the local level model for the Nile data in equations(20)-(21 can be spccificd as follows R> data(Nile", package ="datasets") R> 1. spir < spir:: SS(Fmat function(tt, x, phi)t return(matrix(1)) 1, gmat function(tt,x, phi) i return(matrix(1)) F, Vmat function(tt, x, phi)i return(matrix(exp (phi[1]))) f, Wmat function(tt, x, phi) t return(matrix(exp (phi [2]))) + h, y=as matrix(Nile, ncol =1)) R> str(m1 spir List of 19 y num[1:100,1]1120116096312101160116081312301370 list of 1 .s x: logi NA s Fmat function (tt, x, phi) S gmat function(t七,x,phi) s Vmat function (tt, x, phi) S Wmat function (tt,x, phi) :num[1,1]0 S CO :num[1,1]100 S phi num[1:2]11 nt100 int 1 P logi NA s iteration num O Sm logi NA logi NA s mu logi NA Journal of statistical softwaro s likelihood: logi NA logi NA attr(*, class"=chr " SS The constructor yields(in m1 spir) an object with the passed arguments filled in and some other components(which might also have been specified in the call to sS)set to defaults Among othcrs, we scc y, containing the timc scrics, mO and co(initial mcan and covariance natrix of the stale vector), phi(vector of parameters), etc. The only two paraneters on and a2 have been parameterized as exp(1) and exp(o2)to ensure non negativity. A rich set accessor functions allow to got sct individual componcnts in an SS object. For instance, we can set the initial covariance matrix of the state veclor(Co) and the model paraneters(phi R> CO(m1 spir)<-matrix(10"7, 1, 1) R> phi(m1 spir)<-c(9, 7) and then inspect their values, as in R> phi(m1 spir) [1]97 Functions kfilter and smoother perform Kalman filtering and smoothing on the ss object m1. sspir, returning a new SS object with the original information, the filtered (or smoothed) estimates of the state, and their covariance matrices in slots m and C respectively R> m1 spir.f <-kfilter(m1 spir) Aside from the use of function ss to set up a state space model, spir offers a function ssm bascd on the familiar formula notation of functions like lm, glm, ctc. For instance, the local level model for the Nile data could be specified as follows R> mod2 spir <-ssm(Nilo tvar(1), famil aussian co diag (1)*107, phi exp(c(9, 7))) R> class(mod2. sapir) [1]"ssm""1 The returned object is of class ssm. Marker tvar specifies that a right hand side element is time varying(in the example above, the intercept ). The family argument can take values gaussian, binomial, Gamma, poisson, quasibinomial, quasipoisson and quasi, each with several choices of link functions. II gaussian is specilied, the function kfs runs the stall dard Kalman smoother. For a non-Gaussian family, kfs runs an iterated extended Kalman smoother: it iterates ovcr a scqucncc of approximating Gaussian modcls Function ssm provides a very convenient user interface, even if it can only be used for uni- variate time series. The use of ssm is further enhanced by ancillary functions like polytrig polytime, season and sumseason, which help to define polynomial trends and trigonometric polynomials as well as seasonal dummies Kalman filtering in r sapir includes other useful functions for: simulation from a(Gaussian) state-space model (recursion), simulation smoothing(simulate) and EM estimation(EMalgo Regarding the algorithmic approach followed, sspir implements a covariance filter, with the faster( 8)rather than the more stable( 9) update for the state covariance. Somewhat surpris- ingly, missing values are not a. llowed All in all, the design of the package emphasizes usability and flexibility: the fact that system matrices can(or need) be specified as functions of time, covariates and parameters affords unparallelled fexibility. This flexibility has a price, however, as every single use of a matrix forces its recreation from scratch by the underlying function. As we will see in Section 4, this coupled with the all-R condition of the package, severely impacts performance 3.3. Package dlm Vcrsion 0.7-1 of package dIm was submitted to Cran in August 2006. The version we revicw is 1. 1-1(see Petris 2010). Other than the manual pages, an informative vignette accompanies the package. A recent book, Petris et al. (2009), provides an introduction to the underlying thcory and a wealth of cxamples on the usc of the package Like spir. the state-space model considered is a simplified version of (1)-(2), with no in tercepts Ct, dt and no Rt matrix. The emphasis is on Bayesian analysis of dynamic linear models(DI Ms), but the package can also be used for Kalman filtering/ smoothing and maxi- mum likelihood estimation. The design objectives as stated in the concluding remarks of the vignette Petris(2010)are Flexibility and numcrical stability of the filtering, smoothing, and likelihood evaluation algorithns. These two goals are somewhat related, since naive imple- mentations of the Kalman filter are known to suffer from numerical instabilit for gcncral DLMs. Thcrcforc, in an cnvironment whcrc the uscr is frcc to specify virtually any lype of DLM, it was inportant to try to avoid as nuch as possible the aforementioned instability problems In order to cope with numerical instability, the author of the package has chosen to implement a forin of square root lilter described in Zhang and Li (1996), which propagates factors of the singular value decomposition(SVD) of Pt. It is claimed that the algorithm is more robust and general than standard square root filters. The computation is done in C, with extensive use of LAPACK routines. IniTial conditions aol-1 and Pol-1 have to be specilied (or the defaults accepted ) There is no provision for exact initial diffuse conditions Asidc from the heavy duty compiled functions, the package provides a rich sct of interface functions in R. lo help with the specification of slale-space nodels. A line-invariant Inodel an object of (S3)class "dlm". A general constructor is provided, not unlike SS in package sapir, which returns a dlm object whon invoked with arguments V, W(covariance matrices of the measurement and state equations, respectively ). FF and GG (measurement equation matrix and transition matrix respectively ), and mo, co(prior mean and covariance matrix of the state vector). For instance R>m1.dm<-dm(FF=1,V=0.8,GG=1,W=0.1,m0=0,C0=100) creates a dlm object representing the model in(20)(21)with a2=0.8 and 02=0.1.The class and structure of the object is 【实例截图】
【核心代码】