Introducción

Análisis de series temporales no estacionarias, con raíces unitarias.

Sabemos que un proceso ARMA (p,q) es un proceso de memoria corta, porque los coeficientes \(\psi\) en la representación como un proceso lineal infinito

\[ x_t = \sum_{j=0}^{\infty}\psi_{j} w_{t-j} \]

que se obtienen al resolver \(\phi_{(z)}\psi_{(z)} = \theta_{(z)}\) decaen exponencialmente rápido. Cuando esto ocurre, se cumple que la ACF de un proceso de memoria corta satisface que \(\rho_{(h)} \to 0\)  exponencialmente rápido, cuando \(h \to \infty\).

Cuando ocurre lo contrario, es decir que la ACF muestral decae lentamente, Shumway et al. (2000) recomienda diferenciar la serie de tiempo para que parezca estacionaria.

Pero hay una consideración que hay que tener, la primera diferencia \[ \nabla x_{t} = x_{t} - x_{t-1} = (1- B)x_{t} \] puede ser una modificación demasiado severa, en el sentido de que el modelo no estacionario podría representar una sobre diferenciación del proceso original.

Desarrollo del problema

Para analizar lo que ocurre, supongamos el proceso estacionario AR(1) \[ x_{t} = \phi x_{t-1} + w_{t} \] \(x_{t}\) es estacionaria

\(\phi\) es una constante

\(w_{t}\) es wn(0,\(\sigma_{w}^2\))

Apliquemos la primera diferencia \(\nabla x_{t}\) \[ x_{t} - x_{t-1}= \phi x_{t-1} + w_{t} - \phi x_{t-2} - w_{t-1} \] \[ x_{t} - x_{t-1}= \phi (x_{t-1} - x_{t-2}) + w_{t} - w_{t-1} \] \[ y_{t} = \phi y_{t-1} + w_{t} - w_{t-1} \] donde \(\nabla x_{t} = y_{t}\) y \(\nabla x_{t-1} = y_{t-1}\).

Al realizar esta diferenciación, se presentarán problemas de correlación e invertibilidad. Sabemos que el proceso AR(1) \(x_{t}\) es causal porque lo podemos expresar como

\[ x_{t}= \sum_{j=0}^{\infty}\phi^{j} w_{t-j}, \mbox{ } \left|{\phi}\right| < 1 \]

pero \(y_{t}\) que es el proceso diferenciado de \(x_{t}\), es ahora un proceso no invertible, veamos por qué \[ y_{t} - \phi y_{t-1} = w_{t} - w_{t-1} \] \[ (1- \phi B)y_{t} = (1 - B)w_{t} \] en términos de Z \[ (1- \phi z)y_{t} = (1 - z)w_{t} \] \(\textbf{Causalidad}\) \[ \Phi(z) = 1 - \phi z \] podemos ver que z depende de \(\phi\) y el valor que hace que el polinomio sea cero (raiz) es \[ z = \frac{1}{\phi} \] \(y_{t}\) será causal si dependiendo del valor de \(\phi\) hace que el \(\left|{z}\right| > 1\) y esté fuera del circulo unitario.

\(\textbf{Invertibilidad}\)

\[ \Theta (z) = 1 -z \]

podemos ver que z no depende de \(\phi\)  y el valor que hace que el polinomio sea cero (raiz) es \(z = 1\),  lo que llamamos una raiz unitaria.

\(y_{t}\) no es invertible porque no se cumple que \(\left|{z}\right| > 1\)  lo cual hace que no esté fuera del círculo unitario.

La condición de invertibilidad de un modelo ARMA (que el valor absoluto de las raíces de \(\Theta(z)\) estén fuera del círculo unitario), garantiza que los pesos \(\pi_0,\pi_1,\pi_2,....\) del polinomio de orden infinito \(\Pi(B)\), satisfacen la condición que \(\sum_{j=0}^{\infty}\left|{\pi_j}\right| < \infty\) lo que implica que cuando se escribe

\[ \Pi(B)y_{t} = w_{t} \]

