# Written by Jin-Lung Lin on May 1, 2012
# This program finds the quantiles of ADF with simulation
# First, simulates the random walk process x_t = x_(t-1) + e_t
# regress y_t on y_(t-1} to obtain OLS coefficient, b[i]
# Repeat the two steps above for say, 5000 times
# Find the critical values  for b and t-stat  

rm(list=ls())
library(TSA)
nrep=5000
nobs=5000
nobs2=nobs*nobs
x=matrix(0,nobs,nrep)
xtx=numeric(nrep)
xty=numeric(nrep)
xte=numeric(nrep)
Tb1=numeric(nrep)
bhat =numeric(nrep)
sighat=numeric(nrep)
t_bhat=numeric(nrep)

set.seed(123456789)
for (i in 1:nrep)
{
x[,i]=cumsum(rnorm(nobs))
y0=x[2:nobs,i]
y1=x[1:(nobs-1),i]

xtx[i]=t(y1)%*%y1
xty[i]=t(y1)%*%y0
bhat[i]=xty[i]/xtx[i]
Tb1[i]=nobs*(bhat[i]-1)
sighat[i]=sum( (y0-bhat[i]*y1)*(y0-bhat[i]*y1) )/nobs
t_bhat[i]=sqrt(xtx[i])*(bhat[i]-1)/( sqrt(sighat[i]) )  
}
p <- c(0.1,0.25,0.5,1,2,5,10,50)/100
quantile(t_bhat,p)
quantile(Tb1,p)
"mean of T(b-1) =";  mean(Tb1)
"mean of t_stat ="; mean(t_bhat)
skewness(Tb1)
skewness(t_bhat)


op=par(mfrow=c(2,1))
hist(t_bhat,main="Histogram of t-stat for ADF test",breaks=100,prob=TRUE)
abline(v=0,col="red")
hist(Tb1,main="Histogram of T*(b-1) for ADF test",breaks=100,prob=TRUE)
abline(v=0,col="red")
par(op)