Capítulo 3 Procesos Estacionarios y Modelos Univariados

3.1 Definición de ergodicidad y estacionariedad

Las definiciones de ergodicidad y estacionariedad de esta sección, la función de autocorrelación de la siguiente y las pruebas de autocorrelación y normalidad que las acompañan siguen la exposición de Kirchgässner, Wolters y Hassler (2013, sec. 1.4 y 1.5). El tratamiento de los procesos $AR$, $MA$ y $ARMA$ que se desarrolla más adelante sigue principalmente a Enders (2015, cap. 2).

A partir de esta sección introduciremos mayor formalidad matemática al análisis de las series de tiempo. Por ello cambiaremos un poco la notación y ocuparemos a $X_t$ en lugar de $Z_t$ como objeto de nuestro análisis. Con $X_t$ denotaremos a una serie de tiempo, ya que con $Z_t$ denotaremos a una variable, sin que ella fuera necesariamente una serie de tiempo en los términos que a continuación discutimos. Asimismo, iniciaremos por establecer una serie de definiciones.

De esta forma, definiremos a una serie de tiempo como un vector de variables aleatorias de dimensión $T$, dado como:

\[\begin{equation} X_1, X_2, X_3, \ldots ,X_T \tag{3.1} \end{equation}\]

Cada una de las $X_t$ ($t = 1, 2, \ldots, T$) consideradas como una variable aleatoria. Así, también podemos denotar a la serie de tiempo como:

\[\begin{equation} \{ X_t \}^T_{t = 1} \tag{3.2} \end{equation}\]

Es decir, definiremos a una serie de tiempo como una realización de un proceso estocástico –o un Proceso Generador de Datos (PGD). Consideremos una muestra de los múltiples posibles resultados de muestras de tamaño $T$, la colección dada por:

\[\begin{equation} \{X^{(1)}_1, X^{(1)}_2, \ldots, X^{(1)}_T\} \tag{3.3} \end{equation}\]

Digamos que la ecuación (3.3) es una de las tantas posibles resultantes del proceso estocástico o PGD. Eventualmente podríamos estar dispuestos a observar este proceso indefinidamente, de forma tal que estemos interesados en la secuencia dada por $\{ X^{(1)}_t \}^{\infty}_{t = 1}$, lo cual no dejaría de ser sólo una de las tantas realizaciones o secuencias del proceso estocástico original.

Tan solo por poner un ejemplo, podríamos observar las siguientes realizaciones del mismo PGD:

\[\begin{eqnarray*} & \{X^{(2)}_1, X^{(2)}_2, \ldots, X^{(2)}_T\} & \\ & \{X^{(3)}_1, X^{(3)}_2, \ldots, X^{(3)}_T\} & \\ & \{X^{(4)}_1, X^{(4)}_2, \ldots, X^{(4)}_T\} & \\ & \vdots & \\ & \{X^{(j)}_1, X^{(j)}_2, \ldots, X^{(j)}_T\} & \end{eqnarray*}\]

Donde $j \in \mathbb{Z}$. En lo subsecuente, diremos que una serie de tiempo es una realización del proceso estocástico subyacente. Considerando, en consecuencia, al proceso estocástico con todas sus posibilidades de realización.

Para hacer más sencilla la notación no distinguiremos entre el proceso en sí mismo y una de sus realizaciones, es decir, siempre escribiremos a una serie de tiempo como la secuencia mostrada en la ecuación (3.2), o más precisamente como la siguiente realización: \[\begin{equation} \{ X_1, X_2, \ldots, X_T \} \tag{3.4} \end{equation}\]

O simplemente:

\[\begin{equation} X_1, X_2, \ldots, X_T \tag{3.5} \end{equation}\]

El proceso estocástico de dimensión $T$ puede ser completamente descrito por su función de distribución multivariada de dimensión $T$. No obstante, esto no resulta práctico para el análisis que desarrollaremos más adelante. Por ello, en este libro –como se hace en casi todos los textos– sólo nos enfocaremos en sus primer y segundo momentos. Es decir, en su media o valor esperado:

\[\begin{equation*} \mathbb{E}[X_t] \end{equation*}\]

Para $t = 1, 2, \ldots, T$; o:

\[\begin{equation*} \left[ \begin{array}{c} \mathbb{E}[X_1] \\ \mathbb{E}[X_2] \\ \vdots \\ \mathbb{E}[X_T] \end{array} \right] \end{equation*}\]

\[\begin{equation*} \left[ \begin{array}{c} \mathbb{E}[X_1], \mathbb{E}[X_2], \ldots, \mathbb{E}[X_T] \end{array} \right] \end{equation*}\]

Y en su varianza:

\[\begin{equation*} Var[X_t] = \mathbb{E}[(X_t - \mathbb{E}[X_t])^2] \end{equation*}\]

Para $t = 1, 2, \ldots, T$, y de sus $T(T-1)/2$ covarianzas: \[\begin{equation*} Cov[X_t,X_s] = \mathbb{E}[(X_t - \mathbb{E}[X_t])(X_s - \mathbb{E}[X_s])] \end{equation*}\]

Para $t < s$. Por lo tanto, en la forma matricial podemos escribir lo siguiente:

\[\begin{equation*} \left[ \begin{array}{c c c c} Var[X_1] & Cov[X_1,X_2] & \cdots & Cov[X_1,X_T] \\ Cov[X_2,X_1] & Var[X_2] & \cdots & Cov[X_2,X_T] \\ \vdots & \vdots & \ddots & \vdots \\ Cov[X_T,X_1] & Cov[X_T,X_2] & \cdots & Var[X_T] \\ \end{array} \right] \end{equation*}\]

\[\begin{equation} = \left[ \begin{array}{c c c c} \sigma_1^2 & \sigma_{12} & \cdots & \sigma_{1T} \\ \sigma_{21} & \sigma_2^2 & \cdots & \sigma_{2T} \\ \vdots & \vdots & \ddots & \vdots \\ \sigma_{T1} & \sigma_{T2} & \cdots & \sigma_T^2 \\ \end{array} \right] \tag{3.6} \end{equation}\]

Donde $\sigma^2_t = Var[X_t]$ y $\sigma_{ts} = Cov[X_t, X_s]$. Es claro que en la matriz de la ecuación (3.6) existen $T(T-1)/2$ covarianzas distintas, ya que se cumple que $Cov[X_t,X_s] = Cov[X_s,X_t]$, para $t \neq s$.

El prefijo auto en el nombre de estas covarianzas señala una diferencia importante respecto del caso transversal: aquí no se está midiendo la covariación entre dos variables distintas, sino la de una misma variable consigo misma en dos fechas separadas. Es esa estructura de dependencia interna –y no la relación con otras variables– la que el análisis de series de tiempo busca describir y aprovechar.

Restringirse a los dos primeros momentos supone, desde luego, una pérdida de información. La excepción es el caso normal multivariado: si la distribución conjunta del proceso es normal, la media y la matriz de la ecuación (3.6) la determinan por completo, y entonces trabajar con dos momentos no cuesta nada. Fuera de ese caso, el análisis que sigue debe entenderse como una caracterización parcial del proceso, suficiente para los fines de modelación y pronóstico que nos ocupan.

Ahora introduciremos el concepto de ergodicidad, el cual indica que los momentos muestrales –calculados con base en una serie de tiempo con un número finito de observaciones– tienden, en la medida en que $T \rightarrow \infty$, a los verdaderos valores poblacionales, los cuales definiremos como $\mu$, para la media, y $\sigma^2_X$ para la varianza.

Este concepto sólo es cierto si asumimos que, por ejemplo, el valor esperado y la varianza son como se dice a continuación para todo $t = 1, 2, \ldots, T$:

\[\begin{eqnarray} \mathbb{E}[X_t] = \mu_t = \mu \tag{3.7} \end{eqnarray}\]

\[\begin{eqnarray} Var[X_t] = \sigma^2_X \tag{3.8} \end{eqnarray}\]

Más formalmente, se dice que el PGD o el proceso estocástico es ergódico en la media si:

\[\begin{equation} \displaystyle\lim_{T \to \infty}{\mathbb{E} \left[ \left( \frac{1}{T} \sum^{T}_{t = 1} (X_t - \mu) \right) ^2 \right]} = 0 \tag{3.9} \end{equation}\]

y ergódico en la varianza si:

\[\begin{equation} \displaystyle\lim_{T \to \infty}{\mathbb{E} \left[ \left( \frac{1}{T} \sum^{T}_{t = 1} (X_t - \mu) ^2 - \sigma^2_X \right) ^2 \right]} = 0 \tag{3.10} \end{equation}\]

A estas condiciones se les conoce como propiedades de consistencia para las variables aleatorias. Sin embargo, éstas no pueden ser probadas con una única realización del proceso, por lo que se tratan como un supuesto de trabajo. Lo relevante para lo que sigue es que para que los momentos muestrales de una serie de tiempo aproximen a los momentos poblacionales del proceso, se requiere que éste sea estacionario y ergódico.

Podemos distinguir dos tipos de estacionariedad. Si asumimos que la función común de distribución del proceso estocástico no cambia a lo largo del tiempo, se dice que el proceso es estrictamente estacionario. Como este concepto es difícil de aplicar en la práctica, solo consideraremos a la estacionariedad débil o estacionariedad en sus momentos.

Definiremos a la estacionariedad por sus momentos del correspondiente proceso estocástico dado por $\{X_t\}$:

Estacionariedad en media: Un proceso estocástico es estacionario en media si $E[X_t] = \mu_t = \mu$ es constante para todo $t$.
Estacionariedad en varianza: Un proceso estocástico es estacionario en varianza si $Var[X_t] = \mathbb{E}[(X_t - \mu_t)^2] = \sigma^2_X = \gamma(0)$ es constante y finita para todo $t$.
Estacionariedad en covarianza: Un proceso estocástico es estacionario en covarianza si $Cov[X_t,X_s] = \mathbb{E}[(X_t - \mu_t)(X_s - \mu_s)] = \gamma(|s-t|)$ es sólo una función de la distancia $|s-t|$ entre las dos variables aleatorias, y no del momento $t$ en que se observa el proceso.
Estacionariedad débil: Como la estacionariedad en varianza resulta de forma inmediata de la estacionariedad en covarianza cuando se asume que $s = t$, un proceso estocástico es débilmente estacionario cuando es estacionario en media y covarianza.

Puesto que resulta poco factible asumir una estacionariedad diferente a la débil, en adelante, siempre que digamos que un proceso es estacionario, nos referiremos al caso débil, sin el apelativo de débil.

Ejemplo. Supongamos una serie de tiempo denotada por: $\{U_t\}^T_{t = 0}$. Decimos que el proceso estocástico $\{U_t\}$ es un proceso estocástico puramente aleatorio o simplemente un proceso de ruido blanco, si este tiene las siguientes propiedades:

$\mathbb{E}[U_t] = 0$, $\forall t$;
$Var[U_t] = \mathbb{E}[(U_t - \mu_t)^2] = \mathbb{E}[(U_t - \mu)^2] = \mathbb{E}[(U_t)^2] = \sigma^2$, $\forall t$, y
$Cov[U_t,U_s] = \mathbb{E}[(U_t - \mu_t)(U_s - \mu_s)] = \mathbb{E}[(U_t - \mu)(U_s - \mu)] = \mathbb{E}[U_t U_s] = 0$, $\forall t \neq s$.

En otras palabras, un proceso $U_t$ es un ruido blanco si su valor esperado (promedio) es cero (0), tiene una varianza finita y constante, y además no le importa la historia pasada. Así, su valor presente no se ve influenciado por sus valores pasados, no importando respecto de qué período se tome referencia.

En apariencia, por sus propiedades, este proceso es débilmente estacionario –o simplemente, estacionario–. Todas las variables aleatorias tienen una media de cero, una varianza $\sigma^2$ y no existe correlación entre ellas.

Ahora, supongamos que definimos un nuevo proceso estocástico $\{X_t\}$ como:

