Errores estándar de los valores ajustados de una regresión de series temporales

Question

Realmente quiero entender cómo funcionan las matemáticas aquí. Estoy tratando de obtener el error estándar de los valores ajustados para un modelo de regresión de series de tiempo. En la regresión de series no temporales, sé que puedo tomar la transposición de los datos multiplicados por la matriz de varianza - covarianza de los coeficientes del modelo y luego multiplicar por los valores de los datos de nuevo para obtener los errores estándar de los valores ajustados.

Pero no estoy seguro de cómo hacer esto cuando estoy incluyendo un término autorregresivo.

require(forecast)
require(tserieS)

Variable de respuesta

Sablects <- rnorm(10)

Covariables

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

En mis datos reales, los valores están normalizados, por lo que aquí establezco el intercepto igual a cero.

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

La función predecir me dará los errores estándar de mi predicción dentro de la muestra (los valores ajustados de mi modelo).

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

my.se

Ahora, para comparar la salida de mi.se, quiero hacerlo matemáticamente pero no sé qué usar para los valores del término ar2. Uso 1's como marcador de posición para demostrar que mi salida no es igual a los valores de my.se por encima de

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

C

Creo que este valor debería ser igual al primer valor de mi.se, pero no está produciendo el mismo valor que mi.se

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

Además, no soy muy bueno con la multiplicación de matrices, pero aquí está mi trabajo para obtener todos los valores de se.

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

Especificar que el número máximo de i es igual al número de filas de 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, ])  
    }

se.output

Answer 1

Se pueden obtener expresiones de forma cerrada para la varianza de los errores de predicción en un modelo ARMA. Evaluando $Var(y_{T+i} - E(y_{T+i})) = Var(e_{T+i})$ para la expresión de un modelo ARMA definido para la serie $y_t$ , $t=1,2,...,T$ , encontramos la varianza del $i$ -Error de predicción por pasos.

Por ejemplo, se puede comprobar que en un modelo AR(1): \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}

En un 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}

En un 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}

donde $\sigma^2_\epsilon$ es la varianza del término de perturbación, $\phi$ es el coeficiente AR y $\theta$ es el coeficiente MA.

El filtro de Kalman puede utilizarse para realizar estas operaciones para un ARMA general. Para simplificar, tomemos una serie simulada a partir de un modelo AR(2) y ajustémosla a un modelo AR(2) (sin regresores externos y parámetros fijos):

set.seed(123)
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)

Las previsiones y sus errores estándar son los siguientes:

res1 <- predict(fit, n.ahead = 10)
res1
#$pred
#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
#$se
#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

Para un modelo ARMA general, los errores estándar de las previsiones se obtienen como la multiplicación de la varianza de los errores de predicción devuelta por KalmanForecast por la varianza estimada del término de perturbación (y sacando la raíz cuadrada):

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

La mayor parte de las matemáticas que le interesan son realizadas por KalmanForecast . Para mostrar estas operaciones, primero necesitamos definir las matrices de la representación del espacio de estados del modelo ARMA (objeto ss ) y ejecutar KalmanSmooth para obtener el vector de estado y su matriz de covarianza:

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)

A continuación, las predicciones y las varianzas se obtienen de la siguiente manera (basándose en la función predict.stsmSS en el paquete 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]))
}

Se obtienen las mismas previsiones y errores estándar que las anteriores:

pred
# [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)
res3
# [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

Como alternativa, se puede ejecutar el filtro de Kalman utilizando la función KF en el paquete KFKSDS en lugar de KalmanSmooth :

require("KFKSDS")
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]

La función KFKSDS::KF se escribe en R por lo que puede ser más fácil seguir el código de esta función en lugar del código de la función KalmanSmooth . Puede leer esta función junto con las ecuaciones del filtro de Kalman definidas, por ejemplo, en Durbin y Koopman (2001) citado en ?KalmanRun .