Pronósticos con modelos ARMA
Raúl Godínez
12/12/2020
Una de las funciones principales de los modelos de series de tiempo es el pronóstico. En macroeconomía, una variable de interés para su pronóstico es el PIB, en especial, el PIB trimestral que permite dar seguimiento a la coyuntura económica.
Uno de los métodos ampliamente utilizados es el de los modelos ARMA, el cual, será mostrado a en el presente artículo. Pero antes, se muestra un breve introducción a la teoría y la matemática detrás de estos modelos.
Procesos estocásticos estacionarios:
Se dice que un proceso estocástico es estacionario si su media y su varianza son constantes en el tiempo y si el valor de la covarianza entre dos periodos depende sólo de la distancia o rezago entre estos dos periodos (Gujarati & Porter, 2012). \[\begin{align} \label{CORR1} \mathrm{Corr}[y_t,y_{t-s}]=\frac{\mathrm{cov}[y_t,y_{t-s}]}{\sqrt{\mathrm{var}[y_t]\mathrm{var}[y_{t-s}]}} = \frac{\gamma_s}{\gamma_0} =\rho_s \qquad \mathrm{para}\quad s=1,2,3,... \end{align}\]
Usualmente, cuando se cumplen estas condiciones se le conoce como proceso estocástico débilmente estacionario, estacionario covariante,estacionario de segundo orden o proceso estocástico en amplio sentido.
Propiedades de estacionariedad débil, sea \(Y_t\) una serie de tiempo: \[\begin{align} &\textit{Media:}\qquad \quad E(Y_t)=\mu \\ &\textit{Varianza:}\qquad E(Y_t) =E(Y_t-\mu)^2=\gamma_0\\ &\textit{Covarianza:}\quad \gamma_{k} = E[(Y_t-\mu)(Y_{t+k}-\mu)] \end{align}\]
Porque si una serie de tiempo es no estacionaria, sólo podemos estudiar su comportamiento durante el periodo en consideración. Por tanto, cada conjunto de datos perteneciente a la serie de tiempo corresponder a un episodio particular. En consecuencia, no es posible generalizar para otros periodos. Así, para propósitos de pronóstico, tales series de tiempo (no estacionarias) tienen poco valor práctico (Gujarati & Porter, 2012).
Proceso puramente aleatorio o de ruido blanco La serie \({\varepsilon_t}\) es un proceso de ruido blanco si cada valor en la serie tiene una media de cero, una varianza constante y no está correlacionado con todas las demás realizaciones. Formalmente, si la notación \(E(x)\) denota el valor medio teórico de \(x\), la serie \({\varepsilon_t}\) es un proceso de ruido blanco si para cada período de tiempo \(t\) (Enders, 1948). \[\begin{align} &E(\varepsilon_t)=E(\varepsilon_{t-1})=\dots =0\\ &E(\varepsilon^2_{t})=E(\varepsilon^2_{t-1})=\dots =\sigma^2\qquad \quad [\mathtt{var}(\varepsilon_t)=\mathtt{var} (\varepsilon_{t-1})=\dots=\sigma^2]\\ &E(\varepsilon_t \quad\varepsilon_{t-s})=E(\varepsilon_{t-j}\quad \varepsilon_{t-j-s})=0 \end{align}\]
Procesos estocásticos no estacionarios
El ejemplo clásico de una serie de tiempo no estacionaria es el modelo de caminata aleatoria (MCA).
la cual se clasifica en dos tipos: Caminata aleatoria sin deriva y la caminata aleatoria con deriva (deriva=intercepto)
Caminata aleatoria sin deriva
Suponga que \(u_t\) es un término de error de ruido blanco, con media 0 y varianza constante \(\sigma^2\). Entonces decimos que la serie \(Y_t\) es una caminata aleatoria si \[\begin{align} \label{CASD} &Y_t=Y_{t-1}+u_t \end{align}\] Por tanto, es un modelo AR(1), del modelo anterior podemos escribir
\[\begin{align*} &Y_1 =Y_0+u_1\\ &Y_2 =Y_1+u_1 = Y_0+u_1+u_2\\ &Y_3=Y_2+u_3 =Y_0+u_1+u_2+u_3\\ \end{align*}\] De forma general, su representación cuando inicia en el tiempo 0 es: \[\begin{align} \label{sum} Y_t=Y_0+\sum u_t \end{align}\] Por tanto, \[\begin{align} E(Y_t)=E(Y_0+\sum u_t)=Y_0 \end{align}\] de igual forma se demuestra que \[\begin{align} \mathrm{var}(Y_t)=t\sigma^2 \end{align}\]
Como revelan las expresiones anteriores, la media de \(Y\) es igual a su valor inicial (constante), pero conforme se incrementa t, su varianza aumenya de manera indefinida, lo que viola una condición de estacionariedad. En resumen, el MCA sin deriva es un proceso estocástico no estacionario.
Caminata aleatoria con deriva
Modificando de la siguiente forma \[\begin{align} \label{CACD} Y_t=\delta+Y_{t-1}+u_t \end{align}\]
donde \(\delta\) se conoce como el parámetro de deriva. El término deriva proviene del hecho de que si escribimos la ecuación anterior en su primera diferencia \[\begin{align} \label{PD} Y_t-Y_{t-1}=\Delta Y_t=\delta+u_t \end{align}\]
se demuestra que \(Y_t\) se deriva o desvía hacia arriba o hacia abajo, según \(\delta\) sea positiva o negativa.Observe que el modelo de caminata aleatoria con deriva es un modelo AR(1).
Usando la misma demostración del modelo sin deriva, se comprueba que el modelo con deriva también es no estacionario en el sentido débil y que su varianza aumenta a media que aumenta \(t\).
Por último, es importante saber que el modelo de caminata aleatoria es un ejemplo de lo que se conoce como Proceso de raíz unitaria
Proceso estocástico de raíz unitaria
Escribiendo el MCA como:
\[\begin{align} \label{RU} Y_t=\rho Y_{t-1}+u_t \qquad \quad -1\leq \rho \leq 1 \end{align}\]
Si \(\rho=1\), se convierte en MCA (sin deriva). Si \(\rho\) es en efecto 1, se tiene lo que se conoce como problema de raíz unitaria; es decir, enfrentamos una situación de no estacionariedad.
Por tanto, los términos: caminata aleatoria, raíz unitaria, no estacionariedad y tendencia estocástica se consideran sinónimos.
Sin embargo, si el valor absoluto de \(\rho\) es menor a 1, se puede demostrar que la serie de tiempo es estacionaria.
Procesos estocásticos estacionarios en tendencia (ET) y estacionarios en diferencias (ED)
En términos generales, si la tendencia de una serie de tiempo es del todo predecible y no variable, se le llama tendencia determinista; si no es predecible, se le llama tendencia estocástica. Para formalizar la definición, considere el siguiente modelo de la serie de tiempo \(Y_t\); \[\begin{align} \label{E} Y_t=\beta_1+\beta_2 t +\beta_3 Y_{t-1} u_t \end{align}\]
donde \(u_t\) es un término de error de ruido blanco. Ahora tenemos las siguientes probabilidades:\ Si en el modelo anterior, \(\beta_1=0,\quad \beta_2=0, \quad \beta_3=1\), se obtiene
\[\begin{align} \label{CAP} Y_t=Y_{t-1}+u_t \end{align}\]
que como indentificamos, es un MCA sin deriva y por tanto no estacionario. Pero si expresamos este modelo como: \[\begin{align} \label{PDCAP} \Delta Y_t =(Y_t-Y_{t-1})=u_t \end{align}\] se convierte en estacionaria. Por tanto, un MCA sin deriva es un proceso estacionario en diferencias (PED).
caminata aleatoria con deriva si en modelo \(Y_t=\beta_1+\beta_2 t +\beta_3 Y_{t-1}\) , \(\beta_1\neq 0,\quad \beta_2=0, \quad \beta_3=1\) se obtiene
\[\begin{align} \label{LA} Y_t=\beta_1+Y_{t-1}+u_t \end{align}\]
Que es una caminata aleatoria con deriva y en consecuencia es no estacionaria. Si se expresa como
\[\begin{align} \label{CACD2} Y_t-Y_{t-1}=\Delta Y_t = \beta_1 + u_t \end{align}\] esto significa que \(Y_t\) mostrará una tendencia positiva (\(\beta_1>0\)) o negativa (\(\beta_1<0\)). Tal tendencia se llama tendencia estocástica. La ecuación anterior es un PED porque la no estacionaridad en \(Y_t\) se elimina al tomar primeras diferencias.
: Si en \(Y_t=\beta_1+\beta_2 t +\beta_3 Y_{t-1}\), \(\beta_1\neq 0,\quad \beta_2\neq 0, \quad \beta_3=0\) se obtiene \[\begin{align} \label{5A} Y_t=\beta_1+\beta_2 t+u_t \end{align}\] la cual se llama proceso estacionario en tendencia (PET).
Procesos estocásticos integrados
El modelo de caminata aleatoria no es más que un caso específico de una clase más general de procesos estocásticos conocidos como procesos integrados. Recuerde que el MCA sin deriva es no estacionario, pero su serie de primeras diferencias, es estacionaria. Por tanto, el MCA sin deriva se llama proceso integrado de orden 1 y se denota como I(1). De manera similar, si una serie de tiempo tiene que diferenciarse dos veces (es decir, se toman primeras diferencias de la serie de primeras diferencias) para hacerla estacionaria, esa serie de tiempo se denomina integrada de orden 2. En general, si una serie de tiempo (no estacionaria) debe diferenciarse d veces para hacerla estacionaria, decimos que la serie es . Una serie de tiempo \(Y_t\) integrada de orden d se denota como \(Y_t \sim I(d)\). Si una serie de tiempo es estacionaria desde el principio, es de orden 0.
Propiedades de las series integradas
Sean las series de tiempo \(X_t\), \(Y_t\) y \(Z_t\):
Si \(X_t\sim I(0)\), \(Z_t=(X_t+Y_t)=I(1)\); es decir, una combinación lineal o suma de series de tiempo estacionaria y no estacionaria es no estacionaria.
Si \(X_t\sim I(d)\), \(Z_t=(a+bX_t)=I(d)\), donde \(a\) y \(b\) son constantes. Es decir, una combinación lineal de una serie \(I(d)\) es también es \(I(d)\)
Si \(X_t\sim I(d_1)\) y \(I(d_2)\), \(Z_t=(aX_t+bY_t)\sim I(d_2)\), donde \(d_1<d_2\).
Si \(X_t\sim I(d)\) y \(Y_t\sim I(d)\), \(Z_t=(aX_t+bY_t)\sim I(d^*)\), \(d^*\) es por lo general igual a \(d\), pero en algunos casos \(d^*<d\) \end{itemize}
Modelos ARMA
Procesos autoregresivos (AR)
Un proceso autoregresivo de orden \(p\) o proceso AR(p) es un procesos \(\{Y_t\}^{\infty}_{-\infty}\) para el que se satisface: \[\begin{align} \label{AR} Y_t=a_0+a_1Y_{t-1}+a_2Y_{t-2}+\dots + a_iY_{t-p}+\varepsilon_t \end{align}\] Siendo \(a_0\), \(a_1\), \(a_2\) \(\dots, a_i\) parámetros y \(\{\varepsilon_t\}^{\infty}_{-\infty}\) ruido blanco. De manera resumida se puede expresar como: \[\begin{align} \label{SAR} Y_t=a_0+\sum^p_{i=1}a_iY_{t-i}+\varepsilon_t \end{align}\]
Restricciones de estacionariedad para un AR(1)
Considere un proceso AR(1) \[\begin{align} \label{AR1} Y_t=a_0+a_1Y_{t-1}+\varepsilon_t \end{align}\] Suponga que este proceso empieza en el periodo cero, así que \(Y_0\) es una condición inicial determinista. La solución de esta ecuación es: \[\begin{align} \label{SOAR1} Y_t=a_0\sum_{i=0}^{t-1}a_1^{i}+a_1^{t}Y_0+\sum_{i=0}^{t-1}a_1^{i}\varepsilon_{t-i} \end{align}\] Tomando el valor esperado de la ecuación anterior, se obtiene \[\begin{align} \label{VESOAR} E(Y_t)=a_0\sum_{i=0}^{t-1}a_1^{i}+a^t_{1}Y_0 \end{align}\] y para \(s\) periodos: \[\begin{align} \label{SOAR1S} E(Y_{t+s})=a_0\sum_{i=0}^{t+s-1}a_1^{i}+a^{t+s}_{1}Y_0 \end{align}\] Comparando estas dos últimas ecuaciones, es claro que ambas medias son dependientes del tiempo. Dado que \(E(Y_t)\) no es igual a \(E(Y_{t+s})\) la serie no puede ser estacionaria. Sin embargo, si t es grande, podemos considerar limitar el valor de \(Y_t\). Si \(|a_1|<1\), la expresión \((a_1)^tY_0\) converge a cero cuando \(t\to \infty\) y la suma \(a_0[1+a_1+(a_1)^2+(a_1)^3+\dots ]\) converge a \(a_0/(1-a_1)\) así, cuando \(t\to \infty\) y si \(|a_1|<1\): \[\begin{align} \label{LIMAR1} \lim_{t \to \infty} Y_t= \frac{a_0}{1-a_1} + \sum^{\infty}_{i=0}a_1^{i}\varepsilon_{t-i} \end{align}\]
ahora en el valor esperado, como t es suficientemente grande, \(E(Y_t)=a_0/(1-a_1)\). Así, el valor esperado de \(Y_t\) es finito y es independiente del tiempo, por lo tanto, $E(Y_t)=E(Y_{t-s})=a_0/(1-a_1) $ para todo s.
Mientras la varianza se expresa de la siguiente forma \[\begin{align} \label{VAR1} E(Y_t-\mu)^2 &= E[(\varepsilon_t + a_1\varepsilon_{t-1}+(a_1)^2\varepsilon_{t-2}+\dots )^2]\\ &=\sigma^2[1+(a_1)^2+(a_1)^2+\dots] =\frac{\sigma^2}{1-(a_1)^2} \end{align}\] la cual también es finita y es independiente del tiempo, finalmente es fácil demostrar que el límite de los valores de todas las autocovarianzas son finitas e independientes del tiempo. \[\begin{align} \label{ACARS} E[(Y_t-\mu)(Y_{t-s}-\mu)]&=E\{[\varepsilon_t+a_1\varepsilon_{t-1}+(a_1)^2\varepsilon_{t-2}+\dots ][\varepsilon_{t-s}+a_1\varepsilon_{t-s-1}+(a_1)^2\varepsilon_{t-s-2}+\dots]\}\\ &=\sigma^2(a_1)^s[1+(a_1)^2+(a_1)^4+\dots]\\ &=\frac{\sigma^2(a_1)^s}{[1-(a_1)^2]} \end{align}\]
Si no tuviéramos la condición inicial, poco cambiaría, la suma de la solución homogénea y particular para \(\{Y_t\}\) es \[\begin{align} \label{NLAR1} Y_t=\frac{a_0}{(1-a_1)}+\sum^{\infty}_{i=0}a^i_{1}\varepsilon_{t-i}+A(a_1)^t \end{align}\] Si tomamos el valor esperado es evidente que la serie \(Y_t\) no puede ser estacionaria a no ser que \(A(a_1)^t\) sea igual a cero, o bien la serie debe haber comenzado hace mucho tiempo (así que \(a^t_{1}=0\)) o la constante arbitraria \(A\) es cero.
La constante \(A\) se interpreta como la desviación del equilibrio de largo plazo. Recordando las condiciones de estabilidad, estás pueden ser suficiente.
Procesos de media móvil
Un proceso de medias móviles de orden \(q\) (o proceso MA(\(q\))) es un proceso \(\{Y_t\}^{\infty}_{-\infty}\) para el que se satisface: \[\begin{align} \label{MA} Y_t=\mu +\varepsilon_t +\beta_i\varepsilon_{t-1}+\beta_2\varepsilon_{t-2}+\dots+\beta_q\varepsilon_{t-q} \end{align}\] de manera resumida, un proceos MA(\(q\)) puede ser expresado de la siguiente forma \[\begin{align} \label{MAG} X_t=\sum^q_{i=0} \beta_i\varepsilon_{t-i} \end{align}\] En particular un proceso MA(1) puede escribirse como: \[\begin{align*} Y_t&=\mu +\varepsilon_t+\beta_1\varepsilon_{t-1}\\ E[Y_t]&=E(\mu +\varepsilon_t+\theta\varepsilon_{t-1})=\mu\\ \gamma_0 &= \mathrm{Var}[Y_t] = E[(Y_t-\mu)^2] = E[(\varepsilon_t+\beta_1\varepsilon_{t-1})^2] =E(\varepsilon^2_{1})+\beta_1^{2}E(\varepsilon^2_{t-1})+2\beta_1E(\varepsilon_t\varepsilon_{t-1})= (1+\beta_1^{2})\sigma^2\\ \gamma_1 &= \beta_1\sigma^2\\ \gamma_s&=E[(\varepsilon_t+\beta_1\varepsilon_{t-1})(\varepsilon_{t-s}+\beta_1\varepsilon_{t-s-1})] = 0, \qquad \forall s\geq 2\\ \end{align*}\] Puesto que la media, la varianza y las autocovarianzas no dependen de \(t\) el proceso es estacionario.
Es posible combinar en una ecuación el proceso de medias móviles (MA) y diferenciación lineal (AR) y obtener un modelo autoregresivo de medias móviles (ARMA);
\[\begin{align} \label{ARMA1} Y_t=a_0+a_1Y_{t-1}+\dots a_{p} Y_{t-p}+\varepsilon_{t} - \beta_1 \varepsilon_{t-1} - \dots -\beta_q \varepsilon_{t-q} \end{align}\]
siendo \(a_0,a_1,...a_p, \beta_1,...\beta_q\) parámetros y \(\{\varepsilon\}^{\infty}_{-\infty}\) es ruido blanco, y si expresamos un modelo ARMA en su forma general \[\begin{align} \label{ARMA2} Y_t=a_0+\sum^{p}_{i=1}a_iY_{t-i}+\sum^q_{i=1}\beta_i\varepsilon_{t-i} \end{align}\] \[\begin{align} \label{ARMA3} Y_t=a_0+\sum^{p}_{i=1}a_iY_{t-i}+X_t \end{align}\] Si las raíces características de un modelo ARMA están todas en el circulo unitario, \(\{Y_t\}\) es llamado un modelo ARMA para \(Y_t\).
En un modelo ARMA, es perfectamente permisible permitir que \(p\) y/o \(q\) tiendan a infinito. Sin embargo, si una o más raíces características son mayor a una unidad, la serie \(Y_t\) se adapta a un proceso integrado y es un modelo medias móvil autoregresivo integrado (ARIMA).
Tratando el modelo como una ecuación en diferencias, sugiere que podemos resolver para \(Y_t\) en términos de la serie \(\{\varepsilon_t\}\). La solución de un modelo ARMA expresando \(Y_t\) en términos de la serie \(\{\varepsilon_t\}\) es la representación de medias móvil de \(Y_t\).
Para el modelo AR(1); \(Y_t=a_0+a_1Y_{t-1}+\varepsilon_{t}\), la representación de medias móvil es: \[\begin{align*} Y_t=\frac{a_0}{(1-a_1)} +\sum_{i=0}^{\infty}a^i_{1}\varepsilon_{t-i} \end{align*}\]
Función de autocorrelación (FAC)
Proceso MA(1)
\[\begin{align} Y_t= \varepsilon_t +\beta\varepsilon_{t-1} \end{align}\]
Multiplicando ambos lados de la ecuación por \(Y_{t-s}\) calculando su valor esperado:
\[\begin{align*} \gamma_0&=E(Y_tY_t)=E[(\varepsilon_t+\beta\varepsilon_{t-1})(\varepsilon_t+\beta\varepsilon_{t-1})] =(1+\beta^2)\sigma^2\\ \gamma_1 &=E(Y_tY_{t-1}) =E[(\varepsilon_t+\beta\varepsilon_{t-1})(\varepsilon_{t-1}+\beta\varepsilon_{t-2})] =\beta_1\sigma^2\\ \gamma_s &=E(Y_tY_{t-s}) =E[(\varepsilon_t+\beta\varepsilon_{t-1})(\varepsilon_{t-s}+\beta\varepsilon_{t-s-1})] =0 \quad \forall s>q \end{align*}\]
Ahora aplicando las propiedades de Yule Walker, que son las siguientes: \[\begin{align} E(Y_{t},Y_{t-s})&=E(Y_{t-s},Y_{t})=E(Y_{t-k},Y_{t-k-s})=\gamma_s\\ E(\varepsilon_tY_t)&=\sigma^2\\ E(\varepsilon_tY_{t-s})&=0 \end{align}\]
Se obtiene así la función de autocorrelación para procesos MA(1) \[\begin{align} \rho_0 & = \gamma_0/\gamma_0 = 1\\ \rho_1 & = \gamma_1/\gamma_0= \frac{\beta} {(1+\beta^2)}\\ \rho_s & = 0 \qquad s>q \end{align}\]
Proceso AR(1)
\[\begin{align} Y_t=a_0 +a_1Y_{t-1}+a_2Y_{t-2}+...a_pY_{t-p}+\varepsilon_t \end{align}\]
Función de autocorrelación parcial
El coeficiente de autocorrelación parcial en el retardo s denotado por \(\phi_{s,s}\) es la correlación entre \(Y_t\) e \(Y_{t-s}\) después de extraer la influenza de los retados intermedios.
El cálculo de autocorrelaciones parciales puede basarse en el modelo de regresión múltiple en desviaciones respecto a las medias de esta. \[\begin{align} Y_t =\phi_{1,s}\overline{Y}_{t-1}+\dots+\phi_{s,s}\overline{Y}_{t-s}+u_t \end{align}\] Multiplicando por \(Y_{t-s}\) y tomando su valor esperado \[\begin{align} E(Y_t,Y_{t-s})=\phi_{1,s}E(\overline{Y}_{t-1},Y_{t-s})+\dots+\phi_{s,s}E(Y^2_{t-s})+E(u_t,Y_{t-s}) \end{align}\] Aplicando las reglas de Yule Walker: \[\begin{align} \gamma_s=\phi_{1,s}\gamma_{s-1}+\dots+\phi_{s,s}\gamma_{t-s} \end{align}\]
dividiendo por \(\gamma_0\) podemos especificar el sistema de ecuaciones de Yule Walker
\[\begin{align} \begin{pmatrix} \rho_1\\ \rho_2\\ \vdots\\ \rho_s \end{pmatrix} = \begin{pmatrix} \rho & \rho_1 & \dots & \rho_{s-1}\\ \rho_1 & \rho_0 &\dots &\rho_{s-2}\\ \vdots & \dots & \dots & \vdots\\ \rho_{s-1} & \rho_{s-2} &\dots & \rho_0 \end{pmatrix} \begin{pmatrix} \phi_{1,s}\\ \phi_{2,s}\\ \vdots\\ \phi_{s,s} \end{pmatrix} \end{align}\] Aplicando la regla de cramer
\[\begin{align} \phi_{s,s} = \frac{\begin{vmatrix} 1 & \rho_1 & \dots & \rho_1\\ \rho_1 & 1 & \dots & \rho_2\\ \vdots & \vdots & \dots & \vdots\\ \rho_{s-1} & \rho_{s-2 }& \dots & \rho_0\\ \end{vmatrix} }{\begin{vmatrix} 1 & \rho_1 & \dots & \rho_{s-1}\\ \rho_1 & 1 & \dots & \rho_{s-2}\\ \vdots & \vdots & \dots & \vdots\\ \rho_{s-1} & \rho_{s-2 }& \dots & 1\\ \end{vmatrix}} \end{align}\]
Cuyo resultado son las ecuaciones de Durbin para todo \(s >2\):
\[\begin{align} \phi_{s,s}= \frac{\rho_s -\sum_{j=1}^{s-1} \phi_{s-1,j} \rho_{s-j}}{1-\sum_{j=1}^{s-1} \phi_{s-1,j} \rho_j} \end{align}\]
\[\begin{align} \phi_{11} &= \rho_1\\ \phi_{22} &= \frac{\rho_2-\rho_{1}^{2}}{1-\rho^2_{1}} \end{align}\]
Metodología Box-Jenkins para la selección de modelos
Las estimaciones de los modelos AR (1), ARMA (1, 1) y AR (2) en la sección anterior ilustran la estrategia de Box-Jenkins (1976) para la selección apropiada del modelo. Box y Jenkins popularizaron un método de tres etapas destinado a seleccionar un modelo apropiado con el fin de estimar y pronosticar una serie de tiempo univariante. En la etapa de identificación, el investigador examina visualmente el gráfico de tiempo de la serie, la función de autocorrelación y la función de correlación parcial. El trazado de la ruta de tiempo de la secuencia \(\{Y_t\}\) proporciona información útil sobre valores atípicos, valores perdidos y saltos estructurales en los datos.
Las variables no estacionarias pueden tener una tendencia pronunciada o parecer serpenteantes sin una media o varianza constante a largo plazo. Los valores perdidos y los valores atípicos se pueden corregir en este punto. En un momento, la práctica estándar era diferenciar primero cualquier serie se considera no estacionario.
Actualmente, existe una gran cantidad de literatura sobre procedimientos formales para verificar la no estacionariedad. Una comparación de la muestra FAC y FACP a los de varios procesos teóricos de ARMA pueden sugerir varios modelos plausibles. En la etapa de estimación, cada uno de los modelos tentativos se ajusta y se examinan los diversos coeficientes \(a_i\) y \(\beta_i\). En esta segunda etapa, el objetivo es seleccionar un modelo estacionario y parsimonioso que tenga un buen ajuste. La tercera etapa implica la verificación de diagnóstico para garantizar que los residuos del modelo estimado imiten un proceso de ruido blanco.
Criterios de selección de un modelo ARMA, ARIMA
- Criterio de Akaike: AIC
- Criterio de Schwarz: BIC
- Criterio de Hannan-Quinn:HQC
Se seleccionará el modelo que tenga el menor valor en los criterios de selección.
En la práctica, esta teoría juega como papel principal el entendimiento e interpretación del modelo, ya que, la mayoría de cálculos y pasos necesarios para generar el modelo, están automatizados.
Pronóstico ARMA usando R
Para ejemplificar el método, realizaremos un prónostico del PIB trimestral de Nicaragua. La obtención, limpieza y estructuración de los datos, es el primer paso.
Los datos utilizados corresponden a los publicados por el Banco Central de Nicaragua, los cuales, los puede descargar aquí. Al descargarlos, nos encontramos con el primer inconveniente, la disposición de los datos. Si desea conocer más sobre la importancia de la estructura de los datos, les recomiendo el siguiente post del señor Jilber Urbina.
Figure 1: Datos publicados por el Banco Central
Como se muestra en la figura 1, se resalta en rojo las filas y columnas que no contienen información, siendo estas un problemas al importar la tabla. Luego están los cuadros verdes, estos ocupan el lugar del nombre de las variables y realmente son datos. Ante estos problemas, debemos dar un tratamiento a la disposición de los datos. (Nota: Existen mútiples formas de hacerlo, usando paquetes como dplyr o tidyr; sin embargo, en este ejemplo se utilizó un método paso a paso.)
#Librerias a utilizar
library(readxl)
library(tidyr)
library(dplyr)
library(readxl)
PIBT <- read_excel("Cuadros_de_salida_CNT.xls",
sheet = 3, range="A35:BV38")| Componente | 2006 | …3 | …4 | …6 | 2007 |
|---|---|---|---|---|---|
| NA | I | II | III | NA | I |
| NA | NA | NA | NA | NA | NA |
| Producto interno bruto | 29028.05904817687 | 28612.949339097715 | 29738.797709521561 | NA | 30079.698699867582 |
Seleccionamos la hoja y el rango de los datos a utilizar, en este caso, solo utilizamos el producto interno bruto en millones de córdobas de 2006. A continuación, se detallan los pasoa a seguir para ordenar los datos.
PIBT = PIBT[-2,]#Eliminamos la 2da fila, ya que esta está vacía
PIBT= t(PIBT) #Transponemos los datos
colnames(PIBT) = c("Trimestre","PIB") #Agregamos los nombres a las variables
PIBT = as.data.frame(PIBT) #Convertimos de los datos de nuevo en data frames
PIBT = PIBT[!is.na(PIBT$Trimestre),] #Eliminamos las filas con NAs
row.names(PIBT) = c(1:59) #renombramos las filas con el número de la observación
Año = c(rep(2006:2020, each=4)) #Creamos una variable de año
Año = as.data.frame(Año) #la convertimos en un data frame
PIBT = cbind(Año[1:59,],PIBT) #Unimos las variables
colnames(PIBT) = c("Año","Trimestre","PIBT") #renombramos el conjunto de datos| Año | Trimestre | PIBT |
|---|---|---|
| 2006 | I | 29028.05904817687 |
| 2006 | II | 28612.949339097715 |
| 2006 | III | 29738.797709521561 |
| 2006 | IV | 31457.904109438787 |
| 2007 | I | 30079.698699867582 |
| 2007 | II | 29987.867520322616 |
| 2007 | III | 31305.334838206541 |
| 2007 | IV | 33497.429634839078 |
| 2008 | I | 31789.945730553136 |
| 2008 | II | 32174.945761999381 |
| 2008 | III | 32599.411655709617 |
| 2008 | IV | 32596.218618897215 |
Ahora sí, el conjunto de datos se puede trabajar y entender de una mejor manera, tal como se muestra en la tabla anterior. Graficamos la variable para comprobar que no existan errores.
library(tseries)
library(forecast)
library(TSstudio)
library(plotly)
library(ggplot2)
ts_plot(ts(PIBT$PIBT, start = c(2006,1,1), frequency = 4), title = "PIB Trimestral",
Xgrid = TRUE, Ygrid = TRUE)Una vez tenemos la variable, el siguiente paso recomendado es eliminar el componente estacional (es decir, aquellos eventos que ocurren cada año de forma estacional, son los patrones observados en la serie.) Para ello, seguimos los siguientes pasos:
library(seasonal)## Warning: package 'seasonal' was built under R version 4.1.2#Aplicamos logaritmo a la serie
PIBT$PIBT= as.numeric(PIBT$PIBT)
log_pib = log(PIBT$PIBT)
log_pib = ts(log_pib, start = c(2006,1,1), frequency = 4)
#Desestacionalizar la serie
log_pib_sa = seas(log_pib, x11="")
log_pib_sa =log_pib_sa$data[,3]
log_pib_sa = as.data.frame(log_pib_sa)
#La convertimos nuevamente en una serie de tiempo
log_pib_sa = ts(log_pib_sa$x, start = c(2006,1,1), frequency = 4)
ts_plot(log_pib_sa, Xgrid = TRUE, Ygrid = TRUE, color="red",
title = "Log(PIBT) desestacionalizado")Nuestra serie ya está casi lista para modelarse, pero antes, debemos asegurarnos que la serie sea estacionaria. Para este ejemplo, modelaremos los residuos de una regresión que tiene por objetivo extraer el ciclo económico o proceso estocástico de la variable.
library(mFilter)## Warning: package 'mFilter' was built under R version 4.1.2#Creamos unas variables dummys para contrarestar los efectos de las crisis de 2008 y 2018
#El shock de 2008 lo modelamos como un shock transitorio
d_2008 = c(rep(0, times=9),rep(1, times = 7),rep(0, times = 43))
d_2008 = ts(d_2008, start = c(2006,1,1), frequency = 4)
#El de 2018 como un shock permanente para el resto de la serie
d_2018 = c(rep(0,times=49), rep(1, times=10))
d_2018 = ts(d_2018, start = c(2006,1,1), frequency = 4)
#Extraemos la tendencia para incluirla en el modelo
pibt_hp = hpfilter(log_pib_sa,freq=1600,type="lambda")
#Modelo
detren_pib = lm(log_pib_sa~d_2008+d_2018+pibt_hp$trend+pibt_hp$trend*d_2008)
(summary(detren_pib))##
## Call:
## lm(formula = log_pib_sa ~ d_2008 + d_2018 + pibt_hp$trend + pibt_hp$trend *
## d_2008)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.066572 -0.016131 -0.003036 0.013683 0.047842
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.679919 0.263255 -2.583 0.01254 *
## d_2008 21.061034 6.139233 3.431 0.00116 **
## d_2018 -0.038284 0.008886 -4.308 7e-05 ***
## pibt_hp$trend 1.065243 0.025013 42.588 < 2e-16 ***
## d_2008:pibt_hp$trend -2.031204 0.591881 -3.432 0.00116 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02251 on 54 degrees of freedom
## Multiple R-squared: 0.9793, Adjusted R-squared: 0.9778
## F-statistic: 639.7 on 4 and 54 DF, p-value: < 2.2e-16#Extraemos los residuos como nuestra variable a modelar
cyl_pibt = ts(detren_pib$residuals, start = c(2006,1,1), frequency = 4)
ts_plot(cyl_pibt,Xgrid = TRUE, Ygrid = TRUE, color= "#36B5D9",
title = "Ciclo del PIBT")Prueba de estacionariedad
El siguiente paso es verificar si la serie a modelar es estacionaria, para ello, nos auxiliamos de las pruebas de Dickey-Fuller.
library(aTSA)
adf_test_pibt = adf.test(cyl_pibt)## Augmented Dickey-Fuller Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag ADF p.value
## [1,] 0 -4.32 0.0100
## [2,] 1 -3.16 0.0100
## [3,] 2 -2.54 0.0130
## [4,] 3 -2.53 0.0135
## Type 2: with drift no trend
## lag ADF p.value
## [1,] 0 -4.29 0.0100
## [2,] 1 -3.13 0.0331
## [3,] 2 -2.52 0.1318
## [4,] 3 -2.50 0.1369
## Type 3: with drift and trend
## lag ADF p.value
## [1,] 0 -4.25 0.010
## [2,] 1 -3.08 0.138
## [3,] 2 -2.48 0.375
## [4,] 3 -2.47 0.378
## ----
## Note: in fact, p.value = 0.01 means p.value <= 0.01Con estos resultados, podemos inferir que la serie es estacionaria. El siguiente paso es identificar la estructura ARMA. Esta estructura puede determinarse mediante el uso de correlogramas, estimar modelos y definir un modelo ganador en base a los criterios de información.
Como el objetivo de estos modelos es predecir, realizamos una predicción seudo fuera de muestra. Nuestro conjunto de entrenamiento está comprendido entre 2006.1 y 2014.4 y el set de prueba será de 2017.1 a 2020.3.
#Definimos el conjunto de entrenamiento
PIBT_train = cyl_pibt
#CORRELOGRAMAS
grid <- matrix(c(1,2), nrow=2,
ncol=2, byrow = TRUE)
layout(grid)
acf(PIBT_train, main="")
pacf(PIBT_train, main="")
Los correlogramas nos sugieren dos posibles modelos, un ARMA(1,0) o ARMA(0,1). Por lo cual, procedemos a estimar ambos modelos.
modelo_1 = arma(PIBT_train, order = c(1,0))
modelo_2 = arma(PIBT_train, order = c(0,1))
summary(modelo_1)##
## Call:
## arma(x = PIBT_train, order = c(1, 0))
##
## Model:
## ARMA(1,0)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0801954 -0.0102504 0.0005497 0.0132220 0.0348675
##
## Coefficient(s):
## Estimate Std. Error t value Pr(>|t|)
## ar1 0.5078933 0.1118700 4.540 5.62e-06 ***
## intercept -0.0002207 0.0024293 -0.091 0.928
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Fit:
## sigma^2 estimated as 0.0003543, Conditional Sum-of-Squares = 0.02, AIC = -297.34summary(modelo_2)##
## Call:
## arma(x = PIBT_train, order = c(0, 1))
##
## Model:
## ARMA(0,1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.620e-02 -1.211e-02 -3.912e-05 1.352e-02 3.684e-02
##
## Coefficient(s):
## Estimate Std. Error t value Pr(>|t|)
## ma1 4.885e-01 1.246e-01 3.922 8.77e-05 ***
## intercept 3.852e-05 3.683e-03 0.010 0.992
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Fit:
## sigma^2 estimated as 0.0003731, Conditional Sum-of-Squares = 0.02, AIC = -294.29Ambos modelos son significativos y basandonos en los criterios de información, se elige el modelo ARMA(1,0).
Realización del pronóstico
En la práctica, la realización del pronóstico es más sencilla de lo que parece. Gracias a las librerías y demás funciones disponibles en la web. Sin embargo, cabe aclarar que esto no garantiza que cumplan con la rigurosidad estadística necesaria, es por ello que en este post se trata de explicar brevemente la matemática detrás de estos modelos y los pasos a seguir usando la metodología Box Jenkins.
library(tsutils)
PIB_f = stlf(PIBT_train, method = "arima",h=4)Una vez realizado el pronóstico, debemos regresar la variable a su estado original, es decir, revertir las transformaciones realizadas. Por ejemplo, si se le aplicó logaritmo, usamos \(exp(variable)\)
ts_plot(PIB_f)Y así concluye este ejercicio sobre como realizar pronósticos con modelos ARIMA, aclarar que este ejemplo está en su nivel más básico, por lo tanto, el pronóstico puede mejorar si se utilizaran técnicas como la proyección con ventana móvil o quizás al aplicar la validación de supuestos se determine otro tipo de modelo ARIMA.
Bibliografía
Gujarati, D. N., & Porter, D. C. (2011). Econometria Básica. 5 Ed. McGraw Hill
McLeod, A. I., Yu, H., & Mahdi, E. (2012). Time series analysis with R. In Handbook of statistics (Vol. 30, pp. 661-712). Elsevier.
Urbina, J. (2020, 12 junio). RPubs - Estructuración de datos. RPubs - Jilber Urbina. https://www.rpubs.com/Jilber_Urbina/Estructuracion_de_datos