I really want to understand how the math is working here. I am trying to get the standard error of the fitted values for a time series regression model. In the non-time series regression, I know I can take the transpose of the data multiplied by the variance - covariance matrix of the model coefficients and then multiply by the data values again to get the standard errors of the fitted values.

But I'm not sure how to do this when I am including an autoregressive term.


Response variable

Sablects <- rnorm(10)


my.xreg <- cbind(rnorm(10),rbinom(10,1,0.5))

In my actual data, values are normalized so I set the intercept equal to zero here.

m4<-arima(Sablects, order=c(2,0,0),fixed=c(0,NA,0,NA,NA),xreg=my.xreg) 

The predict function will give me standard errors on my in-sample prediction (the fitted values of my model).

my.se <- predict(m4, newxreg = my.xreg, n.ahead = 10)$se         


Now to compare the output of my.se, I want to do this mathematically but I don't know what to use for the values of the ar2 term. I use 1's as a placeholder to demonstrate that my output does not equal the values from my.se above

C <- cbind(rep(1, nrow(my.xreg)), my.xreg[, 1], my.xreg[, 2])


I think this value should equal the first value in my.se, but is not producing the same value as my.se

sqrt(t(C[1, ]) %*% vcov(m4) %*% C[1, ])

Also, I'm not so great with matrix multiplication but here is my work around for getting all of the se values.

se.output <- matrix(nrow=nrow(C))

Specify that the max number of i is equal to number of rows of C.

  for(i in 1:nrow(C)){

    # Loop through your multiplication for each row (i) of `C`. For each iteration, save the new data into the new row of se.output

    se.output[i] <- sqrt(t(C[i, ]) %*% vcov(m4) %*% C[i, ])  

  • 2
    $\begingroup$ This doesn't really sound like a programming question. If you're not sure how certain statistical values are calculated, that's really a general statistical question and not specific to R or any programming language. So perhaps you could edit your question to make it a more specific programming question, or consider posting to Cross Validated instead. $\endgroup$
    – MrFlick
    Sep 7, 2014 at 16:15
  • $\begingroup$ you'll need to dig into forecast:::predict.Arima, and from there into KalmanForecast, and from there into the underlying C code in KalmanFore at github.com/wch/r-source/blob/… $\endgroup$
    – Ben Bolker
    Sep 7, 2014 at 16:24

1 Answer 1


Closed-form expressions can be obtained for the variance of the prediction errors in an ARMA model. Evaluating $Var(y_{T+i} - E(y_{T+i})) = Var(e_{T+i})$ for the expression of an ARMA model defined for the series $y_t$, $t=1,2,...,T$, we find the variance of the $i$-steps-ahead prediction error.

For example, it can be checked that in an AR(1) model: \begin{eqnarray} \begin{array}{ll} Var(e_{T+1})) = \sigma^2_\epsilon \\ Var(e_{T+2})) = \sigma^2_\epsilon (1 + \phi^2) \end{array} \end{eqnarray}

In a MA(1): \begin{eqnarray} \begin{array}{ll} Var(e_{T+1})) = \sigma^2_\epsilon \\ Var(e_{T+2})) = \sigma^2_\epsilon (1 + \theta^2) \end{array} \end{eqnarray}

In an ARMA(1,1): \begin{eqnarray} \begin{array}{ll} Var(e_{T+1})) = \sigma^2_\epsilon \\ Var(e_{T+2})) = \sigma^2_\epsilon (1 + (\phi - \theta)^2) \end{array} \end{eqnarray}

where $\sigma^2_\epsilon$ is the variance of the disturbance term, $\phi$ is the AR coefficient and $\theta$ is the MA coefficient.

The Kalman filter can be used to do this operations for a general ARMA. For simplicity let's take a series simulated from an AR(2) model and fit an AR(2) model (without external regressors and fixed parameters):

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

The forecasts and their standard errors are the following:

res1 <- predict(fit, n.ahead = 10)
#Time Series:
#Start = 121 
#End = 130 
#Frequency = 1 
# [1]  0.98712621  0.63595091 -0.02690017 -0.26972965 -0.14517678  0.02388100
# [7]  0.07183011  0.03197952 -0.01022221 -0.01869105
#Time Series:
#Start = 121 
#End = 130 
#Frequency = 1 
# [1] 0.9006326 1.0402947 1.0419657 1.0697430 1.0760609 1.0764649 1.0783412
# [8] 1.0786303 1.0786862 1.0788097

For a general ARMA model, the standard errors of the forecasts are obtained as the multiplication of the variance of the prediction errors returned by KalmanForecast by the estimated variance of the disturbance term (and taking the square root):

KFvar <- KalmanForecast(n.ahead = 10, fit$model)[[2]]
res2 <- sqrt(KFvar * fit$sigma2)
# [1] 0.9006326 1.0402947 1.0419657 1.0697430 1.0760609 1.0764649 1.0783412
# [8] 1.0786303 1.0786862 1.0788097
all.equal(as.vector(res1$se), res2)
#[1] TRUE

Most of the maths that you are interested in are carried out by KalmanForecast. In order to show those operations, we first need to define the matrices of the state space representation of the ARMA model (object ss) and run KalmanSmooth to obtain the state vector and its covariance matrix:

n.ahead <- 10
ss <- list(Z = rbind(fit$model$Z), T = fit$model$T, 
  H = fit$model$h, Q = fit$model$V, 
  a0 = fit$model$a, P0 = fit$model$P)
kf <- KalmanSmooth(y, fit$model)

Then the predictions and the variances are obtained as follows (based on function predict.stsmSS in package KFKSDS:

m <- ncol(ss$Z)
a <- matrix(nrow = n.ahead + 1, ncol = m)
P <- array(NA, dim = c(m, m, n.ahead + 1))
n <- length(y)
a[1,] <- if (m == 1) kf$smooth[n] else kf$smooth[n,]
P[,,1] <- if (m == 1) kf$var[n] else kf$var[n,,]
pred <- var <- rep(0, n.ahead)
for (i in seq_len(n.ahead))
  ip1 <- i + 1
  a[ip1,] <- ss$T %*% a[i,]
  P[,,ip1] <- ss$T %*% P[,,i] %*% t(ss$T) + ss$Q
  pred[i] <- ss$Z %*% a[ip1,]
  var[i] <- ss$H + sum(diag(crossprod(ss$Z) %*% P[,,ip1]))

The same forecasts and standard errors as above are obtained:

# [1]  0.98712621  0.63595091 -0.02690017 -0.26972965 -0.14517678  0.02388100
# [7]  0.07183011  0.03197952 -0.01022221 -0.01869105
all.equal(as.vector(res1$pred), pred)
#[1] TRUE
res3 <- sqrt(var * fit$sigma2)
# [1] 0.9006326 1.0402947 1.0419657 1.0697430 1.0760609 1.0764649 1.0783412
# [8] 1.0786303 1.0786862 1.0788097
all.equal(res2, res3)
#[1] TRUE

Alternatively, you could run the Kalman filter using the function KF in package KFKSDS instead of KalmanSmooth:

kf <- KFKSDS::KF(y, ss)
a[1,] <- if (m == 1) kf$a.upd[n] else kf$a.upd[n,]
P[,,1] <- if (m == 1) kf$P.upd[n] else kf$P.upd[,,n]

The function KFKSDS::KF is written in R so it may be easier to follow the code of this function instead of the code of function KalmanSmooth. You may read this function along with the Kalman filter equations as defined for example in Durbin and Koopman (2001) cited in ?KalmanRun.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.