\[\begin{equation} X_t = \left\{ \begin{array}{l} U_0 \mbox{ para } t = 0 \\ X_{t-1} + U_t \mbox{ para } t = 1, 2, 3, \ldots \end{array}\right. \tag{3.11} \end{equation}\]

Donde $\{ U_t \}$ es un proceso puramente aleatorio. Este proceso estocástico, conocido como caminata aleatoria sin deriva (drift), puede ser reescrito como:

\[\begin{equation} X_t = \sum^t_{j = 0} U_j \tag{3.12} \end{equation}\]

Tratemos de dar más claridad al ejemplo, para ello asumamos que generamos a $\{U_t\}$ por medio del lanzamiento de una moneda. Donde obtenemos una cara (águila) con una probabilidad de $0.5$, en cuyo caso decimos que la variable aleatoria $U_t$ tomará el valor de $+1$, y una cruz (sol) con una probabilidad de $0.5$, en cuyo caso decimos que la variable aleatoria $U_t$ toma el valor de $-1$.

Este planteamiento cumple con las propiedades enunciadas ya que:

$\mathbb{E}[U_t] = 0.5 \times -1 + 0.5 \times 1 = 0$, $\forall t$
$Var[U_t] = \mathbb{E}[(U_t - 0)^2] = \frac{1}{2}((-1)^2) + \frac{1}{2}((1)^2) = 1$, $\forall t$
$Cov[U_t,U_s] = \mathbb{E}[(U_t - 0)(U_s - 0)] = \mathbb{E}[U_t \cdot U_s] = 0$, $\forall t \neq s$.

Retomando a nuestro proceso $X_t$, supondremos que $X_0 = 0$ en $t = 0$. Si verificamos cuáles son sus primeros y segundos momentos de $\{X_t\}$, tenemos:

\[\begin{equation} \mathbb{E}[X_t] = \mathbb{E}\left[ \sum^t_{j=1} U_j \right] = \sum^t_{j=1} \mathbb{E}[U_j] = 0 \tag{3.13} \end{equation}\]

En cuanto a la varianza:

\[\begin{eqnarray} Var[X_t] & = & Var \left[ \sum^t_{j=1} U_j \right] \nonumber \\ & = & \sum^t_{j=1} Var[U_j] + 2 \sum_{j < k} Cov[U_j,U_k] \nonumber \\ & = & \sum^t_{j=1} 1 \nonumber \\ & = & t \tag{3.14} \end{eqnarray}\]

Lo anterior, dado que hemos supuesto que en la caminata aleatoria todas las variables aleatorias son independientes, es decir, $Cov[U_t,U_s] = E[U_t \cdot U_s] = 0$. Por su parte, la covarianza del proceso estocástico se puede ver como:

\[\begin{eqnarray*} Cov[X_t,X_s] & = & \mathbb{E} \left[ \left( \sum^t_{j=1} U_j - 0 \right) \left( \sum^s_{i=1} U_i - 0 \right) \right] \\ & = & \mathbb{E}[(U_1 + U_2 + \ldots + U_t)(U_1 + U_2 + \ldots + U_s)] \\ & = & \sum^t_{j=1} \sum^s_{i=1} \mathbb{E}[U_j U_i] \\ & = & \mathbb{E}[U^2_1] + \mathbb{E}[U^2_2] + \ldots + \mathbb{E}[U^2_{min(t,s)}] \\ & = & min(t,s) \cdot \sigma^2 \\ & = & min(t,s) \end{eqnarray*}\]

Donde la última igualdad usa que en este ejemplo $\sigma^2 = 1$; sólo sobreviven los términos cruzados con $i = j$, que son exactamente $min(t,s)$.

Así, el proceso estocástico dado por la caminata aleatoria sin deriva es estacionario en media, pero no en varianza ni en covarianza y, consecuentemente, no es estacionario –condición contraria al caso del proceso simple descrito en $U_t$–.

Es fácil ver que muchas de las posibilidades de realización de este proceso estocástico (series de tiempo) pueden tomar cualquiera de las rutas consideradas en la Figura 3.1.

set.seed(1234)
# Utilizaremos una función guardada en un archivo aparte
# Llamamos a la función:
source("Caminata.R")

# Definimos argumentos de la función
Opciones <- c(-1, 1)
#
Soporte <- 10000

# Número de trayectorias (realizaciones) del mismo proceso estocástico
Caminos <- 10

# Cada llamada a Caminata() devuelve una lista con el vector de tiempo y la
# suma acumulada de los choques; guardamos cada realización en una columna.
Trayectorias <- sapply(seq_len(Caminos),
                       function(i) Caminata(Opciones, Soporte)[[2]])

Tiempo <- seq_len(Soporte)

matplot(Tiempo, Trayectorias, type = "l", lty = 1,
        col = rainbow(Caminos, v = 0.8),
        xlab = "Tiempo", ylab = "Ganancias acumuladas")
abline(h = 0, lty = 2, col = "gray40")

# Bandas de +/- 2 desviaciones estándar: Var[X_t] = t
lines(Tiempo,  2*sqrt(Tiempo), lty = 2, lwd = 2, col = "red")
lines(Tiempo, -2*sqrt(Tiempo), lty = 2, lwd = 2, col = "red")

$Diez trayectorias simuladas de una caminata aleatoria en la que sólo son posibles cambios de $+1$ y $-1$. Las líneas rojas punteadas son las bandas $\pm 2 \sqrt{t}$, es decir, $\pm 2$ desviaciones estándar del proceso: ilustran cómo la dispersión crece con el tiempo$

Figura 3.1: Diez trayectorias simuladas de una caminata aleatoria en la que sólo son posibles cambios de $+1$ y $-1$. Las líneas rojas punteadas son las bandas $\pm 2 \sqrt{t}$, es decir, $\pm 2$ desviaciones estándar del proceso: ilustran cómo la dispersión crece con el tiempo

3.2 Función de autocorrelación

Para ampliar la discusión, es posible calcular la fuerza o intensidad de la dependencia de las variables aleatorias dentro de un proceso estocástico, ello mediante el uso de las autocovarianzas. Cuando las covarianzas son normalizadas respecto de la varianza, el resultado es un término que es independiente de las unidades de medida aplicadas, y se conoce como la función de autocorrelación.

Para procesos estacionarios, dicha función de autocorrelación está dada por:

\[\begin{equation} \rho(\tau) = \frac{\mathbb{E}[(X_t - \mu)(X_{t+\tau} - \mu)]}{\mathbb{E}[(X_t - \mu)^2]} = \frac{\gamma(\tau)}{\gamma(0)} \tag{3.15} \end{equation}\]

Donde $\tau = \ldots, -2, -1, 0, 1, 2, \ldots$. Dicha función tiene las siguientes propiedades:

$\rho(0) = 1$. Es fácil demostrar que la función $\rho(0)$ es:

\[\begin{equation} \rho(0) = \frac{\mathbb{E}[(X_t - \mu)(X_{t + 0} - \mu)]}{\mathbb{E}[(X_t - \mu)^2]} = \frac{\mathbb{E}[(X_t - \mu)^2]}{\mathbb{E}[(X_t - \mu)^2]} = 1 \end{equation}\]

$\rho(\tau) = \rho(-\tau)$. Partiendo de la definición de $\rho(\tau)$ podemos ver que la distancia que existe entre $t$ y $t + \tau$ es $\tau$, de esta forma la autocorrelación de la variable $X$ entre los periodos antes señalados debería ser la misma para el caso en que $\rho(-\tau)$. Partamos de la ecuación para ver más claramente:

\[\begin{equation} \rho(\tau) = \frac{\mathbb{E}[(X_t - \mu)(X_{t + \tau} - \mu)]}{\mathbb{E}[(X_t - \mu)^2]} = \frac{\mathbb{E}[(X_t - \mu)(X_{t - \tau} - \mu)]}{\mathbb{E}[(X_t - \mu)^2]} = \rho(-\tau) \end{equation}\]

$\lvert\rho(\tau)\lvert \leq 1$, para todo $\tau$.

Derivado de las propiedades 1 y 2 antes descritas se puede concluir que sólo es necesario conocer la función de autocorrelación para el caso de $\tau = 1, 2, 3, \ldots$, ya que de estos casos podemos derivar los valores de la función de autocorrelación complementarios de $\tau = \ldots, -3, -2, -1$.

Partiendo de los supuestos de ergodicidad en relación con la media, varianza y covarianzas de un proceso estacionario, podemos estimar dichos parámetros con las siguientes formulaciones o propuestas de estimadores puntuales:

\[\begin{equation} \hat{\mu} = \frac{1}{T} \sum^T_{t=1} X_t \tag{3.16} \end{equation}\]

\[\begin{equation} \hat{\gamma}(0) = \frac{1}{T} \sum^T_{t=1} (X_t - \hat{\mu})^2 = \hat{\sigma}^2 \tag{3.17} \end{equation}\] \[\begin{equation} \hat{\gamma}(\tau) = \frac{1}{T} \sum^{T - \tau}_{t=1} (X_t - \hat{\mu})(X_{t+\tau} - \hat{\mu}) \mbox{, para } \tau = 1, 2, \ldots, T-1 \tag{3.18} \end{equation}\]

No hacemos la demostración en este libro –sería deseable que el alumno revisara la siguiente afirmación–, pero estos últimos son estimadores consistentes de $\mu$, $\gamma(0)$ y $\gamma(\tau)$. Por su parte, un estimador consistente de la función de autocorrelación estará dado por: \[\begin{equation} \hat{\rho}(\tau) = \frac{\sum^{T - \tau}_{t=1} (X_t - \hat{\mu})(X_{t+\tau} - \hat{\mu})}{\sum^T_{t=1} (X_t - \hat{\mu})^2} = \frac{\hat{\gamma}(\tau)}{\hat{\gamma}(0)} \tag{3.19} \end{equation}\]

El estimador de la ecuación (3.19) es asintóticamente insesgado, y su comportamiento en muestras grandes es lo que permite usarlo para contrastar hipótesis. Conviene detenerse en el caso de referencia: si la serie fuera ruido blanco, todas las autocorrelaciones poblacionales serían cero y $\hat{\rho}(\tau)$ se distribuiría aproximadamente como una normal de media cero y varianza $1/T$. De ahí se obtiene la regla práctica que usaremos al leer cualquier correlograma: bajo la hipótesis de ruido blanco, alrededor del 95% de los coeficientes estimados debería caer dentro de la banda $\pm 2/\sqrt{T}$.

Dos advertencias sobre esa banda, que en la práctica se pasan por alto con frecuencia. La primera es que se construye bajo la hipótesis nula de ausencia de autocorrelación; si la serie realmente está autocorrelacionada, la varianza de $\hat{\rho}(\tau)$ es mayor y la banda subestima la incertidumbre. La segunda es que se trata de un intervalo puntual, coeficiente por coeficiente: al inspeccionar simultáneamente veinte o cuarenta rezagos, es esperable que algunos crucen la banda por azar. Ése es precisamente el motivo por el que se recurre a las pruebas conjuntas que veremos enseguida, en vez de decidir a partir de la inspección visual del correlograma.

El uso más frecuente de estas herramientas no es describir la serie original, sino diagnosticar un modelo ya estimado. La lógica es la siguiente: si un modelo capturó toda la dinámica sistemática de la serie, lo que queda en sus residuales no debería tener estructura alguna; cualquier autocorrelación remanente es señal de que quedó información sin modelar. Como acabamos de ver que examinar los coeficientes uno por uno tiene el inconveniente de las comparaciones múltiples, lo que se contrasta es la hipótesis conjunta de que los primeros $m$ coeficientes son simultáneamente nulos: \[\begin{equation} H_0 : \rho(\tau) = 0 \mbox{, para todo } \tau = 1, 2, \ldots, m \mbox{ y } m < T \tag{3.20} \end{equation}\]

Esta expresión se puede interpretar como una prueba respecto de si la correlación entre la información de periodos atrás es cero con la información contemporánea. Para hacer una prueba global de la hipótesis de si un número $m$ de coeficientes de autocovarianzas son cero, Box y Pierce (1970) desarrollaron la siguiente estadística:

\[\begin{equation} Q^* = T \sum_{j = 1}^{m} \hat{\rho} (j)^2 \tag{3.21} \end{equation}\]

Bajo la hipótesis nula esta estadística se distribuye asintóticamente como una chi cuadrado ($\chi^2$) con $m-k$ grados de libertad y con $k$ que representa al número de parámetros estimados.

La distribución de esta estadística sólo se sostiene asintóticamente. Por ello, Ljung y Box (1978) propusieron la siguiente modificación de la estadística para muestras pequeñas:

\[\begin{equation} Q = T(T + 2) \sum_{j = 1}^{m} \frac{\hat{\rho} (j)^2}{T - j} \tag{3.22} \end{equation}\]

La cual también se distribuye asintóticamente como $\chi^2$ con $m-k$ grados de libertad.

También es intuitivamente claro que la hipótesis nula de no autocorrelación de residuales debería ser rechazada si alguno de los valores $\hat{\rho} (j)$ es muy grande, es decir, si $Q$ o $Q^*$ es muy grande. O más precisamente, si estas estadísticas son más grandes que los correspondientes valores críticos de la distribución $\chi^2$ con $m-k$ grados de libertad a algún grado dado de significancia.

Una alternativa para esta prueba es una del tipo Multiplicadores de Lagrange (o LM) desarrollada por Breusch (1978) y Godfrey (1978), en la cual, al igual que en las estadísticas $Q$ y $Q^*$, la hipótesis nula está dada por:

$H_0$: Los residuales no están autocorrelacionados.

$H_a$: Los residuales muestran alguna autocorrelación de forma autorregresiva o de medias móviles.

La prueba consiste en realizar una regresión auxiliar en la cual los residuales se estiman en función de las variables explicativas del modelo original y en los residuales mismos pero rezagados hasta el término $m$ (regresión auxiliar). La prueba resulta en una estadística con una distribución $\chi^2$ con $m$ grados de libertad la cual está dada por la expresión:

\[\begin{equation} LM = T \times R^2 \tag{3.23} \end{equation}\]

Donde $R^2$ es el resultante de la regresión auxiliar y $T$ es el número de observaciones totales.

En comparación con una prueba Durbin - Watson que es comúnmente usada en la econometría tradicional, para probar autocorrelación de los residuales, las estadísticas $Q$, $Q^*$ y $LM$ tienen las siguientes ventajas:

Su alcance no se limita al primer rezago. La Durbin - Watson está construida específicamente contra una alternativa $AR(1)$, de modo que puede no detectar patrones que aparecen en rezagos más largos –típicamente los estacionales, como el rezago 12 en datos mensuales–, mientras que $Q$, $Q^*$ y $LM$ examinan de una vez el bloque completo de los primeros $m$ rezagos;
Siguen siendo válidas cuando entre los regresores hay rezagos de la variable dependiente, que es la situación normal en series de tiempo. En ese caso la Durbin - Watson está sesgada hacia el valor 2 –es decir, hacia no detectar autocorrelación– y deja de ser utilizable, y
No tienen una región de indecisión: la prueba Durbin - Watson se contrasta contra dos valores críticos ($d_L$ y $d_U$) que delimitan una zona en la que la prueba no permite concluir, mientras que las estadísticas $Q$, $Q^*$ y $LM$ se comparan contra un único valor crítico de la distribución $\chi^2$.

El hecho de que los residuales no estén autocorrelacionados no implica que estos sean independientes y normalmente distribuidos: la ausencia de autocorrelación sólo implica independencia estocástica cuando las variables se distribuyen de forma normal.

A menudo se asume que estos residuales están distribuidos normalmente, ya que la mayoría de las pruebas estadísticas tienen este supuesto detrás. No obstante, ello también depende de los otros momentos de la distribución, específicamente del tercer y cuarto momento. Los cuales se expresan como:

\[\begin{equation*} \mathbb{E}[(X_t - \mathbb{E}[X_t])^i] \mbox{, } i = 3, 4 \end{equation*}\]

El tercer momento es necesario para determinar el sesgo, el cual está dado como:

\[\begin{equation} \hat{S} = \frac{1}{T} \frac{\sum_{t = 1}^{T} (X_t - \hat{\mu})^3}{\sqrt{\hat{\gamma}(0)^3}} \tag{3.24} \end{equation}\]

Para distribuciones simétricas (como en el caso de la distribución normal) el valor teórico para el sesgo es cero.

La curtosis, la cual está dada en función del cuarto momento, se puede expresar como:

\[\begin{equation} \hat{K} = \frac{1}{T} \frac{\sum_{t = 1}^{T} (X_t - \hat{\mu})^4}{\hat{\gamma}(0)^2} \tag{3.25} \end{equation}\]

Para el caso de una distribución normal, esta estadística toma el valor de 3. Valores más grandes que 3 indican que la distribución tiene colas anchas. En tales casos se ubican los datos financieros.

Usando el valor de las estadísticas para medir el sesgo y la curtosis, $S$ y $K$, respectivamente, Jarque y Bera (1980) propusieron una prueba de normalidad, la cual puede ser aplicada a series de tiempo en niveles o en diferencias indistintamente. Dicha prueba se expresa como: \[\begin{equation} JB = \frac{T}{6} \left(\hat{S}^2 + \frac{1}{4} (\hat{K} - 3)^2 \right) \tag{3.26} \end{equation}\]

La cual tiene una distribución $\chi^2$ con $2$ grados de libertad y donde $T$ es el tamaño de la muestra. La hipótesis de que las observaciones están distribuidas de forma normal se rechaza si el valor de la estadística de prueba es más grande que los correspondientes valores críticos en tablas.

library(ggplot2)
library(dplyr)
library(readxl)

Datos <- read_excel("BD/Base_Transporte.xlsx", 
                    sheet = "Datos", col_names = TRUE)

## Rótulos del periodo muestral, calculados de los datos para que el texto,
## los cuadros y los títulos de las figuras no queden desfasados cuando se
## actualice la base.
meses_es  <- c("enero", "febrero", "marzo", "abril", "mayo", "junio",
               "julio", "agosto", "septiembre", "octubre", "noviembre",
               "diciembre")
meses_abr <- c("Ene", "Feb", "Mar", "Abr", "May", "Jun",
               "Jul", "Ago", "Sep", "Oct", "Nov", "Dic")

T_trans   <- nrow(Datos)
fin_trans <- as.Date(max(Datos$Periodo))
fin_mes   <- as.integer(format(fin_trans, "%m"))
fin_anio  <- as.integer(format(fin_trans, "%Y"))

fin_txt   <- paste(meses_es[fin_mes], "de", fin_anio)
per_txt   <- paste0("Ene-2000 a ", meses_abr[fin_mes], "-", fin_anio)

## Horizonte y arranque del pronóstico de la última sección del capítulo: se
## toman del archivo de variables dicotómicas futuras, que debe empezar el mes
## siguiente al último dato observado.
Predict_Datos <- read_excel("BD/Predict_Base_Transporte_ARIMA.xlsx",
                            sheet = "Datos", col_names = TRUE)

h_f   <- nrow(Predict_Datos)
ini_f <- as.Date(min(Predict_Datos$Periodo))
ini_f <- c(as.integer(format(ini_f, "%Y")), as.integer(format(ini_f, "%m")))

# Verificación: el archivo de dummies futuras debe empezar justo después del
# último dato observado.
stopifnot(all(ini_f == c(fin_anio + (fin_mes == 12), fin_mes %% 12 + 1)))

Ejemplo. Veamos un ejemplo para ilustrar el uso de la función de autocorrelación. Tomemos como variable al número de pasajeros transportados por el sistema de transporte del metro de la CDMX.² Los datos empleados fueron tomados del INEGI y son una serie de tiempo en el período que va de enero de 2000 a mayo de 2026, es decir, 317 observaciones. Como se puede apreciar en la Figura 3.2, el número de pasajeros por mes ha oscilado significativamente a lo largo del tiempo. Incluso podemos observar un cambio estructural de la serie entre 2011 y 2012. Asimismo, podemos ubicar una caída atípica que ocurrió en septiembre de 2017. Pero lo más relevante es la caída asociada a la pandemia de COVID-19 de 2020.

ggplot(data = Datos, aes(x = Periodo, y = Pax_Metro)) + 
  geom_line(linewidth = 0.5, color = "darkblue") +
  #geom_point(size = 1.0, color = "darkblue") + 
  #theme_bw() + 
  xlab("Tiempo") + 
  ylab("Millones de pasajeros") + 
  theme(plot.title = element_text(size = 11, face = "bold", hjust = 0)) + 
  theme(plot.subtitle = element_text(size = 10, hjust = 0)) + 
  theme(plot.caption = element_text(size = 10, hjust = 0)) +
  theme(plot.margin = unit(c(1,1,1,1), "cm")) +
  labs(
    title = "Pasajeros Transportados en el Metro de la CDMX",
    subtitle = paste0("(", per_txt, ")"),
    caption = "Fuente: Elaboración propia con información del INEGI, \nhttps://www.inegi.org.mx/app/indicadores/?tm=0&t=1090"
  )

Figura 3.2: Evolución del número de pasajeros en el Metro de la CDMX, enero de 2000 a mayo de 2026

A esta serie de tiempo le calculamos los principales estadísticos hasta ahora estudiados y obtenemos el Cuadro 3.1. En dicho cuadro se destaca que se muestra la función de autocorrelación para los tres primeros rezagos. Para mayor detalle, en la Figura 3.3 se muestra la función de autocorrelación, en donde las bandas descritas por las líneas azules son el intervalo de confianza dentro de las cuales no se puede rechazar la hipótesis nula de que $H_0: \rho(\tau) = 0$, para todo $\tau = 1, 2, \ldots, T-1$.

Cuadro 3.1: Estadísticas descriptivas del número de pasajeros en el Metro de la CDMX, enero de 2000 a mayo de 2026
Estadística	Valor
$\hat{\mu} = \frac{1}{T} \sum^T_{t=1} X_t$	115.6156
$\hat{\gamma}(0) = \frac{1}{T} \sum^T_{t=1} (X_t - \hat{\mu})^2$	418.1772
$\hat{\gamma}(1) = \frac{1}{T} \sum^{T - 1}_{t=1} (X_t - \hat{\mu})(X_{t+1} - \hat{\mu})$	372.7523
$\hat{\gamma}(2) = \frac{1}{T} \sum^{T - 2}_{t=1} (X_t - \hat{\mu})(X_{t+2} - \hat{\mu})$	363.3531
$\hat{\gamma}(3) = \frac{1}{T} \sum^{T - 3}_{t=1} (X_t - \hat{\mu})(X_{t+3} - \hat{\mu})$	341.5693
$\hat{\rho}(1) = \hat{\gamma}(1)/\hat{\gamma}(0)$	0.8914
$\hat{\rho}(2) = \hat{\gamma}(2)/\hat{\gamma}(0)$	0.8689
$\hat{\rho}(3) = \hat{\gamma}(3)/\hat{\gamma}(0)$	0.8168
$Q^* = T \sum_{j = 1}^{1} \hat{\rho}(j)^2$	251.8716
$Q^* = T \sum_{j = 1}^{2} \hat{\rho}(j)^2$	491.2010

Pax_Metro <- ts(Datos$Pax_Metro,
                start = 2000,
                freq = 12)

acf(Pax_Metro,
    lag.max = 150,
    xlab = 'Rezagos k en meses',
    main = "Función de Autocorrelación del número de pasajeros del metro")

Figura 3.3: Función de Autocorrelación: 150 rezagos del número de pasajeros en el Metro de la CDMX, enero de 2000 a mayo de 2026

3.3 Procesos estacionarios univariados

En este capítulo analizaremos el método o metodología de análisis de series de tiempo propuesto por Box y Jenkins (1970). Los modelos propuestos dentro de esta metodología o conjunto de métodos se han vuelto indispensables para efectos de realizar pronósticos de corto plazo.

En este sentido, se analizarán los métodos más importantes en series de tiempo: procesos autorregresivos (AR, por sus siglas en inglés) y procesos de medias móviles (MA, por sus siglas en inglés). Asimismo, se realizará un análisis de los procesos que resultan de la combinación de ambos, conocida como ARMA, los cuales son más comúnmente usados para realizar pronósticos.

3.4 Procesos Autorregresivos (AR)

Los procesos autorregresivos tienen su origen en los trabajos de Yule (1927) y Slutsky, quienes mostraron que series aparentemente cíclicas pueden generarse acumulando choques puramente aleatorios. Su incorporación al análisis de regresión se debe a Cochrane y Orcutt (1949), quienes modelaron los residuales de una regresión clásica como un proceso autorregresivo con el fin de corregir la autocorrelación. En este libro asumimos conocido el modelo de regresión clásica (véase la sección de prerrequisitos de la Introducción).

3.4.1 AR(1)

Como primer caso analizaremos al proceso autorregresivo de primer orden, $AR(1)$, el cual podemos definir como una Ecuación Lineal en Diferencia de Primer Orden Estocástica. Diremos que una Ecuación Lineal en Diferencia de Primer Orden es estocástica si en su representación analítica considera un componente estocástico como en la ecuación (3.27) descrita a continuación:

\[\begin{equation} X_t = a_0 + a_1 X_{t-1} + U_t \tag{3.27} \end{equation}\]

Donde $a_0$ es un término constante, $U_t$ es un proceso estacionario, con media cero (0), una varianza finita y constante ($\sigma^2$) y una covarianza que depende de la distancia entre $t$ y cualquier $t-s$ ($\gamma_s$)–que no depende de los valores pasados o futuros de la variable–, $X_0$ es el valor inicial del proceso $X_t$. No obstante, en ocasiones vamos a asumir que la covarianza será cero (0), por lo que en esos casos tendremos un proceso puramente aleatorio. Considerando la ecuación (3.27) y un proceso de sustitución sucesivo podemos establecer lo siguiente, empezando con $X_1$:

\[\begin{eqnarray*} X_{1} & = & a_0 + a_1 X_{0} + U_{1} \end{eqnarray*}\]

Para $X_2$:

\[\begin{eqnarray*} X_{2} & = & a_0 + a_1 X_{1} + U_{2} \\ & = & a_0 + a_1 (a_0 + a_1 X_{0} + U_{1}) + U_{2} \\ & = & a_0 + a_1 a_0 + a_1^2 X_{0} + a_1 U_{1} + U_{2} \end{eqnarray*}\]

Para $X_3$:

\[\begin{eqnarray*} X_{3} & = & a_0 + a_1 X_{2} + U_{3} \\ & = & a_0 + a_1 (a_0 + a_1 a_0 + a_1^2 X_{0} + a_1 U_{1} + U_{2}) + U_{3} \\ & = & a_0 + a_1 a_0 + a_1^2 a_0 + a_1^3 X_{0} + a_1^2 U_{1} + a_1 U_{2} + U_{3} \end{eqnarray*}\]

Así, para cualquier $X_t$, $t = 1, 2, 3, \ldots$, obtendríamos: \[\begin{eqnarray} X_{t} & = & a_0 + a_1 X_{t - 1} + U_{t} \nonumber \\ & = & a_0 + a_1 (a_0 + a_1 a_0 + a_1^2 a_0 + \ldots + a_1^{t-2} a_0 + a_1^{t-1} X_{0} \nonumber \\ & & + a_1^{t-2} U_{1} + \ldots + a_1 U_{t - 2} + U_{t - 1}) + U_{t} \nonumber \\ & = & a_0 + a_1 a_0 + a_1^2 a_0 + a_1^3 a_0 + \ldots + a_1^{t-1} a_0 + a_1^{t} X_{0} \nonumber \\ & & + a_1^{t-1} U_{1} + \ldots a_1^2 U_{t - 2} + a_1 U_{t - 1} + U_{t} \nonumber \\ & = & (1 + a_1 + a_1^2 + a_1^3 + \ldots + a_1^{t-1}) a_0 + a_1^{t} X_{0} \nonumber \\ & & + a_1^{t-1} U_{1} + \ldots + a_1^2 U_{t - 2} + a_1 U_{t - 1} + U_{t} \nonumber\\ & = & \frac{1 - a_1^t}{1 - a_1} a_0 + a_1^{t} X_{0} + \sum^{t-1}_{j = 0} a_1^{j} U_{t - j} \tag{3.28} \end{eqnarray}\]

De esta forma en la ecuación (3.28) observamos un proceso que es explicado por dos partes: una que depende del tiempo y otra que depende de un proceso estocástico. Asimismo, debe notarse que la condición de convergencia es idéntica al caso de ecuaciones en diferencia estudiadas en el capítulo previo: $|a_1| < 1$, por lo que cuando $t \to \infty$, la expresión (3.28) será la siguiente: \[\begin{equation} X_t = a_0 \frac{1}{1 - a_1} + \sum^{\infty}_{j = 0} a_1^{j} U_{t - j} \tag{3.29} \end{equation}\]

Así, desaparece la parte dependiente del tiempo y únicamente prevalece la parte que es dependiente del proceso estocástico. Esta es la solución de largo plazo del proceso $AR(1)$, la cual depende del proceso estocástico. Notemos, además, que esta solución implica que la variable o la serie de tiempo $X_t$ es también un proceso estocástico que hereda las propiedades de $U_t$. Así, $X_t$ es también un proceso estocástico estacionario, como demostraremos más adelante.

Observemos que la ecuación (3.29) se puede reescribir si consideramos la formulación que en la literatura se denomina como la descomposición de Wold, en la cual se define que es posible asumir que $\psi_j = a_1^j$ y se considera el caso en el cual $|a_1| < 1$. De esta forma tendremos, por ejemplo, que:

\[\begin{equation*} \sum^{\infty}_{j = 0} \psi^2_j = \sum^{\infty}_{j = 0} a_1^{2j} = \frac{1}{1 - a_1^2} \end{equation*}\]

Alternativamente y de forma similar a las ecuaciones en diferencia estudiadas previamente, podemos escribir el proceso $AR(1)$ mediante el uso del operador rezago como:

\[\begin{eqnarray} X_t & = & a_0 + a_1 L X_t + U_t \nonumber \\ X_t - a_1 L X_t & = & a_0 + U_t \nonumber \\ (1 - a_1 L) X_t & = & a_0 + U_t \nonumber \\ X_t & = & \frac{a_0}{1 - a_1 L} + \frac{1}{1 - a_1 L} U_t \tag{3.30} \end{eqnarray}\]

En esta última ecuación retomamos el siguiente término para reescribirlo como:

\[\begin{equation} \frac{1}{1 - a_1 L} = 1 + a_1 L + a_1^2 L^2 + a_1^3 L^3 + \ldots \end{equation}\]

Tomando este resultado para sustituirlo en la ecuación (3.30), obtenemos la siguiente expresión:

\[\begin{eqnarray} X_t & = & (1 + a_1 L + a_1^2 L^2 + a_1^3 L^3 + \ldots) a_0 + (1 + a_1 L + a_1^2 L^2 + a_1^3 L^3 + \ldots) U_t \nonumber \\ & = & (1 + a_1 + a_1^2 + a_1^3 + \ldots) a_0 + U_t + a_1 U_{t-1} + a_1^2 U_{t-2} + a_1^3 U_{t-3} + \ldots \nonumber \\ & = & \frac{a_0}{1 - a_1} + \sum^{\infty}_{j = 0} a_1^j U_{t-j} \tag{3.31} \end{eqnarray}\]

Donde la condición de convergencia y estabilidad del proceso descrito en esta ecuación es que $|a_1| < 1$. Por lo que hemos demostrado que, mediante el uso del operador de rezago, es posible llegar al mismo resultado que obtuvimos mediante el procedimiento de sustituciones iterativas.

La ecuación (3.31) se puede interpretar como sigue. La solución o trayectoria de equilibrio de un AR(1) se divide en dos partes. La primera es una constante que depende de los valores de $a_0$ y $a_1$. La segunda parte es la suma ponderada de las desviaciones o errores observados y acumulados en el tiempo hasta el momento $t$.

Ahora obtendremos los momentos que describen a la serie de tiempo cuando se trata de un proceso $AR(1)$. Para ello debemos obtener la media, la varianza y las covarianzas de $X_t$. Para los siguientes resultados debemos recordar y tener en mente que si $U_t$ es un proceso puramente aleatorio, entonces:

$\mathbb{E}[U_t] = 0$ para todo $t$
$Var[U_t] = \sigma^2$ para todo $t$
$Cov[U_t, U_s] = 0$ para todo $t \neq s$

Dicho lo anterior y partiendo de la ecuación (3.31), el primer momento o valor esperado de la serie de tiempo será el siguiente: \[\begin{eqnarray} \mathbb{E}[X_t] & = & \mathbb{E} \left[ \frac{a_0}{1 - a_1} + \sum^{\infty}_{j = 0} a_1^j U_{t-j} \right] \nonumber \\ & = & \frac{a_0}{1 - a_1} + \sum^{\infty}_{j = 0} a_1^j \mathbb{E}[U_{t-j}] \nonumber \\ & = & \frac{a_0}{1 - a_1} = \mu \tag{3.32} \end{eqnarray}\]

Respecto de la varianza podemos escribir la siguiente expresión a partir de la ecuación (3.31):

\[\begin{eqnarray} Var[X_t] & = & \mathbb{E}[(X_t - \mu)^2] \nonumber \\ & = & \mathbb{E} \left[ \left( \frac{a_0}{1 - a_1} + \sum^{\infty}_{j = 0} a_1^j U_{t-j} - \frac{a_0}{1 - a_1} \right)^2 \right] \nonumber \\ & = & \mathbb{E}[(U_{t} + a_1 U_{t-1} + a_1^2 U_{t-2} + a_1^3 U_{t-3} + \ldots)^2] \nonumber \\ & = & \mathbb{E}[U^2_{t} + a_1^2 U^2_{t-1} + a_1^4 U^2_{t-2} + a_1^6 U^2_{t-3} + \ldots \nonumber \\ & & + 2 a_1 U_t U_{t-1} + 2 a_1^2 U_t U_{t-2} + \ldots] \nonumber \\ & = & \mathbb{E}[U^2_{t}] + a_1^2 \mathbb{E}[U^2_{t-1}] + a_1^4 \mathbb{E}[U^2_{t-2}] + a_1^6 \mathbb{E}[U^2_{t-3}] + \ldots \nonumber \\ & = & \sigma^2 + a_1^2 \sigma^2 + a_1^4 \sigma^2 + a_1^6 \sigma^2 + \ldots \nonumber \\ & = & \sigma^2 (1 + a_1^2 + a_1^4 + a_1^6 + \ldots) \nonumber \\ & = & \sigma^2 \frac{1}{1 - a_1^2} = \gamma(0) \tag{3.33} \end{eqnarray}\]

Previo a analizar la covarianza de la serie de tiempo, recordemos que para el proceso puramente aleatorio $U_t$ su varianza y covarianza pueden verse como $\mathbb{E}[U_t U_s] = \sigma^2$, para $t = s$, y $\mathbb{E}[U_t U_s] = 0$, para cualquier otro caso, respectivamente.

Dicho lo anterior, partiendo de la ecuación (3.31) la covarianza de la serie estará dada por:

\[\begin{eqnarray} Cov(X_t, X_{t-\tau}) & = & \mathbb{E}[(X_t - \mu)(X_{t-\tau} - \mu)] \nonumber \\ & = & \mathbb{E} \left[ \left( \frac{a_0}{1 - a_1} + \sum^{\infty}_{j = 0} a_1^j U_{t-j} - \frac{a_0}{1 - a_1} \right) \right. \nonumber \\ & & \left. \times \left( \frac{a_0}{1 - a_1} + \sum^{\infty}_{j = 0} a_1^j U_{t-\tau-j} - \frac{a_0}{1 - a_1} \right) \right] \nonumber \\ & = & a_1^{\tau} \mathbb{E}[U^2_{t-\tau} + a_1^2 U^2_{t-\tau-1} + a_1^4 U^2_{t-\tau-2} + a_1^6 U^2_{t-\tau-3} + \ldots] \nonumber \\ & = & a_1^{\tau} \sigma^2 \frac{1}{1 - a_1^2} = \gamma(\tau) \tag{3.34} \end{eqnarray}\]

Nótese que con estos resultados en las ecuaciones (3.33) y (3.34) podemos construir la función de autocorrelación teórica como sigue:

\[\begin{eqnarray} \rho(\tau) & = & \frac{\gamma(\tau)}{\gamma(0)} \nonumber \\ & = & a_1^\tau \end{eqnarray}\]

Donde $\tau = 1, 2, 3, \ldots$ y $|a_1| < 1$. Este último resultado significa que cuando el proceso autorregresivo es de orden 1 (es decir, AR(1)) la función de autocorrelación teóricamente es igual al parámetro $a_1$ elevado al número de rezagos considerados. No obstante, esto no significa que la autocorrelación estimada en una muestra coincida exactamente con la teórica: en muestras finitas, la autocorrelación observada será ligeramente distinta. Ahora veamos algunos ejemplos.

Ejemplo. En el primer ejemplo simularemos una serie y mostraremos el análisis de un proceso construido considerando un proceso puramente aleatorio como componente $U_t$. Por su parte, en un segundo ejemplo aplicaremos el análisis a una serie de tiempo de una variable económica observada.

Para el primer ejemplo, consideremos un proceso dado por la forma de un $AR(1)$ como en la ecuación (3.30) cuya solución está dada por la ecuación (3.31). En específico, supongamos que el término o componente estocástico $U_t$ es una serie generada a partir de números aleatorios de una función normal con media $0$ y desviación estándar $4$. Los detalles del proceso simulado se muestran en las siguientes gráficas.

La Figura 3.4 ilustra la trayectoria que resulta de construir la serie de forma iterativa a partir de la ecuación (3.27). Por su parte, la Figura 3.5 contrasta esa misma realización con los momentos teóricos que derivamos arriba: la media $\mu = a_0/(1 - a_1)$ de la ecuación (3.32) y las bandas de $\pm 2$ desviaciones estándar no condicionales, $\pm 2 \sqrt{\gamma(0)}$, con $\gamma(0)$ dada por la ecuación (3.33). Finalmente, las Figuras 3.6 y 3.7 comparan la función de autocorrelación estimada sobre la serie simulada con la función de autocorrelación teórica $\rho(\tau) = a_1^{\tau}$.

library(ggplot2)
library(dplyr)
library(latex2exp)

a0 <- 5; a1 <- 0.9; sigma <- 4; X_0 <- (a0/(1 - a1)); T <- 1000

X_t <- data.frame(Tiempo = c(0:T))


set.seed(12345)

# Agregamos el término estocástico (ruido blanco) al data frame

X_t$U_t <- rnorm(T + 1, mean = 0, sd = sigma)

# Columna para la serie simulada; inicia en el valor inicial X_0
X_t$XR_t <- NA
X_t$XR_t[1] <- X_0

# Columna para la función de autocorrelación teórica: rho(tau) = a1^tau
X_t$rho <- NA

for (i in 2:(T + 1)) {
  # Recursión del AR(1):
  X_t$XR_t[i] = a0 + a1*X_t$XR_t[i-1] + X_t$U_t[i-1]

  # Autocorrelación teórica del rezago (i-1):
  X_t$rho[i-1] = a1^(i-1)
}

# Momentos teóricos del proceso (ecuaciones de la media y la varianza)
mu_teo    <- a0/(1 - a1)
sd_teo    <- sqrt(sigma^2/(1 - a1^2))

ggplot(data = X_t, aes(x = Tiempo, y = XR_t)) + 
  geom_line(linewidth = 0.5, color = "darkred") +
  #theme_bw() + 
  xlab("Tiempo") + 
  ylab(TeX("$X_t$")) + 
  theme(plot.title = element_text(size = 11, face = "bold", 
                                  hjust = 0)) + 
  theme(plot.subtitle = element_text(size = 10, hjust = 0)) + 
  theme(plot.caption = element_text(size = 10, hjust = 0)) +
  theme(plot.margin = unit(c(1,1,1,1), "cm")) +
  labs(
    title = "Trayectoria simulada de un AR(1) con a0 = 5 y a1 = 0.9",
    subtitle = "Con un error con Distribución Normal (media = 0, desviación estándar = 4)",
    caption = "Fuente: Elaboración propia."
  )

Figura 3.4: Trayectoria simulada de un proceso $AR(1)$ con $a_0 = 5$ y $a_1 = 0.9$

ggplot(data = X_t, aes(x = Tiempo, y = XR_t)) + 
  geom_line(linewidth = 0.4, color = "gray50") +
  geom_hline(yintercept = mu_teo, color = "darkblue", linewidth = 0.8) +
  geom_hline(yintercept = mu_teo + 2*sd_teo, color = "darkblue",
             linetype = "dashed") +
  geom_hline(yintercept = mu_teo - 2*sd_teo, color = "darkblue",
             linetype = "dashed") +
  xlab("Tiempo") + 
  ylab(TeX("$X_t$")) + 
  theme(plot.title = element_text(size = 11, face = "bold", 
                                  hjust = 0)) + 
  theme(plot.subtitle = element_text(size = 10, hjust = 0)) + 
  theme(plot.caption = element_text(size = 10, hjust = 0)) +
  theme(plot.margin = unit(c(1,1,1,1), "cm")) +
  labs(
    title = "Serie simulada y momentos teóricos del AR(1)",
    subtitle = paste0("Media teórica = ", round(mu_teo, 2),
                      "; desviación estándar teórica = ", round(sd_teo, 2)),
    caption = "Fuente: Elaboración propia."
  )

$Serie simulada frente a sus momentos teóricos: media $\mu = a_0/(1-a_1)$ (línea azul) y bandas de $\pm 2$ desviaciones estándar no condicionales (líneas punteadas)$

Figura 3.5: Serie simulada frente a sus momentos teóricos: media $\mu = a_0/(1-a_1)$ (línea azul) y bandas de $\pm 2$ desviaciones estándar no condicionales (líneas punteadas)

acf(X_t$XR_t, lag.max = 30, col = "blue", 
    ylab = "Autocorrelación",
    xlab="Rezagos", 
    main="Función de Autocorrelación estimada")

Figura 3.6: Función de autocorrelación estimada sobre la serie simulada

barplot(X_t$rho[1:30], names.arg = c(1:30), col = "blue", 
        border="blue", density = c(10,20), 
        ylab = "Autocorrelación", 
        xlab="Rezagos", 
        main="Función de Autocorrelación Teórica")

$Función de autocorrelación teórica del $AR(1)$: $\rho(\tau) = a_1^{\tau}$$

Figura 3.7: Función de autocorrelación teórica del $AR(1)$: $\rho(\tau) = a_1^{\tau}$

Recordemos que la solución de un $AR(1)$ es la de la ecuación (3.31): el nivel de la serie es un promedio ponderado de los choques pasados con pesos $a_1^j$ que decaen geométricamente. En consecuencia, el valor del parámetro $a_1$ determina qué tan persistente es el efecto de un choque: mientras más cercano a uno, más lentamente se disipa. La Figura 3.8 ilustra esta idea comparando dos procesos $AR(1)$ construidos con los mismos choques $U_t$ y la misma media teórica, pero con $a_1 = 0.9$ y $a_1 = 0.5$: el primero se aleja de su media durante episodios prolongados, mientras que el segundo regresa a ella con rapidez.

# Segundo proceso: mismo ruido blanco y misma media teórica (mu = 50),
# pero con menor persistencia (a1 = 0.5 implica a0 = 25)
a1_b <- 0.5
a0_b <- mu_teo*(1 - a1_b)

X_t$XB_t <- NA
X_t$XB_t[1] <- mu_teo

for (i in 2:(T + 1)) {
  X_t$XB_t[i] = a0_b + a1_b*X_t$XB_t[i-1] + X_t$U_t[i-1]
}

ggplot(data = X_t, aes(x = Tiempo)) +
  geom_line(aes(y = XR_t, color = "a1 = 0.9"), linewidth = 0.5) +
  geom_line(aes(y = XB_t, color = "a1 = 0.5"), linewidth = 0.5) +
  geom_hline(yintercept = mu_teo, linetype = "dashed",
             color = "gray30") +
  scale_color_manual(values = c("a1 = 0.9" = "darkred",
                                "a1 = 0.5" = "darkblue")) +
  xlab("Tiempo") + 
  ylab(TeX("$X_t$")) + 
  theme(legend.position = "bottom", legend.title = element_blank()) +
  theme(plot.title = element_text(size = 11, face = "bold", 
                                  hjust = 0)) + 
  theme(plot.subtitle = element_text(size = 10, hjust = 0)) + 
  theme(plot.caption = element_text(size = 10, hjust = 0)) +
  theme(plot.margin = unit(c(1,1,1,1), "cm")) +
  labs(
    title = "Efecto de la persistencia en un AR(1)",
    subtitle = "Mismos choques y misma media teórica; sólo cambia a1",
    caption = "Fuente: Elaboración propia."
  )

$Persistencia de los choques en un $AR(1)$: dos procesos construidos con los mismos choques y la misma media teórica ($\mu = 50$), uno con $a_1 = 0.9$ y otro con $a_1 = 0.5$$

Figura 3.8: Persistencia de los choques en un $AR(1)$: dos procesos construidos con los mismos choques y la misma media teórica ($\mu = 50$), uno con $a_1 = 0.9$ y otro con $a_1 = 0.5$

Ejemplo. Para el segundo ejemplo consideremos una aplicación a una serie de tiempo en específico: Pasajeros transportados mensualmente en el Sistema de Transporte Colectivo Metro (pasajeros medidos en millones).³

A la serie se le aplicará una metodología de estimación dada por el método de Máxima Verosimilitud (ML, por sus siglas en inglés). Antes de realizar el proceso de estimación, consideremos una transformación de diferencias logarítmicas con el objeto de obtener una serie de tiempo expresada en tasas de crecimiento⁴ y con un comportamiento parecido a un proceso estacionario.

Así, para cada una de las series que analicemos en diferencias logarítmicas respecto del momento $k$ las expresaremos bajo la siguiente transformación:

\[\begin{equation*} DLX_t = log(X_t) - log(X_{t-k}) \end{equation*}\]

Donde $k = 1, 2, 3, \ldots$ y $log(.)$ es la función logaritmo natural. Esta expresión se puede interpretar como una tasa de crecimiento, puesto que asumimos variaciones pequeñas para las cuales se cumple que: $log(X_t) - log(X_{t-k}) \approx \frac{X_t - X_{t-k}}{X_{t-k}}$.

Primero, para realizar el análisis de una serie de tiempo deberemos decidir si éste se realizará para la serie en niveles o en diferencias. Por convención, decimos que la serie está en niveles si ésta se analiza sin hacerle ninguna transformación o si se analiza aplicando solo logaritmos. Cuando la serie se analiza en diferencias significa que la diferencia se hace sin aplicar logaritmos o aplicando logaritmos. Sin embargo, la convención es hacer un análisis en diferencias logarítmicas.

Para decidir cómo analizar la serie de pasajeros en el metro de la CDMX en la Figura 3.9 se muestra la gráfica de la serie en niveles (sin transformación logarítmica y con transformación logarítmica) y en diferencias logarítmicas mensuales (es decir, con $k = 1$).

library(ggplot2)
library(dplyr)
library(readxl)
library(stats)

Datos <- read_excel("BD/Base_Transporte.xlsx", 
                    sheet = "Datos", col_names = TRUE)

# En Niveles
Pax_Metro <- ts(Datos$Pax_Metro, start = c(2000, 1), 
                freq = 12)

# En Logaritmos:
Pax_LMetro <- ts(log(Datos$Pax_Metro), start = c(2000, 1), 
                freq = 12)

# Diferencias mensuales:
Pax_DLMetro <- ts( log(Datos$Pax_Metro) - 
                   dplyr::lag( log(Datos$Pax_Metro), 1 ),
                 start = c(2000, 1), freq = 12)

#
par(mfrow = c(3,1))

plot(Pax_Metro, xlab = "Tiempo", 
     main = "Pasajeros transportados (Millones) en el Metro CDMX",
     col = "darkgreen")

plot(Pax_LMetro, xlab = "Tiempo", 
     main = "LN Pasajeros transportados (Millones) en el Metro CDMX",
     col = "darkblue")

plot(Pax_DLMetro, xlab = "Tiempo", 
     main = "Diff LN Pasajeros transportados (Millones) en el Metro CDMX", 
     col = "darkred")

Figura 3.9: Pasajeros transportados (Millones) en el Metro de la CDMX en niveles y en diferencias logarítmicas

par(mfrow=c(1,1))

A continuación, estimaremos una $AR(1)$ para la serie en niveles bajo la transformación logarítmica ($PaxLMetro_t$) y en diferencias logarítmicas ($PaxDLMetro_t$). Los resultados se muestran en los Cuadros 3.2 y 3.3.

source("arroots.R")
source("plot.armaroots.R")
AR_LMetro  <- arima(Pax_LMetro,  order = c(1, 0, 0), method = "ML")
AR_DLMetro <- arima(Pax_DLMetro, order = c(1, 0, 0), method = "ML")

arima_kable <- function(mod, cap) {
  cf  <- mod$coef
  se  <- sqrt(diag(mod$var.coef))
  tv  <- round(cf / se, 3)
  pv  <- 2 * pnorm(-abs(tv))
  sig <- ifelse(pv < 0.001, "***", ifelse(pv < 0.01, "**", ifelse(pv < 0.05, "*", ifelse(pv < 0.1, ".", ""))))
  df  <- data.frame(Estimado = round(cf, 4), SE = round(se, 4), t = tv, Signif = sig)
  # Con escape = FALSE, los guiones bajos de nombres de coeficientes
  # (p. ej. dummies como D_Sep2017_MA) deben escaparse solo para LaTeX
  if (knitr::is_latex_output()) rownames(df) <- gsub("_", "\\\\_", rownames(df))
  knitr::kable(df, col.names = c("Estimado", "Error Est.", "Estad. $t$", "Signif."),
               caption = paste0(cap, " ($\\hat{\\sigma}^2 = ", round(mod$sigma2, 6),
                                "$; AIC $=$ ", round(AIC(mod), 2), "$)$."),
               align = c("r","r","r","c"), booktabs = TRUE, escape = FALSE)
}

arima_kable(AR_LMetro, "AR(1) para la variable $PaxLMetro_t$")

Cuadro 3.2: AR(1) para la variable $PaxLMetro_t$ ($\hat{\sigma}^2 = 0.009442$; AIC $=$ -570.86$)$.
	Estimado	Error Est.	Estad. $t$	Signif.
ar1	0.8907	0.0249	35.730	***
intercept	4.7280	0.0487	97.061	***

arima_kable(AR_DLMetro, "AR(1) para la variable $PaxDLMetro_t$")

Cuadro 3.3: AR(1) para la variable $PaxDLMetro_t$ ($\hat{\sigma}^2 = 0.009591$; AIC $=$ -565.61$)$.
	Estimado	Error Est.	Estad. $t$	Signif.
ar1	-0.2025	0.0550	-3.681	***
intercept	-0.0003	0.0046	-0.057

En ambos casos observamos que el parámetro asociado al componente AR es significativo y cumple con la restricción de ser en valor absoluto menor a 1, por lo que la solución asociada al proceso será convergente. También en ambos casos se reporta la estadística o Criterio de Información de Akaike (AIC, por sus siglas en inglés), cuya importancia y aplicación discutiremos más adelante.

3.4.2 AR(2)

Una vez analizado el caso de $AR(1)$ analizaremos el caso del $AR(2)$. La ecuación generalizada del proceso autorregresivo de orden 2 (denotado como $AR(2)$) puede ser escrita como:

\[\begin{equation} X_t = a_0 + a_1 X_{t-1} + a_2 X_{t-2} + U_t \tag{3.35} \end{equation}\]

Donde $U_t$ denota un proceso puramente aleatorio con media cero ($0$), varianza constante ($\sigma^2$) y autocovarianza cero ($Cov(U_t, U_s) = 0$, con $t \neq s$), y un parámetro $a_2 \neq 0$. Así, utilizando el operador rezago podemos reescribir la ecuación (3.35) como:

\[\begin{eqnarray*} X_t - a_1 X_{t-1} - a_2 X_{t-2} & = & a_0 + U_t \\ (1 - a_1 L^1 - a_2 L^2) X_t & = & a_0 + U_t \end{eqnarray*}\]

Donde, vamos a denotar a $\alpha (L) = (1 - a_1 L^1 - a_2 L^2)$, y lo denotaremos como un polinomio que depende del operador rezago y que es distinto de cero. De esta forma podemos reescribir la ecuación (3.35) como:

\[\begin{equation} \alpha(L) X_t = a_0 + U_t \end{equation}\]

Ahora, supongamos que existe el inverso multiplicativo del polinomio $\alpha(L)$, el cual será denotado como: $\alpha^{-1}(L)$ y cumple con que:

\[\begin{equation} \alpha^{-1}(L) \alpha(L) = 1 \end{equation}\]

Así, podemos escribir la solución a la ecuación (3.35) como: \[\begin{equation*} X_t = \alpha^{-1}(L) \delta + \alpha^{-1}(L) U_t \end{equation*}\]

Si utilizamos el hecho que $\alpha^{-1}(L)$ se puede descomponer a través del procedimiento de Wold en un polinomio de forma similar al caso de $AR(1)$, tenemos que:

\[\begin{equation} \alpha^{-1}(L) = \psi_0 + \psi_1 L + \psi_2 L^2 + \ldots \end{equation}\]

Por lo tanto, el inverso multiplicativo $\alpha^{-1}(L)$ se puede ver como:

\[\begin{equation} 1 = (1 - a_1 L^1 - a_2 L^2) (\psi_0 + \psi_1 L + \psi_2 L^2 + \ldots) \tag{3.36} \end{equation}\]

Desarrollando la ecuación (3.36) tenemos la siguiente expresión:

\[\begin{eqnarray*} 1 & = & \psi_0 + \psi_1 L + \psi_2 L^2 + \psi_3 L^3 + \ldots \\ & & - a_1 \psi_0 L - a_1 \psi_1 L^2 - a_1 \psi_2 L^3 - \ldots \\ & & - a_2 \psi_0 L^2 - a_2 \psi_1 L^3 - \ldots \end{eqnarray*}\]

Ahora, podemos agrupar todos los términos en función del exponente asociado al operador rezago $L$. La siguiente es una solución particular y es una de las múltiples que podrían existir que cumpla con la ecuación (3.36). Sin embargo, para efectos del análisis, sólo necesitamos una de esas soluciones. Utilizaremos las siguientes condiciones que deben cumplirse en una de las posibles soluciones: \[\begin{eqnarray*} L^0 & : & \Rightarrow \psi_0 = 1 \\ L^1 & : & \psi_1 - a_1 \psi_0 = 0 \Rightarrow \psi_1 = a_1 \\ L^2 & : & \psi_2 - a_1 \psi_1 - a_2 \psi_0 = 0 \Rightarrow \psi_2 = a^2_1 + a_2 \\ L^3 & : & \psi_3 - a_1 \psi_2 - a_2 \psi_1 = 0 \Rightarrow \psi_3 = a^3_1 + 2 a_1 a_2 \\ & \vdots & \end{eqnarray*}\]

De esta forma podemos observar que en el límite siempre obtendremos una ecuación del tipo $\psi_j - a_1 \psi_{j-1} - a_2 \psi_{j-2} = 0$ asociada a cada uno de los casos en que exista un $L^j$, donde $j \neq 0, 1$, y la cual siempre podremos resolver conociendo que las condiciones iniciales son: $\psi_0 = 1$ y $\psi_1 = a_1$.

Así, de las relaciones antes mencionadas y considerando que $\alpha^{-1} (L)$ aplicada a una constante como $a_0$ tendrá como resultado otra constante, podemos escribir la solución del proceso AR(2) de la ecuación (3.35) mediante una expresión como sigue: \[\begin{equation} X_t = \frac{a_0}{1 - a_1 - a_2} + \sum^{\infty}_{j = 0} \psi_j U_{t - j} \tag{3.37} \end{equation}\]

Donde todos los parámetros $\psi_i$ están determinados por los parámetros $a_0$, $a_1$ y $a_2$. En particular, $\psi_0 = 1$ y $\psi_1 = a_1$, como describimos anteriormente. Al igual que en el caso del $AR(1)$, en la ecuación (3.37) las condiciones de estabilidad estarán dadas por las soluciones del siguiente polinomio característico:⁵ \[\begin{equation} \lambda^2 - \lambda a_1 - a_2 = 0 \end{equation}\]

Así, la condición de estabilidad de la trayectoria es que $|\lambda_i| < 1$, para $i = 1, 2$. Es decir, es necesario que cada una de las raíces sea, en valor absoluto, siempre menor que la unidad. Estas son las condiciones de estabilidad para el proceso $AR(2)$.

Finalmente, al igual que en un $AR(1)$, a continuación determinamos los momentos de una serie que sigue un proceso $AR(2)$. Iniciamos con la determinación de la media de la serie:

\[\begin{equation} \mathbb{E}[X_t] = \mu = \frac{a_0}{1 - a_1 - a_2} \end{equation}\]

Lo anterior es cierto puesto que $\mathbb{E}[U_{t - i}] = 0$, para todo $i = 0, 1, 2, \ldots$. Para determinar la varianza utilizaremos las siguientes relaciones basadas en el uso del valor esperado, varianza y covarianza de la serie. Adicionalmente, para simplificar el trabajo asumamos que $a_0 = 0$, lo cual implica que $\mu = 0$. Dicho lo anterior, partamos de:

\[\begin{eqnarray*} \mathbb{E}[X_t X_{t - \tau}] & = & \mathbb{E}[(a_1 X_{t-1} + a_2 X_{t-2} + U_t) X_{t - \tau}]\\ & = & a_1 \mathbb{E}[X_{t - 1} X_{t - \tau}] + a_2 \mathbb{E}[X_{t - 2} X_{t - \tau}] + \mathbb{E}[U_{t} X_{t - \tau}] \end{eqnarray*}\]

Donde $\tau = 0, 1, 2, 3, \ldots$ y $\mathbb{E}[U_{t} X_{t - \tau}] = 0$ para todo $\tau \neq 0$.⁶ Dicho esto, podemos derivar el valor del valor esperado para diferentes valores de $\tau$:

\[\begin{eqnarray*} \tau = 0 & : & \gamma(0) = a_1 \gamma(1) + a_2 \gamma(2) + \sigma^2 \\ \tau = 1 & : & \gamma(1) = a_1 \gamma(0) + a_2 \gamma(1) \\ \tau = 2 & : & \gamma(2) = a_1 \gamma(1) + a_2 \gamma(0) \\ & \vdots & \end{eqnarray*}\]

Donde debe ser claro que $\mathbb{E}[(X_{t} - \mu)(X_{t - \tau} - \mu)] = \mathbb{E}[X_{t} X_{t - \tau}] = \gamma(\tau)$. Así, en general cuando $\tau \neq 0$:

\[\begin{equation} \gamma(\tau) = a_1 \gamma(\tau - 1) + a_2 \gamma(\tau - 2) \end{equation}\]

Realizando la sustitución recursiva y solucionando el sistema respectivo obtenemos que la varianza y las covarianzas estarán determinadas por: \[\begin{equation} Var[X_t] = \gamma(0) = \frac{1 - a_2}{(1 + a_2)[(1 - a_2)^2 - a^2_1]} \sigma^2 \end{equation}\]

\[\begin{equation} \gamma(1) = \frac{a_1}{(1 + a_2)[(1 - a_2)^2 - a^2_1]} \sigma^2 \end{equation}\] \[\begin{equation} \gamma(2) = \frac{a^2_1 + a_2 - a^2_2}{(1 + a_2)[(1 - a_2)^2 - a^2_1]} \sigma^2 \end{equation}\]

Recordemos que las funciones de autocorrelación se obtienen de la división de cada una de las funciones de covarianza ($\gamma(\tau)$) por la varianza ($\gamma(0)$). Así, podemos construir la siguiente expresión:

\[\begin{equation} \rho(\tau) - a_1 \rho(\tau - 1) - a_2 \rho(\tau - 2) = 0 \end{equation}\]

Ejemplo. Utilizaremos la serie de Pasajeros en vuelos nacionales (en vuelos de salidas) para estimar un $AR(2)$ mediante el método de máxima verosimilitud (ML, por sus siglas en inglés). Antes de realizar el proceso de estimación, consideremos una transformación de la serie en logaritmos y una más en diferencias logarítmicas; lo anterior con el objeto de obtener una serie de tiempo suavizada y expresada en tasas de crecimiento, con un comportamiento parecido a un proceso estacionario.

Así, para cada una de las series que analicemos en diferencias logarítmicas, las expresaremos bajo la siguiente transformación: \[\begin{equation*} DLX_t = log(X_t) - log(X_{t-k}) \end{equation*}\]

Donde $k = 1, 2, 3, \ldots$ y $log(.)$ es la función logaritmo natural. Por convención, decimos que la serie está en niveles si esta se analiza sin hacerle ninguna transformación o se analiza en logaritmos. Cuando la serie se analiza en diferencias significa que la diferencia se hace sin aplicar logaritmos. Y cuando la serie analizada está en diferencias logarítmicas también diremos que está en diferencias. Sin embargo, lo común es hacer un análisis en logaritmos y en diferencias logarítmicas.

Primero, para decidir si se realizará un AR(2) para la serie en niveles o en diferencias, analizaremos su gráfica. La serie en niveles, en niveles bajo una transformación logarítmica y en diferencias logarítmicas mensuales de los pasajeros en vuelos nacionales se muestra en la Figura 3.10.

library(ggplot2)
library(dplyr)
library(readxl)
library(stats)

Datos <- read_excel("BD/Base_Transporte.xlsx", 
                    sheet = "Datos", col_names = TRUE)

# En Niveles
Pax_Nal <- ts(Datos$Pax_Nal, 
              start = c(2000, 1),
              freq = 12)

# Logaritmos:
LPax_Nal <- ts(log(Datos$Pax_Nal), 
               start = c(2000, 1), 
               freq = 12)

# Diferencias mensuales:
DLPax_Nal <- ts(log(Datos$Pax_Nal) - 
                dplyr::lag(log(Datos$Pax_Nal), 1),
                start = c(2000, 1), freq = 12)

#
par(mfrow=c(3,1))

plot(Pax_Nal, xlab = "Tiempo", ylab = "Pasajeros",
     main = "Pasajeros en vuelos nacionales de salida",
     col = "darkgreen")

plot(LPax_Nal, xlab = "Tiempo", ylab = "LN Pasajeros",
     main = "LN Pasajeros en vuelos nacionales de salida",
     col = "darkblue")

plot(DLPax_Nal, xlab = "Tiempo", ylab = "DLN Pasajeros",
     main = "Diff LN Pasajeros en vuelos nacionales de salida", 
     col = "darkred")

Figura 3.10: Pasajeros en vuelos de salidas nacionales en niveles y en diferencias logarítmicas

par(mfrow=c(1,1))

A continuación, estimaremos un $AR(2)$ para la serie en niveles bajo la transformación logarítmica ($LPaxNal_t$) y en diferencias logarítmicas ($DLPaxNal_t$). Los resultados se muestran en los Cuadros 3.4 y 3.5.

AR_LPax_Nal  <- arima(LPax_Nal,  order = c(2, 0, 0), method = "ML")
AR_DLPax_Nal <- arima(DLPax_Nal, order = c(2, 0, 0), method = "ML")
arima_kable(AR_LPax_Nal, "AR(2) para la variable $LPaxNal_t$")

Cuadro 3.4: AR(2) para la variable $LPaxNal_t$ ($\hat{\sigma}^2 = 0.027178$; AIC $=$ -233.21$)$.
	Estimado	Error Est.	Estad. $t$	Signif.
ar1	1.0358	0.0557	18.589	***
ar2	-0.1094	0.0560	-1.953	.
intercept	14.7970	0.1215	121.775	***

arima_kable(AR_DLPax_Nal, "AR(2) para la variable $DLPaxNal_t$")

Cuadro 3.5: AR(2) para la variable $DLPaxNal_t$ ($\hat{\sigma}^2 = 0.025978$; AIC $=$ -248.62$)$.
	Estimado	Error Est.	Estad. $t$	Signif.
ar1	0.0911	0.0539	1.689	.
ar2	-0.2819	0.0539	-5.234	***
intercept	0.0039	0.0076	0.517

Para ambos casos se reportan los errores estándar y el criterio de Akaike (AIC). Finalmente, en la Figura 3.11 mostramos las raíces del polinomio característico para verificar la convergencia de las soluciones. Las funciones arroots(), maroots() y plot.armaroots() que usamos para graficarlas están tomadas del código publicado por Rob J. Hyndman en su blog Hyndsight, que es también la base de la función plot.Arima() del paquete forecast.⁷

par(mfrow=c(1,2))

plot.armaroots(arroots(AR_LPax_Nal), 
               main="Inverse AR roots of \nAR(2): LN Pax Nal")

#
plot.armaroots(arroots(AR_DLPax_Nal), 
               main="Inverse AR roots of \nAR(2): Diff LN Pax Nal")

Figura 3.11: Inverso de las Raíces del polinomio característico

par(mfrow=c(1,1))

3.4.3 AR(p)

Veremos ahora una generalización de los procesos autorregresivos (AR). Esta generalización es conocida como un proceso $AR(p)$ y que puede ser descrito por la siguiente ecuación en diferencia estocástica: \[\begin{equation} X_t = a_0 + a_1 X_{t-1} + a_2 X_{t-2} + a_3 X_{t-3} + \ldots + a_p X_{t-p} + U_t \tag{3.38} \end{equation}\]

Donde $a_p \neq 0$, y $U_t$ es un proceso puramente aleatorio con media cero (0), varianza constante ($\sigma^2$) y covarianza cero (0). Usando el operador rezago, $L^k$, para $k = 0, 1, 2, \ldots, p$, obtenemos la siguiente expresión de la ecuación (3.38):

\[\begin{equation} (1 - a_1 L - a_2 L^2 - a_3 L^3 - \ldots - a_p L^p) X_t = a_0 + U_t \end{equation}\]

Definamos el polinomio $\alpha(L)$ como:

\[\begin{equation} \alpha(L) = 1 - a_1 L - a_2 L^2 - a_3 L^3 - \ldots - a_p L^p \tag{3.39} \end{equation}\]

De forma similar que en los procesos $AR(1)$ y $AR(2)$, las condiciones de estabilidad del proceso $AR(p)$ estarán dadas por la solución de la ecuación característica:

\[\begin{equation} \lambda^p - a_1 \lambda^{p-1} - a_2 \lambda^{p-2} - a_3 \lambda^{p-3} - \ldots - a_p = 0 \end{equation}\]

Así, solo si el polinomio anterior tiene raíces cuyo valor absoluto sea menor a uno ($|\lambda_i| < 1$) podremos decir que el proceso es convergente y estable. Lo anterior significa que la ecuación (3.39) puede expresarse en términos de la descomposición de Wold o como la suma infinita de términos como: \[\begin{equation} \frac{1}{1 - a_1 L - a_2 L^2 - a_3 L^3 - \ldots - a_p L^p} = \psi_0 + \psi_1 L + \psi_2 L^2 + \psi_3 L^3 + \ldots \end{equation}\]

Donde, por construcción de $\alpha(L) \alpha^{-1}(L) = 1$ implica que $\psi_0 = 1$. De forma similar a los procesos AR(1) y AR(2), es posible determinar el valor de los coeficientes $\psi_j$ en términos de los coeficientes $a_i$. Así, la solución del proceso $AR(p)$ estará dada por:

\[\begin{equation} X_t = \frac{a_0}{1 - a_1 - a_2 - a_3 - \ldots - a_p} + \sum^{\infty}_{j = 0} \psi_j U_{t-j} \tag{3.40} \end{equation}\]

Considerando la solución de la ecuación (3.38) expresada en la ecuación (3.40) podemos determinar los momentos del proceso y que estarán dados por una media como:

\[\begin{equation} \mathbb{E}[X_t] = \mu = \frac{a_0}{1 - a_1 - a_2 - a_3 - \ldots - a_p} \end{equation}\]

Lo anterior, considerado que $\mathbb{E}[U_t] = 0$, para todo $t$. Para determinar la varianza del proceso, sin pérdida de generalidad, podemos definir una ecuación: $\gamma(\tau) = \mathbb{E}[X_{t - \tau} X_t]$, la cual (omitiendo la constante, ya que la correlación de una constante con cualquier variable aleatoria que depende del tiempo es cero (0)) puede ser escrita como:

\[\begin{equation} \gamma(\tau) = \mathbb{E}[(X_{t - \tau}) \cdot (a_1 X_{t-1} + a_2 X_{t-2} + a_3 X_{t-3} + \ldots + a_p X_{t-p} + U_t)] \end{equation}\]

Donde $\tau = 0, 1, 2, \ldots, p$ y $a_0 = 0$, lo que implica que $\mu = 0$. De lo anterior obtenemos el siguiente conjunto de ecuaciones mediante sustituciones de los valores de $\tau$:

\[\begin{eqnarray} \gamma(0) & = & a_1 \gamma(1) + a_2 \gamma(2) + \ldots + a_p \gamma(p) + \sigma^2 \nonumber \\ \gamma(1) & = & a_1 \gamma(0) + a_2 \gamma(1) + \ldots + a_p \gamma(p-1) \nonumber \\ \vdots \nonumber \\ \gamma(p) & = & a_1 \gamma(p-1) + a_2 \gamma(p-2) + \ldots + a_p \gamma(0) \nonumber \end{eqnarray}\]

De esta forma, es fácil observar que la ecuación general para $p > 0$ estará dada por:

\[\begin{equation} \gamma(\tau) - a_1 \gamma(\tau - 1) - a_2 \gamma(\tau - 2) - \ldots - a_p \gamma(\tau - p) = 0 \tag{3.41} \end{equation}\]

Dividiendo la ecuación (3.41) por $\gamma(0)$, se obtiene la siguiente ecuación:

\[\begin{equation} \rho(\tau) - a_1 \rho(\tau - 1) - a_2 \rho(\tau - 2) - \ldots - a_p \rho(\tau - p) = 0 \end{equation}\]

Así, podemos escribir el siguiente sistema de ecuaciones: \[\begin{eqnarray} \rho(1) & = & a_1 + a_2 \rho(1) + a_3 \rho(2) + \ldots + a_p \rho(p-1) \nonumber \\ \rho(2) & = & a_1 \rho(1) + a_2 + a_3 \rho(1) + \ldots + a_p \rho(p-2) \nonumber \\ & \vdots & \nonumber \\ \rho(p) & = & a_1 \rho(p-1) + a_2 \rho(p-2) + \ldots + a_p \nonumber \end{eqnarray}\]

Lo anterior se puede expresar como un conjunto de vectores y matrices de la siguiente forma:

\[\begin{equation} \left[ \begin{array}{c} \rho(1) \\ \rho(2) \\ \vdots \\ \rho(p) \end{array} \right] = \left[ \begin{array}{c c c c} 1 & \rho(1) & \ldots & \rho(p - 1) \\ \rho(1) & 1 & \ldots & \rho(p - 2) \\ \rho(2) & \rho(1) & \ldots & \rho(p - 3) \\ \vdots & \vdots & \ldots & \vdots \\ \rho(p - 1) & \rho(p - 2) & \ldots & 1 \\ \end{array} \right] \left[ \begin{array}{c} a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_p \\ \end{array} \right] \end{equation}\]

De lo anterior podemos escribir la siguiente ecuación que es una forma alternativa para expresar los valores de los coeficientes $a_i$ de la solución del proceso $AR(p)$:

\[\begin{equation} \boldsymbol{\rho} = \mathbf{R} \mathbf{a} \end{equation}\]

Es decir, podemos obtener la siguiente expresión:

\[\begin{equation} \mathbf{a} = \mathbf{R}^{-1} \boldsymbol{\rho} \end{equation}\]

Ejemplo. Utilizaremos la serie de Pasajeros en vuelos internacionales de salida para estimar un $AR(p)$ mediante el método de máxima verosimilitud (ML). Antes de realizar el proceso de estimación, consideremos una transformación de la serie en logaritmos y una más en diferencias logarítmicas; lo anterior con el objeto de obtener una serie de tiempo suavizada y expresada en tasas de crecimiento, con un comportamiento parecido a un proceso estacionario.

Primero, para decidir si se realizará un $AR(p)$ para la serie en niveles o en diferencias, analizaremos su gráfica. La serie de Pasajeros en vuelos internacionales de salidas se muestra en la Figura 3.12. En esta se muestra la gráfica de la serie en niveles (sin transformación logarítmica y con transformación logarítmica) y en diferencias logarítmicas mensuales (es decir, con diferencia respecto del mes inmediato anterior).

library(ggplot2)
library(dplyr)
library(readxl)
library(stats)

Datos <- read_excel("BD/Base_Transporte.xlsx", 
                    sheet = "Datos", col_names = TRUE)

# En Niveles
Pax_Int <- ts(Datos$Pax_Int, 
              start = c(2000, 1), 
              freq = 12)

# Logaritmos:
LPax_Int <- ts(log(Datos$Pax_Int), 
               start = c(2000, 1), 
               freq = 12)

# Diferencias mensuales:
DLPax_Int <- ts(log(Datos$Pax_Int) - dplyr::lag(log(Datos$Pax_Int), 1),
                start = c(2000, 1), 
                freq = 12)

#
par(mfrow=c(3,1))

plot(Pax_Int, xlab = "Tiempo", ylab = "Pasajeros",
     main = "Pasajeros en vuelos internacionales de salida",
     col = "darkgreen")

plot(LPax_Int, xlab = "Tiempo", ylab = "LN Pasajeros",
     main = "LN Pasajeros en vuelos internacionales de salida",
     col = "darkblue")

plot(DLPax_Int, xlab = "Tiempo", ylab = "DLN Pasajeros",
     main = "Diff LN Pasajeros en vuelos internacionales de salida", 
     col = "darkred")

Figura 3.12: Pasajeros en vuelos internacionales de salida en niveles y en diferencias logarítmicas

par(mfrow=c(1,1))

De la gráfica en la Figura 3.12 observamos que quizá la mejor forma de estimar un $AR(p)$ es mediante la serie en diferencias, ya que ésta es la que parece ser una serie estacionaria. A continuación, estimaremos una AR(4) para la serie en diferencias logarítmicas ($DLPaxInt_t$):

source("arroots.R")
source("plot.armaroots.R")
AR_LPax_Int  <- arima(LPax_Int,  order = c(4, 0, 0), method = "ML")
AR_DLPax_Int <- arima(DLPax_Int, order = c(4, 0, 0), method = "ML")

arima_kable(AR_DLPax_Int, "AR(4) para la variable $DLPaxInt_t$")

Cuadro 3.6: AR(4) para la variable $DLPaxInt_t$ ($\hat{\sigma}^2 = 0.059659$; AIC $=$ 18.2$)$.
	Estimado	Error Est.	Estad. $t$	Signif.
ar1	0.0767	0.0562	1.363
ar2	-0.2926	0.0557	-5.250	***
ar3	-0.1409	0.0557	-2.531
ar4	-0.0236	0.0560	-0.421
intercept	0.0040	0.0100	0.397

Los resultados se muestran en el Cuadro 3.6. Finalmente, en la Figura 3.13 mostramos las raíces del polinomio característico; la inspección visual confirma que el AR(4) tiene solución convergente y estable.

par(mfrow=c(1,2))

plot.armaroots(arroots(AR_LPax_Int),
               main="Inverse AR roots of \nAR(4): LN Pax Int")

#
plot.armaroots(arroots(AR_DLPax_Int),
               main="Inverse AR roots of \nAR(4): Diff LN Pax Int")

Figura 3.13: Inverso de las Raíces del polinomio característico

par(mfrow=c(1,1))

3.5 Procesos de Medias Móviles (MA)

3.5.1 MA(1)

Una vez planteado el proceso generalizado de $AR(p)$, iniciamos el planteamiento de los procesos de medias móviles, denotados como $MA(q)$. Iniciemos con el planteamiento del proceso $MA(1)$, que se puede escribir como una ecuación como la siguiente:

\[\begin{equation} X_t = \mu + U_t - b_1 U_{t-1} \tag{3.42} \end{equation}\]

O como:

\[\begin{equation} X_t - \mu = (1 - b_1 L) U_{t} \end{equation}\]

Donde $U_t$ es un proceso puramente aleatorio, es decir, con $\mathbb{E}[U_t] = 0$, $Var[U_t] = \sigma^2$, y $Cov[U_t, U_s] = 0$. Así, un proceso $MA(1)$ puede verse como un proceso AR con una descomposición de Wold en la que $\psi_0 = 1$, $\psi_1 = - b_1$ y $\psi_j = 0$ para todo $j > 1$.

Al igual que los procesos autorregresivos, determinaremos los momentos de un proceso $MA(1)$. En el caso de la media, observamos que será: \[\begin{eqnarray} \mathbb{E}[X_t] & = & \mu + \mathbb{E}[U_t] - b_1 \mathbb{E}[U_{t - 1}] \nonumber \\ & = & \mu \end{eqnarray}\]

Por su parte, la varianza estará dada por:

\[\begin{eqnarray} Var[X_t] & = & \mathbb{E}[(X_t - \mu)^2] \nonumber \\ & = & \mathbb{E}[(U_t - b_1 U_{t-1})^2] \nonumber \\ & = & \mathbb{E}[U_t^2 - 2 b_1 U_t U_{t-1} + b_1^2 U_{t - 1}^2] \nonumber \\ & = &\mathbb{E}[U_t^2] - 2 b_1 \mathbb{E}[U_t U_{t-1}] + b_1^2 \mathbb{E}[U_{t - 1}^2] \nonumber \\ & = & \sigma^2 + b_1^2 \sigma^2 \nonumber \\ & = & (1 + b_1^2) \sigma^2 = \gamma(0) \end{eqnarray}\]

De esta forma, la varianza del proceso es constante en cualquier periodo $t$. Para determinar la covarianza utilizaremos la siguiente ecuación: \[\begin{eqnarray} \mathbb{E}[(X_t - \mu)(X_{t + \tau} - \mu)] & = & \mathbb{E}[(U_t - b_1 U_{t-1})(U_{t + \tau} - b_1 U_{t + \tau - 1})] \nonumber \\ & = & \mathbb{E}[U_t U_{t + \tau} - b_1 U_t U_{t + \tau - 1} - b_1 U_{t - 1} U_{t + \tau} \nonumber \\ & & + b_1^2 U_{t - 1} U_{t + \tau - 1}] \nonumber \\ & = & \mathbb{E}[U_t U_{t + \tau}] - b_1 \mathbb{E}[U_t U_{t + \tau - 1}] \nonumber \\ & & - b_1 \mathbb{E}[U_{t - 1} U_{t + \tau}] + b_1^2 \mathbb{E}[U_{t - 1} U_{t + \tau - 1}] \tag{3.43} \end{eqnarray}\]

Si hacemos sustituciones de diferentes valores de $\tau$ en la ecuación (3.43) notaremos que la covarianza será distinta de cero únicamente para el caso de $\tau = 1, -1$. En ambos casos tendremos como resultado:

\[\begin{eqnarray} \mathbb{E}[(X_t - \mu)(X_{t + 1} - \mu)] & = & \mathbb{E}[(X_t - \mu)(X_{t - 1} - \mu)] \nonumber \\ & = & - b_1 \mathbb{E}[U_t U_{t}] \nonumber \\ & = & - b_1 \mathbb{E}[U_{t - 1} U_{t - 1}] \nonumber \\ & = & - b_1 \sigma^2 = \gamma(1) \end{eqnarray}\]

De esta forma tendremos que las funciones de autocorrelación estarán dadas por los siguientes casos:

\[\begin{eqnarray} \rho(0) & = & 1 \nonumber \\ \rho(1) & = & \frac{- b_1}{1 + b_1^2} \nonumber \\ \rho(\tau) & = & 0 \text{ para todo } \tau > 1 \nonumber \end{eqnarray}\]

Ahora regresando a la ecuación (3.42), su solución la podemos expresar como:

\[\begin{eqnarray} U_ t & = & - \frac{\mu}{1 - b_1} + \frac{1}{1 - b_1 L} X_t \nonumber \\ & = & - \frac{\mu}{1 - b_1} + X_t + b_1 X_{t-1} + b_1^2 X_{t-2} + \ldots \nonumber \end{eqnarray}\]

Donde la condición para que se cumpla esta ecuación es que $|b_1| < 1$. La manera de interpretar esta condición es como una condición de estabilidad de la solución y como una condición de invertibilidad. Notemos que un $MA(1)$ (y en general un $MA(q)$) es equivalente a un $AR(\infty)$, es decir, cuando se invierte un MA se genera un AR con infinitos rezagos.

En esta sección no desarrollaremos un ejemplo, primero explicaremos en qué consiste una modelación del tipo $MA(q)$ y después plantearemos un ejemplo en concreto.

3.5.2 MA(q)

En general, el proceso de medias móviles de orden $q$, $MA(q)$, puede ser escrito como:

\[\begin{equation} X_t = \mu + U_t - b_1 U_{t-1} - b_2 U_{t-2} - \ldots - b_q U_{t-q} \tag{3.44} \end{equation}\]

Podemos reescribir la ecuación (3.44) utilizando el operador rezago. Así, tendremos el proceso de $MA(q)$ como:

\[\begin{eqnarray} X_t - \mu & = & (1 - b_1 L - b_2 L^2 - \ldots - b_q L^q) U_{t} \nonumber \\ X_t - \mu & = & \beta(L) U_t \tag{3.45} \end{eqnarray}\]

Donde $U_t$ es un proceso puramente aleatorio con $\mathbb{E}[U_t] = 0$, $Var[U_t] = \mathbb{E}[U_t^2] = \sigma^2$ y $Cov[U_t, U_s] = \mathbb{E}[U_t U_s] = 0$, y $\beta(L) = 1 - b_1 L - b_2 L^2 - \ldots - b_q L^q$ es un polinomio del operador rezago $L$. La ecuación (3.45) puede ser interpretada como un polinomio de rezagos de orden $q$ aplicado a la serie $U_t$.

Ahora determinemos los momentos de un proceso $MA(q)$:

\[\begin{eqnarray} \mathbb{E}[X_t] & = & \mathbb{E}[\mu + U_t - b_1 U_{t-1} - b_2 U_{t-2} - \ldots - b_q U_{t-q}] \nonumber \\ & = & \mu + \mathbb{E}[U_t] - b_1 \mathbb{E}[U_{t-1}] - b_2 \mathbb{E}[U_{t-2}] - \ldots - b_q \mathbb{E}[U_{t-q}] \nonumber \\ & = & \mu \end{eqnarray}\]

En el caso de la varianza tenemos que se puede expresar como: \[\begin{eqnarray} Var[X_t] & = & \mathbb{E}[(X_t - \mu)^2] \nonumber \\ & = & \mathbb{E}[(U_t - b_1 U_{t-1} - b_2 U_{t-2} - \ldots - b_q U_{t-q})^2] \nonumber \\ & = & \mathbb{E}[U_t^2 + b_1^2 U_{t-1}^2 + b_2^2 U_{t-2}^2 + \ldots + b_q^2 U_{t-q}^2 \nonumber \\ & & - 2 b_1 U_t U_{t - 1} - \ldots - 2 b_{q - 1} b_q U_{t - q + 1} U_{t - q}] \nonumber \\ & = & \mathbb{E}[U_t^2] + b_1^2 \mathbb{E}[U_{t-1}^2] + b_2^2 \mathbb{E}[U_{t-2}^2] + \ldots + b_q^2 \mathbb{E}[U_{t-q}^2] \nonumber \\ & & - 2 b_1 \mathbb{E}[U_t U_{t - 1}] - \ldots - 2 b_{q - 1} b_q \mathbb{E}[U_{t - q + 1} U_{t - q}] \nonumber \\ & = & \sigma^2 + b^2_1 \sigma^2 + b^2_2 \sigma^2 + \ldots + b^2_q \sigma^2 \nonumber \\ & = & (1 + b^2_1 + b^2_2 + \ldots + b^2_q) \sigma^2 \end{eqnarray}\]

En el caso de las covarianzas podemos utilizar una idea similar al caso del $AR(p)$, construir una expresión general para cualquier rezago $\tau$:

\[\begin{eqnarray} Cov[X_t, X_{t + \tau}] & = & \mathbb{E}[(X_t - \mu)(X_{t + \tau} - \mu)] \nonumber \\ & = & \mathbb{E}[(U_t - b_1 U_{t-1} - b_2 U_{t-2} - \ldots - b_q U_{t-q}) \nonumber \\ & & (U_{t + \tau} - b_1 U_{t + \tau -1} - b_2 U_{t + \tau -2} - \ldots - b_q U_{t + \tau - q})] \nonumber \end{eqnarray}\]

La expresión anterior se puede desarrollar para múltiples casos de $\tau = 1, 2, \ldots, q$. De esta forma tenemos el siguiente sistema: \[\begin{eqnarray} \tau = 1 & : & \gamma(1) = (- b_1 + b_1 b_2 + \ldots + b_{q-1} b_q) \sigma^2 \nonumber \\ \tau = 2 & : & \gamma(2) = (- b_2 + b_1 b_3 + \ldots + b_{q-2} b_q) \sigma^2 \nonumber \\ & \vdots & \nonumber \\ \tau = q & : & \gamma(q) = - b_q \sigma^2 \nonumber \end{eqnarray}\]

Donde $\gamma(\tau) = 0$ para todo $\tau > q$. Es decir, todas las autocovarianzas y autocorrelaciones con órdenes superiores a $q$ son cero (0). De esta forma, esta característica teórica permite identificar el orden de $MA(q)$ visualizando la función de autocorrelación y verificando a partir de cuál valor de rezago la autocorrelación es no significativa.

Regresando al problema original que es el de determinar una solución para la ecuación (3.44), tenemos que dicha solución estará dada por un $AR(\infty)$ en términos de $U_t$:

\[\begin{eqnarray} U_t & = & - \frac{\mu}{1 - b_1 - b_2 - \ldots - b_q} + \beta(L)^{-1} X_t \nonumber \\ & = & - \frac{\mu}{1 - b_1 - b_2 - \ldots - b_q} + \sum_{j = 0}^{\infty} c_j X_{t-j} \tag{3.46} \end{eqnarray}\]

Donde $c_0 = 1$ y se cumple que: $1 = (1 - b_1 L^1 - b_2 L^2 - \ldots - b_q L^q)(c_0 + c_1 L + c_2 L^2 + \ldots)$, de modo que los coeficientes $c_j$ se pueden determinar por el método de coeficientes indeterminados en términos de los valores $b_i$. De igual forma que en el caso de la ecuación (3.38), en la ecuación (3.46) se deben cumplir condiciones de estabilidad asociadas con las raíces del polinomio característico dado por:

\[\begin{equation} 1 - b_1 x - b_2 x^2 - \ldots - b_q x^q = 0 \end{equation}\]

El cual debe cumplir que todas sus raíces se ubiquen fuera del círculo unitario ($|x_i| > 1$) o, equivalentemente, que las raíces del polinomio $\lambda^q - b_1 \lambda^{q-1} - \ldots - b_q = 0$ cumplan $|\lambda_i| < 1$ –la condición de invertibilidad del proceso–.

Ejemplo. Ahora veamos un ejemplo del proceso $MA(q)$, para lo cual retomaremos la serie de Pasajeros transportados en el metro de la CDMX ($PaxMetro$). Estimaremos el $MA(q)$ mediante el método de máxima verosimilitud (ML). Antes de realizar el proceso de estimación, consideremos una transformación de la serie en logaritmos y una más en diferencias logarítmicas; lo anterior con el objeto de obtener una serie de tiempo suavizada y expresada en tasas de crecimiento, con un comportamiento parecido a un proceso estacionario.

La serie de Pasajeros transportados en el metro de la CDMX se muestra en la Figura 3.9. En esta, se muestra la gráfica de la serie en niveles (sin transformación logarítmica y con transformación logarítmica) y en diferencias logarítmicas mensuales (es decir, con una diferencia respecto del mes inmediato anterior). Utilizaremos la serie en diferencias, ya que es la que parece ser estacionaria. Esta serie tiene la peculiaridad de que tiene un salto a la baja y uno al alza entre septiembre de 2017 y octubre de 2017. Para controlar ese efecto, en nuestro modelo $MA(q)$ incluiremos dos variables dummies para dichos meses.

A continuación, estimaremos una $MA(4)$ para la serie en diferencias. Los resultados se muestran en el Cuadro 3.7.

library(readxl)
Datos_MA <- read_excel("BD/Base_Transporte_ARIMA.xlsx",
                       sheet = "Datos", col_names = TRUE)
DLPax_Metro_MA <- ts(log(Datos_MA$Pax_Metro) - dplyr::lag(log(Datos_MA$Pax_Metro), 1),
                     start = c(2000, 1), freq = 12)
D_Sep2017_MA <- ts(Datos_MA$D_Sep2017, start = c(2000, 1), freq = 12)
t_idx <- time(DLPax_Metro_MA)
D_Oct2017_MA <- as.numeric(abs(t_idx - (2017 + 9/12)) < 0.01)

MA_DLPax_Metro <- arima(DLPax_Metro_MA, order = c(0, 0, 4),
                         xreg = cbind(D_Sep2017_MA, D_Oct2017_MA),
                         method = "ML")
arima_kable(MA_DLPax_Metro, "MA(4) para la variable $DLPaxMetro_t$")

Cuadro 3.7: MA(4) para la variable $DLPaxMetro_t$ ($\hat{\sigma}^2 = 0.008268$; AIC $=$ -602.3$)$.
	Estimado	Error Est.	Estad. $t$	Signif.
ma1	-0.1785	0.0638	-2.796	**
ma2	-0.0468	0.0761	-0.616
ma3	-0.1843	0.0976	-1.888	.
ma4	-0.0967	0.0791	-1.223
intercept	-0.0003	0.0026	-0.107
D_Sep2017_MA	-0.4019	0.0887	-4.533	***
D_Oct2017_MA	0.3940	0.0888	4.436	***

El cuadro reporta, además del coeficiente estimado, su error estándar, la estadística $t$ y el criterio de Akaike (AIC). De los cuatro coeficientes de medias móviles, sólo el primero resulta significativo al 1% y el tercero al 10%, lo que sugiere que un orden menor podría ser suficiente –más adelante formalizamos esta decisión con los criterios de información–. En cambio, las dos variables dicotómicas son altamente significativas y tienen signos opuestos y magnitudes similares, tal como corresponde a la caída de septiembre de 2017 y al rebote de octubre del mismo año. Finalmente, podemos determinar si la solución es invertible; para ello, en la Figura 3.14 mostramos los inversos de las raíces del polinomio de medias móviles. De la inspección visual –todas las raíces inversas se ubican dentro del círculo unitario– podemos concluir que el MA(4) estimado es invertible.

source("maroots.R")

source("plot.armaroots.R")

plot.armaroots(maroots(MA_DLPax_Metro),
               main = "Inverso de las raíces MA del \nMA(4): Diff LN Pax Metro")

Figura 3.14: Inverso de las raíces del polinomio de medias móviles del $MA(4)$ estimado para $DLPaxMetro_t$

3.6 Procesos ARMA(p, q) y ARIMA(p, d, q)

Hemos establecido algunas relaciones de los procesos AR y los procesos MA, es decir, cómo un $MA(q)$ de la serie $X_t$ puede ser reexpresado como un $AR(\infty)$ de la serie $U_t$ y, viceversa, un $AR(p)$ de la serie $X_t$ puede ser reexpresado como un $MA(\infty)$.

En este sentido, para cerrar esta sección veamos el caso de la especificación que conjunta ambos modelos en un modelo general conocido como $ARMA(p, q)$ o $ARIMA(p, d, q)$. La diferencia entre el primero y el segundo es el número de veces que se tuvo que diferenciar la serie analizada, registro que se lleva en el índice $d$ de los parámetros dentro del concepto $ARIMA(p, d, q)$. No obstante, en general nos referiremos al modelo como $ARMA(p, q)$ y dependerá del analista si modela la serie en niveles (por ejemplo, en logaritmos) o en diferencias logarítmicas (o diferencias sin logaritmos).

3.6.1 ARMA(1, 1)

Dicho lo anterior, vamos a empezar con el análisis de un $ARMA(1, 1)$. Un proceso $ARMA(1, 1)$ puede verse como:

\[\begin{equation} X_t = \delta + a_1 X_{t - 1} + U_t - b_1 U_{t - 1} \tag{3.47} \end{equation}\]

Aplicando el operador rezago podemos rescribir la ecuación (3.47) como:

\[\begin{equation} (1 - a_1 L) X_t = \delta + (1 - b_1 L) U_t \end{equation}\]

Donde $U_t$ es un proceso puramente aleatorio como en los casos de $AR(p)$ y $MA(q)$, y $X_t$ puede ser una serie en niveles o en diferencias (ambas, en términos logarítmicos).

Así, el modelo $ARMA(p, q)$ también tiene una representación de Wold que estará dada por las siguientes expresiones:

\[\begin{equation} X_t = \frac{\delta}{1 - a_1} + \frac{1 - b_1 L}{1 - a_1 L} U_t \tag{3.48} \end{equation}\]

Donde $a_1 \neq b_1$, puesto que en caso contrario $X_t$ sería un proceso puramente aleatorio con una media $\mu = \frac{\delta}{1 - a_1}$. Así, podemos reescribir la descomposición de Wold a partir del componente de la ecuación (3.48):

\[\begin{equation} \frac{1 - b_1 L}{1 - a_1 L} = \psi_0 + \psi_1 L + \psi_2 L^2 + \psi_3 L^3 + \ldots \tag{3.49} \end{equation}\]

Esta ecuación es equivalente a la expresión:

\[\begin{eqnarray} (1 - b_1 L) & = & (1 - a_1 L)(\psi_0 + \psi_1 L + \psi_2 L^2 + \psi_3 L^3 + \ldots) \nonumber \\ & = & \psi_0 + \psi_1 L + \psi_2 L^2 + \psi_3 L^3 + \ldots \nonumber \\ & & - a_1 \psi_0 L - a_1 \psi_1 L^2 - a_1 \psi_2 L^3 - a_1 \psi_3 L^4 - \ldots \nonumber \end{eqnarray}\]

De esta forma podemos establecer el siguiente sistema de coeficientes indeterminados:

\[\begin{eqnarray*} L^0 & : & \Rightarrow \psi_0 = 1 \\ L^1 & : & \psi_1 - a_1 \psi_0 = - b_1 \Rightarrow \psi_1 = a_1 - b_1 \\ L^2 & : & \psi_2 - a_1 \psi_1 = 0 \Rightarrow \psi_2 = a_1(a_1 - b_1) \\ L^3 & : & \psi_3 - a_1 \psi_2 = 0 \Rightarrow \psi_3 = a^2_1(a_1 - b_1) \\ & \vdots & \\ L^j & : & \psi_j - a_1 \psi_{j - 1} = 0 \Rightarrow \psi_j = a^{j - 1}_1(a_1 - b_1) \end{eqnarray*}\]

Así, la solución a la ecuación (3.47) estará dada por la siguiente generalización:

\[\begin{equation} X_t = \frac{\delta}{1 - a_1} + U_t + (a_1 - b_1) U_{t - 1} + a_1(a_1 - b_1) U_{t - 2} + a_1^2(a_1 - b_1) U_{t - 3} + \ldots \tag{3.50} \end{equation}\]

En la ecuación (3.50) las condiciones de estabilidad y de invertibilidad del sistema (de un MA a un AR, y viceversa) estarán dadas por: $|a_1| < 1$ y $|b_1| < 1$. Adicionalmente, la ecuación (3.50) expresa cómo una serie que tiene un comportamiento $ARMA(1, 1)$ es equivalente a una serie modelada bajo un $MA(\infty)$.

Al igual que en los demás modelos, ahora determinaremos los momentos del proceso $ARMA(1, 1)$. La media estará dada por:

\[\begin{eqnarray} \mathbb{E}[X_t] & = & \mathbb{E}[\delta + a_1 X_{t-1} + U_t - b_1 U_{t-1}] \nonumber \\ & = & \delta + a_1 \mathbb{E}[X_{t-1}] \nonumber \\ & = & \frac{\delta}{1 - a_1} \nonumber \\ & = & \mu \end{eqnarray}\]

Donde hemos utilizado que $\mathbb{E}[X_t] = \mathbb{E}[X_{t-1}] = \mu$. Es decir, la media de un $ARMA(1, 1)$ es idéntica a la de un $AR(1)$.

Para determinar la varianza tomaremos una estrategia similar a los casos de $AR(p)$ y $MA(q)$. Por lo que para todo $\tau \geq 0$, y suponiendo por simplicidad que $\delta = 0$ (lo que implica que $\mu = 0$) tendremos:

\[\begin{eqnarray} \mathbb{E}[X_{t-\tau} X_t] & = & \mathbb{E}[(X_{t-\tau}) \cdot (a_1 X_{t-1} + U_t - b_1 U_{t-1})] \nonumber \\ & = & a_1 \mathbb{E}[X_{t-\tau} X_{t-1}] + \mathbb{E}[X_{t-\tau} U_t] - b_1 \mathbb{E}[X_{t-\tau} U_{t-1}] \tag{3.51} \end{eqnarray}\]

De la ecuación (3.51) podemos determinar una expresión para el caso de $\tau = 0$:

\[\begin{eqnarray} \mathbb{E}[X_{t} X_t] & = & \gamma(0) \nonumber \\ & = & a_1 \gamma(1) + \mathbb{E}[U_t X_t] - b_1 \mathbb{E}[X_t U_{t-1}] \nonumber \\ & = & a_1 \gamma(1) + \sigma^2 - b_1 \mathbb{E}[U_{t-1} (a_1 X_{t-1} + U_t - b_1 U_{t-1})] \nonumber \\ & = & a_1 \gamma(1) + \sigma^2 - b_1 a_1 \sigma^2 + b_1^2 \sigma^2 \nonumber \\ & = & a_1 \gamma(1) + (1 - b_1 (a_1 - b_1)) \sigma^2 \end{eqnarray}\]

Para el caso en que $\tau = 1$:

\[\begin{eqnarray} \mathbb{E}[X_{t-1} X_t] & = & \gamma(1) \nonumber \\ & = & a_1 \gamma(0) + \mathbb{E}[X_{t-1} U_t] - b_1 \mathbb{E}[X_{t-1} U_{t-1}] \nonumber \\ & = & a_1 \gamma(0) - b_1 \sigma^2 \end{eqnarray}\]

Estas últimas expresiones podemos resolverlas como sistema para determinar los siguientes valores:

\[\begin{eqnarray} \gamma(0) & = & \frac{1 + b_1^2 - 2 a_1 b_1}{1 - a_1^2} \sigma^2 \\ \gamma(1) & = & \frac{(a_1 - b_1)(1 - a_1 b_1)}{1 - a_1^2} \sigma^2 \end{eqnarray}\]

En general, para cualquier valor $\tau \geq 2$ tenemos que la autocovarianza y la función de autocorrelación serán:

\[\begin{eqnarray} \gamma(\tau) = a_1 \gamma(\tau - 1) \\ \rho(\tau) = a_1 \rho(\tau - 1) \end{eqnarray}\]

Por ejemplo, para el caso de $\tau = 1$ tendríamos:

\[\begin{equation} \rho(1) = \frac{(a_1 - b_1)(1 - a_1 b_1)}{1 + b_1^2 - 2 a_1 b_1} \end{equation}\]

De esta forma, la función de autocorrelación de un $ARMA(1, 1)$ arranca en un valor $\rho(1)$ que depende de $a_1$ y $b_1$ y, a partir de ahí, decae geométricamente a razón de $a_1$ –si $a_1 < 0$, el decaimiento es oscilante–. Este patrón es justamente el que dificulta distinguir un $ARMA(1,1)$ de un $AR(1)$ mediante la sola inspección del correlograma.

3.6.2 ARMA(p, q)

La especificación general de un $ARMA(p, q)$ (donde $p, q \in \mathbb{N}$) puede ser descrita por la siguiente ecuación: \[\begin{eqnarray} X_t & = & \delta + a_1 X_{t - 1} + a_2 X_{t - 2} + \ldots + a_p X_{t - p} \nonumber \\ & & + U_t - b_1 U_{t - 1} - b_2 U_{t - 2} - \ldots - b_q U_{t - q} \tag{3.52} \end{eqnarray}\]

Donde $U_t$ es un proceso puramente aleatorio, y $X_t$ puede ser modelada en niveles o en diferencias (ya sea en logaritmos o sin transformación logarítmica).

Mediante el uso del operador rezago se puede escribir la ecuación (3.52) como:

\[\begin{equation} (1 - a_1 L - a_2 L^2 - \ldots - a_p L^p) X_t = \delta + (1 - b_1 L - b_2 L^2 - \ldots - b_q L^q) U_t \tag{3.53} \end{equation}\]

En la ecuación (3.53) definamos dos polinomios: $\alpha(L) = (1 - a_1 L - a_2 L^2 - \ldots - a_p L^p)$ y $\beta(L) = (1 - b_1 L - b_2 L^2 - \ldots - b_q L^q)$. Así, podemos reescribir la ecuación (3.53) como:

\[\begin{equation} \alpha(L) X_t = \delta + \beta(L) U_t \end{equation}\]

Asumiendo que existe el polinomio inverso tal que: $\alpha(L)^{-1}\alpha(L) = 1$. La solución entonces puede ser escrita como:

\[\begin{eqnarray} X_t & = & \alpha(L)^{-1} \delta + \alpha(L)^{-1} \beta(L) U_t \nonumber \\ & = & \frac{\delta}{1 - a_1 - a_2 - \ldots - a_p} + \frac{\beta(L)}{\alpha(L)} U_t \nonumber \\ & = & \frac{\delta}{1 - a_1 - a_2 - \ldots - a_p} + U_t + \psi_1 L U_t + \psi_2 L^2 U_t + \ldots \tag{3.54} \end{eqnarray}\]

Donde la ecuación (3.54) nos permite interpretar que un ARMA(p, q) se puede reexpresar e interpretar como un $MA(\infty)$ y donde las condiciones para la estabilidad de la solución y la invertibilidad son que las raíces de las ecuaciones características asociadas a $\alpha(L)$ y $\beta(L)$ sean, en valor absoluto, menores a 1 –equivalentemente, que las raíces de los polinomios en $L$ se ubiquen fuera del círculo unitario–.

Adicionalmente, la fracción en la ecuación (3.54) se puede descomponer como en la forma de Wold:

\[\begin{equation} \frac{\beta(L)}{\alpha(L)} = 1 + \psi_1 L + \psi_2 L^2 + \ldots \end{equation}\]

Bajo los supuestos de estacionariedad del componente $U_t$, los valores de la media y varianza de un proceso $ARMA(p, q)$ serán como describimos ahora. Para el caso de la media, podemos partir de la ecuación (3.54) para generar:

\[\begin{eqnarray} \mathbb{E}[X_t] & = & \mathbb{E}\left[ \frac{\delta}{1 - a_1 - a_2 - \ldots - a_p} + U_t + \psi_1 U_{t-1} + \psi_2 U_{t-2} + \ldots \right] \nonumber \\ & = & \frac{\delta}{1 - a_1 - a_2 - \ldots - a_p} \nonumber \\ & = & \mu \end{eqnarray}\]

Esta expresión indica que en general un proceso $ARMA(p, q)$ converge a una media idéntica a la de un proceso $AR(p)$. Para determinar la varianza utilizaremos la misma estrategia que hemos utilizado para otros modelos $AR(p)$ y $MA(q)$.

Sin pérdida de generalidad podemos asumir que $\delta = 0$, lo que implica que $\mu = 0$, de lo que podemos establecer una expresión de autocovarianzas para cualquier valor $\tau = 0, 1, 2, \ldots$: \[\begin{eqnarray} \gamma(\tau) & = & \mathbb{E}[X_{t-\tau} X_t] \nonumber \\ & = & \mathbb{E}[X_{t-\tau} (\delta + a_1 X_{t - 1} + a_2 X_{t - 2} + \ldots + a_p X_{t - p} \nonumber \\ & & + U_t - b_1 U_{t - 1} - b_2 U_{t - 2} - \ldots - b_q U_{t - q})] \nonumber \\ & = & a_1 \gamma(\tau - 1) + a_2 \gamma(\tau - 2) + \ldots + a_p \gamma(\tau - p) \nonumber \\ & & + \mathbb{E}[X_{t-\tau} U_{t}] - b_1 \mathbb{E}[X_{t-\tau} U_{t-1}] - \ldots - b_q \mathbb{E}[X_{t-\tau} U_{t-q}] \end{eqnarray}\]

Ahora veamos un ejemplo. Utilizaremos la serie de Pasajeros en vuelos nacionales de salida para estimar un $ARMA(p, q)$ mediante el método de máxima verosimilitud (ML). Antes de realizar el proceso de estimación, consideremos una transformación de la serie en diferencias logarítmicas, ya que, según la gráfica en la Figura 3.10, esa es la que puede ser estacionaria.

A continuación, estimaremos una $ARMA(1, 1)$ para la serie en diferencias logarítmicas ($DLPaxNal_t$). También incorporaremos al análisis variables exógenas tales como dummies de estacionalidad. En particular, utilizaremos los meses de enero, febrero, julio y diciembre. No debe pasar desapercibido que un análisis de estacionalidad más formal debería considerar todos los meses para separar del término de error la parte que puede ser explicada por los ciclos estacionales.

Los resultados se muestran en el Cuadro 3.8.

D_Ene_Nal <- as.numeric(cycle(DLPax_Nal) == 1)
D_Feb_Nal <- as.numeric(cycle(DLPax_Nal) == 2)
D_Jul_Nal <- as.numeric(cycle(DLPax_Nal) == 7)
D_Dic_Nal <- as.numeric(cycle(DLPax_Nal) == 12)

ARMA_DLPax_Nal <- arima(DLPax_Nal, order = c(1, 0, 1),
                         xreg = cbind(D_Ene_Nal, D_Feb_Nal,
                                      D_Jul_Nal, D_Dic_Nal),
                         method = "ML")
arima_kable(ARMA_DLPax_Nal, "ARMA(1, 1) para la variable $DLPaxNal_t$")

Cuadro 3.8: ARMA(1, 1) para la variable $DLPaxNal_t$ ($\hat{\sigma}^2 = 0.022281$; AIC $=$ -289.08$)$.
	Estimado	Error Est.	Estad. $t$	Signif.
ar1	-0.5903	0.0958	-6.163	***
ma1	0.8197	0.0621	13.206	***
intercept	0.0035	0.0116	0.305
D_Ene_Nal	-0.0910	0.0320	-2.843	**
D_Feb_Nal	-0.1280	0.0328	-3.903	***
D_Jul_Nal	0.1701	0.0295	5.764	***
D_Dic_Nal	0.0603	0.0313	1.924	.

El cuadro reporta el coeficiente estimado, su error estándar y la estadística $t$ de cada parámetro, además del criterio de Akaike (AIC). Finalmente, podemos determinar si las soluciones serán convergentes, para ello en la Figura 3.15 mostramos las raíces asociadas a cada uno de los polinomios. De la inspección visual, podemos concluir que tenemos una solución convergente y estable. Por su parte, la Figura 3.16 muestra los residuales de la estimación del $ARMA(1, 1)$.

source("maroots.R")

par(mfrow=c(1,2))

plot.armaroots(arroots(ARMA_DLPax_Nal),
               main="Inverse AR roots of \nARMA(1,1): Diff LN Pax Nal")

plot.armaroots(maroots(ARMA_DLPax_Nal),
               main="Inverse MA roots of \nARMA(1,1): Diff LN Pax Nal")

Figura 3.15: Inverso de las Raíces de los polinomios característicos de un ARMA(1,1)

par(mfrow=c(1,1))

plot(residuals(ARMA_DLPax_Nal), col = "darkred",
     ylab = "Residuales", xlab = "Tiempo",
     main = "Residuales del ARMA(1,1) de Diff LN Pax Nal")

Figura 3.16: Residuales de un $ARMA(1, 1)$ de la serie $DLPaxNal_t$

En lo que resta de este capítulo, utilizaremos la serie en diferencias logarítmicas de los pasajeros en vuelos nacionales de salida, $DLPaxNal_t$, para discutir los ejemplos que ilustran cada uno de los puntos teóricos que a continuación exponemos.

3.7 Función de Autocorrelación Parcial

Ahora introduciremos el concepto de Función de Autocorrelación Parcial (PACF, por sus siglas en inglés). Primero, dadas las condiciones de estabilidad y de convergencia, si suponemos que los procesos AR, MA, ARMA o ARIMA tienen toda la información de los rezagos de la serie en conjunto y toda la información de los promedios móviles del término de error, resulta importante construir una métrica para distinguir el efecto de $X_{t - \tau}$ o el efecto de $U_{t - \tau}$ (para cualquier $\tau$) sobre $X_t$ de forma individual.

La idea es construir una métrica de la correlación que existe entre las diferentes variables aleatorias, si para tal efecto se ha controlado el efecto del resto de la información. Así, podemos definir la ecuación que puede responder a este planteamiento como:

\[\begin{equation} X_t = \phi_{k1} X_{t-1} + \phi_{k2} X_{t-2} + \ldots + \phi_{kk} X_{t-k} + U_t \tag{3.55} \end{equation}\]

Donde $\phi_{ki}$ es el coeficiente de la variable dada con el rezago $i$ si el proceso tiene un orden $k$. Así, los coeficientes $\phi_{kk}$ son los coeficientes de la autocorrelación parcial (considerando un proceso AR(k)). Observemos que la autocorrelación parcial mide la correlación entre $X_t$ y $X_{t-k}$ que se mantiene cuando el efecto de las variables $X_{t-1}$, $X_{t-2}$, $\ldots$ y $X_{t-(k-1)}$ en $X_{t}$ y $X_{t-k}$ ha sido eliminado.

Dada la expresión considerada en la ecuación (3.55), podemos resolver el problema de establecer el valor de cada $\phi_{ki}$ mediante la solución del sistema que se representa en lo siguiente: \[\begin{equation} \left[ \begin{array}{c} \rho(1) \\ \rho(2) \\ \vdots \\ \rho(k) \end{array} \right] = \left[ \begin{array}{c c c c} 1 & \rho(1) & \ldots & \rho(k - 1)\\ \rho(1) & 1 & \ldots & \rho(k - 2)\\ \rho(2) & \rho(1) & \ldots & \rho(k - 3)\\ \vdots & \vdots & \ldots & \vdots\\ \rho(k - 1) & \rho(k - 2) & \ldots & 1\\ \end{array} \right] \left[ \begin{array}{c} \phi_{k1} \\ \phi_{k2} \\ \phi_{k3} \\ \vdots \\ \phi_{kk} \\ \end{array} \right] \end{equation}\]

Del cual se puede derivar una solución, resolviendo por cualquier método que consideremos y que permita calcular la solución de sistemas de ecuaciones.

Posterior al planteamiento analítico, plantearemos un enfoque para interpretar las funciones de autocorrelación y autocorrelación parcial. Este enfoque pretende aportar al principio de parsimonia, en el cual podemos identificar el número de parámetros que posiblemente puede describir mejor a la serie en un modelo ARMA(p, q).

En el siguiente cuadro mostramos un resumen de las características que debemos observar para determinar el número de parámetros de cada uno de los componentes AR y MA. Lo anterior por observación de las funciones de autocorrelación y autocorrelación parcial. Este enfoque no es el más formal, más adelante implementaremos uno más formal y que puede ser más claro de cómo determinar el número de parámetros.

Cuadro 3.9: Relación entre la Función de Autocorrelación y la Función de Autocorrelación Parcial de una serie $X_t$.
Modelo	Función de Autocorrelación	Función de Autocorrelación Parcial
MA(q)	Se corta después del rezago $q$	Decae gradualmente
AR(p)	Decae gradualmente	Se corta después del rezago $p$

Ejemplo. Continuando con el ejemplo en la Figura 3.17 mostramos tanto la Función de Autocorrelación como la Función de Autocorrelación Parcial. En esta identificamos que ambas gráficas muestran que el modelo que explica a la variable $DLPaxNal_t$ tiene tanto componentes AR como MA. Sin embargo, dado lo errático del comportamiento de ambas funciones, resulta complicado determinar cuál sería un buen número de parámetros $p$ y $q$ a considerar en el $ARMA(p,q)$. Por esta razón, a continuación, plantearemos algunas pruebas más formales para determinar dichos parámetros.

par(mfrow=c(2,1))

acf(na.omit(DLPax_Nal), lag.max = 48,
    xlab = "Rezagos k en meses",
    main = "Función de Autocorrelación de Diff LN Pax Nal")

pacf(na.omit(DLPax_Nal), lag.max = 48,
     xlab = "Rezagos k en meses",
     main = "Función de Autocorrelación Parcial de Diff LN Pax Nal")

Figura 3.17: Función de Autocorrelación y la Función de Autocorrelación Parcial de una serie $DLPaxNal_t$

par(mfrow=c(1,1))

3.8 Selección de las constantes p, q, d en un AR(p), un MA(q), un ARMA(p, q) o un ARIMA(p, d, q)

Respecto de cómo estimar un proceso ARMA(p, q) –en general utilizaremos este modelo para discutir, pero lo planteado en esta sección es igualmente aplicable en cualquier otro caso como aquellos modelos que incluyen variables exógenas– existen diversas formas de estimar los parámetros $a_i$ y $b_i$: i) por máxima verosimilitud y ii) por mínimos cuadrados ordinarios. El primer caso requiere que conozcamos la distribución del proceso aleatorio $U_t$. El segundo, por el contrario, no requiere el mismo supuesto. No obstante, en este libro utilizaremos el método de máxima verosimilitud.

