El modelo de regresión lineal en un contexto de series de tiempo es:

\[ \begin{align} y_t &= \beta_1 + \beta_2 x_t + \varepsilon_t;\\ \varepsilon_t &\sim WN(0,\sigma^2) \end{align} \]

¿Cómo usamos este modelo para hacer predicciones?

Podemos utilizar el modelo de regresión lineal para hacer un pronóstico de \(y_t\) condicional a los supuestos sobre la variable \(x_t\). Por ejemplo:

\[ \mathbb{E}\left(y_{T+h}|x_{T+h}\right)^\ast = \beta_1 + \beta_2 x_{T+h}^\ast \]

donde hemos supuesto que el valor de \(x_t\) en el periodo \(T+h\) es \(x_{T+h}^*\). Este tipo de análisis se conoce en la práctica como análisis de escenarios o análisis de contingencia.

Por supuesto, si suponemos que \(\varepsilon_t\) es un ruido blanco gaussiano, entonces la distribución condicional del pronóstico es:

\[ y_{T+h}|x_{T+h}^\ast \sim \mathcal{N}(y_{T+h}|x_{T+h}^\ast,\sigma^2) \]

Esto nos permite calcular intervalos de confianza para la predicción.

Naturalmente, el pronóstico que acabamos de realizar no está exento de incertidumbre. Para ver esto, estimamos el modelo por MCO y lo evaluamos en \(x_{T+h}^*\). En este caso, el pronóstico para el periodo \(T+h\) es:

\[ \hat{y}_{T+h,T}|x_{T+h}^\ast = b_1 + b_2 x_{T+h}^\ast, \\ \] y el residual:

\[ \begin{align} e_{T+h,T} &= y_{T+h} - \hat{y}_{T+h,T}|x_{T+h}^\ast\\ & = (\beta_1 - b_1) + (\beta_2 - b_2) x_{T+h}^\ast + \varepsilon_{T+h} \end{align} \]

Si calculamos la varianza del residual:

\[\mathbb{V}ar(e_{T+h,T}) = \mathbb{V}ar(b_1) + x_{T+h}^{\ast 2}\mathbb{V}ar(b_2) + \sigma^2\] donde los primeros dos términos reflejan la incertidumbre en los parámetros del modelo. Entonces, para construir un intervalo de confianza para el pronóstico que refleje esta incertidumbre:

\[\hat{y}_{T+h}|x_{T+h}^\ast \sim \mathcal{N}(\hat{y}_{T+h}|x_{T+h}^\ast,\mathbb{V}ar(b_1) + x_{T+h}^{\ast 2}\mathbb{V}ar(b_2) + \sigma^2)\]

En ocasiones deseamos encontrar el mejor pronóstico posible para \(y_t\). Esto es, un pronóstico no condicional. Para hacer esto, es necesario insertar un pronóstico óptimo para \(x_t\), \(x_{T+h,T}\), en la ecuación.

\[y_{T+h,T} = \beta_1 + \beta_2 x_{T+h,T}\]

Naturalmente, no siempre contamos con este valor futuro de \(x_t\) (tenemos un modelo para \(y_t\), no para \(x_t\)).

Una solución consiste en estimar \(y_t\) usando como regresores los valores pasados de \(x_t\): \(x_{t-h}, x_{t-h-1}, x_{t-h-2},\dots\). Por ejemplo, si deseamos pronosticar un periodo hacia adelante, podemos usar el modelo:

\[y_t = \beta_1 + \beta_2 x_{t-1} + \varepsilon_t\]

Por supuesto, es posible incluir más retardos de \(x_t\). También, es importante resaltar que si el objetivo es pronosticar \(y_t\) \(h\) periodos hacia adelante, las variables \(x_t\) deben retardarse al menos \(h\) periodos.

El modelo anterior puede generalizarse como (nota el cambio de notación):

\[ y_t = \phi_0 + \sum_{i=1}^{N_x}\gamma_i x_{t-i} + \varepsilon_t \]

En este caso decimos que \(y_t\) depende de los rezagos distribuidos de \(x_t\). Pero, ¿cuánto vale \(N_x\)?

Nota. ¿Por qué este modelo tendrá este nombre?

Para elegir el valor de \(N_x\) tenemos varias opciones:

• Usar pruebas de significancia estadística.
• Usar los criterios de información (AIC y BIC).
• Revisar el comportamiento de los residuales.

Los modelos con rezagos distribuidos (si las series son estacionarias) pueden estimarse por MCO o por MCG.

Los modelos con rezagos distribuidos pueden generalizarse aun más:

\[ y_t = \phi_0 + \sum_{i=1}^{N_y}\phi_i y_{t-i} + \sum_{j=1}^{N_x}\gamma_j x_{t-j} + \varepsilon_t \]

En este modelo hemos agregado retardos de la variable dependiente. Este es un modelo AR con rezagos distribuidos (ARDL) de orden \(N_y, N_x\). Este puede estimarse (si las series son estacionarias y no cointegradas) por Máxima Verosimilitud, MCO o MCG.

Nota

Se puede mostrar que los modelos ARDL pueden expresarse como modelos con errores ARMA. Para ver esto, considera el modelo \(ARDL(N_y, N_x)\):

\[ y_t = \phi_0 + \sum_{i=1}^{N_y}\phi_i y_{t-i} + \sum_{j=1}^{N_x}\gamma_i x_{t-j} + \varepsilon_t \] El operador de rezagos nos permite reescribir el modelo como:

\[ y_t = \frac{\phi_0}{\Phi(L)} + \frac{\Gamma(L)}{\Phi(L)} x_{t} + \frac{1}{\Phi(L)}\varepsilon_t; \quad \Phi(L)=1-\sum_{i=1}^{N_y}\phi_i L^i;\quad \Gamma(L)=\sum_{j=0}^{N_x}\gamma_j L^j. \]

Nota que el término\(\frac{1}{\Phi(L)}\varepsilon_t\) induce un \(MA(\infty)\) en el error, lo cual es equivalente a un modelo ARMA bajo invertibilidad.

Consideremos el siguiente modelo:

\[ u_t = \phi_0 + \phi_1 u_{t-1} + \phi_2 u_{t-2} + \gamma_1 g_{t-1} + \varepsilon_t \]

