# Monte Carlo study of the behavior of the Dickey-Fuller ratio 
set.seed(2135)
n = 100
sig2 = 1
M = 10000
DF = rep(0,M)
for (j in 1:M){
a = rnorm(n,0,sqrt(sig2))
y = rep(0,1)
for (t in 2:n)
  y[t] = y[t-1]+a[t]

phi.mle = sum(y[2:n]*y[1:(n-1)])/sum(y[1:(n-1)]*y[1:(n-1)])
sig2.mle = sum((y[2:n]-phi.mle*y[1:(n-1)])^2)/(n-1)
se.mle = sqrt(sig2.mle/sum(y[1:(n-1)]*y[1:(n-1)]))

DF[j] = (phi.mle-1)/se.mle
}
DF = sort(DF)

plot(density(DF),main="",xlab="DF ratio")
lines(DF,dnorm(DF),col=2)