Otra duda que debe quedar hasta el momento es ¿cómo determinar el orden $p$ y $q$ del proceso ARMA(p, q)? La manera más convencional y formal que existe para tal efecto es utilizar los criterios de información. Así, el orden se elige de acuerdo con aquel criterio de información que resulta ser el mínimo. En el caso de $d$ se selecciona revisando la gráfica que parezca más estacionaria–más adelante mostraremos un proceso más formal para su selección.

Los criterios de información más comunes son los siguientes:

FPE (Final Prediction Error): \[\begin{equation} FPE = \frac{T+m}{T-m}\frac{1}{T}\sum_{t=1}^{T} \left( \hat{U}_t^{(p)} \right) ^2 \end{equation}\]
Akaike: \[\begin{equation} AIC = ln \left[ \frac{1}{T} \sum_{t=1}^{T} \left( \hat{U}_t^{(p)} \right) ^2 \right] + m \frac{2}{T} \end{equation}\]
Schwarz: \[\begin{equation} SC = ln \left[ \frac{1}{T} \sum_{t=1}^{T} \left( \hat{U}_t^{(p)} \right) ^2 \right] + m \frac{ln(T)}{T} \end{equation}\]
Hannan - Quinn: \[\begin{equation} HQ = ln \left[ \frac{1}{T} \sum_{t=1}^{T} \left( \hat{U}_t^{(p)} \right) ^2 \right] + m \frac{2 ln(ln(T))}{T} \end{equation}\]