\(y_{t}\) es un proceso estacionario tal que la correlación parcial de  \(y_{t}\) y  \(y_{t-k}\) tienda a cero a medida que el retardo \(k\) aumenta. Como no es invertible, lo argumentado antes, no se cumple.

En general, la no invertibilidad indica un uso incorrecto por exceso del operador \(\nabla\) (sobre diferenciación), lo que puede llevar a caracterizar erróneamente la tendencia de un proceso.

El tema ahora es determinar si \(\phi=1\), que haría que \(x_{t}\)  sea un proceso AR(1) no estacionario conocido como caminata aleatoria (su función de autocovarianza depende de t  \(Cov(x_s,x_t) = min (s,t)\sigma_{w}^2\),  y su diferencia de orden (1)  \(y_{t}\) será un proceso estacionario iid\((0,\sigma_{w}^2)\),  en cuyo caso el proceso \(\nabla x_{t} = y_{t}\)  será ruido blanco de acuerdo a Mauricio (2007).

Cuando un proceso tiene raíces unitarias, los estimadores (MCO, MLE, etc) y los estadísticos de contraste habituales (t, F, etc) no siguen distribuciones estándar ni siquiera asintóticamente. Por lo cual en estos modelos con raíces unitarias, cualquier inferencia realizada sobre la base de dichas distribuciones típicas (t, F, normal, etc) es incorrecta y puede llevar a conclusiones erróneas.

Unit Root Test

Unit root nos permite definir un estadístico que me dé una buena aproximación para hacer inferencias correctas en presencia de procesos con raíces unitarias. Veamos un ejemplo:

Una caminata aleatoria constituye el caso más sencillo de lo que suele denominarse un proceso integrado de orden 1 ARIMA (0,1,0) (es decir, un proceso no estacionario R.W cuya diferencia regular de orden 1, es un proceso estacionario).

Unit root me permite probar las siguientes hipótesis:

\[ \mbox{Si } \mbox{ } H_0: \phi =1 \mbox{ } \mbox{ , $x_{t}$ no es causal y es caminata aleatoría pura} \] \[ \mbox{o } \mbox{ } H_1: \left|{\phi}\right| < 1 \mbox{ } \mbox{ , $x_{t}$ es causal, es estacionario} \]

es decir

\[ H_0: \phi =1 \mbox{ } \mbox{ versus } \mbox{ } \mbox{ } H_1: \left|{\phi}\right| < 1 \]

Unit root me da un procedimiento para probar \(H_0\) y \(H_1\). Un procedimiento apropiado sería estimar \(\phi\) con \(\hat{\phi}\)  y definir el estadístico \((\hat{\phi} - 1)\),  apropiadamente normalizado, con la esperanza que el estadístico de prueba sea asintóticamente normal.

\(\hat{\phi}\) es un estimador óptimo y se calcula como lo vimos en la sección de Estimación del libro Shumway et al. (2000). El problema es que el procedimiento visto no se puede aplicar en un proceso no estacionario, que sería lo que ocurre si \(H_0:\) no se rechaza, es decir \(\phi=1\). Sin embargo al límite de la estacionaridad, se pueden utilizar la formula hallada para el estimador de \(\phi\).

Para examinar el estadístico propuesto \((\hat{\phi} - 1)\) bajo la \(H_0: \phi = 1\), es decir el modelo es una caminata aleatoria pura, lo cual implica que

\[ x_{t}= \sum_{j=1}^{t}w_{j} \mbox{ } \mbox{ ,o } x_{t} = x_{t-1} + w_{t} \mbox{ } \mbox{ con } x_{0} = 0 \]

vamos a considerar calcular \(\hat{\phi}\) con los estimadores de mínimos cuadrados condicionales, cuya fórmula es

\[ \hat{\phi} = \frac{\sum_{t=2}^{n}(x_{t}-\overline{x_{2}})(x_{t-1}-\overline{x_{1}})}{\sum_{t=2}^{n}(x_{t-1}-\overline{x_{1}})^2} \]

Bajo \(H_0\), la \(E(x_{t}) = \mu_x = 0\), el estimador de mínimos cuadrados quedaría escrito

\[ \hat{\phi} = \frac{\sum_{t=1}^{n}(x_{t})(x_{t-1})}{\sum_{t=1}^{n}(x_{t-1})^2} \]

como sé que \(x_{t} = x_{t-1} + w_{t}\), al reemplazar

\[ \hat{\phi} = \frac{\sum_{t=1}^{n}(x_{t-1} + w_{t})(x_{t-1})}{\sum_{t=1}^{n}(x_{t-1})^2} = \frac{\sum_{t=1}^{n}(x_{t-1})^2 + \sum_{t=1}^{n}w_{t}x_{t-1}}{\sum_{t=1}^{n}(x_{t-1})^2} \] \[ \hat{\phi} = 1 - \frac{\sum_{t=1}^{n}w_{t}x_{t-1}}{\sum_{t=1}^{n}(x_{t-1})^2} \] \[ \hat{\phi} - 1 = \frac{\frac{1}{n\sigma_w^2}\sum_{t=1}^{n}w_{t}x_{t-1}}{\frac{1}{n\sigma_w^2}\sum_{t=1}^{n}(x_{t-1})^2} \]

vamos a desarrollar primero el numerador de la ecuación anterior. Como sé que \(x_{t} = x_{t-1} + w_{t}\), elevo ambos lados de la ecuación al cuadrado y obtengo \(x_{t}^2 = x_{t-1}^2 + 2x_{t-1}w_{t} + w_{t}^2\) así

\[ w_{t}x_{t-1} = \frac{1}{2}(x_{t}^2 - x_{t-1}^2 - w_{t}^2) \]

el numerador quedaría

\[ \frac{1}{n\sigma_w^2}\sum_{t=1}^{n}w_{t}x_{t-1} = \frac{1}{2n\sigma_w^2}\sum_{t=1}^{n}(x_{t}^2 - x_{t-1}^2) - \frac{1}{2n\sigma_w^2}\sum_{t=1}^{n}w_{t}^2 \] \[ = \frac{1}{2}\left(\frac{x_n^2}{n\sigma_w^2} - \frac{\sum_{t=1}^{n}w_{t}^2}{n\sigma_w^2}\right) \]

• Como \(x_n = \sum_{1}^{n}w_{t}\), tenemos que \(x_n \sim N(0,n\sigma^2)\).
• Así que \(\chi_1^2 = \frac{1}{n\sigma_w^2}x_n^2\) tiene una distribución chi - cuadrado con 1 grado de libertad.
• Además como \(w_{t}\) es ruido blanco gaussiano, entonces por la propiedad de consistencia, la varianza muestral \(S_{n}^2 = \frac{1}{n}\sum_{1}^{n}(w_{t} - \overline{w})^2\), sabiendo que la media \(= 0\), podemos decir que \(S_{n}^2\) converge en probabilidad a \(\sigma_w^2\). Entonces \(\frac{1}{n}\sum_{1}^{n}w_{t}^2 \xrightarrow {\;\; P \;\; } \sigma_w^2\) o lo que es lo mismo \(\frac{1}{n\sigma_w^2}\sum_{1}^{n}w_{t}^2 \xrightarrow {\;\; P \;\; } 1\).
• Consecuentemente cuando \(n \to \infty\),

\[ \frac{1}{n\sigma_w^2}\sum_{t=1}^{n}w_{t}x_{t-1} \xrightarrow {\;\; D \;\; } \frac{1}{2}(\chi_1^2 - 1). \]

Seguidamente analizaremos el denominador de la expresión hallada para el estadístico propuesto, con el fin de poder determinar la forma que podría asintóticamente distribuir. Para eso vamos a revisar primero la siguiente definición de Movimiento Browniano estándar.

Movimiento Browniano estándar

Un proceso en tiempo continuo \(\left\{W(t); t \geq 0\right\}\) se llama Movimiento Browniano estándar si satisface las siguientes \(3\) condiciones:

• [i] \(W(0) = 0\) y \(\sigma_w^2 =1\).
• [ii] \(\left\{W(t_2) - W(t_1), W(t_3) - W(t_2),.........W(t_n) - W(t_{n-1}) \right\}\) son independientes para cualquier colección de puntos, \(0 \leq t_1 \leq t_2 \cdot \cdot \cdot < t_n\) y el entero \(n > 2\). Esto significa que el proceso tiene incrementos independientes.
• [iii] \(W(t + \varDelta t) - W(t) \sim N(0,\varDelta t)\) para \(\varDelta t > 0\). Esto significa que para cualesquiera tiempos, el incremento tiene distribución normal de media cero y varianza igual al incremento. Además se desprende que los incrementos del proceso son estacionarios.

Una caminata aleatoria es el proceso en tiempo discreto correspondiente a un movimiento Browniano en tiempo continuo, en otras palabras un movimiento Browniano es una caminata aleatoria en tiempo continuo.

De acuerdo al Teorema Donsker’s, que a menudo se denomina el principio de invarianza o Teorema del límite funcional, y dado que el límite es el movimiento browniano en este caso, se denomina un teorema del límite central funcional de acuerdo a Billingsley (2013). El teorema plantea que si tengo \(n\) variables aleatorias iid con media cero y varianza \(\sigma^2\) y \(X^{n}\) es una función aleatoria definida de una forma particular, entonces \(X^{n} \to W\).

El denominador de la expresión que estamos analizando, utiliza este teorema, que en términos mas explícitos dice si \(\xi_1, \xi_2, \xi_3, ......\xi_n\) es una secuencia de variables aleatorias iid con media cero y varianza \(1\), entonces para \(0 \leq t \leq 1\), el proceso en tiempo continuo se denota por

\[ S_n(t) = \frac{1}{\sqrt{n}}\sum_{j=1}^{\|nt\|}\xi_j \xrightarrow {\;\; D \;\; } W(t) \]

podemos ver a \(\xi_j\) como los desplazamientos independientes en tiempos sucesivos, con \(n \to \infty\), donde \(\|\mbox{ } \|\) es la mayor función entera y \(W(t)\) es un movimiento Browniano estándar entre \([0,1]\).

Recordemos que bajo la hipótesis nula \(H_0: \phi = 1\), el proceso es una caminata aleatoria \(x_{t}= \sum_{j=1}^{t}w_{j}\), entonces en un valor \(s\),  \(x_{s} = w_1 + w_2 + w_3 + ......w_s\)  distribuirá   \(\sim N(0,s\sigma_w^2)\) y basados en la distribución del proceso de tiempo continuo definido anteriormente, tenemos que \(\frac{x_s}{\sigma_w\sqrt{n}} \xrightarrow {\;\; D \;\; } W(s)\).

A partir de este hecho se puede demostrar que el denominador

\[ \frac{1}{n\sigma_w^2}\sum_{t=1}^{n}(x_{t-1})^2 = \sum_{t=1}^{n}\left(\frac{x_{t-1}}{\sigma_w\sqrt{n}}\right)^2 \]

cuando \((n \to \infty)\)

\[ \sum_{t=1}^{n}\left(\frac{x_{t-1}}{\sigma_w\sqrt{n}}\right)^2\frac{1}{n} \xrightarrow {\;\; D \;\; } \int_{0}^{1}W^{2}(t)dt \]

debo ajustar el estimador en un factor \(n^{-1}\), para lo cual multiplico a ambos lados de la ecuación por \(n\), quedando

\[ n(\hat{\phi} - 1) = n\frac{\frac{1}{n\sigma_w^2}\sum_{t=1}^{n}w_{t}x_{t-1}}{\frac{1}{n\sigma_w^2}\sum_{t=1}^{n}(x_{t-1})^2}= \frac{\frac{1}{n\sigma_w^2}\sum_{t=1}^{n}w_{t}x_{t-1}}{\frac{1}{n^2\sigma_w^2}\sum_{t=1}^{n}(x_{t-1})^2} \]

\[ n(\hat{\phi} - 1) \xrightarrow {\;\; D \;\; } \frac{\frac{1}{2}(\chi_1^2 - 1)}{\int_{0}^{1}W^{2}(t)dt} \]

Al estadístico de prueba \(n(\hat{\phi} - 1)\) se le conoce como estadístico de raiz unitaria o estadístico Dickey-Fuller DF, Fuller (2009). Aunque el actual estadístico de prueba DF se normaliza un poco diferente.

Debido a que la distribución del estadístico de prueba no tiene una forma cerrada, los cuantiles de la distribución se deben calcular por aproximación numérica o simulación.

Así hemos visto como el DF se estableció para un modelo básico de la forma \(x_{t} = \phi x_{t-1} + w_{t}\), que al diferenciarlo \(\nabla x_{t} = (\phi-1)x_{t-1} + w_{t} = \gamma x_{t-1} + w_{t}\) y las hipótesis se podrían plantear en términos de \(\gamma\) así:

\[ H_0: \gamma =0 \mbox{ } \mbox{ versus } \mbox{ } \mbox{ } H_1: \gamma < 0 \]

probando \(H_0\) haciendo una regresión de \(\nabla x_{t}\) en \(x_{t-1}\). De aquí sacaron un estadístico de WALD y derivaron su distribución límite.

La prueba se amplió para adaptarse a modelos \(AR(p)\), donde se prueba \(H_0: \gamma =0\) estimando \(\gamma\) en la regresión de \(\nabla x_{t}\) en \(x_{t-1},\nabla x_{t-1},.......\nabla x_{t-p+1}\) y se formó una prueba de WALD basada en \(t_y = \hat{\gamma}/se(\hat{\gamma})\). Está prueba conduce a la llamada prueba Dickey-Fuller aumentada ADF. Escoger el \(p\) es crucial.

Para modelos ARMA(p,q) el test ADF asume que \(p\) es lo suficientemente grande para capturar la estructura de correlación esencial. Existe la prueba de Phillips-Perron (PP), que se diferencia de las pruebas ADF principalmente en cómo tratan la correlación serial y la heterocedasticidad en los errores.

El modelo se puede extender al incluirle una constante o una tendencia no estocástica. Por ejemplo consideremos el modelo

\[ x_t = \beta_0 + \beta_{1}t + \phi x_{t-1}+ w_{t} \]

si asumimos que \(\beta_1 = 0\) bajo la hipótesis nula que \(\phi = 1\), el proceso será una caminata aleatoria con deriva \(\beta_0\). Bajo la hipótesis alternativa el proceso es AR(1) causal con media \(\beta_0(1-\phi)\).

Si no podemos asumir que \(\beta_1 = 0\) entonces el interés es probar la nula de manera conjunta \((\beta_1,\phi) = (0,1)\), versus la alternativa de que \(\beta_1 \neq 0\) y \(\left|{\phi}\right| < 1\).

La hipótesis nula es que el proceso es una caminata aleatoria con deriva, versus la alternativa de que es un proceso con tendencia estacionaria.

Ejemplo en R

Se hará la prueba a la serie VARVE del paquete astsa

• Varves (sedimentary deposits) can be used as proxies for paleoclimatic parameters, such as temperature.

tsplot(varve, main="varve", ylab="")

y la transformación logaritmica

tsplot(log(varve), main="log(varve)", ylab="" )

adf.test(log(varve), k=0) # DF test

## Warning in adf.test(log(varve), k = 0): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  log(varve)
## Dickey-Fuller = -12.857, Lag order = 0, p-value = 0.01
## alternative hypothesis: stationary

adf.test(log(varve)) # ADF test

## 
##  Augmented Dickey-Fuller Test
## 
## data:  log(varve)
## Dickey-Fuller = -3.5166, Lag order = 8, p-value = 0.04071
## alternative hypothesis: stationary

pp.test(log(varve)) # PP test

## Warning in pp.test(log(varve)): p-value smaller than printed p-value

## 
##  Phillips-Perron Unit Root Test
## 
## data:  log(varve)
## Dickey-Fuller Z(alpha) = -304.54, Truncation lag parameter = 6, p-value
## = 0.01
## alternative hypothesis: stationary

En cada prueba se rechaza la hipótesis nula que la transformación logarítmica de la serie varve (log(varve)) tiene raices unitarias.