Start up values for the Kalman filter

Question

I am trying to understand how the start up values (initialisation) are calculated in the Kalman filter. As an example, I simulated the MA(2) model below.

y <- arima.sim(n = 120, model = list(order = c(0,0,2), ma = c(0.6,0.4)))
fit <- arima(y, order = c(0,0,2), include.mean = FALSE)

Extracted matrices needed for the Kalman Filter.

sigma2 <- fit$sigma2
T <- fit$model$T
R <- c(1,fit$coef)
Z <- t(fit$model$Z)
I <- diag(9)

And calculated start up values for predicted state vector a_pred1 and variance P_pred1 as below,

A <- solve(I-kronecker(T, T)) 
vecc <- A%*%vec(R%*%t(R))
P_pred1 <- as.matrix(cbind(vecc[1:3],vec[4:6],vec[7:9]))
a_pred1 <- c(0,0,0)

But I do not get the same start up values as produced by KalmanRun.

F <- as.numeric(Z%*%P_pred1%*%t(Z))
v <- as.numeric(y[i]-Z%*%a_pred1)
a_filt1 <- a_pred1 + P_pred1%*%t(Z)*(F^-1)*v   
kfrun <- KalmanRun(y, fit$model)
cbind(kfrun$states[1,],a_filt1)
          [,1]      [,2]
[1,] 1.1964116 1.1964116
[2,] 0.7583928 0.6427514
[3,] 0.5353670 0.2493952

Can any one explain what I am doing wrong in my calculations?

Answer 1

The initialisation is done by makeARIMA and sets the initial state vector equal to zeros. Do not get confuse with fit$mod$a, which returns the contemporaneous state vector at the last iteration, not the first one.

It seems that you are obtaining the prediction error, v, before running the prediction step, a_filt1.

---

The Kalman filter recursions that reproduce the output from KalmanRun can be implemented as follows (based on source code of KalmanRun and package KFKSDS).

set.seed(123)
y <- arima.sim(n = 120, model = list(order = c(0,0,2), ma = c(0.6,0.4)))
fit <- arima(y, order = c(0,0,2), include.mean = FALSE)
# matrices of the state space representation of the model
ss <- fit$model
# storage for the contemporaneous state estimates and residuals
a.upd <- matrix(nrow = length(y), ncol = length(ss$a))
v <- rep(NA, length(y))
# initial state vector and corresponding covariance matrix
a.upd <- rbind(rep(0, length(ss$a)), a.upd)
P.upd <- ss$P
# Kalman filter recursions
for (i in seq_along(y))
{
  # prediction
  a.pred <- ss$T %*% a.upd[i,]
  P.pred <- ss$T %*% P.upd %*% t(ss$T) + ss$V
  # prediction error
  v[i] <- y[i] - ss$Z %*% a.pred
  # variance of prediction error, 'f'
  M <- crossprod(P.pred, ss$Z)
  f <- drop(ss$Z %*% M + ss$h)
  # update of state vector and its covariance matrix
  a.upd[i+1,] <- a.pred + M * v[i] / f
  P.upd <- P.pred - tcrossprod(M) / f
}
# remove the initial state vector containing zeros
a.upd <- a.upd[-1,]

Now compare with KalmanRun (remember setting the initial state vector to zeros):

ss$a <- rep(0, length(ss$a))
kf <- KalmanRun(y, ss)
head(cbind(kf$states[1,], a.upd[1,]))
#           [,1]      [,2]
# [1,] 1.1964116 1.1964116
# [2,] 0.7150462 0.7150462
# [3,] 0.3612108 0.3612108
all.equal(kf$states, a.upd)
# [1] TRUE

The residuals also match those returned by KalmanRun:

all.equal(kf$resid, v)
# [1] TRUE

However, this does not match the residuals from the fitted model:

head(cbind(residuals(fit), v))
# [1,]      0.9941325  1.1964116
# [2,]      0.2668942  0.1986161
# [3,]      0.3804713  0.3151606
# [4,]      1.5094031  1.5725179
# [5,]      0.5247281  0.5066900
# [6,]     -1.0697975 -1.0800751

Some further details remain to be checked to fully reproduce the fitted model, but I think this will help you to reproduce KalmanRun.

---

For completeness I show the equations of the Kalman filter discussed above (Durbin and Koopman, 2001 Time Series Analysis by State Space Methods Section 4.2):

\begin{eqnarray} \begin{array}{lll} a_{t+1} = T a_{t|t} & P_{t+1} = T P_{t|t} T^\top + V \\ v_t = y_t - Z a_t & M_t = P_t Z^\top & f_t = Z M_t + H \\ a_{t|t} = a_t + M_t v_t / f_t & P_{t|t} = P_t - M_t M_t^\top / f_t \\ t = 1,2,\dots,n. \end{array} \end{eqnarray}

Answer 2

I would have added this answer as a comment to javlacalle's excellent answer above but I do not have enough reputation.

javlacalle's answer is mostly correct, however, there are a few details that are missing.

The reason why the residuals from R's arima fit and javalacalle's fit do not match is because the residuals R uses are, from Durbin & Koopman, 2012, 2.12.1, equation 2.65, "standardised one-step ahead forecast errors"

$$e_{t} = \frac{v_{t}}{\sqrt{F_{t}}},\qquad t = 1,\dots,n$$

The unscaled first residual, $v_{1}$, is always equal to $y_{1}$ since $a_{1}$ is always set equal to a vector of zeros, $Z$ is nothing more than a selection vector containing a one in the first position and zeros everywhere else, and, as javalacalle points out, $v_{1} = y_{1} - Za_{1}$.

Further, in the R code provided lines 5 and 11 should be replaced with,

ss <- makeARIMA(phi = numeric(), theta = c(0.6, 0.4), Delta = numeric())
...
P.upd <- ss$Pn

respectively, and,

v[1] <- v[1] / sqrt(f)

should be inserted after line 22. The reason for the first change is that R uses ss\$Pn from a call to makeARIMA to calculate f while fit\$model represents the final state as opposed to the initial state of the model.

For additional information, the scaling of the residuals by $F_{t}$ can be seen on line 719 (as of July 2016) of arima.c located here, while the construction of the ARIMA state space model can be found here.