Donde $\hat{U}_t^{(p)}$ son los residuales estimados mediante un proceso ARIMA y $m$ es el número de parámetros estimados: $m = p + q + 1$ (los $p$ coeficientes autorregresivos, los $q$ de medias móviles y el término constante; asumimos $d = 0$). Una propiedad que no se debe perder de vista es que los criterios de información cumplen la siguiente relación: \[\begin{equation} orden(SC) \leq orden(HQ) \leq orden(AIC) \end{equation}\]

Es decir, el criterio de Schwarz selecciona órdenes menores o iguales que los del criterio de Hannan-Quinn y éste, a su vez, menores o iguales que los del criterio de Akaike (el criterio SC también se conoce como BIC, por Bayesian Information Criterion). En este libro utilizaremos principalmente el criterio de Akaike, ya que en el trabajo aplicado con fines de pronóstico suele ser más costoso subparametrizar la dinámica –dejar autocorrelación en los residuales– que estimar algunos parámetros de más. Conviene advertir, no obstante, que el criterio de Schwarz es consistente en la selección del orden verdadero cuando éste es finito, mientras que el de Akaike tiende a sobreestimarlo; por ello, la práctica recomendable es reportar ambos y verificar que el modelo elegido supere las pruebas de diagnóstico de los residuales.

