Start up values for the Kalman filter

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?

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)
#  TRUE


The residuals also match those returned by KalmanRun:

all.equal(kf$resid, v) #  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} 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 <- v / 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.