Could you please explain in simple terms how the table of critical values for the augmented Dickey–Fuller (ADF) test is created?
1 Answer
I am not sure an easy answer is possible here.
As can be found in many textbooks, the limiting null distribution of the "Dickey-Fuller t-statistic" is that of a nonstandard random variable which may be expressed through a functional of a Brownian motion $W$.
Denote by $\hat{\rho}_T$ the OLS estimate of a regression of $y_t$ on $y_{t-1}$, and by $t_T$ the standard t-ratio for the null that $\rho=1$.
In the simplest case without constant or trend in the test regression, we have
\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*} This random variable has no easy expression for its density or cdf, but can be simulated, noting that, for a suitable distribution of $u$ like the standard normal, $1/\sqrt{T}\sum_{t=1}^{[sT]}u_t$ will behave like $W(s)$ for $T$ "large", where $[sT]$ denotes the integer part of $sT$.
Hence, the DF distribution can be simulated as follows:
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)))
}
The resulting (magenta) simulated (kernel-density estimated) distribution is given here, with the green standard normal for comparison:
plot(density(DFstats),lwd=2,col=c("deeppink2"))
xax <- seq(-4,4,by=.1)
lines(xax,dnorm(xax),lwd=2,col=c("chartreuse4"))
We see that the DF distribution is shifted to the left relative to the standard normal, and skewed.
The critical values are then, as usual, the quantiles of the null distribution of the test statistic, and can be retrieved via (note that we - typically - perform the DF test against the left-tailed alternative of stationarity that $|\rho|<1$)
CriticalValues <- sort(DFstats)[c(0.01,0.05,0.1)*reps]
> CriticalValues
[1] -2.571179 -1.943025 -1.611253
That is, these are the values such that 1%, 5% or 10% of the probability mass is to the left (=the rejection region) of them, producing the desired rejection rates.
Graphically (zooming in on the more relevant left part of the distribution):
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)
Given that these are simulated critical values, these may thus differ slightly from those reported in published tables.
EDIT: In response to the comment below, for the case with constant we would modify the code as follows, following Dickey-Fuller unit root test with no trend and supressed constant in 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)))
}
In particular, the line W_mu <- W - mean(W)
creates the demeaned Brownian motion $W^\mu(r)=W(r)-\int W(s)ds$. Its first element W_mu[1]^2
corresponds to the initial entry of that demeaned Brownian motion, $W^\mu(0)^2$, from the link.
The link also deals with the expression for the trend case, which may therefore be simulated via
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])
-
1$\begingroup$ @Michal, OK, see the edit. $\endgroup$ Commented Sep 21, 2016 at 6:38
-
1$\begingroup$ I am sorry, the edit was of course for the case with constant and no trend, not with trend. I corrected that. I also tried to further explain the code. $\endgroup$ Commented Sep 21, 2016 at 11:11
-
1$\begingroup$ You are right about the right critical values, but for me it does produce these values:
> (CriticalValues <- sort(DFstats)[c(0.01,0.05,0.1)*reps])
yields[1] -3.412332 -2.868623 -2.569988
$\endgroup$ Commented Sep 21, 2016 at 12:53 -
1
-
1$\begingroup$ @confused, unlike the t-distribution, the DF distribution is not charcterized by any degree of freedom parameter. $\endgroup$ Commented Jul 28, 2020 at 4:16