Ejemplo. Ahora veamos un ejemplo de estimación del número de rezagos óptimo de un $ARMA(p, q)$. Retomemos la serie en diferencias logarítmicas de los pasajeros en vuelos nacionales de salidas, pero ahora incluiremos las variables exógenas de dummies estacionales: enero, febrero, julio y diciembre.

Como mencionamos, las gráficas de las funciones de autocorrelación permiten observar el valor de la correlación existente entre la variable en el momento $t$ con cada uno de los rezagos. Incluso la Función de Autocorrelación Parcial puede ayudar a determinar el número máximo de rezagos que se debe incluir en el proceso $AR(p)$. No obstante, una métrica más formal es el uso de los criterios de información. En nuestro caso, dado lo discutido, sólo utilizaremos el criterio de Akaike.

Al respecto, en el siguiente cuadro reportamos el criterio de Akaike que resulta de aplicar dicho criterio a los residuales resultantes de cada combinación de procesos $ARMA(p, q)$. La forma de escoger será aquel modelo que reporta el criterio de Akaike menor. En la última columna del cuadro se marca con un asterisco el modelo cuyo criterio de información resulta ser el mínimo de todos los posibles.

xreg_base <- cbind(D_Ene_Nal, D_Feb_Nal, D_Jul_Nal, D_Dic_Nal)
ciarma_rows <- vector("list", 36)
k <- 0L
for (p_ord in 1:6) {
  for (q_ord in 1:6) {
    k <- k + 1L
    fit_aic <- tryCatch(
      AIC(arima(DLPax_Nal, order = c(p_ord, 0L, q_ord),
                xreg = xreg_base, method = "ML")),
      error = function(e) NA_real_
    )
    ciarma_rows[[k]] <- data.frame(Modelo = k, AR = p_ord, MA = q_ord,
                                    AIC = round(fit_aic, 4))
  }
}
df_ciarma <- do.call(rbind, ciarma_rows)
opt_idx <- which.min(df_ciarma$AIC)
opt_p <- df_ciarma$AR[opt_idx]
opt_q <- df_ciarma$MA[opt_idx]
df_ciarma$Optimo <- ifelse(seq_len(nrow(df_ciarma)) == opt_idx, "*", "")
knitr::kable(df_ciarma,
             col.names = c("Modelo", "$AR(p)$", "$MA(q)$", "Akaike (AIC)", "Óptimo"),
             caption = "Criterio de Akaike para diferentes modelos $ARMA(p, q)$ de la serie $DLPaxNal_t$.",
             align = c("c","c","c","r","c"), booktabs = TRUE, escape = FALSE)