Donde:

• \(u_t\) es la tasa de desempleo en tiempo t,
• \(g_t\) es la tasa de crecimiento del PIB en tiempo t.

Usando datos de los Estados Unidos tenemos que:

Podemos hacer pronósticos usando modelos autorregresivos con rezagos distribuidos usando las mismas técnicas aprendidas antes (ver modelos ARMA).

La representación ARMA del modelo ARDL nos permite calcular su función de impulso-respuesta. En el ejemplo tenemos:

Los modelos con vectores autorregresivos (VAR) son una generalización de los modelos autorregresivos uniecuacionales (AR) y de rezagos distribuidos (DL). Estos sirven para capturar la evolución y las interdependencias entre múltiples series de tiempo.

Como ejemplo, consideramos la relación entre el inicio y el fin de construcción de nuevas viviendas en los Estados Unidos:

Un VAR involucra \(k\) variables. En un vector autorregresivo de orden \(p\) con \(k\) variables, VAR(\(p\)), tenemos que estimar \(k\) diferentes ecuaciones. En cada ecuación tenemos como regresores \(p\) rezagos de la variable dependiente y \(p\) rezagos de todas las demás variables del sistema.

Por ejemplo, para un VAR(\(1\)) con dos variables, \(y_{1,t}\) y \(y_{2,t}\), tenemos:

\[ \begin{align} y_{1,t} &= \phi^1_{1,1}y_{1,t-1} + \phi^1_{2,1}y_{2,t-1} + \varepsilon_{1,t} \\ y_{2,t} &= \phi_{1,1}^2y_{1,t-1} + \phi_{2,1}^2y_{2,t-1} + \varepsilon_{2,t} \\ \varepsilon_{1,t} & \sim WN(0,\sigma^2_1) \\ \varepsilon_{2,t} & \sim WN(0,\sigma^2_2) \\ \end{align} \] y:

\[ \begin{align} \mathbb{C}ov(\varepsilon_{1,t},\varepsilon_{2,\tau}) & = \begin{cases} \sigma_{12}; & \ t=\tau \\ 0; & \ t \neq \tau \\ \end{cases} \end{align} \]

Un VAR(\(1\)) puede escribirse en una forma más compacta:

\[ \mathbf{y}_t = \boldsymbol{\phi}_1 \mathbf{y}_{t-1} + \boldsymbol{\varepsilon}_t \]

donde:

\[ \mathbf{y}_t = \begin{pmatrix} y_{1,t} \\ y_{2,t}\end{pmatrix}; \qquad \boldsymbol{\phi}_1 = \begin{pmatrix} \phi_{1,1}^1 & \phi_{2,1}^1 \\ \phi^2_{1,1} & \phi^2_{2,1}\end{pmatrix}; \qquad \boldsymbol{\varepsilon}_t = \begin{pmatrix} \varepsilon_{1,t} \\ \varepsilon_{2,t}\end{pmatrix} \]

Ejercicio

Considera el siguiente modelo VAR(\(2\)):

\[ \begin{align} y_{1,t} &= \phi^1_{1,1}y_{1,t-1} + \phi_{2,1}^1y_{2,t-1} + \varepsilon_{1,t} \\ y_{2,t} &= \phi_{1,1}^2y_{1,t-1} + \phi_{2,2}^2y_{2,t-2} + \varepsilon_{2,t} \\ \varepsilon_{1,t-1} & \sim WN(0,\sigma^2_1) \\ \varepsilon_{2,t-1} & \sim WN(0,\sigma^2_2) \\ \mathbb{C}ov(\varepsilon_{1,t},\varepsilon_{2,\tau}) & = \begin{cases} \sigma_{12}; & \ t=\tau \\ 0; & \ t \neq \tau \\ \end{cases} \end{align} \]

Escribe el modelo en forma compacta.

En general, un VAR(\(p\)) puede escribirse como:

\[ \mathbf{y}_t = \boldsymbol{\phi}_1 \mathbf{y}_{t-1} + \boldsymbol{\phi}_2 \mathbf{y}_{t-2} + \dots + \boldsymbol{\phi}_p \mathbf{y}_{t-p} + \boldsymbol{\varepsilon}_t \]

O como:

\[ \boldsymbol{\Phi}(L)\mathbf{y}_t=\boldsymbol{\varepsilon}_t \]

Donde \(\boldsymbol{\Phi}(L)\) es un polinomio de orden \(p\):

\[ \boldsymbol{\Phi}(L) = \mathbf{I}_k - \boldsymbol{\phi}_1 L - \boldsymbol{\phi}_2 L^2 - \boldsymbol{\phi}_p L^p \]

Además, es fácil deducir que, para un VAR(\(p\)):

