¿Cómo se calcula la tabla de valores críticos de la prueba Dickey-Fuller aumentada (ADF)?

Question

¿Podría explicar de forma sencilla cómo se crea la tabla de valores críticos para la prueba Dickey-Fuller aumentada (ADF)?

Answer 1

No estoy seguro de que sea posible una respuesta fácil en este caso.

Como se puede encontrar en muchos libros de texto, la distribución nula limitante del "estadístico t de Dickey-Fuller" es la de una variable aleatoria no estándar que se puede expresar a través de un funcional de un movimiento browniano $W$ .

Denota por $\hat{\rho}_T$ la estimación OLS de una regresión de $y_t$ en $y_{t-1}$ y por $t_T$ el cociente t estándar para la hipótesis nula de que $\rho=1$ .

En el caso más sencillo, sin constante ni tendencia en la regresión de prueba, tenemos

\begin{eqnarray*} t_T&=&\frac{\hat{\rho}_T-1}{s.e.(\hat{\rho}_T)}\\&=&\frac{T(\hat{\rho}_T-1)}{\{s^2_T\}^{1/2}}\left\{T^{-2}\sum_{t=1}^Ty_{t-1}^2\right\}^{1/2}\\ &\Rightarrow&\frac{1/2\{W(1)^2-1\}}{\int_0^1W(r)^2d r}\frac{1}{\sigma}\left\{\sigma^2\int_0^1W(r)^2dr\right\}^{1/2}\\ &=&\frac{W(1)^2-1}{2 \left\{\int_0^1W(r)^2dr\right\}^{1/2}} \end{eqnarray*} Esta variable aleatoria no tiene una expresión fácil para su densidad o cdf, pero se puede simular, observando que, para una distribución adecuada de $u$ como el estándar normal, $1/\sqrt{T}\sum_{t=1}^{[sT]}u_t$ se comportará como $W(s)$ para $T$ "grande", donde $[sT]$ denota la parte entera de $sT$ .

Por lo tanto, la distribución DF puede simularse como sigue:

T <- 5000
reps <- 50000

DFstats <- rep(NA,reps)

for (i in 1:reps){
  u <- rnorm(T)
  W <- 1/sqrt(T)*cumsum(u)
  DFstats[i] <- (W[T]^2-1)/(2*sqrt(mean(W^2)))
}

La distribución simulada (magenta) resultante (densidad del núcleo estimada) se da aquí, con la normal estándar verde para la comparación:

plot(density(DFstats),lwd=2,col=c("deeppink2"))
xax <- seq(-4,4,by=.1)
lines(xax,dnorm(xax),lwd=2,col=c("chartreuse4"))

Vemos que la distribución DF está desplazada a la izquierda respecto a la normal estándar, y sesgada.

Los valores críticos son entonces, como es habitual, los cuantiles de la distribución nula del estadístico de la prueba, y pueden obtenerse a través de (nótese que -típicamente- realizamos la prueba DF contra la alternativa de cola izquierda de estacionariedad que $|\rho|<1$ )

CriticalValues <- sort(DFstats)[c(0.01,0.05,0.1)*reps]
> CriticalValues
[1] -2.571179 -1.943025 -1.611253

Es decir, son los valores tales que el 1%, el 5% o el 10% de la masa de probabilidad está a la izquierda (=la región de rechazo) de ellos, produciendo los índices de rechazo deseados.

Gráficamente (ampliando la parte izquierda más relevante de la distribución):

plot(DFdensity,lwd=2,col=c("deeppink2"), xlim=c(-4,0))

xshade1 <- DFdensity$x[DFdensity$x <= CriticalValues[1]]
yshade1 <- DFdensity$y[DFdensity$x <= CriticalValues[1]]
polygon(c(xshade1[1],xshade1,CriticalValues[1]),c(0,yshade1,0),col="deeppink2", border = "deeppink2", lwd=2)

xshade2 <- DFdensity$x[DFdensity$x > CriticalValues[1] & DFdensity$x <= CriticalValues[2]]
yshade2 <- DFdensity$y[DFdensity$x > CriticalValues[1] & DFdensity$x <= CriticalValues[2]]
polygon(c(xshade2[1],xshade2,CriticalValues[2]),c(0,yshade2,0),col="deeppink4", border = "deeppink4", lwd=2)

xshade3 <- DFdensity$x[DFdensity$x > CriticalValues[2] & DFdensity$x <= CriticalValues[3]]
yshade3 <- DFdensity$y[DFdensity$x > CriticalValues[2] & DFdensity$x <= CriticalValues[3]]
polygon(c(xshade3[1],xshade3,CriticalValues[3]),c(0,yshade3,0),col="darkorchid4", border = "darkorchid4", lwd=2) 

Dado que se trata de valores críticos simulados, éstos pueden diferir ligeramente de los indicados en las tablas publicadas.

EDIT: En respuesta al comentario de abajo, para el caso con constante modificaríamos el código de la siguiente manera, siguiendo Prueba de raíz unitaria de Dickey-Fuller sin tendencia y con constante suprimida en Stata

for (i in 1:reps){
  u <- rnorm(T)
  W <- 1/sqrt(T)*cumsum(u)
  W_mu <- W - mean(W)
  DFstats[i] <- (W_mu[T]^2-W_mu[1]^2-1)/(2*sqrt(mean(W_mu^2)))
}

En particular, la línea W_mu <- W - mean(W) crea el movimiento browniano degradado $W^\mu(r)=W(r)-\int W(s)ds$ . Su primer elemento W_mu[1]^2 corresponde a la entrada inicial de ese movimiento browniano degradado, $W^\mu(0)^2$ , desde el enlace.

El enlace también se ocupa de la expresión para el caso de la tendencia, que por lo tanto se puede simular a través de

s <- seq(0,1,length.out = T)

for (i in 1:reps){
     u <- rnorm(T)
     W <- 1/sqrt(T)*cumsum(u)
     W_tau <- W - (4-6*s)*mean(W) - (12*s-6)*mean(s*W)
     DFstats[i] <- (W_tau[T]^2-W_tau[1]^2-1)/(2*sqrt(mean(W_tau^2)))
}
(CriticalValues <- sort(DFstats)[c(0.01,0.05,0.1)*reps])