Cuadro 3.10: Criterio de Akaike para diferentes modelos $ARMA(p, q)$ de la serie $DLPaxNal_t$.
Modelo	$AR(p)$	$MA(q)$	Akaike (AIC)
1	1	1	-289.0818
2	1	2	-320.3322
3	1	3	-319.2098
4	1	4	-317.7313
5	1	5	-315.7919
6	1	6	-331.8578
7	2	1	-314.6954
8	2	2	-320.9990
9	2	3	-321.5749
10	2	4	-320.3133
11	2	5	-324.5436
12	2	6	-329.9303
13	3	1	-318.3860
14	3	2	-321.6786
15	3	3	-323.4454
16	3	4	-389.8864
17	3	5	-388.5931
18	3	6	-387.5791
19	4	1	-290.3528
20	4	2	-320.0947
21	4	3	-347.0402
22	4	4	-348.1351
23	4	5	-349.3588
24	4	6	-383.8267
25	5	1	-315.6990
26	5	2	-319.8355
27	5	3	-387.9538
28	5	4	-350.6625
29	5	5	-348.6939
30	5	6	-385.8841
31	6	1	-331.6773
32	6	2	-329.9577
33	6	3	-388.3933
34	6	4	-373.6583
35	6	5	-384.0884
36	6	6	-383.7144

El cuadro anterior reporta un resumen de los resultados para 36 diferentes modelos, todos incluyen variables exógenas. Como resultado del análisis concluimos que el modelo más adecuado es el 16, el cual considera un $ARMA(3, 4)$, con variables dummies para controlar la estacionalidad de los meses de enero, febrero, julio y diciembre. Más adelante mostramos los resultados del modelo.

No obstante, una inspección de los residuales del modelo nos permite sospechar que se requiere incluir un par de variables dicotómicas más, ambas asociadas con la caída del transporte aéreo de mayo y junio de 2009: en esos meses coincidieron la epidemia de influenza A(H1N1) en México –que llevó a la suspensión de actividades y a restricciones de viaje– y la crisis financiera mundial. La Figura 3.18 muestra los residuales mencionados.

mod1 <- arima(DLPax_Nal, order = c(opt_p, 0, opt_q),
              xreg = xreg_base, method = "ML")

plot(residuals(mod1), col = "darkred",
     ylab = "Residuales", xlab = "Tiempo",
     main = paste0("Residuales del ARMA(", opt_p, ",", opt_q,
                   ") de Diff LN Pax Nal"))

Figura 3.18: Residuales del ARMA(3, 4) de una serie $DLPaxNal_t$

Una vez incluidas dos dummies más para mayo y junio de 2009, reestimamos el modelo con el orden óptimo antes seleccionado. El siguiente cuadro muestra los resultados de ambas especificaciones. No lo mostramos en esta sección, pero ambos modelos reportados tienen raíces de sus respectivos polinomios característicos menores a 1 en valor absoluto. En la Figura 3.19 mostramos los residuales ahora ajustados por las dummies de mayo y junio de 2009.

t_idx_nal <- time(DLPax_Nal)
D_May2009 <- as.numeric(abs(t_idx_nal - (2009 + 4/12)) < 0.01)
D_Jun2009 <- as.numeric(abs(t_idx_nal - (2009 + 5/12)) < 0.01)

mod2 <- arima(DLPax_Nal, order = c(opt_p, 0, opt_q),
              xreg = cbind(xreg_base, D_May2009, D_Jun2009), method = "ML")

lbl_map <- c(setNames(paste0("$DLPaxNal_{t-", 1:6, "}$"), paste0("ar", 1:6)),
             setNames(paste0("$\\hat{U}_{t-", 1:6, "}$"), paste0("ma", 1:6)),
             D_Ene_Nal="$DEne_{t}$", D_Feb_Nal="$DFeb_{t}$",
             D_Jul_Nal="$DJul_{t}$", D_Dic_Nal="$DDic_{t}$",
             intercept="Constante",
             D_May2009="$DMay2009_{t}$", D_Jun2009="$DJun2009_{t}$")

make_col <- function(mod) {
  cf <- mod$coef; se <- sqrt(diag(mod$var.coef))
  data.frame(nm = names(cf), cf = round(cf, 4), se = round(se, 4))
}
d1 <- make_col(mod1); d2 <- make_col(mod2)
d_all <- merge(d1, d2, by = "nm", all = TRUE)
d_all$lbl <- lbl_map[d_all$nm]
d_out <- d_all[, c("lbl", "cf.x", "se.x", "cf.y", "se.y")]
d_out <- rbind(d_out,
               data.frame(lbl = c("$\\hat{\\sigma}^2$", "AIC"),
                          cf.x = c(round(mod1$sigma2, 6), round(AIC(mod1), 2)),
                          se.x = NA_real_, cf.y = c(round(mod2$sigma2, 6), round(AIC(mod2), 2)),
                          se.y = NA_real_))
d_out$se.x <- ifelse(is.na(d_out$se.x), "", as.character(d_out$se.x))
d_out$se.y <- ifelse(is.na(d_out$se.y), "", as.character(d_out$se.y))
d_out$cf.x <- ifelse(is.na(d_out$cf.x), "", as.character(d_out$cf.x))
d_out$cf.y <- ifelse(is.na(d_out$cf.y), "", as.character(d_out$cf.y))
knitr::kable(d_out,
             col.names = c("Variable", "Coef. M1", "SE M1", "Coef. M2", "SE M2"),
             caption = "Modelo $ARMA(p, q)$ de la serie $DLPaxNal_t$.",
             align = c("l","r","r","r","r"), booktabs = TRUE)

Cuadro 3.11: Modelo $ARMA(p, q)$ de la serie $DLPaxNal_t$.
Variable	Coef. M1	SE M1	Coef. M2	SE M2
$DLPaxNal_{t-1}$	-1.0184	0.0403	-1.0322	0.0413
$DLPaxNal_{t-2}$	0.2354	0.0698	0.2108	0.0717
$DLPaxNal_{t-3}$	0.7133	0.0403	0.6992	0.0414
$DDic_{t}$	0.0536	0.0365	0.0352	0.0358
$DEne_{t}$	-0.1249	0.0343	-0.1207	0.034
$DFeb_{t}$	-0.0253	0.0366	-0.037	0.0359
$DJul_{t}$	0.2644	0.0446	0.2582	0.0447
$DJun2009_{t}$			0.2031	0.0868
$DMay2009_{t}$			-0.372	0.0866
Constante	-0.0098	0.0084	-0.0064	0.0084
$\hat{U}_{t-1}$	1.2669	0.0023	1.3318	0.0079
$\hat{U}_{t-2}$	-0.3332	0.0214	-0.2852	0.0156
$\hat{U}_{t-3}$	-1.3883	0.0207	-1.4342	0.008
$\hat{U}_{t-4}$	-0.5395		-0.6012	0.0119
$\hat{\sigma}^2$	0.015149		0.014124
AIC	-389.89		-406.67

plot(residuals(mod2), col = "darkblue",
     ylab = "Residuales", xlab = "Tiempo",
     main = paste0("Residuales del ARMA(", opt_p, ",", opt_q,
                   ") con dummies de mayo y junio de 2009"))

Figura 3.19: Residuales del ARMA(3, 4) con dummies de 2009 de una serie $DLPaxNal_t$

3.9 Pronósticos

Para pronosticar el valor de la serie es necesario determinar cuál es el valor esperado de la serie en un momento $t + \tau$ condicional en que ésta se comporta como un $AR(p)$, un $MA(q)$ o un $ARMA(p, q)$ y a que los valores hasta $t$ están dados. Denotemos con $\mathbb{E}_t[\cdot]$ a la esperanza condicional en la información disponible en $t$, es decir, $\mathbb{E}_t[\cdot] = \mathbb{E}[\cdot \, | \, X_t, X_{t-1}, \ldots]$. Aplicando este operador a la ecuación (3.52) obtenemos la regla recursiva de pronóstico:

\[\begin{eqnarray} \mathbb{E}_t[X_{t+\tau}] & = & \delta + a_1 \mathbb{E}_t[X_{t+\tau-1}] + a_2 \mathbb{E}_t[X_{t+\tau-2}] + \ldots + a_p \mathbb{E}_t[X_{t+\tau-p}] \nonumber \\ & & - b_1 \mathbb{E}_t[U_{t+\tau-1}] - b_2 \mathbb{E}_t[U_{t+\tau-2}] - \ldots - b_q \mathbb{E}_t[U_{t+\tau-q}] \tag{3.56} \end{eqnarray}\]

Donde, para aplicar la recursión con $\tau = 1, 2, \ldots$, se usan las dos reglas siguientes:

$\mathbb{E}_t[X_{t+j}] = X_{t+j}$ si $j \leq 0$ (el valor ya se observó) y $\mathbb{E}_t[X_{t+j}]$ es el pronóstico calculado en el paso previo si $j > 0$.
$\mathbb{E}_t[U_{t+j}] = \hat{U}_{t+j}$ si $j \leq 0$ (el residual ya es conocido) y $\mathbb{E}_t[U_{t+j}] = 0$ si $j > 0$.

Nótese que los términos de medias móviles no desaparecen del pronóstico de manera automática: sólo se anulan los errores futuros. En consecuencia, los componentes $MA(q)$ dejan de contribuir al pronóstico únicamente a partir del horizonte $\tau > q$, mientras que para $\tau \leq q$ los residuales ya observados siguen apareciendo en la ecuación (3.56).

De lo anterior se derivan dos resultados importantes. Primero, el pronóstico así construido es óptimo en el sentido de que minimiza el error cuadrático medio de predicción, ya que la esperanza condicional es la mejor aproximación en ese sentido. Segundo, usando la representación de Wold de la ecuación (3.54), el error de pronóstico a $\tau$ periodos es una suma de los choques ocurridos entre $t+1$ y $t+\tau$:

\[\begin{equation} e_t(\tau) = X_{t+\tau} - \mathbb{E}_t[X_{t+\tau}] = U_{t+\tau} + \psi_1 U_{t+\tau-1} + \ldots + \psi_{\tau-1} U_{t+1} \tag{3.57} \end{equation}\]

Por lo que su varianza es:

\[\begin{equation} Var[e_t(\tau)] = \sigma^2 \left( 1 + \psi_1^2 + \psi_2^2 + \ldots + \psi_{\tau-1}^2 \right) \tag{3.58} \end{equation}\]

La ecuación (3.58) es creciente en $\tau$: la incertidumbre del pronóstico se acumula con el horizonte. Además, en un proceso estacionario la suma converge a $\gamma(0)$, de modo que para horizontes largos el intervalo de predicción se estabiliza en torno a la media $\mu$ con una amplitud dada por la desviación estándar no condicional de la serie –a diferencia de un proceso con raíz unitaria, en el que la varianza del error crece sin límite (véase el Capítulo 4)–. Suponiendo normalidad de los errores, el intervalo de predicción al 95% se construye como:

\[\begin{equation} \mathbb{E}_t[X_{t+\tau}] \pm 1.96 \sqrt{Var[e_t(\tau)]} \end{equation}\]

En la práctica, $\sigma^2$ y los $\psi_j$ se sustituyen por sus estimaciones, por lo que el intervalo resultante subestima ligeramente la incertidumbre total: no incorpora el error de estimación de los parámetros ni la incertidumbre asociada a la elección del modelo. La función predict() de R devuelve tanto el pronóstico puntual ($pred) como los errores estándar de la ecuación (3.58) ($se), con los que se construyen estos intervalos.

Continuando con el ejemplo, en la Figura 3.20 mostramos el resultado del pronóstico de la serie a 24 meses. Para este ejercicio estimamos un $ARMA(6, 6)$ sobre la serie en diferencias logarítmicas que, además de las dummies estacionales, incorpora variables dummy para los meses atípicos de la pandemia de COVID-19 (marzo, abril, junio y julio de 2020, y marzo de 2021); el pronóstico de la tasa de crecimiento se acumula después sobre el último nivel observado de la serie.

library(ggplot2)
library(dplyr)
library(readxl)
library(stats)

Datos <- read_excel("BD/Base_Transporte_ARIMA.xlsx", 
                    sheet = "Datos", col_names = TRUE)

Pax_Nal <- ts(Datos$Pax_Nal, 
              start = c(2000, 1), 
              freq = 12)

LPax_Nal <- ts(log(Datos$Pax_Nal), 
               start = c(2000, 1), 
               freq = 12)

DLPax_Nal <- ts(log(Datos$Pax_Nal) - dplyr::lag(log(Datos$Pax_Nal), 1),
                start = c(2000, 1), 
                freq = 12)

D_Mar2020   <- ts(Datos$D_Mar2020, 
                start = c(2000, 1), 
                freq = 12)

D_Abr2020   <- ts(Datos$D_Abr2020, 
                start = c(2000, 1), 
                freq = 12)

D_Jun2020   <- ts(Datos$D_Jun2020, 
                start = c(2000, 1), 
                freq = 12)

