#########################################################################
# Exploration of the Dickey-Fuller Test for Non-stationary AR(1) Model
##########################################################################


#####################
#To get started, let's write a simple function to compute the 
#Dickey-Fuller statistic
###############################

DF<-function(x){
ar1.out<-arima.mle(x, model=list(ar=0));
temp<-(ar1.out$model$ar-1)/sqrt(ar1.out$var.coef);
return(as.numeric(temp))
}

##################################
# The function as.numeric() is used to get a pure number out of DF()
# Otherwise DF would return a 1-by-1 matrix
# This is not a big deal, but it makes the out-put look simpler.
#
# Since arima.mle is a bit "fussy" it is useful to note
# error messages from arima.mle will flow through to 
# the user of DF.
# 
# Here by "fussy" I mean that the optimizer may not 
# converge, the assumptions of the esitimator may be violated, etc.
##################################

###############################
# To see if DF works. We first simulate a Gaussian Random walk.
# This is the nicest unit-root model one can have.
# With luch DF should "fail to reject" the unit root hypothesis.
###############################

set.seed(101)
GRW<-cumsum(rnorm(100))
tsplot(GRW)

#Now get the value of the test statistic

DF.out<-DF(GRW)

# Now use punitroot() to look at the p-value
# You may want to review what the text says about punitroot()
# For the moment: "nc" means no constant and "t" means Dickey=Fuller

punitroot(DF.out,trend="nc", statistic="t")

#########################################
# As suspected, we fail to reject.
# What if we  give a longer series?
##########################################


set.seed(101)
GRW1000<-cumsum(rnorm(1000))
tsplot(GRW1000)

DF.out<-DF(GRW1000)

punitroot(DF.out,trend="nc", statistic="t")

####################################################################
# Again we fail to reject. This is just as it should be.
# What if we introduce some mild heteroscedasticity?
####################################

w<-1:200
w<-(1.02)^w  #weights that grow by 2% each step
set.seed(101)
ERW<-cumsum(w*rnorm(200))
tsplot(ERW)
punitroot(DF(ERW), trend="nc", statistic="t")

#####################################################
# DF is still on track. It is not fooled by a 
# modest violation of the assumptions.
# Let's try another violation --- an AR(1) with
# a phi BIGGER than one.
# First lets try arima.sim()
######################################################

ARns<-arima.sim(model=list(ar=1.02), n=200)

##############################################
# As expected this fails.  We need to use bare hands.
################################################

tempARns<-rep(0,300) #We'll generate 300 observations and just use the last 200
set.seed(101)
for(j in 1:300){tempARns[j]<-(1.002)*tempARns[j-1]+rnorm(1)}
ARns<-tempARns[101:300]

tsplot(ARns)

punitroot(DF(ARns), trend="nc", statistic="t")

############################################
# Again, we fail to reject the "unit root" hypothesis.
# The result this time is a little more curious, but it
# reflects the one-sided nature of the test.
# In essence, the null hypothesis of DF is that
# phi is equal to OR BIGGER THAN 1.
############################################

tempARns<-rep(0,300) #We'll generate 300 observations and just use the last 200
set.seed(101)
for(j in 1:300){tempARns[j]<-(1.01)*tempARns[j-1]+rnorm(1)}

#That is a HUGE 1% perday kicker.

ARns1<-tempARns[101:300]

tsplot(ARns1)

punitroot(DF(ARns1), trend="nc", statistic="t")

###################
# OOPS, arima.mle broke. A good implementation of the 
# DF test will need to work around this problem.
# This is done in the S-Plus code unitroot()
# which gives an implementation of the AUGMENTED DF test.  
#####################################################################