Lección 8. Predicción de modelos ARIMA

Cargamos las funciones auxiliares (véase la carpeta src/)

import warnings
warnings.filterwarnings("ignore", category=UserWarning)
%run -i ./src/analisis_armas.py

Consideremos una serie temporal \(\boldsymbol{y}_1^T=(y_1,\ldots,y_T)\). En relación con la predicción de los valores futuros de la serie, hay tres aspectos a considerar :

• Necesitamos un método para hacer predicciones basadas en un modelo ajustado a las observaciones.

• Necesitamos evaluar la incertidumbre asociada a estas predicciones. Sin ello la utilidad de los pronósticos se ve comprometida.

  Es común expresar la incertidumbre con un intervalo de confianza para la predicción; un rango que con alta probabilidad (como 0.95 o 0.99) incluye el valor futuro pronosticado.

• Además, la previsión es un proceso dinámico: al recibir nuevos datos, es necesario actualizar los pronósticos con la nueva información.

Los modelos lineales facilitan el cálculo de estos tres aspectos.

Queremos predecir \(Y_{t+\ell}\) usando una función \(g\) de variables aleatorias con índice \(n\leq t\), es decir, una función del conjunto \(\boldsymbol{Y}_{|:t} = \{Y_n \mid n \leq t\}\).

Nos referiremos al conjunto \(\boldsymbol{Y}_{|:t}\) como ``el pasado'' de \(\boldsymbol{Y}\).

Para evaluar diferentes funciones de predicción, estableceremos como criterio minimizar el error cuadrático medio (MSE): \[ MSE\big(Y_{t+\ell},\,g_t(\ell)\big) = E\Big(Y_{t+\ell}-g_t(\ell)\Big)^2. \]

Preferiremos aquel método \(g_t(\ell)\) que minimiza la esperanza del error cuadrático \(Y_{t+\ell}-g_t(\ell)\); donde \(g_t(\ell)\) es la función de predicción de \({Y}_{t+\ell}\), \(t\) es el índice de la última variable aleatoria considerada y \(\ell\) representa el horizonte de predicción.

La función \(g_t(\ell)\) que minimiza dicho criterio MSE es la esperanza condicional de \(Y_{t+\ell}\): \[ E\big(Y_{t+\ell}\mid \boldsymbol{Y}_{|:t}\big). \]

Las predicciones resultantes se llaman predicciones de error cuadrático medio mínimo.

El problema de la esperanza condicional radica en que, según el caso, su cálculo puede ser muy difícil… o sencillamente inviable.

Como solución alternativa podemos limitarnos a buscar predictores entre:

• las funciones de \(\boldsymbol{Y}_{|:t}\) (i.e., funciones ``del pasado'') que son lineales;
• para seleccionar aquella con menor error cuadrático medio (MSEL).

Al adoptar esta estrategia, aceptamos cometer mayores errores de predicción, \(MSEL \geq MSE\), pues al excluir funciones no lineales, podríamos estar omitiendo la que genera los errores más pequeños.

Pregunta: ¿Cómo tiene que ser la esperanza condicional para que \(MSEL\) y \(MSE\) sean iguales? ¿En qué caso ocurrirá?

Sea \( \boldsymbol{Y} \) un proceso ARIMA\((p,d,q)\times(P,D,Q)_S\) \[ \boldsymbol{\phi}_p(\mathsf{B})\,\boldsymbol{\Phi}_P(\mathsf{B}^S)\,\nabla^d\,\nabla_{_S}^D\, Y_t = c +\boldsymbol{\theta}_q(\mathsf{B})\,\boldsymbol{\Theta}_q(\mathsf{B}^S)\, U_t; \quad t\in\mathbb{N}. \] donde \(c\) puede ser distinto de cero para permitir procesos con esperanza no nula.

Denotemos con \(\boldsymbol{\varphi}\) al polinomio de grado \(p+P+d+D\) resultante del producto convolución de los polinomios de la parte AR: \[ \boldsymbol{\varphi}(z)=\boldsymbol{\phi}_p(z)*\boldsymbol{\Phi}_P(z^S)*\nabla^d*\nabla_{_S}^D. \] Y denotemos con \(\boldsymbol{\vartheta}\) al polinomio de grado \(q+Q\) resultante del producto convolución de los polinomios de la parte MA: \[\boldsymbol{\vartheta}(z)=\boldsymbol{\theta}_q(z)*\boldsymbol{\Theta}_q(z^S).\] Entonces podemos reescribir el modelo ARIMA de una manera mucho más compacta: \[ \boldsymbol{\varphi}(\mathsf{B})\, Y_t = c + \boldsymbol{\vartheta}(\mathsf{B})\, U_t. \]

Debido a la simple estructura de los modelos ARIMA, tenemos que

De \eqref{eqn:ARMA-en-t}, para \( \ell = 0, \pm 1, \pm 2, \dots \) y \( t \) arbitrario