\[ \begin{align} \mathbb{E}(\boldsymbol{\varepsilon}_t) &= \mathbf{0} \\ \mathbb{E}(\boldsymbol{\varepsilon}_t \boldsymbol{\varepsilon}_\tau') &= \begin{cases} \boldsymbol{\Omega} & \ t=\tau \\ \mathbf{0} & \ t \neq \tau \end{cases} \end{align} \] donde \(\boldsymbol{\Omega}\) es una matriz cuadrada de tamaño \(k\), simétrica y definida positiva.

Igual que para los procesos AR(\(p\)), un proceso VAR(\(p\)) es estacionario en tanto las raíces del polinomio:

\[ |\mathbf{I}_k - \boldsymbol{\phi}_1 L - \boldsymbol{\phi}_2 L^2 - \boldsymbol{\phi}_p L^p|=0 \]

se encuentren fuera del círculo unitario.

Ejemplo

Consideremos el proceso:

\[ \mathbf y_t = \begin{pmatrix} 0.7 & 0.2\\ 0.1 & 0.6 \end{pmatrix} \mathbf y_{t-1} + \boldsymbol\varepsilon_t, \qquad \boldsymbol\varepsilon_t \sim WN(\mathbf 0, \boldsymbol\Omega) \]

El polinomio característico del proceso es:

\[ |\mathbf I_2 - \boldsymbol{\Phi}_1 L| = 0 \]

En este caso, es fácil mostrar que raíces de este polinomio son:

\[ L_1 = \frac{1}{0.79}>1, \qquad L_2 = \frac{1}{0.51}>1 \]

Nota

La condición de estacionariedad se puede encontrar a partir de los valores propios de una matriz cuadrada. Para ver esto, considera un proceso VAR(\(p\)):

\[ \mathbf{y}_t = \boldsymbol{\phi}_1 \mathbf{y}_{t-1} + \boldsymbol{\phi}_2 \mathbf{y}_{t-2} + \dots + \boldsymbol{\phi}_p \mathbf{y}_{t-p} + \boldsymbol{\varepsilon}_t \]

Definimos la matriz compañera como:

\[ \mathbf{A} = \begin{pmatrix} \boldsymbol{\phi}_1 & \boldsymbol{\phi}_2 & \dots & \boldsymbol{\phi}_{p-1} & \boldsymbol{\phi}_{p}\\ \mathbf{I}_k & \mathbf{0} & \dots & \mathbf{0} & \mathbf{0}\\ \mathbf{0} & \mathbf{I}_k & \dots & \mathbf{0} & \mathbf{0}\\ \vdots & \vdots & \ddots & \vdots & \vdots \\ \mathbf{0} & \mathbf{0} & \dots & \mathbf{I}_k & \mathbf{0} \end{pmatrix} \] Para que este proceso sea estacionario, es necesario que los valores propios de la matriz compañera se encuentren dentro del círculo unitario.

Ejemplo

Consideremos el proceso:

\[ \mathbf y_t = \begin{pmatrix} 0.7 & 0.2\\ 0.1 & 0.6 \end{pmatrix} \mathbf y_{t-1} + \boldsymbol\varepsilon_t, \qquad \boldsymbol\varepsilon_t \sim WN(\mathbf 0, \boldsymbol\Omega) \]

Muestra que los valores propios de la matriz \(\boldsymbol{\Phi}_1\) son:

\[ \lambda_1 = 0.79<1, \qquad \lambda_2 = 0.51<1 \]

Ejercicio

Considera el proceso:

\[ \mathbf y_t=\begin{pmatrix} 0.50 & 0\\ 0 & 0.40 \end{pmatrix}\,\mathbf y_{t-1}+\begin{pmatrix} 0.20 & 0\\ 0 & 0.10 \end{pmatrix}\,\mathbf y_{t-2}+\boldsymbol\varepsilon_t,\quad \boldsymbol\varepsilon_t\sim WN(\mathbf 0,\boldsymbol\Omega), \]

1. Encuentra las raíces del polinomio \(|\mathbf I-\Phi_1 L-\Phi_2 L^2|=0\).

2. Encuentra la matriz compañera.

3. Encuentra los valores propios de la matriz compañera.

4. ¿Es este proceso estacionario?

Considera un proceso VAR(\(p\)):

\[ \mathbf y_t = \sum_{i=1}^{p} \boldsymbol\Phi_i\,\mathbf y_{t-i} + \boldsymbol\varepsilon_t, \qquad \mathbb E(\boldsymbol\varepsilon_t)=\mathbf 0, \qquad \mathbb E(\boldsymbol\varepsilon_t\boldsymbol\varepsilon_\tau’)=\begin{cases}\boldsymbol\Omega,& t=\tau\\ \mathbf 0,& t\ne \tau\end{cases} \]

La función de autocovarianzas puede calcularse como:

\[ \boldsymbol\Gamma(h) \;=\; \mathbb E\big[(\mathbf y_t-\boldsymbol\mu)(\mathbf y_{t-h}-\boldsymbol\mu)’\big],\quad h=0,1,2,\dots \]

Nota. Al igual que en el caso de los procesos AR(\(p\)), tenemos estacionariedad si \(\boldsymbol\Gamma(h)\) no depende de \(t\).

Por simplicidad, supongamos que \(\boldsymbol\mu = \mathbf{0}\). En este caso:

\[ \boldsymbol\Gamma(h) \;=\; \mathbb E\big[\mathbf y_t\mathbf y_{t-h}’\big],\quad h=0,1,2,\dots \]

Nota

Nota que ahora:

\[ \begin{align} \boldsymbol\Gamma(h) \;&=\; \mathbb E\big[\mathbf y_t\mathbf y_{t-h}’\big]\\ &=\; \mathbb E\big[\mathbf y_{t+h}\mathbf y_{t+h-h}’\big]\\ &=\; \mathbb E\big[\mathbf y_{t+h}\mathbf y_{t}’\big]\\ &\neq\; \boldsymbol\Gamma(-h) \end{align} \]

Sin embargo:

\[ \begin{align} \boldsymbol\Gamma(h)' &=\; \mathbb E\big[\mathbf y_{t+h}\mathbf y_{t}’\big]'\\ \;&=\; \mathbb E\big[\mathbf y_t \mathbf y_{t+h}’\big]\\ &=\; \boldsymbol\Gamma(-h) \end{align} \]

Las ecuaciones de Yule-Walker se calculan de manera similar al caso de los procesos AR(\(p\)). Estas nos ayudan a sistematizar el cálculo de la matriz de las autocovarianzas de un VAR(\(p\)).

Paso 1

Multiplicamos el proceso por \(\mathbf y_t’\) y calculamos el valor esperado:

\[ \mathbb E(\mathbf y_t\,\mathbf y_t’) = \sum_{i=1}^{p} \boldsymbol\Phi_i\,\mathbb E(\mathbf y_{t-i}\,\mathbf y_t’) \; +\; \mathbb E(\boldsymbol\varepsilon_t\,\mathbf y_t’)\\ \]

Usando la definición de un VAR(\(p\)) y las propiedades de \(\boldsymbol\varepsilon_t\):

\[ \boldsymbol\Gamma(0) = \sum_{i=1}^{p} \boldsymbol\Phi_i\,\boldsymbol\Gamma(i)’ + \boldsymbol\Omega \]

Paso 2

Ahora multiplicamos el proceso por \(\mathbf y_{t-h}’\) para \(h\geq 1\). Obtenemos que:

\[ \begin{align} \mathbb{E}(\mathbf y_t\,\mathbf y_{t-h}’) &= \sum_{i=1}^{p} \boldsymbol\Phi_i\,\mathbb E(\mathbf y_{t-i}\,\mathbf y_{t-h}’) \; +\;\mathbb E(\boldsymbol\varepsilon_t\,\mathbf y_{t-h}’) \\ &=\sum_{i=1}^{p} \boldsymbol\Phi_i\,\mathbb E(\mathbf y_{t-i}\,\mathbf y_{t-h}’) \end{align} \] Esto es:

\[ \boldsymbol\Gamma(h) = \sum_{i=1}^{p} \boldsymbol\Phi_i\,\boldsymbol\Gamma(h-i) \]

Considera un VAR(\(1\)):

\[ \mathbf{y}_t = \boldsymbol{\phi}_1 \mathbf{y}_{t-1} + \boldsymbol{\varepsilon}_t \]

El paso 1 implica:

\[ \boldsymbol\Gamma(0) = \boldsymbol\Phi_1\,\boldsymbol\Gamma(1)' + \boldsymbol\Omega \]

El paso 2 implica:

\[ \boldsymbol\Gamma(h) =\boldsymbol\Phi_1\,\boldsymbol\Gamma(h-1) \]

Entonces:

\[ \boldsymbol\Gamma(0) = \boldsymbol\Phi_1\,\boldsymbol\Gamma(0)\, \boldsymbol\Phi_1' + \boldsymbol\Omega \]

Esta se conoce como ecuación de Lyapunov discreta y nos permite obtener la matriz de varianzas y covarianzas contemporáneas de un proceso VAR(\(1\)) estacionario.

Nota

La ecuación de Lyapunov discreta se extiende de forma natural a procesos VAR(\(p\)). Para esto, el modelo puede reescribirse en su forma estado-espacio:

\[ \boldsymbol\xi_t = F\,\boldsymbol\xi_{t-1} + \boldsymbol\nu_t, \] La covarianza estacionaria de \(\boldsymbol\xi_t\) satisface:

\[ \Gamma_\xi(0) = F\,\Gamma_\xi(0)\,F' + Q. \]

El bloque superior izquierdo de \(\Gamma_\xi(0)\) corresponde a \(\Gamma(0)\) del VAR original.

Ejemplo

Considera el VAR(1):

\[ \mathbf y_t = \begin{pmatrix} 0.7 & 0.2\\ 0.1 & 0.6 \end{pmatrix} \mathbf y_{t-1} + \boldsymbol\varepsilon_t, \qquad \boldsymbol{\Omega} = \begin{pmatrix} 1 & 0.5\\ 0.5 & 1 \end{pmatrix} \]

Utiliza las ecuaciones de Yule–Walker para encontrar que:

\[ \boldsymbol\Gamma(0) \approx \begin{pmatrix} 2.85 & 1.71\\ 1.71 & 1.81 \end{pmatrix}. \]

A partir del resultado anterior, podemos encontrar que:

\[ \begin{align} \boldsymbol\Gamma(1) &= \Phi_1 \boldsymbol\Gamma(0)\\ &\approx \begin{pmatrix} 2.47 & 1.64\\ 1.28 & 1.31 \end{pmatrix} \end{align} \] Y:

\[ \begin{align} \boldsymbol\Gamma(2) &= \Phi_1 \boldsymbol\Gamma(1)\\ &\approx \begin{pmatrix} 2.01 & 1.39\\ 1.04 & 1.03 \end{pmatrix} \end{align} \]

Supongamos que las raíces del polinomio \(|\boldsymbol\Phi(L)|=0\) se encuentran fuera del círculo unitario. En este caso, el polinomio \(\boldsymbol\Psi(L) = \boldsymbol\Phi(L)^{-1}\) existe, lo que nos permite expresar un proceso VAR(\(p\)) en la forma de un proceso VMA(\(\infty\)):

\[ \begin{align} \mathbf{y}_t&=\boldsymbol{\Phi}(L)^{-1}\boldsymbol{\varepsilon}_t\\ &=\boldsymbol{\Psi}(L)\boldsymbol{\varepsilon}_t\\ &= \sum_{i=0}^\infty \boldsymbol{\Psi}_i\boldsymbol{\varepsilon}_{t-i} \end{align} \] Con:

\[ \begin{align} \boldsymbol{\Psi}_0 &= \mathbf{I}_k \\ \boldsymbol{\Psi}_1 &= \boldsymbol{\phi}_1 \\ \boldsymbol{\Psi}_2 &= \boldsymbol{\phi}_1 \boldsymbol{\Psi}_1 + \boldsymbol{\phi}_2\\ \vdots \\ \boldsymbol{\Psi}_s &= \boldsymbol{\phi}_1 \boldsymbol{\Psi}_{s-1} + \boldsymbol{\phi}_2 \boldsymbol{\Psi}_{s-2} + \dots + \boldsymbol{\phi}_p \boldsymbol{\Psi}_{s-p} \end{align} \]

Ejemplo

Consideremos un VAR(\(1\)) bivariado:

\[ \mathbf y_t=\boldsymbol{\Phi}_1\,\mathbf y_{t-1}+\boldsymbol\varepsilon_t, \qquad \boldsymbol{\Phi}_1=\begin{pmatrix}0.7&0.2\\[2pt]0.1&0.6\end{pmatrix}, \quad \mathbb E(\boldsymbol\varepsilon_t\boldsymbol\varepsilon_t’)=\boldsymbol\Omega =\begin{pmatrix}1&0.5\\[2pt]0.5&1\end{pmatrix}. \]

(Usando las propiedades del operador de rezagos, las cuales se extienden al caso matricial) No es difícil demostrar que la representación VMA(\(\infty\)) es:

\[ \mathbf y_t \;=\; \sum_{h=0}^{\infty} \boldsymbol{\Psi}_h\,\boldsymbol\varepsilon_{t-h},\qquad \boldsymbol{\Psi}_0=\mathbf I_2,\ \ \boldsymbol{\Psi}_h=\boldsymbol{\Phi}_1^{\,h}\ (h\ge1). \]

Ejemplo (continuación)

Específicamente:

\[ \begin{align} \boldsymbol{\Psi}_0&=\begin{pmatrix}1&0\\0&1\end{pmatrix}\\ \boldsymbol{\Psi}_1&=\begin{pmatrix}0.7&0.2\\0.1&0.6\end{pmatrix}\\ \boldsymbol{\Psi}_2&=\begin{pmatrix}0.51&0.26\\0.13&0.38\end{pmatrix}\\ \boldsymbol{\Psi}_3&=\begin{pmatrix}0.383&0.258\\0.129&0.254\end{pmatrix}\\ \boldsymbol{\Psi}_4&\approx\begin{pmatrix}0.294&0.231\\0.116 & 0.178\end{pmatrix}\\ &\vdots \end{align} \]

La matriz de autocovarianzas de un proceso VAR(\(p\)) también puede obtenerse a partir de su representación VMA(\(\infty\)). En particular:

\[ \begin{align} \boldsymbol\Gamma(h) &= \mathbb{E}(\mathbf y_t\mathbf y_{t-h}’) \\ &= \mathbb{E}\left[\left(\sum_{j=0}^\infty \boldsymbol\Psi_j\boldsymbol\varepsilon_{t-j}\right) \left(\sum_{i=0}^\infty \boldsymbol\Psi_i\boldsymbol\varepsilon_{t-h-i}\right)’\right] \\ &= \sum_{j=0}^\infty \boldsymbol\Psi_j\boldsymbol\Omega\boldsymbol\Psi_{j+h}’ \end{align} \]

Nota. Esta expresión es equivalente a la obtenida con las ecuaciones de Yule–Walker, pero resulta más práctica si disponemos de las matrices \(\boldsymbol\Psi_h\).

Ejercicio

Considera un VAR(\(1\)) bivariado:

\[ \mathbf y_t=\boldsymbol{\Phi}_1\,\mathbf y_{t-1}+\boldsymbol\varepsilon_t, \qquad \boldsymbol{\Phi}_1=\begin{pmatrix}0.7&0.2\\[2pt]0.1&0.6\end{pmatrix}, \quad \mathbb E(\boldsymbol\varepsilon_t\boldsymbol\varepsilon_t’)=\boldsymbol\Omega =\begin{pmatrix}1&0.5\\[2pt]0.5&1\end{pmatrix}. \] Utiliza la fórmula de \(\boldsymbol\Gamma(h)\) para mostrar que esta es equivalente a la encontrada usando el procedimiento de Yule-Walker.

Pista. Quizás sea necesario usar tu computadora.

Consideremos un proceso AR(\(1\)):

\[ y_t = \phi_1 y_{t-1} + \varepsilon_t, \qquad |\phi_1|<1. \]

Como hemos visto, este admite la siguiente representación MA(\(\infty\)):

\[ y_t = \varepsilon_t + \phi_1 \varepsilon_{t-1} + \phi_1^2 \varepsilon_{t-2} + \dots \]

La Función de Impulso-Respuesta (FIR) muestra cómo un choque en el término de error, \(\varepsilon_t\), afecta a una serie de tiempo, \(y_t\), en periodos futuros. En el caso de un AR(\(1\)) tenemos:

Considera ahora un VAR(\(1\)):

\[ \mathbf y_t = \boldsymbol{\Phi}_1 \mathbf y_{t-1} + \boldsymbol\varepsilon_t, \qquad \mathbb E(\boldsymbol\varepsilon_t\boldsymbol\varepsilon_t')=\boldsymbol\Omega. \]

Como hemos visto, este puede representarse como:

\[ \begin{align} \mathbf y_t &= \boldsymbol\varepsilon_t + \boldsymbol{\Phi}_1 \boldsymbol\varepsilon_{t-1} + \boldsymbol{\Phi}_1 ^2\boldsymbol\varepsilon_{t-2} + \boldsymbol{\Phi}_1 ^3\boldsymbol\varepsilon_{t-3} + \cdots \end{align} \]

Cada matriz \(\boldsymbol{\Phi}_1 ^h\) describe cómo un choque en (cualquiera de los elementos de) \(\boldsymbol\varepsilon_t\) afecta la trayectoria de \(\mathbf y_t\).

Ejemplo

Consideremos un VAR(\(1\)) bivariado:

\[ \mathbf y_t = \begin{pmatrix} 0.7 & 0.2 \\ 0.1 & 0.6 \end{pmatrix} \mathbf y_{t-1} + \boldsymbol\varepsilon_{t} \]

En este caso:

\[ \mathbf y_t = \begin{pmatrix} 1 & 0\\ 0 & 1 \end{pmatrix} \boldsymbol\varepsilon_t + \begin{pmatrix} 0.7 & 0.2\\ 0.1 & 0.6 \end{pmatrix} \boldsymbol\varepsilon_{t-1} + \begin{pmatrix} 0.51 & 0.26\\ 0.13 & 0.38 \end{pmatrix} \boldsymbol\varepsilon_{t-2} + \begin{pmatrix} 0.383 & 0.258\\ 0.129 & 0.254 \end{pmatrix} \boldsymbol\varepsilon_{t-3} +\;\cdots \]

Para interpretar, consideremos un choque de una unidad en \(\varepsilon_{1,t}\):

Por su parte, el efecto de un choque de una unidad en \(\varepsilon_{2,t}\):

Para un VAR(\(p\)) aplicamos la misma lógica. Tenemos:

\[ \mathbf y_t \;=\; \sum_{i=1}^p \boldsymbol{\Phi}_i\,\mathbf y_{t-i} \;+\; \boldsymbol\varepsilon_t, \qquad \mathbb E(\boldsymbol\varepsilon_t\boldsymbol\varepsilon_t')=\boldsymbol\Omega. \]

Su representación VMA(\(\infty\)) es:

\[ \mathbf y_t \;=\; \sum_{h=0}^{\infty} \boldsymbol\Psi_h\,\boldsymbol\varepsilon_{t-h}, \qquad \boldsymbol\Psi_0=\mathbf I_k,\quad \boldsymbol\Psi_h \;=\; \sum_{i=1}^p \boldsymbol{\Phi}_i\,\boldsymbol\Psi_{h-i}\ \ (h\ge 1), \]

Las matrices \(\boldsymbol\Psi_h\) son las funciones impulso–respuesta (FIR). Cada una indica cómo un choque actual afecta a las variables endógenas en los siguientes periodos.

Ejercicio

Considera el siguiente proceso VAR(\(2\)):

\[ \mathbf y_t=\begin{pmatrix}0.6&0.2\\[2pt]0.1&0.5\end{pmatrix} y_{t-1} + \begin{pmatrix}0.1&0\\[2pt]0&0.1\end{pmatrix} y_{t-2}+\boldsymbol\varepsilon_t,\qquad \boldsymbol\Omega = \mathbf{I}_2 \]

1. Verifica que el proceso es estable.

2. Encuentra la forma VMA(\(\infty\)) del proceso.

3. Construye la FIR para un choque de una unidad en \(\varepsilon_{1,t}\) y en \(\varepsilon_{2,t}\) hasta \(h=3\).

Considera un VAR(\(p\)):

\[ \mathbf y_t = \sum_{i=1}^p \boldsymbol{\Phi}_i \mathbf y_{t-i} + \boldsymbol\varepsilon_t, \qquad \mathbb E(\boldsymbol\varepsilon_t\boldsymbol\varepsilon_t')=\boldsymbol\Omega. \]

Nota que, dado que \(\boldsymbol\Omega\) no es necesariamente diagonal, un un choque unitario en \(\varepsilon_{k,t}\) podría afectar a todas las ecuaciones del sistema. Por lo tanto, las FIR no pueden interpretarse como choques estructurales puros.

Para evitar este problema, podemos construir innovaciones mutuamente ortogonales. Esto es:

\[ \mathbb E(\boldsymbol u_t\boldsymbol u_t')=\mathbf I_k, \]

En este caso, cada impulso \(u_{j,t}\) representa un shock estructural independiente.

Para construir innovaciones ortogonales, podemos utilizar la descomposición de Cholesky:

\[ \boldsymbol\Omega = \mathbf P \mathbf P', \]

donde \(\mathbf P\) es una matriz triangular inferior.

En este caso:

\[ \boldsymbol\varepsilon_t = \mathbf P \boldsymbol u_t, \qquad \Rightarrow \qquad \mathbb E(\boldsymbol u_t\boldsymbol u_t')=\mathbf I_k. \]

Ahora, sustituimos en la forma VMA(\(\infty\)) de un VAR(\(p\)):

\[ \begin{align} \mathbf y_t &= \sum_{h=0}^\infty \boldsymbol\Psi_h \boldsymbol\varepsilon_{t-h}\\ &= \sum_{h=0}^\infty \boldsymbol\Psi_h \mathbf P \boldsymbol u_{t-h}. \end{align} \]

De aquí, podemos construir la FIR ortogonalizada o FIR estructural. Esta es:

\[ \boldsymbol\Theta_h = \boldsymbol\Psi_h \mathbf P. \]

Consideremos nuevamente el VAR(\(1\)) bivariado:

\[ \mathbf y_t = \begin{pmatrix} 0.7 & 0.2 \\ 0.1 & 0.6 \end{pmatrix}\mathbf y_{t-1} + \boldsymbol\varepsilon_t, \qquad \boldsymbol\Omega = \begin{pmatrix} 1 & 0.6\\ 0.6 & 1 \end{pmatrix}. \]

Como primer paso, obtenemos la descomposición de Cholesky de la matriz \(\boldsymbol\Omega\):

\[ \mathbf P = \begin{pmatrix} 1 & 0\\ 0.6 & 0.8 \end{pmatrix}, \qquad \boldsymbol\varepsilon_t = \mathbf P \boldsymbol u_t, \quad \mathbb E(\boldsymbol u_t\boldsymbol u_t') = \mathbf I_2. \]

La FIR ortogonalizada no es más que:

\[ \boldsymbol\Theta_h = \boldsymbol\Psi_h\,\mathbf P. \]

Por ejemplo:

\[ \boldsymbol\Theta_1 = \begin{pmatrix} 0.7&0.2\\[2pt]0.1&0.6 \end{pmatrix} \begin{pmatrix} 1&0\\[2pt]0.6&0.8 \end{pmatrix} = \begin{pmatrix} 0.82&0.16\\[2pt]0.46&0.48 \end{pmatrix}. \]

Continuando con el procedimiento, no es difícil demostrar que la FIR ortogonalizada es:

Choque estructural en \(u_{1,t}\):

Choque estructural en \(u_{2,t}\):

La Descomposición de la varianza del error de pronóstico (FEVD) (Forecast Error Variance Decomposition) mide qué proporción de la varianza del error de pronóstico de cada variable se debe a los diferentes choques estructurales del sistema.

Para calcular esta, partimos de la representación VMA(\(\infty\)) ortogonalizada:

\[ \mathbf y_t = \sum_{h=0}^\infty \boldsymbol\Theta_h \boldsymbol u_{t-h}, \qquad \mathbb E(\boldsymbol u_t \boldsymbol u_t') = \mathbf I_k, \]

donde:

\[ \boldsymbol\Theta_h = \boldsymbol\Psi_h \mathbf P, \qquad \mathbf P \mathbf P' = \boldsymbol\Omega. \]

El error de pronóstico a \(k\) pasos para \(\mathbf y_t\) es:

\[ \mathbf e_{t+k|t} = \mathbf y_{t+k} - \mathbb E_t(\mathbf y_{t+k}) = \sum_{h=0}^{k-1} \boldsymbol\Theta_h \boldsymbol u_{t+k-h}. \]

Su matriz de varianza–covarianza es:

\[ \mathbb E(\mathbf e_{t+k|t}\mathbf e_{t+k|t}') = \sum_{h=0}^{k-1} \boldsymbol\Theta_h \boldsymbol\Theta_h'. \]

Cada elemento diagonal corresponde a la varianza del error de pronóstico de cada variable a \(k\) pasos.

Para cada variable \(y_i\), la varianza de su error de pronóstico a \(h\) pasos es:

\[ \text{Var}(e_{i,t+h|t}) = \sum_{s=0}^{h-1} \sum_{m=1}^{k} (\boldsymbol\Theta_s)_{im}^2. \]

La contribución del choque estructural \(u_j\) a dicha varianza es:

\[ \text{Var}_{(i\leftarrow j)}(e_{i,t+h|t}) = \sum_{s=0}^{h-1} (\boldsymbol\Theta_s)_{ij}^2. \]

Entonces, la proporción de la varianza del error de pronóstico atribuible al error \(u_j\) es:

\[ \text{FEVD}_{ij}(h) = \frac{\displaystyle \sum_{s=0}^{h-1} (\boldsymbol\Theta_s)_{ij}^2} {\displaystyle \sum_{s=0}^{h-1} \sum_{m=1}^{k} (\boldsymbol\Theta_s)_{im}^2}. \] En general, la FEVD:

• Cuantifica qué tanto de la incertidumbre futura de cada variable se debe a los diferentes shocks estructurales.
• Las proporciones \(\text{FEVD}_{ij}(k)\) para un \(i\) fijo suman 1.
• A medida que \(k \to \infty\), la descomposición converge al efecto de largo plazo de cada choque sobre la variable \(i\).

Consideremos el VAR(\(1\)) bivariado:

\[ \mathbf y_t = \begin{pmatrix} 0.7 & 0.2 \\ 0.1 & 0.6 \end{pmatrix}\mathbf y_{t-1} + \boldsymbol\varepsilon_t, \qquad \boldsymbol\Omega = \begin{pmatrix} 1 & 0.6\\ 0.6 & 1 \end{pmatrix}. \] Donde, como hemos visto:

\[ \boldsymbol\Theta_0 = \begin{pmatrix}1&0\\[2pt]0.6&0.8\end{pmatrix},\quad \boldsymbol\Theta_1 = \begin{pmatrix}0.82&0.16\\[2pt]0.46&0.48\end{pmatrix},\quad \boldsymbol\Theta_2 = \begin{pmatrix}0.636&0.208\\[2pt]0.328&0.356\end{pmatrix},\quad \boldsymbol\Theta_3 = \begin{pmatrix}0.505&0.190\\[2pt]0.244&0.277\end{pmatrix}. \]

Para \(h=1\), usando:

\[ \boldsymbol{\Theta}_0 = \begin{pmatrix} 1 & 0\\ 0.6 & 0.8 \end{pmatrix}, \] la FEVD se calcula como:

\[ \text{FEVD}_{ij}(1) = \frac{(\boldsymbol{\Theta}_0)_{ij}^2} {\sum_{m=1}^{k}(\boldsymbol{\Theta}_0)_{im}^2}, \qquad k=2. \]

Entonces:

\[ \begin{aligned} \text{FEVD}_{1\leftarrow 1}(1) &= \frac{1^2}{1^2 + 0^2}\\ &= 1.00 \\ \text{FEVD}_{1\leftarrow 2}(1) &= \frac{0^2}{1^2 + 0^2}\\ &= 0.00\\ \text{FEVD}_{2\leftarrow 1}(1) &= \frac{0.6^2}{0.6^2 + 0.8^2}\\ &= 0.36\\ \text{FEVD}_{2\leftarrow 2}(1) &= \frac{0.8^2}{0.6^2 + 0.8^2}\\ & = 0.64 \end{aligned} \]

Para \(h=2\), usamos las matrices: \[ \boldsymbol{\Theta}_0 = \begin{pmatrix} 1 & 0\\ 0.6 & 0.8 \end{pmatrix},\quad \boldsymbol{\Theta}_1 = \begin{pmatrix} 0.82 & 0.16\\ 0.46 & 0.48 \end{pmatrix}, \] Entonces:

\[ \text{FEVD}_{ij}(2)= \frac{\sum_{s=0}^{1}(\boldsymbol{\Theta}_s)_{ij}^2} {\sum_{s=0}^{1}\sum_{m=1}^{k}(\boldsymbol{\Theta}_s)_{im}^2},\quad k=2. \]

Sustituyendo:

\[ \begin{align} \text{FEVD}_{1\leftarrow 1}(2) &= \dfrac{1.6724}{1.6980}\\ & \approx 0.985 \\ \text{FEVD}_{1\leftarrow 2}(2) &= \dfrac{0^2+0.16^2}{1.6980}\\ & \approx 0.015\\ \text{FEVD}_{2\leftarrow 1}(2) &= \dfrac{0.5716}{1.4420}\\ & \approx 0.397 \\ \text{FEVD}_{2\leftarrow 2}(2) &= \dfrac{0.8^2+0.48^2}{1.4420}\\ & \approx 0.603 \end{align} \]

Ejercicio

Muestra que la FEVD para \(y_1\) es:

Y para \(y_2\):

Consideremos el modelo:

\[ \mathbf y_t = \boldsymbol c + \sum_{i=1}^p \boldsymbol\Phi_i \mathbf y_{t-i} + \boldsymbol\varepsilon_t, \qquad \boldsymbol\varepsilon_t \sim \text{i.i.d. } \mathcal N(\mathbf 0, \boldsymbol\Omega). \]

Para simplificar, definimos:

\[ \mathbf x_t = \begin{bmatrix} 1 & \mathbf y_{t-1} & \cdots & \mathbf y_{t-p} \end{bmatrix}, \qquad \boldsymbol\beta = \begin{bmatrix} \boldsymbol c \\ \boldsymbol\Phi_1 \\ \vdots \\ \boldsymbol\Phi_p \end{bmatrix}. \] Esto nos permite reescribir el modelo como:

\[ \mathbf y_t = \mathbf x_t \boldsymbol\beta + \boldsymbol\varepsilon_t \]

Si suponemos que los errores son normales (multivariados):

\[ \boldsymbol\varepsilon_t \sim \text{i.i.d. } \mathcal N(\mathbf 0, \boldsymbol\Omega), \qquad t = p+1,\dots,T \]

y que las observaciones (\(\mathbf{y}_{-p+1}, \dots, \mathbf{y}_{0}\)) son conocidas, la función de verosimilitud (condicional) es:

\[ L(\boldsymbol\beta, \boldsymbol\Omega) = (2\pi)^{-\frac{nk}{2}} |\boldsymbol\Omega|^{-\frac{n}{2}} \exp\left\{ -\frac{1}{2}\sum_{t=p+1}^T (\mathbf y_t - \mathbf x_t\boldsymbol\beta)' \boldsymbol\Omega^{-1}(\mathbf y_t - \mathbf x_t\boldsymbol\beta) \right\}; \ n=T-p \]

La función de log-verosimilitud (condicional) es:

\[ \ell(\boldsymbol\beta, \boldsymbol\Omega) = -\frac{nk}{2}\ln(2\pi) -\frac{n}{2}\ln|\boldsymbol\Omega| -\frac{1}{2}\sum_{t=p+1}^T (\mathbf y_t - \mathbf x_t\boldsymbol\beta)' \boldsymbol\Omega^{-1}(\mathbf y_t - \mathbf x_t\boldsymbol\beta). \]

A partir de lo anterior, podemos encontrar que el estimador de máxima verosimilitud es:

\[ \widehat{\boldsymbol\beta} = \left(\sum_{t=p+1}^{T}\mathbf x_t’\mathbf x_t\right)^{-1} \left(\sum_{t=p+1}^{T}\mathbf x_t’\mathbf y_t\right). \]

Este estimador coincide con el estimador de Mínimos Cuadrados Ordinarios (MCO) ecuación por ecuación.

Por otra parte, el estimador de máxima verosimilitud para la matriz de varianzas y covarianzas de los errores es:

\[ \widehat{\boldsymbol\Omega} = \frac{1}{n}\sum_{t=p+1}^{T} \widehat{\boldsymbol\varepsilon}_t\widehat{\boldsymbol\varepsilon}_t’, \qquad \widehat{\boldsymbol\varepsilon}_t = \mathbf y_t - \mathbf x_t \widehat{\boldsymbol\beta}. \]

La causalidad de Granger evalúa si los valores pasados de una variable \(y_{2t}\) ayudan a predecir a otra variable \(y_{1t}\), más allá de la información contenida en los rezagos de \(y_{1t}\) misma.

Formalmente, decimos que \(y_{2t}\) no causa a \(y_{1t}\) en el sentido de Granger si:

\[ \mathbb E(y_{1t} \mid \Omega_{t-1}) = \mathbb E(y_{1t} \mid \Omega_{t-1}\setminus {y_{2,t-1}, y_{2,t-2}, \dots}), \]

donde \(\Omega_{t-1}\) es el conjunto de toda la información disponible hasta \(t-1\).

Para un VAR(\(p\)) bivariado tenemos:

\[ y_{1t} = c_1 + \sum_{i=1}^{p} \phi_{11,i} y_{1,t-i} + \sum_{i=1}^{p} \phi_{12,i} y_{2,t-i} + \varepsilon_{1t},\\ y_{2t} = c_2 + \sum_{i=1}^{p} \phi_{21,i} y_{1,t-i} + \sum_{i=1}^{p} \phi_{22,i} y_{2,t-i} + \varepsilon_{2t}. \]

Entonces:

• \(y_{2t}\) no causa a \(y_{1t}\) en el sentido de Granger si y solo si \(\phi_{12,1}=\cdots=\phi_{12,p}=0\).
• \(y_{1t}\) no causa a \(y_{2t}\) en el sentido de Granger si y solo si \(\phi_{21,1}=\cdots=\phi_{21,p}=0\).

Podemos verificar la causalidad de Granger empíricamente. Para esto, planteamos:

\[ H_0: \phi_{ij,1} = \phi_{ij,2} = \cdots = \phi_{ij,p} = 0, \qquad i \ne j. \] Usamos la prueba F:

\[ F = \frac{(SSR_R - SSR_{UR}) / q}{SSR_{UR}/(T - k)}, \] donde:

• \(SSR_R\): suma de residuos al cuadrado del modelo restringido (sin rezagos de \(y_2\))
• \(SSR_{UR}\): suma de residuos del modelo no restringido
• \(q=p\): número de restricciones.

Estimamos un VAR(\(4\)) para el inicio y el fin de construcción de nuevas viviendas en los Estados Unidos.

## 
## Call:
## lm(formula = y ~ -1 + ., data = datamat)
## 
## Coefficients:
## starts.l1   comps.l1  starts.l2   comps.l2  starts.l3   comps.l3  starts.l4  
##  0.655743   0.007276   0.235371  -0.118989   0.138673  -0.021141   0.011877  
##  comps.l4      const  
##  0.002711   0.131220

## 
## Call:
## lm(formula = y ~ -1 + ., data = datamat)
## 
## Coefficients:
## starts.l1   comps.l1  starts.l2   comps.l2  starts.l3   comps.l3  starts.l4  
##   0.06802    0.20486    0.05966    0.21775    0.03358    0.10504    0.09309  
##  comps.l4      const  
##   0.18231    0.04524

Para predecir el futuro usando un modelo VAR(\(p\)), usamos exactamente las mismas técnicas que para los modelos ARMA(\(p\),\(q\)). Usamos la muestra 1968:01-1991:12 para estimar el modelo y pronosticamos el periodo 1992:01-1996:06. Para el inicio de construcción de nuevas viviendas tenemos:

Para el fin de construcción de nuevas viviendas tenemos:

La siguiente tabla muestra el resultado del test de Granger para el inicio de construcción de viviendas:

## $Granger
## 
##  Granger causality H0: comps do not Granger-cause starts
## 
## data:  VAR object var_fit
## F-Test = 2.3066, df1 = 4, df2 = 658, p-value = 0.05685
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: comps and starts
## 
## data:  VAR object var_fit
## Chi-squared = 11.478, df = 1, p-value = 0.0007043

El test de Granger para el fin de construcción de viviendas:

## $Granger
## 
##  Granger causality H0: starts do not Granger-cause comps
## 
## data:  VAR object var_fit
## F-Test = 32.125, df1 = 4, df2 = 658, p-value < 2.2e-16
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: starts and comps
## 
## data:  VAR object var_fit
## Chi-squared = 11.478, df = 1, p-value = 0.0007043

Las funciones de impulso-respuesta para nuetro ejemplo son:

Otra forma de caracterizar la dinámica de un VAR —muy relacionada a las funciones impulso-respuesta, es la descomposición de la varianza. Esta nos permite analizar el error de predicción en cada variable a un determinado horizonte de tiempo. Este error se atribuye a la incertidumbre acerca de la evolución futura en cada una de las variables.