D_Jul2020 <- ts(Datos$D_Jul2020, 
                start = c(2000, 1), 
                freq = 12)

D_Mar2021 <- ts(Datos$D_Mar2021, 
                start = c(2000, 1), 
                freq = 12)

D_Ene <- ts(Datos$D_Ene, 
            start = c(2000, 1), 
            freq = 12)

D_Feb <- ts(Datos$D_Feb, 
            start = c(2000, 1), 
            freq = 12)

D_Jul <- ts(Datos$D_Jul, 
            start = c(2000, 1), 
            freq = 12)

D_Dic <- ts(Datos$D_Dic, 
            start = c(2000, 1), 
            freq = 12)

ARMA_Ex_DLPax_Nal_2 <- arima(DLPax_Nal, order = c(6, 0, 6),
                             xreg = cbind(D_Ene, D_Feb, D_Jul, 
                                          D_Dic, D_Mar2020, 
                                          D_Abr2020, D_Jun2020, 
                                          D_Jul2020, D_Mar2021),
                             method = "ML")

Predict_Datos <- read_excel("BD/Predict_Base_Transporte_ARIMA.xlsx", 
                            sheet = "Datos", col_names = TRUE)

D_Mar2020_f <- ts(Predict_Datos$D_Mar2020, 
                start = ini_f, 
                freq = 12)

D_Abr2020_f <- ts(Predict_Datos$D_Abr2020, 
                start = ini_f, 
                freq = 12)

D_Jun2020_f <- ts(Predict_Datos$D_Jun2020, 
                start = ini_f, 
                freq = 12)

D_Jul2020_f <- ts(Predict_Datos$D_Jul2020, 
                start = ini_f, 
                freq = 12)

D_Mar2021_f <- ts(Predict_Datos$D_Mar2021, 
                start = ini_f, 
                freq = 12)

D_Ene_f <- ts(Predict_Datos$D_Ene, 
            start = ini_f, 
            freq = 12)

D_Feb_f <- ts(Predict_Datos$D_Feb, 
            start = ini_f, 
            freq = 12)

D_Jul_f <- ts(Predict_Datos$D_Jul, 
           start = ini_f, 
            freq = 12)

D_Dic_f <- ts(Predict_Datos$D_Dic, 
            start = ini_f, 
            freq = 12)

DLPax_Nal_f <- predict(ARMA_Ex_DLPax_Nal_2, n.ahead = h_f, 
                      newxreg = cbind(D_Ene_f, D_Feb_f, D_Jul_f, 
                                      D_Dic_f, D_Mar2020_f, 
                                      D_Abr2020_f, D_Jun2020_f, 
                                      D_Jul2020_f, D_Mar2021_f))

Pronostico_Arima <- read_excel("BD/Pronostico_ARIMA.xlsx",
                               sheet = "Datos", col_names = TRUE)

Pronostico_Arima$Pax_Nal_f <- Pronostico_Arima$Pax_Nal

# Último dato observado: el archivo trae NA en los 24 meses a pronosticar,
# por lo que el índice se calcula (y no se fija a mano) para que el código
# siga siendo válido cuando se actualice la base.
n_obs <- sum(!is.na(Pronostico_Arima$Pax_Nal))

# El pronóstico está en diferencias logarítmicas, de modo que el nivel se
# reconstruye con el factor exp(DL) --exacto-- y no con (1 + DL), que es
# sólo la aproximación de primer orden.
for(i in 1:h_f){
  Pronostico_Arima$Pax_Nal_f[n_obs + i] <-
    Pronostico_Arima$Pax_Nal_f[n_obs + i - 1]*exp(DLPax_Nal_f$pred[i])
}

options(scipen = 999) # NO notacion cientifica

ggplot(data = Pronostico_Arima, aes(x = Periodo)) +
  geom_line(aes(y = Pax_Nal_f, color = "Pronóstico (24 meses)")) +
  geom_line(aes(y = Pax_Nal, color = "Observado")) +
  scale_color_manual(values = c("Observado" = "darkblue",
                                "Pronóstico (24 meses)" = "darkred")) +
  #theme_bw() + 
  theme(legend.position = "bottom") +
  theme(legend.title = element_blank()) +
  guides(col = guide_legend(nrow = 1, byrow = TRUE)) + 
  xlab("Tiempo") + 
  ylab("Pasajeros") + 
  theme(plot.title = element_text(size = 11, face = "bold", 
                                  hjust = 0)) + 
  theme(plot.subtitle = element_text(size = 10, hjust = 0)) + 
  theme(plot.caption = element_text(size = 10, hjust = 0)) +
  theme(plot.margin = unit(c(1,1,1,1), "cm")) +
  labs(
    title = "Pasajeros en vuelos nacionales (Salidas)",
    subtitle = paste0("(", per_txt, " y pronóstico a ", h_f, " meses)"),
    caption = "Fuente: Elaboración propia"
  )

Figura 3.20: Pronóstico de la serie $PaxNal_t$ a partir de un ARMA(6,6) en diferencias logarítmicas

3.10 Desestacionalización y filtrado de Series

3.10.1 Motivación

Existen múltiples enfoques para la desestacionalización de series. Algunos modelos, por ejemplo, pueden estar basados en modelos ARIMA como un conjunto de dummies. No obstante, el caso particular que discutiremos en este libro estará basado en un modelo ARIMA de la serie. Este enfoque está basado en el modelo X11 de la oficina del censo de Estados Unidos (Census Bureau) el cual es conocido como el modelo X13-ARIMA-SEATS.⁸ El modelo X13-ARIMA-SEATS es, como su nombre lo indica, la combinación de un modelo ARIMA con componentes estacionales (por la traducción literal de la palabra: seasonal).

Un modelo ARIMA estacional emplea la serie en diferencias y como regresores los valores rezagados de las diferencias de la serie tantas veces como procesos estacionales $s$ existan en ésta, con el objeto de remover los efectos aditivos de la estacionalidad. Sólo para entender qué significa este mecanismo, recordemos que cuando se utiliza la primera diferencia de la serie respecto del periodo inmediato anterior se remueve la tendencia. Por su parte, cuando se incluye la diferencia respecto del mes $s$ se está en el caso en que se modela la serie como una media móvil en términos del rezago $s$.

El modelo ARIMA estacional incluye como componentes autorregresivos y de medias móviles a los valores rezagados de la serie en el periodo $s$ en diferencias. El ARIMA(p, d, q)(P, D, Q) estacional puede ser expresado de la siguiente manera utilizando el operador rezago $L$: \[\begin{equation} \Theta_P(L^s) \theta_p(L) (1 - L^s)^D (1 - L)^d X_t = \Psi_Q(L^s) \psi_q(L) U_t \tag{3.59} \end{equation}\]

Donde $\Theta_P(.)$, $\theta_p(.)$, $\Psi_Q(.)$ y $\psi_q(.)$ son polinomios de $L$ de orden $P$, $p$, $Q$ y $q$ respectivamente. En general, la representación es de una serie no estacionaria; no obstante, si $D = d = 0$ y todas las raíces de los polinomios del lado izquierdo de la ecuación (3.59) son mayores que 1 en valor absoluto, el proceso modelado será estacionario.

Si bien es cierto que existen otras formas de modelar la desestacionalización, como la modelación en diferencias con dummies para identificar ciertos patrones regulares, en los algoritmos disponibles como el X11 o X13-ARIMA-SEATS se emplea la formulación de la ecuación (3.59). A continuación, implementaremos la desestacionalización de una serie.

Como ejemplo utilizaremos la serie del Índice Nacional de Precios al Consumidor (INPC).⁹ Podemos ver que la serie original del INPC y su ajuste estacional bajo una metodología X13-ARIMA-SEATS son como se muestra en la Figura 3.21.

library(ggplot2)
library(dplyr)
library(readxl)
library(stats)
library(seasonal)

## 
## Attaching package: 'seasonal'

## The following object is masked from 'package:tibble':
## 
##     view

Datos <- read_excel("BD/Base_VAR.xlsx", sheet = "Datos", col_names = TRUE)

INPC <- ts(Datos$INPC, 
           start = c(2000, 1), 
           freq = 12)

Seas_INPC <- seas(INPC)

INPC_Ad <- final(Seas_INPC)

plot(Seas_INPC)

Figura 3.21: Índice Nacional de Precios al Consumidor ($INPC_t$) y su serie desestacionalizada utilizando un proceso X13-ARIMA-SEATS

TC <- ts(Datos$TC, 
         start = c(2000, 1),
         freq = 12)

Seas_TC <- seas(TC)

TC_Ad <- final(Seas_TC)

CETE28 <- ts(Datos$CETE28, 
             start = c(2000, 1), 
             freq = 12)

Seas_CETE28 <- seas(CETE28)

CETE28_Ad <- final(Seas_CETE28)

IGAE <- ts(Datos$IGAE, 
           start = c(2000, 1), 
           freq = 12)

Seas_IGAE <- seas(IGAE)

IGAE_Ad <- final(Seas_IGAE)

IPI <- ts(Datos$IPI, 
          start = c(2000, 1), 
          freq = 12)

Seas_IPI <- seas(IPI)

IPI_Ad <- final(Seas_IPI)

Datos_Ad <- data.frame(cbind(INPC_Ad, TC_Ad, CETE28_Ad, IGAE_Ad, IPI_Ad))

Datos_Ad <- cbind(Datos, Datos_Ad)

# Se guarda con el nombre que consumen los capítulos 4 y 5, de modo que al
# actualizar BD/Base_VAR.xlsx la desestacionalización se propague sola y no
# queden dos versiones de la misma base.
save(Datos_Ad, file = "BD/Datos_Ad.RData")

El mismo procesamiento puede ser seguido para todas las series que busquemos analizar. En particular, en adelante, además del INPC que incluimos en la lista, utilizaremos las siguientes series, así como su versión desestacionalizada:

Índice Nacional de Precios al Consumidor (base 2QJul2018 = 100), $INPC_t$.
Tipo de Cambio FIX, $TC_t$
Tasa de rendimiento promedio mensual de los Cetes 28, en por ciento anual, $CETE28_t$
Indicador global de la actividad económica (base 2013 = 100), $IGAE_t$
Industrial Production Index o Índice de Producción Industrial de los Estados Unidos (base 2012 = 100), $IPI_t$

3.10.2 Filtro Hodrick-Prescott

Como último tema de los procesos univariados y que no necesariamente aplican a series estacionarias, a continuación desarrollaremos el procedimiento conocido como filtro de Hodrick y Prescott (1997). El trabajo de estos autores era determinar una técnica de regresión que permitiera utilizar series agregadas o macroeconómicas para separarlas en dos componentes: uno de ciclo de negocios y otro de tendencia. En su trabajo original, Hodrick y Prescott (1997) utilizaron datos trimestrales de algunas series como el Producto Nacional Bruto (GNP, por sus siglas en inglés), los agregados monetarios M1, empleo, etc., de los Estados Unidos que fueron observados posteriormente a la Segunda Guerra Mundial.

El marco conceptual de Hodrick y Prescott (1997) parte de suponer una serie $X_t$ que se puede descomponer en la suma de componente de crecimiento tendencial, $g_t$, y su componente de ciclo de negocios, $c_t$, de esta forma para $t = 1, 2, \ldots, T$ tenemos que: \[\begin{equation} X_t = g_t + c_t \tag{3.60} \end{equation}\]

En la ecuación (3.60) se asume que la suavidad de la ruta seguida por $g_t$ se controla penalizando la suma de los cuadrados de su segunda diferencia. En esa misma ecuación asumiremos que $c_t$ son las desviaciones de $g_t$, las cuales en el largo plazo tienen una media igual a cero (0). Por esta razón, se suele decir que el filtro de Hodrick y Prescott representa una descomposición de la serie en su componente de crecimiento natural y de sus desviaciones transitorias que en promedio son cero, en el largo plazo.

Estas consideraciones que hemos mencionado señalan que el procedimiento de Hodrick y Prescott (1997) implica resolver el siguiente problema de minimización para determinar cada uno de los componentes en que $X_t$ se puede descomponer:

\[\begin{equation} \min_{\{ g_t \}^T_{t = -1} } \left[ \sum^T_{t = 1} c^2_t + \lambda \sum^T_{t = 1} [ \Delta g_t - \Delta g_{t-1}]^2 \right] \end{equation}\]

Donde $\Delta g_t = g_t - g_{t-1}$ y $\Delta g_{t-1} = g_{t-1} - g_{t-2}$; $c_t = X_t - g_t$, y el parámetro $\lambda$ es un número positivo que penaliza la variabilidad en el crecimiento de las series. El valor de $\lambda$ debe fijarse de acuerdo con la periodicidad de las series. Los valores de uso más extendido en el trabajo aplicado son:

100 si la serie es de datos anuales
1,600 si la serie es de datos trimestrales
14,400 si la serie es de datos mensuales

Conviene precisar que el único valor propuesto por Hodrick y Prescott (1997) fue $\lambda = 1{,}600$ para datos trimestrales de la economía estadounidense; los valores para datos anuales y mensuales son convenciones que se obtienen al escalar $\lambda$ con el cuadrado del cambio de frecuencia. Ravn y Uhlig (2002) mostraron que el ajuste correcto es con la cuarta potencia del cambio de frecuencia, lo que implica $\lambda \approx 6.25$ para datos anuales y $\lambda = 129{,}600$ para datos mensuales. En este libro usamos $\lambda = 14{,}400$ por ser el valor que reportan por defecto los paquetes estadísticos más usados, pero el lector debe tener presente que la elección de $\lambda$ no es inocua: un valor mayor produce una tendencia más suave y, por lo tanto, un ciclo de mayor amplitud.

En resumen, podemos decir que el filtro de Hodrick y Prescott (1997) es un algoritmo que minimiza las distancias o variaciones de la trayectoria de largo plazo. De esta forma, determina una trayectoria estable de largo plazo, por lo que las desviaciones respecto de esta trayectoria serán componentes de ciclos de negocio o cambios transitorios (tanto positivos como negativos).

A continuación, ilustraremos el filtro de Hodrick y Prescott (1997) para dos series desestacionalizadas: $INPC_t$ y $TC_t$. Las Figura 3.22 y Figura 3.23 muestran los resultados de la implementación del filtro.

library(ggplot2)
library(dplyr)
library(readxl)
library(stats)
library(mFilter)
library(plm)

## 
## Attaching package: 'plm'

## The following objects are masked from 'package:dplyr':
## 
##     between, lag, lead

load("BD/Datos_Ad.RData")

INPC <- ts(Datos_Ad$INPC_Ad, start = c(2000, 1), freq = 12)

TC <- ts(Datos_Ad$TC_Ad, start = c(2000, 1), freq = 12)

INPC_hpf <- hpfilter(INPC, freq = 14400)

# ask = FALSE dibuja los dos paneles de corrido; sin él, el método de
# mFilter pide un Enter entre gráficas cuando se ejecuta en la consola.
plot(INPC_hpf, ask = FALSE)

Descomposición del Índice Nacional de Precios al Consumidor ($INPC_t$) en su tendencia o trayectoria de largo plazo y su ciclo de corto plazo utilizando un filtro Hodrick y Prescott [-@hodrickprescott1997]

Figura 3.22: Descomposición del Índice Nacional de Precios al Consumidor ($INPC_t$) en su tendencia o trayectoria de largo plazo y su ciclo de corto plazo utilizando un filtro Hodrick y Prescott (1997)

TC_hpf <- hpfilter(TC, freq = 14400)

plot(TC_hpf, ask = FALSE)

Descomposición del Tipo de Cambio FIX ($TC_t$) en su tendencia o trayectoria de largo plazo y su ciclo de corto plazo utilizando un filtro Hodrick y Prescott [-@hodrickprescott1997]

Figura 3.23: Descomposición del Tipo de Cambio FIX ($TC_t$) en su tendencia o trayectoria de largo plazo y su ciclo de corto plazo utilizando un filtro Hodrick y Prescott (1997)

3.10.2.1 Filtro Hodrick-Prescott planteado por St-Amant y van Norden

Existe una técnica adicional que atiende una limitación conocida del filtro HP: el llamado problema de fin de muestra o end-point problem (a veces referido como el problema de las “colas” de la muestra). Al ser un filtro de dos colas, la estimación de la tendencia en las primeras y últimas observaciones utiliza menos información que en el centro de la muestra, por lo que es la parte de la tendencia –y por lo tanto de la brecha– que más se revisa cuando llegan datos nuevos. Justamente el final de la muestra es el que interesa para la política económica. La propuesta de St-Amant y van Norden (1997) consiste en anclar el crecimiento de la tendencia al final de la muestra a una tasa de crecimiento de largo plazo. Analicemos ese caso. El método tradicional de HP consiste en elegir la serie $\{ \tau_t \}_{t=1}^T$ que minimiza:

\[\begin{equation} \sum_{t=1}^T (y_t - \tau_t)^2 + \lambda \sum_{t=2}^{T-1} [(\tau_{t+1} - \tau_{t}) - (\tau_{t} - \tau_{t-1})]^2 \end{equation}\]

Donde $\lambda$ es un parámetro fijo (determinado ex-ante) y $\tau_t$ es un componente de tendencia de $y_t$. Nótese que la penalización sólo puede evaluarse en las $T-2$ segundas diferencias que la muestra permite calcular; en particular, no se imponen valores iniciales a $\tau_0$ ni a $\tau_{-1}$.