Dividiendo los sumatorios de \(\eqref{eqn:ARMA-en-t+l}\) en dos partes: por un lado la suma desde \(k=1\) hasta \(\ell-1\), y por otro la suma del resto (de \(\ell\) en adelante); tenemos que

Indicar la relación de cada elemento con el subespacio \(\mathcal{H}_{Y_t}\) nos ayudará a calcular las previsiones:

Por linealidad de la proyección ortogonal, el predictor \(\ell\)-pasos adelante de \( Y_t \) es

donde lo que era ortogonal a \(\mathcal{H}_{Y_t}\) se anula, y lo que pertenece a \(\mathcal{H}_{Y_t}\) no cambia.

A partir de esta expresión podemos obtener fórmulas recursivas tanto para los predictores como para sus actualizaciones cuando se disponen de nuevas observaciones.

Consideremos primero modelos lineales, causales y estacionarios.

Sea \(Y_t = \psi(B)U_t\) la representación MA del proceso, donde \(\boldsymbol{\psi}=\boldsymbol{\phi}^{-1}*\boldsymbol{\theta}\). Entonces: \[ Y_{t+\ell} = \sum_{i=0}^{\infty} \psi_i U_{t+\ell-i}. \] La predicción lineal óptima será su proyección ortogonal sobre \(\mathcal{H}_{Y_t}\):

Restando \(\widehat{Y_{t+\ell|t}}\) de \(Y_{t+\ell}\) obtenemos el error de predicción: \[ e_t(\ell) \;=\; Y_{t+\ell} - \widehat{Y_{t+\ell|t}} \;=\; U_{t+\ell} + \psi_1 U_{t+\ell-1} + \cdots + \psi_{\ell-1} U_{t+1}, \] cuya esperanza es cero.

Asumiendo que \(\boldsymbol{Y}\) tiene distribución gaussiana, \(\widehat{Y_{t+\ell|t}}=E(Y_{t+\ell}\mid \boldsymbol{Y}_{|:t})\); y por tanto

\[ E\Big[e_t(\ell)^2\Big]= Var(e_t(\ell))= %E\Big[\big(Y_{t+\ell} - E(Y_{t+\ell}\mid \boldsymbol{Y}_{|:t})\big)^2\Big]= \sigma^2 (1 + \psi_1^2 + \cdots + \psi_{\ell-1}^2) \] donde \(1 + \psi_1^2 + \cdots + \psi_{\ell-1}^2\) es la suma del cuadrado de los \(\ell\) primeros coeficientes de la secuencia \(\boldsymbol{\psi}=\boldsymbol{\phi}^{-1}*\boldsymbol{\theta}\), que es una secuencia de cuadrado sumable.

Y ahora consideremos modelos lineales, causales NO estacionarios.

• Suponga que \(\boldsymbol{Y}\) es SARIMA (\(\boldsymbol{\phi}\) tiene raíces en el círculo unidad) pero que simultáneamente \(\boldsymbol{U}\) tiene un comienzo (que \(U_t=0\) para \(t<0\); por tanto \(Y_t = \psi(B)U_n=0\) para \(t<0\)).
• El cálculo sigue siendo el mismo. Aunque ahora \(\boldsymbol{\psi}=\boldsymbol{\phi}^{-\triangleright}*\boldsymbol{\theta}\) no es de cuadrado sumable, en las sumas solo hay una cantidad finita de términos no nulos.
• La incertidumbre en las predicciones difiere entre modelos estacionarios y no estacionarios. En los estacionarios, como \(\psi_k \to 0\) cuando \(k\) aumenta, la varianza de la predicción tiene a un valor constante: la varianza marginal del proceso.

Por ejemplo, en un modelo AR(1), donde \(\psi_k = \phi^k\), cuando \(\ell\) tiende a infinito \[ \lim_{\ell\to\infty}\,Var(e_T(\ell)) = \sigma^2 (1 + \psi_1^2 + \psi_{2}^2 + \cdots ) = \sigma^2 / (1 - \phi^2). \] Si \(|\phi|\) está próximo uno, la varianza marginal será mucho más grande que la varianza condicionada, \(\sigma^2\). Pero la incertidumbre es finita pues la varianza marginal también lo es.

• En modelos no estacionarios, como la serie \(\sum \psi_i^2\) no converge, la incertidumbre en la predicción a largo plazo crece indefinidamente: no se puede anticipar el comportamiento de un proceso no estacionario a largo plazo.

Contamos con predicciones para \(t+1\) hasta \(t+j\) basadas en \(\mathcal{H}_{Y_t}\). Queremos actualizarlas con \(\mathcal{H}_{Y_{t+1}}\). Según \eqref{eqn:predicion-MA-infinito}, la predicción de \(Y_{t+\ell}\) con información hasta \(t\) es: \[ \widehat{Y_{t+\ell|t}} = \psi_\ell U_t + \psi_{\ell+1} U_{t-1} + \ldots. \] Comparando \(\widehat{Y_{t+1|t}}\) con \(Y_{t+1}\) obtenemos el error de predicción \(U_{t+1} = Y_{t+1} - \widehat{Y_{t+1|t}}\).

