3
$\begingroup$

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?

$\endgroup$
3
$\begingroup$

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}

$\endgroup$
0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.