Reuniendo la serie observada y la tendencia en los vectores $Y$ y $G$, ambos de dimensión $T \times 1$, la forma matricial del filtro HP es: \[\begin{equation} (Y - G)'(Y - G) + \lambda G' K' K G \end{equation}\]

La condición de primer orden respecto de $G$ es:

\[\begin{equation} -2 Y + 2 G + 2 \lambda K' K G = 0 \end{equation}\]

Despejando:

\[\begin{equation} G_{hp} = [I_T + \lambda K' K]^{-1} Y \end{equation}\]

Donde $G$ es el vector de tendencia, $Y$ es el vector de la serie de datos, $\lambda$ es la constante tradicional, e $I_T$ es la matriz identidad de orden $T$. La matriz $K$ es la de segundas diferencias, de dimensión $(T-2) \times T$:

\[\begin{equation} K = \begin{pmatrix} 1 & -2 & 1 & 0 & \ldots & 0 & 0 \\ 0 & 1 & -2 & 1 & \ldots & 0 & 0 \\ 0 & 0 & 1 & -2 & \ldots & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & \ldots & -2 & 1 \\ \end{pmatrix} \end{equation}\]

De esta forma, el renglón $t$-ésimo del producto $K G$ es la segunda diferencia $\tau_{t+2} - 2 \tau_{t+1} + \tau_t$, y la forma cuadrática $G' K' K G$ reproduce exactamente la suma de segundas diferencias al cuadrado de la función objetivo. Conviene subrayar que $K$ no es cuadrada: escribirla como una matriz $T \times T$ equivale a fijar $\tau_0 = \tau_{-1} = 0$, con lo que la penalización actúa sobre el nivel de las dos primeras observaciones –y no sobre su curvatura– y la tendencia estimada al inicio de la muestra se colapsa hacia cero.

El método modificado de HP consiste en elegir la serie $\{ \tau_t \}_{t=1}^T$ que minimiza:

\[\begin{equation} \sum_{t=1}^T (y_t - \tau_t)^2 + \lambda \sum_{t=2}^{T-1} [(\tau_{t+1} - \tau_t) - (\tau_t - \tau_{t-1})]^2 + \lambda_{ss} \sum_{t=T-j}^{T} [\Delta \tau_t - u_{ss}]^2 \tag{3.61} \end{equation}\]

Donde $\tau_t$ es el componente de tendencia de $y_t$ y $\lambda$ es el parámetro de suavizamiento de siempre. Los dos primeros sumandos son los del filtro HP tradicional; el tercero es la novedad, y es un castigo por las desviaciones del crecimiento de la tendencia respecto de una tasa de largo plazo $u_{ss}$, aplicado únicamente a los últimos $j+1$ periodos y graduado por el parámetro $\lambda_{ss} \geq 0$. Cuando $\lambda_{ss} = 0$ el tercer término desaparece y se recupera el filtro HP tradicional.

La lógica de esta modificación es la del propio problema que atiende: como el filtro HP es un procedimiento de dos colas, la tendencia del final de la muestra se estima con menos información que la del centro y por eso se revisa tanto cuando llegan datos nuevos. El castigo sustituye la información que falta por un supuesto explícito sobre el crecimiento de largo plazo.

Para escribir el problema (3.61) en forma matricial, definamos la matriz $D_j$ de dimensión $(j+1) \times T$ que extrae las primeras diferencias de la tendencia de los últimos $j+1$ periodos, es decir, $D_j G = (\Delta \tau_{T-j}, \Delta \tau_{T-j+1}, \ldots, \Delta \tau_T)'$, y sea $\boldsymbol{\iota}$ un vector de unos de dimensión $(j+1)$. Con esta notación, la función objetivo del filtro HP-SAVN es: \[\begin{equation} (Y - G)'(Y - G) + \lambda G' K' K G + \lambda_{ss} (D_j G - u_{ss} \boldsymbol{\iota})'(D_j G - u_{ss} \boldsymbol{\iota}) \end{equation}\]

La condición de primer orden respecto de $G$ es:

\[\begin{equation} -2 Y + 2 G + 2 \lambda K' K G + 2 \lambda_{ss} D_j' (D_j G - u_{ss} \boldsymbol{\iota}) = 0 \end{equation}\]

Despejando obtenemos la tendencia del filtro modificado: \[\begin{equation} G_{SAVN} = \left[ I_T + \lambda K' K + \lambda_{ss} D_j' D_j \right]^{-1} \left( Y + \lambda_{ss} u_{ss} D_j' \boldsymbol{\iota} \right) \tag{3.62} \end{equation}\]

La ecuación (3.62) deja ver con claridad el papel de cada componente: cuando $\lambda_{ss} = 0$ se recupera exactamente el filtro HP tradicional, $G_{hp} = [I_T + \lambda K' K]^{-1} Y$; y a medida que $\lambda_{ss}$ crece, la tendencia de los últimos $j+1$ periodos se fuerza a crecer a la tasa $u_{ss}$, lo que estabiliza la estimación del final de la muestra.

Antes de llevar la ecuación (3.62) a los datos hay que fijar sus tres parámetros –$\lambda$, $u_{ss}$ y $\lambda_{ss}$–, y conviene tratarlos por separado porque sólo el primero admite un criterio sistemático.

El parámetro de suavizamiento $\lambda$. Marcet y Ravn (2004) observan que los valores convencionales de $\lambda$ se calibraron con datos de Estados Unidos y que no hay razón para trasladarlos sin más a economías cuyos ciclos tienen otra amplitud. Su propuesta –la regla 1– consiste en resolver el problema del filtro sujeto a una restricción explícita:

\[\begin{equation} \min_{\{\tau_t\}} \sum_{t=1}^T (y_t - \tau_t)^2 \quad \text{sujeto a} \quad \frac{\sum_{t=2}^{T-1} (\Delta^2 \tau_{t+1})^2}{\sum_{t=1}^T (y_t - \tau_t)^2} \leq V \tag{3.63} \end{equation}\]

donde $\Delta^2 \tau_{t+1} = (\tau_{t+1} - \tau_t) - (\tau_t - \tau_{t-1})$ es la segunda diferencia de la tendencia y $V$ mide la variabilidad de la aceleración de la tendencia relativa a la variabilidad del ciclo. Si $\gamma$ es el multiplicador de Lagrange de esa restricción, el problema restringido resulta equivalente al filtro HP cuando

\[\begin{equation} \lambda = \frac{\gamma}{1 - \gamma V} \tag{3.64} \end{equation}\]

de modo que fijar $\lambda$ o fijar $V$ son dos maneras de plantear el mismo problema. La ventaja de razonar en términos de $V$ es que se trata de una cantidad interpretable y comparable entre series, mientras que $\lambda$ no lo es: el valor convencional de 1,600 puede leerse entonces como el $\lambda$ que corresponde al $V$ de la economía estadounidense.

El procedimiento tiene dos pasos. Primero se calcula $V$ sobre una serie de referencia con el $\lambda$ convencional; Marcet y Ravn (2004) usan el PIB trimestral de Estados Unidos con $\lambda = 1600$ y obtienen $V \approx 1.81 \times 10^{-4}$. Después, para la serie de interés se define

\[\begin{equation} F(\lambda) = \frac{\sum_{t=2}^{T-1} [(\tau_{t+1}(\lambda) - \tau_t(\lambda)) - (\tau_t(\lambda) - \tau_{t-1}(\lambda))]^2}{\sum_{t=1}^T [y_t - \tau_t(\lambda)]^2} = \frac{G' K' K G}{(Y - G)'(Y - G)} \tag{3.65} \end{equation}\]

donde $\tau_t(\lambda)$ –o $G = G_{hp}(\lambda)$ en forma matricial– es la tendencia que produce el filtro HP con ese valor de $\lambda$, y se resuelve numéricamente $F(\lambda) = V$. Como $F$ es decreciente en $\lambda$ –más suavizamiento significa una tendencia más plana y, por construcción, un ciclo más amplio–, la solución es única y se obtiene con cualquier método de bisección.

Nótese que la condición no es circular: el valor objetivo $V$ se calcula sobre la serie de referencia, y $F$ se evalúa sobre la serie que se quiere filtrar. Nótese también que $V$ depende de la frecuencia de los datos, de modo que ambas series deben estar en la misma.

Al aplicar esta regla al PIB trimestral de México, Antón Sarabia (2010) obtiene $\lambda = 1096$ en lugar de 1,600: el ciclo mexicano es más volátil que el estadounidense, así que la misma restricción sobre $V$ admite una tendencia menos suave.

La tasa de estado estacionario $u_{ss}$. No se estima: la fija el investigador a partir de lo que sabe de la serie. Antón (2010) adopta $u_{ss} = 0$ para las series que en teoría son estacionarias y, para las demás, la tasa de crecimiento promedio de la serie durante el periodo de referencia.

El castigo $\lambda_{ss}$. Aquí conviene ser explícito: no existe una referencia a priori sobre el valor de $\lambda_{ss}$. Antón (2010) resuelve el punto adoptando $\lambda_{ss} = \lambda = 1096$, es decir, castigando las desviaciones respecto del crecimiento de largo plazo con el mismo peso con que se castiga la curvatura de la tendencia. Es una convención defendible –y la que se usa en el ejemplo de esta sección– pero es una convención, no un resultado. Lo que sí corresponde al usar el método es reportar el valor empleado y mostrar qué tan sensible es la brecha estimada a esa elección.

A diferencia del filtro tradicional, la variante de St-Amant y van Norden no está implementada en mFilter ni en los demás paquetes que usamos en este libro, de modo que hay que programarla. No es un obstáculo: la ecuación (3.62) se traduce de forma directa a código, porque basta construir $K$ y $D_j$, armar la matriz del sistema y resolverlo. La función HP_tail() del archivo HP_tail.R hace exactamente eso y recibe la serie, los dos parámetros de suavizamiento, la tasa de crecimiento de estado estacionario g_ss –declarada en por ciento anual– y el número de periodos anclados. Como el filtro trabaja sobre la serie en logaritmos, las diferencias que penaliza $D_j$ son log-diferencias, de modo que la tasa por periodo es la logarítmica, $u_{ss} = \ln(1 + g_{ss})/s$, con $s$ la frecuencia de los datos.

Como serie de trabajo usaremos el IGAE desestacionalizado, una de las cinco que construimos con X-13 unas páginas atrás y que forma parte de Datos_Ad. La elección no es indiferente: el filtro modificado se diseñó para estimar la brecha del producto –es el problema que motiva a St-Amant y van Norden (1997) y a Antón (2010)–, de modo que un indicador de actividad económica es el caso natural, mientras que hablar de la “brecha” del INPC o del tipo de cambio tendría poco sentido económico. Es además la misma serie sobre la que el Capítulo 5 vuelve a estimar una brecha, ahí con el filtro de Kalman.

Conviene empezar por comprobar que la implementación es correcta, y para eso sirve la propiedad que ya señalamos: con $\lambda_{ss} = 0$ la ecuación (3.62) colapsa al filtro HP tradicional, así que HP_tail() debe reproducir lo que devuelve hpfilter() sobre la misma serie.

source("HP_tail.R")

IGAE_ds <- ts(Datos_Ad$IGAE_Ad, start = c(2000, 1), freq = 12)

# Con lambda_ss = 0 no hay ancla de fin de muestra: debe coincidir con mFilter.
propia    <- HP_tail(IGAE_ds, lambda = 14400, lambda_ss = 0)$tendencia
referencia <- mFilter::hpfilter(log(IGAE_ds), freq = 14400)$trend

max(abs(propia - referencia))

## [1] 0.000000000009663381

La diferencia es del orden de $10^{-12}$: error de redondeo. Con eso verificado, lo que sigue es elegir los parámetros y activar el ancla.

El mismo archivo incluye calibrar_lambda(), que resuelve $F(\lambda) = V$ de la regla 1 de Marcet y Ravn (2004). Como $V$ depende de la frecuencia, el ejemplo agrega el IGAE mensual a frecuencia trimestral, que es la de la serie de referencia:

# El IGAE mensual se agrega a trimestres (promedio) para que la comparacion
# con el valor de referencia de EUA sea a la misma frecuencia.
n_trim    <- floor(nrow(Datos_Ad) / 3) * 3
IGAE_trim <- ts(colMeans(matrix(Datos_Ad$IGAE_Ad[1:n_trim], nrow = 3)),
                start = c(2000, 1), freq = 4)

lambda_mx <- calibrar_lambda(IGAE_trim, V = 1.81e-4)

knitr::kable(
  data.frame(
    Concepto = c("$F$ del IGAE con el $\\lambda$ convencional (1,600)",
                 "$V$ de referencia (PIB de EUA, $\\lambda = 1{,}600$)",
                 "$\\lambda$ calibrado para el IGAE"),
    Valor    = c(sprintf("%.2e", F_MR(IGAE_trim, 1600)),
                 sprintf("%.2e", 1.81e-4),
                 sprintf("%.1f", lambda_mx))),
  escape = FALSE, booktabs = TRUE, row.names = FALSE,
  caption = "Calibración del parámetro de suavizamiento del IGAE trimestral por la regla 1 de Marcet y Ravn")

Cuadro 3.12: Calibración del parámetro de suavizamiento del IGAE trimestral por la regla 1 de Marcet y Ravn
Concepto	Valor
$F$ del IGAE con el $\lambda$ convencional (1,600)	8.77e-05
$V$ de referencia (PIB de EUA, $\lambda = 1{,}600$)	1.81e-04
$\lambda$ calibrado para el IGAE	920.5

El resultado, $\lambda \approx 920$, queda por debajo del convencional 1,600 y es del mismo orden que el 1,096 que reporta Antón (2010) con el PIB y con un periodo distinto. La lectura económica coincide en ambos casos: el ciclo mexicano es más volátil que el estadounidense, de modo que imponer la misma restricción sobre $V$ admite una tendencia menos suave.

Volvamos ahora al IGAE mensual. Conservamos el $\lambda = 14{,}400$ convencional –el $\lambda$ calibrado del cuadro anterior corresponde a la frecuencia trimestral y no puede trasladarse sin más a datos mensuales– y seguimos la convención de Antón (2010) de fijar $\lambda_{ss} = \lambda$. El ancla se pone en 2% anual, del orden del crecimiento promedio del IGAE en el periodo, y se aplica a los últimos $j+1 = 8$ meses. La Figura 3.24 compara la brecha que resulta de ambos filtros.

# HP es el filtro tradicional (lambda_ss = 0, ya verificado arriba); SAVN
# activa el ancla de fin de muestra con lambda_ss = lambda.
HP   <- HP_tail(IGAE_ds, lambda = 14400, lambda_ss = 0)
SAVN <- HP_tail(IGAE_ds, lambda = 14400, lambda_ss = 14400, g_ss = 2,
                j = 7, frecuencia = 12)

plot(100 * HP$ciclo, type = "l", col = "black", lwd = 1,
     xlab = "Tiempo", ylab = "Brecha (%)")
lines(100 * SAVN$ciclo, col = "red", lwd = 1)
abline(h = 0, lty = 3)
legend("bottomleft", c("HP tradicional", "HP-SAVN"),
       lty = 1, col = c("black", "red"), cex = 0.8, bty = "n")

$Brecha del IGAE desestacionalizado estimada con el filtro Hodrick-Prescott tradicional y con la variante de St-Amant y van Norden [-@stamantvannorden1997], que ancla el crecimiento de la tendencia de los últimos ocho meses a una tasa de largo plazo de 2\% anual$

Figura 3.24: Brecha del IGAE desestacionalizado estimada con el filtro Hodrick-Prescott tradicional y con la variante de St-Amant y van Norden (1997), que ancla el crecimiento de la tendencia de los últimos ocho meses a una tasa de largo plazo de 2% anual

El Cuadro 3.13 muestra el efecto sobre el tramo que interesa. La tendencia del filtro tradicional se ajusta a las últimas observaciones y absorbe parte de lo que debería ser ciclo; la del filtro modificado crece a una tasa que converge al ancla 0.165% mensual, de modo que la brecha final deja de estar mecánicamente cerca de cero.

n_igae <- length(IGAE_ds)
ult    <- (n_igae - 5):n_igae

# Los rotulos se derivan del propio objeto ts, para que sigan siendo correctos
# cuando se actualice la base.
fechas_igae <- as.Date(paste(floor(time(IGAE_ds)),
                             round(12 * (time(IGAE_ds) %% 1)) + 1, 1, sep = "-"))

comparacion_hp <- data.frame(
  Periodo          = format(fechas_igae[ult], "%b %Y"),
  `Brecha HP`      = round(100 * as.numeric(HP$ciclo)[ult], 3),
  `Brecha HP-SAVN` = round(100 * as.numeric(SAVN$ciclo)[ult], 3),
  check.names = FALSE)

knitr::kable(comparacion_hp, row.names = FALSE, booktabs = TRUE,
             caption = "Brecha del IGAE (\\%) en las últimas seis observaciones de la muestra: filtro HP tradicional frente a la variante de St-Amant y van Norden")

Cuadro 3.13: Brecha del IGAE (%) en las últimas seis observaciones de la muestra: filtro HP tradicional frente a la variante de St-Amant y van Norden
Periodo	Brecha HP	Brecha HP-SAVN
dic 2025	0.069	-0.440
ene 2026	-0.862	-1.456
feb 2026	-0.760	-1.440
mar 2026	-0.271	-1.036
abr 2026	0.679	-0.171
may 2026	0.400	-0.536

Conviene cerrar con una advertencia: el resultado depende de los tres parámetros que el método añade, y ninguno se estima de los datos. La elección de $j$ fija cuántos periodos se anclan, $u_{ss}$ traslada al filtro un supuesto sobre el crecimiento potencial, y $\lambda_{ss}$ gradúa la fuerza del ancla. El filtro modificado no elimina la incertidumbre del final de la muestra: la sustituye por un supuesto explícito, que es justamente su ventaja frente al filtro tradicional, donde ese supuesto queda implícito en la aritmética del suavizamiento.

3.11 Resumen del capítulo

En este capítulo estudiamos los fundamentos de los procesos estocásticos estacionarios y los principales modelos univariados de series de tiempo.

Estacionariedad y ergodicidad. Un proceso es débilmente estacionario si su media, varianza y estructura de covarianzas no dependen del tiempo. La ergodicidad garantiza que los promedios temporales convergen a los promedios del conjunto. El proceso de ruido blanco, $U_t \sim (0, \sigma^2)$, es el bloque elemental a partir del cual se construyen todos los modelos de este capítulo.

Herramientas de diagnóstico. La función de autocorrelación (ACF) y la función de autocorrelación parcial (PACF) son el puente entre la teoría y la práctica: permiten identificar el orden de los componentes AR y MA, respectivamente. Las pruebas Box-Pierce ($Q^*$), Ljung-Box ($Q$) y Breusch-Godfrey (LM) contrastan formalmente si los residuales de un modelo presentan autocorrelación; la prueba de Jarque-Bera evalúa normalidad.

Modelos AR, MA y ARMA. Los procesos autorregresivos $AR(p)$ modelan $X_t$ en función de sus propios rezagos; son estacionarios si todas las raíces $\lambda_i$ del polinomio característico $\lambda^p - a_1 \lambda^{p-1} - \ldots - a_p = 0$ caen dentro del círculo unitario o, equivalentemente, si las raíces del polinomio en el operador rezago $\alpha(L) = 0$ caen fuera de él. Los procesos de medias móviles $MA(q)$ expresan $X_t$ como combinación lineal de errores pasados; son siempre estacionarios y son invertibles bajo la condición análoga sobre el polinomio $\beta(L)$. El modelo $ARMA(p,q)$ combina ambos componentes y, bajo diferenciación $d$ veces, da lugar al $ARIMA(p,d,q)$.

Selección de modelos. Los criterios AIC y SC –este último también llamado BIC– equilibran el ajuste del modelo con el principio de parsimonia. El patrón de la ACF y la PACF —corte abrupto vs. decaimiento gradual— orienta la selección inicial del orden $(p,q)$.

Pronósticos. En el horizonte $h$, el pronóstico óptimo (mínimo error cuadrático medio) es la esperanza condicional $\mathbb{E}[X_{t+h} | \mathcal{F}_t]$. El intervalo de predicción se amplía con $h$ porque la incertidumbre se acumula.

Desestacionalización y filtrado. El filtro Hodrick-Prescott separa la tendencia de largo plazo ($\tau_t$) del ciclo de corto plazo; la elección del parámetro $\lambda$ depende de la frecuencia de los datos ($\lambda = 1600$ para datos trimestrales, $\lambda = 14400$ para datos mensuales). La variante St-Amant y van Norden atenúa el problema de fin de muestra al anclar el crecimiento de la tendencia en las últimas observaciones a una tasa de largo plazo.

El siguiente capítulo extiende este marco a series no estacionarias y presenta las pruebas de raíz unitaria que determinan si es necesario diferenciar una serie antes de modelarla con las herramientas aquí desarrolladas.

3.12 Ejercicios

Los siguientes ejercicios combinan preguntas teóricas con aplicaciones en R. Los marcados como (Teórico) se resuelven analíticamente; los marcados como (Computacional) utilizan los datos y las funciones empleadas en este capítulo.

(Teórico) Considere el proceso AR(1) dado por $X_t = a_0 + a_1 X_{t-1} + U_t$, donde $|a_1| < 1$ y $U_t$ es un proceso puramente aleatorio con media cero y varianza $\sigma^2$.
1. Demuestre que la media del proceso es $\mu = a_0 / (1 - a_1)$.
2. Demuestre que la varianza es $\gamma(0) = \sigma^2 / (1 - a_1^2)$.
3. Demuestre que la función de autocorrelación es $\rho(\tau) = a_1^{\tau}$, $\tau = 1, 2, \ldots$, y compare la velocidad de decaimiento para $a_1 = 0.5$ y $a_1 = 0.9$ calculando los primeros cinco valores en cada caso.
4. ¿Qué ocurre con $\gamma(0)$ y con $\rho(\tau)$ cuando $a_1 \to 1$? Relacione su respuesta con la discusión de procesos no estacionarios del Capítulo 4.
(Teórico) Considere el proceso MA(2) dado por $X_t = \mu + U_t - b_1 U_{t-1} - b_2 U_{t-2}$, con $U_t$ puramente aleatorio con media cero y varianza $\sigma^2$.
1. Calcule la media y la varianza $\gamma(0)$ del proceso.
2. Demuestre que $\gamma(1) = (-b_1 + b_1 b_2) \sigma^2$, que $\gamma(2) = -b_2 \sigma^2$ y que $\gamma(\tau) = 0$ para todo $\tau > 2$.
3. Con base en el inciso anterior, ¿qué patrón debe exhibir la función de autocorrelación muestral de un MA(2)? ¿En qué difiere del patrón de un AR(2)?
4. Explique por qué el proceso es estacionario sin necesidad de imponer restricciones sobre $b_1$ y $b_2$, mientras que la invertibilidad sí las requiere. Enuncie la condición de invertibilidad para el caso MA(1).
(Teórico) Considere el proceso ARMA(1,1) dado por $X_t = a_1 X_{t-1} + U_t - b_1 U_{t-1}$, estacionario e invertible.
1. Demuestre que para $\tau \geq 2$ las autocovarianzas satisfacen la recursión $\gamma(\tau) = a_1 \gamma(\tau - 1)$.
2. Concluya que la ACF decae geométricamente a partir de $\tau = 1$, de modo que ni la ACF ni la PACF se cortan de manera abrupta.
3. Con base en lo anterior y en el Cuadro 3.9, explique cómo distinguiría en la práctica un ARMA(1,1) de un AR(1) puro y de un MA(1) puro.
(Computacional) Utilice la función arima.sim() de R para simular $T = 500$ observaciones de los siguientes procesos: (i) un AR(2) con $a_1 = 0.6$ y $a_2 = 0.2$; (ii) un MA(2) con $b_1 = 0.7$ y $b_2 = 0.3$; y (iii) un ARMA(1,1) con $a_1 = 0.8$ y $b_1 = 0.4$. Tenga presente que la convención de signos de arima.sim() para el componente MA es contraria a la utilizada en este capítulo (el argumento ma corresponde a $-b_i$).
1. Grafique cada serie simulada y sus funciones ACF y PACF con acf() y pacf().
2. Verifique que los patrones muestrales coinciden con los patrones teóricos resumidos en el Cuadro 3.9.
3. Repita el ejercicio con $T = 100$ y comente cómo afecta el tamaño de muestra a la claridad de los patrones.
(Computacional) El archivo BD/Base_Transporte.xlsx incluye la serie mensual de pasajeros transportados en vuelos nacionales (columna Pax_Nal), de enero de 2000 al último mes disponible. Replique para esta serie la metodología aplicada al Metro de la CDMX en este capítulo:
1. Grafique la serie en niveles y en logaritmos. ¿Considera necesaria la transformación logarítmica?
2. Calcule la diferencia logarítmica mensual y grafique su ACF y PACF. ¿Qué órdenes $p$ y $q$ sugieren los correlogramas?
3. Estime al menos cuatro especificaciones ARMA($p$, $q$) alternativas con arima() y compare sus criterios AIC y BIC. Considere incluir una variable dummy para la caída asociada a la pandemia por COVID-19 (abril y mayo de 2020), construyéndola usted mismo como se hizo en este capítulo con D_Abr2020 y D_May2020.
4. Para el modelo seleccionado, aplique a los residuales la prueba de Ljung-Box y la prueba de normalidad de Jarque-Bera, e interprete los resultados.
(Computacional) Continúe con el modelo seleccionado en el ejercicio anterior.
1. Genere un pronóstico a 24 meses con predict() y grafíquelo junto con la serie observada, incluyendo el intervalo de predicción de $\pm 2$ errores estándar.
2. Explique por qué el intervalo de predicción se amplía conforme crece el horizonte $h$.
3. Por separado, aplique el filtro Hodrick-Prescott a la serie en logaritmos con la función hpfilter() del paquete mFilter, usando $\lambda = 14400$ por tratarse de datos mensuales. Grafique la tendencia y el ciclo resultantes.
4. Comente qué limitación presenta el filtro HP en los extremos de la muestra y cómo la aborda la variante de St-Amant y van Norden discutida en este capítulo.

Los datos y el algoritmo están disponibles en el repositorio de GitHub del proyecto.↩︎
Fuente: INEGI, https://www.inegi.org.mx/app/indicadores/?tm=0&t=1090.↩︎
Estas tasas no son porcentuales, para hacerlas porcentuales faltaría multiplicar por 100 cada valor de la serie.↩︎
Note que estas raíces son los recíprocos de las del polinomio en el operador rezago $1 - a_1 x - a_2 x^2 = 0$, que reordenado es $a_2 x^2 + a_1 x - 1 = 0$ –el mismo polinomio característico del Capítulo 2–.↩︎
Es fácil demostrar esta afirmación, sólo requiere desarrollar la expresión y utilizar el hecho de que $U_t$ es un proceso puramente aleatorio, por lo que la covarianza es cero (0).↩︎
Disponible en https://robjhyndman.com/hyndsight/arma-roots/. Los archivos arroots.R, maroots.R y plot.armaroots.R del repositorio conservan esa atribución en su encabezado.↩︎
La información y material respecto del modelo está disponible en la dirección https://www.census.gov/srd/www/x13as/.↩︎
Las definiciones de las variables empleadas son: INPC: Índice Nacional de Precios al Consumidor (2QJul2018 = 100); TC: Tipo de Cambio FIX; CETE28: Tasa de rendimiento promedio mensual de los Cetes 28, en por ciento anual; IGAE: Indicador Global de la Actividad Económica (2013 = 100); IPI: Industrial Production Index (2012 = 100).↩︎

Estadística	Valor
\(\hat{\mu} = \frac{1}{T} \sum^T_{t=1} X_t\)	115.6156
\(\hat{\gamma}(0) = \frac{1}{T} \sum^T_{t=1} (X_t - \hat{\mu})^2\)	418.1772
\(\hat{\gamma}(1) = \frac{1}{T} \sum^{T - 1}_{t=1} (X_t - \hat{\mu})(X_{t+1} - \hat{\mu})\)	372.7523
\(\hat{\gamma}(2) = \frac{1}{T} \sum^{T - 2}_{t=1} (X_t - \hat{\mu})(X_{t+2} - \hat{\mu})\)	363.3531
\(\hat{\gamma}(3) = \frac{1}{T} \sum^{T - 3}_{t=1} (X_t - \hat{\mu})(X_{t+3} - \hat{\mu})\)	341.5693
\(\hat{\rho}(1) = \hat{\gamma}(1)/\hat{\gamma}(0)\)	0.8914
\(\hat{\rho}(2) = \hat{\gamma}(2)/\hat{\gamma}(0)\)	0.8689
\(\hat{\rho}(3) = \hat{\gamma}(3)/\hat{\gamma}(0)\)	0.8168
\(Q^* = T \sum_{j = 1}^{1} \hat{\rho}(j)^2\)	251.8716
\(Q^* = T \sum_{j = 1}^{2} \hat{\rho}(j)^2\)	491.2010

Modelo	Función de Autocorrelación	Función de Autocorrelación Parcial
MA(q)	Se corta después del rezago \(q\)	Decae gradualmente
AR(p)	Decae gradualmente	Se corta después del rezago \(p\)

Variable	Coef. M1	SE M1	Coef. M2	SE M2
\(DLPaxNal_{t-1}\)	-1.0184	0.0403	-1.0322	0.0413
\(DLPaxNal_{t-2}\)	0.2354	0.0698	0.2108	0.0717
\(DLPaxNal_{t-3}\)	0.7133	0.0403	0.6992	0.0414
\(DDic_{t}\)	0.0536	0.0365	0.0352	0.0358
\(DEne_{t}\)	-0.1249	0.0343	-0.1207	0.034
\(DFeb_{t}\)	-0.0253	0.0366	-0.037	0.0359
\(DJul_{t}\)	0.2644	0.0446	0.2582	0.0447
\(DJun2009_{t}\)			0.2031	0.0868
\(DMay2009_{t}\)			-0.372	0.0866
Constante	-0.0098	0.0084	-0.0064	0.0084
\(\hat{U}_{t-1}\)	1.2669	0.0023	1.3318	0.0079
\(\hat{U}_{t-2}\)	-0.3332	0.0214	-0.2852	0.0156
\(\hat{U}_{t-3}\)	-1.3883	0.0207	-1.4342	0.008
\(\hat{U}_{t-4}\)	-0.5395		-0.6012	0.0119
\(\hat{\sigma}^2\)	0.015149		0.014124
AIC	-389.89		-406.67

Concepto	Valor
\(F\) del IGAE con el \(\lambda\) convencional (1,600)	8.77e-05
\(V\) de referencia (PIB de EUA, \(\lambda = 1{,}600\))	1.81e-04
\(\lambda\) calibrado para el IGAE	920.5