La nueva predicción para \(Y_{t+\ell}\), incorporando \(Y_{t+1}\) como información adicional, será: \[ \widehat{Y_{t+\ell-1|t+1}} = \psi_{\ell-1} U_t + \psi_\ell U_{t-1} + \ldots \] Restando las dos previsiones, tenemos que: \[ \widehat{Y_{t+\ell-1|t+1}} - \widehat{Y_{t+\ell|t}} = \psi_{\ell-1} U_{t+1}. \] Al observar \(Y_{t+1}\) y obtener el error de previsión \(U_{t+1}\), podremos revisar las previsiones: \[ \widehat{Y_{t+\ell-1|t+1}} = \widehat{Y_{t+\ell|t}} + \psi_{\ell-1} U_{t+1} %, \tag{8.23} \] (las predicciones se revisan sumando un múltiplo del último error de predicción).

Si no hay error de previsión un periodo adelante \((U_{t+1} = 0)\), no hay revisión del resto.

Si se tomaron logaritmos para estabilizar la varianza, es necesario recuperar las previsiones de la serie original a partir de las de la serie transformada. Si \[ Y_{t+\ell} = \ln(X_{t+\ell}) \] entonces la previsión puntual de la serie original es \[ \widehat{X_{t+\ell|t}} = \exp\left( \widehat{Y_{t+\ell|t}} + \frac{1}{2} Var\big(e_t(\ell)\big) \right) \] Y el Intervalo de confianza \[ I_{1-\alpha}(X_{t+\ell}) = (a, b) \] donde

que es no simétrico alrededor de la previsión.

Sea un proceso estocástico lineal causal y estacionario de segundo orden definido por \[ Y_t = \sum_{k=0}^{\infty} \psi_k\, U_{t-k}, \] donde:

• \(\sum_{k=0}^{\infty} \psi_k^2 < \infty\),
• \(\boldsymbol{U}\) es un proceso de ruido blanco con \(E[U_t]=0\) y \(E[U_t^2]=\sigma^2 < \infty\),
• el proceso es invertible, es decir, existe una secuencia \(\{\pi_j\}_{j\ge 0}\) con \(\sum_{j=0}^{\infty} \pi_j^2 < \infty\) tal que \[ U_t = \sum_{j=0}^{\infty} \pi_j\, Y_{t-j}. \]

Bajo las hipótesis anteriores se cumple que los subespacios cerrados de \(L^2\) generados por el pasado de \(Y_t\) y por el pasado de las innovaciones \(U_t\) son iguales: \[ \mathcal{H}_{Y_t} \;=\; \mathcal{H}_{U_t} \] donde \(\mathcal{H}_{Z_n} = \overline{\operatorname{sp}}\{1,Z_n,Z_{n-1},\ldots\}\) denota la clausura en la norma de \(L^2\) del subespacio engendrado por las variables aleatorias \(\{1,Z_n,Z_{n-1},Z_{n-2},\ldots\}\). He incluido la constante \(1\) para poder considerar procesos con esperanza no nula. La ``H'' de la notación \(\mathcal{H}_{Z_n}\) tiene como objetivo recordar que se trata de un espacio de Hilbert y que dicho espacio está engendrado por las variables aleatorias que constituyen la ``Historia'' del proceso desde \(Z_n\) hacia atrás.

Demostración:

1. Primera inclusión \(\subseteq\):

   Cada \(Y_{t}\) se puede expresar como combinación lineal de las innovaciones pasadas: \[ Y_{t} = \sum_{k=0}^{\infty} \psi_k\, U_{t-k}. \] Por tanto, todo \(Y_{t}\) pertenece al subespacio cerrado generado por \(\{1,U_t,U_{t-1}, U_{t-2}, \dots\}\), y al incluir la clausura en \(L^2\) obtenemos \[ \mathcal{H}_{Y_t} % \overline{\operatorname{sp}}\{\boldsymbol{Y}_{|:(t-1)}\} \subseteq \mathcal{H}_{U_t} % \overline{\operatorname{sp}}\{\boldsymbol{U}_{|:(t-1)}\} \]

2. Segunda inclusión \(\supseteq\):

   Por invertibilidad, cada \(U_{t}\) se puede escribir como combinación de las \(Y\) pasadas: \[ U_{t} = \sum_{m=0}^{\infty} \pi_m\, Y_{t-m}. \] Así, cada \(U_{t}\) pertenece al subespacio cerrado generado por \(\{1,Y_t,Y_{t-1},Y_{t-2},\dots\}\), y al incluir la clausura en \(L^2\) obtenemos la inclusión opuesta. \[ \mathcal{H}_{Y_t} % \overline{\operatorname{sp}}\{\boldsymbol{Y}_{|:(t-1)}\} \supseteq \mathcal{H}_{U_t} % \overline{\operatorname{sp}}\{\boldsymbol{U}_{|:(t-1)}\} \]