Dickey Fuller simulation in R

Question

so here's my problem, I have to simulate a random walk without drift:

$$ y_t = y_{t-1} + \epsilon_t $$ (with $\epsilon_t$ being a Gaussian white noise). Then I have to estimate this equation: $$ y_t = \alpha + \rho y_{t-1} + u_t $$ and finally test the null hypothesis (unit root test) $ \quad H_0 : \rho = 1 \quad $ vs $ \quad H_1 : \rho < 1$. I have to repeat this operation 10,000 times to get the critical values of the Dickey Fuller distribution.

Now, I did it in the following way (in R), and the resulting critical values that I get are not correct.

> dfcrit <- function(nobs){ e <- rnorm(nobs, 0, 1) yt <- as.matrix(cumsum(e)) x <- cbind(1:1, yt) LS.est <- (solve( t(x) %*% x )) %*% t(x) %*% yt residuals <- yt - x %*% LS.est s.squared <- (t(residuals) %*% residuals)/(nobs - (ncol(x))) se.matrix <- (solve(t(x) %*% x)) * s.squared[1] t.test <- (LS.est[2] - 1)/sqrt((se.matrix[2,2])) t.test }

> DF.100 <- rep(NA, 10000)

> for(i in 1:10000){ DF.100[i] <- dfcrit(100) } > sort(DF.100)

> quantile(DF.100, c(0.01, 0.05), na.rm = TRUE)

So what I think my problem is, in the computation of my LS.estimators (so $\alpha$ and $\rho$ I didn't use $y_{t-1}$ but $y_t$. But the problem is I don't know how to use $y_{t-1}$ in R since I'm still at a beginner... I tried by taking yt[-1] but then the matrices have different sizes and I cannot multiply them anymore.

I need your help for this! Thank you in advance, and sorry for the bad formating of my code :)

Answer 1

Is this what you're after?

install.packages(urca)
library(urca)
nobs = 100
DF.100 <- rep(NA, 10000)
for (i in 1:10000){
yt <- cumsum(rnorm(nobs))
DF.100[i] <- ur.df(yt,lags = 0,type = "drift")@testreg$coefficients[2,3]
}

quantile(DF.100,probs = c(0.01,0.05))

The results will be slightly different each time because you're only drawing 10000 samples from the DF distribution, but should be pretty close to each other and the theoretical quantiles.

To do an OLS fit yourself:

nobs = 100
DF.100 <- rep(NA, 10000)
for (i in 1:10000){

yt <- cumsum(rnorm(nobs))
X <- cbind(1,yt[-nobs]) # remove last element of yt and form predictor 
# matrix

beta <- solve(t(X) %*% X) %*% t(X) %*% yt[-1] # get ols estimate, note 
# observation vector is yt with first element removed.

ssr <- sum((yt[-1] - X %*% beta)^2)

se <- sqrt(ssr / (nobs-3)) / sqrt(sum((X[,2] - mean(X[,2]))^2)) # note
# (nobs - 3) instead of (nobs - 2) as we're only using (nobs - 1) data 
# points 

DF.100[i] <- (beta[2] - 1) / se
}

quantile(DF.100,probs = c(0.01,0.05))