July 2009 - these now work with the newest version of R. I've marked the main scripts at the top with a version number, which is the date (month.day.year) the script is posted. The latest version of everything is 7.21.2009.
There are three levels (0,1,2) of filtering and smoothing, and they are contained in ss0.R, ss1.R, ss2.R. You can download them individually below, or get them in one file: ssall.R. Keep reading for details and examples.
Each ss*.R is a script for providing the Kalman filter and smoother, the innovations and the corresponding variance-covariance matrices, and the value of the innovations likelihood at the location of the parameter values passed to the script. MLE is then accomplished by calling the script that runs the filter. The model is specified by passing the model parameters. Code for running the examples in the text are provided below.
Level 0 is for the case of a fixed measurement matrix and no inputs; i.e., if At = A for all t, and there are no inputs, then use the code at level 0. If the measurement matrices are time varying or there are inputs, use the code at a higher level (1 or 2). Many of the examples in the text can be done at level 0. Level 1 allows for time varying measurement matrices and inputs, and level 2 adds the possibility of correlated noise processes. The models for each case are (below, x is state, y is observation, and t = 1, ..., n):
Finally, at the bottom of this page you'll find the code for the material that uses switching: §6.8 DLM with Switching and §6.10 Stochastic Volatility Models.
To start, download and source ss0.R. Once you source ss0, you'll have Kfilter0 and Ksmooth0, for running the filter and smoother. The call to Kfilter0 is, Kfilter0(n,y,A,mu0,Sigma0,Phi,cQ,cR) in fairly obvious notation, except that cQ=chol(Q) and cR=chol(R) are the Cholesky decompositions of Q and R. The call to Ksmooth0 is similar. I didn't write the code to include inputs because I just want to keep it simple at level 0.
Now some examples... note that the names of the files are the same as the old versions, but the files are different in that they're using ss0.
- Code to duplicate Example 6.5 [§6.2]: ex65.txt
- Code to duplicate Example 6.6 [§6.3]: ex66.txt
- Code to duplicate Example 6.7 [§6.3]: ex67.txt
- Code to do Example 6.10 via BFGS [§6.5]: ex610.txt
The example in the text uses the EM algorithm. The code here uses optim to do quasi-MLE with slightly different starting values. The model is the same and the results are similar.
For the EM algorithm with measurement matrices that don't change (which means no observations are missing), you can use EM0.R [§6.3], which has to be sourced. EM0 uses the code in ss0. Standard errors can be obtained by computing the Hessian at the location of convergence; this uses the package nlme which has to be downloaded and installed. The call is EM0(n,y,A,mu0,Sigma0,Phi,cQ,cR,max.iter,tol) and it isn't fancy... it just estimates all the parameter matrices whether or not they're fixed. So, for example, if you want to hold a Phi parameter fixed, or hold the initial conditions fixed, then you'll have to edit EM0. Here's an example:
- Code to duplicate Example 6.8 [§6.3]: ex68.txt
First, download and source ss1.R. Once you source ss1, you'll have Kfilter1 and Ksmooth1, for running the filter and smoother. The call to the filter is Kfilter1(n,y,A,mu0,Sigma0,Phi,Ups,Gam,cQ,cR,input), where A is an array with dim=c(q,p,n), Ups is Υ [p × r], Gam is Γ [q × r], and input is the matrix of inputs that has the same row dimension as y (which is n × q, so input is n × r; the state dimension is p). The call to Ksmooth1 is similar. You can set Ups or Gam or input to 0 (zero) if you don't want to use them. Here are some examples... be sure to source ss1 first.
To use the EM algorithm with missing observations or with inputs, you'll need EM1.R, which has to be sourced. The call to EM1 is EM1(n,y,A,mu0,Sigma0,Phi,Ups,Gam,cQ,cR,input,max.iter,tol). Here's an example... be sure to source ss1 too.
- Code to duplicate Example 6.9 [§6.4]: ex69.txt
This is the filter that can be used in bootstrapping state space models. Download and source ss2.R, after which you'll have Kfilter2 and Ksmooth2 for running the filter and smoother. Kfilter2 runs the filter given in Property 6.5. The call to the function is Kfilter2(n,y,A,mu1,Sigma1,Phi,Ups,Gam,cQ,cR,S,input), which is similar to Kfilter1, but where mu1 = E(x1), Sigma1 = var(x1). And the examples:
- Code to duplicate Example 6.12 [§6.7]: ex612.txtNB: The example in the text was done in Gauss (which is similar to Matlab). The results using R are close to, but not the same as, the results in the text. The R optimization routine, optim, is a time killer, so I set the tolerance in the bootstrap estimation at .1 [that is, optim(..., control=list(reltol=.1)) at line 65 of the example]. If you run optim at its default tolerance, the analysis will take hours to complete.
NOTE: ex612.txt was updated on 10/31/2008 due to a typo noticed by Lasse Bork of the Aarhus School of Business, Denmark. The line y.star[4:50] = z[4:50]*xp[4:50] + ... was changed to y.star[4:50] = z[4:50]*xp.star[4:50] + ....
- Code to bootstrap an AR [§6.6-6.7]: ex333.txtHere's an example using Property 6.6 to bootstrap an AR(1) via the state space model; see Example 3.33 for details about the data and the setup. This example is not in the text.
Dynamic Linear Models with Switching: Here's the code for the analysis of the flu data, Example 6.16. Please note that the example in the text was not done in R, so the results are similar, but not the same because different optimization routines were used.
Stochastic Volatility Models: You'll need SVfilter.R, which has to be sourced. The call to the function looks like this: SVfilter(n,y,phi0,phi1,sQ,alpha,sR0,mu1,sR1), where phi0, phi1, sQ (= σw) are state parameters, and alpha, sR0 (= σ0), mu1, sR1 (= σ1) are observation equation parameters as presented in Section 6.10. This is similar to Kfilter0, but it does the special case switching filter for the SV model. SVfilter is a stand alone program. And, of course, an example:
- Code for Example 6.21 [§6.10]: ex621.txt