# Dickey Fuller simulation in R

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 :)

• Are these two different problems? How is $u_t$ defined? – Michael R. Chernick Jun 3 '17 at 20:27
• The first equation specifies how $y_t$ is generated, while the second equation is regression equation that I must estimate, and $u_t$ is the residuals from that regression. – Olivier Jun 3 '17 at 20:50
• Isn't nobs = ncol(x)? – Will Jun 3 '17 at 21:56
• Also your function doesn't return anything... – Will Jun 3 '17 at 22:04
• No nobs it's the number of time I simulate $y_t$. And for the return, I forgot (I'm not comfortable with R yet), I added t.test in the end. – Olivier Jun 4 '17 at 6:44

install.packages(urca)
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))  • Right, the thing is that I'm trying to achieve this by estimating the least squares my self (in matrix form). In general if my regression is$y_t = \alpha + \beta x_t + \epsilon_t $The OLS regression is$ b = (X'X)^{-1} X' y$where$X$is a matrix with a column of 1's for the intercept and with$y_t\$ as the second column (the observations). – Olivier Jun 4 '17 at 6:40