# unitroot.ssc					script file for examples used in 
#									Unit Root chapter
#
# created: November 28, 2001
# revised: March 30, 2002
#
####################################################################

# simulate TS and DS data
set.seed(201)
e = rnorm(250)
z.DS = cumsum(e)
z.TS = arima.sim(model=list(ar=0.75),
		 innov=e)
y.DS = 5 + 0.1*seq(250)+ z.DS
y.TS = 5+ 0.1*seq(250) + z.TS
tsplot(y.DS,y.TS)
legend(25,30,legend=c("I(1)","I(0)"),lty=c(1,3))

#
# simulate functions of Wiener processes
#

wiener = function(nobs) {
	e = rnorm(nobs)
	y = cumsum(e)
	ym1 = y[1:(nobs-1)]
	intW2 = nobs^(-2) * sum(ym1^2)
	intWdW = nobs^(-1) * sum(ym1*e[2:nobs])
	ans = list(intW2=intW2,
				intWdW=intWdW)
	ans
}

#
# simulate DF distributions
#

set.seed(378)
nobs = 1000
nsim = 1000
NB = rep(0,nsim)
DF = rep(0,nsim)
for (i in 1:nsim) {
	BN.moments = wiener(nobs)
	NB[i] = BN.moments$intWdW/BN.moments$intW2
	DF[i] = BN.moments$intWdW/sqrt(BN.moments$intW2)
}

#
# compute empirical quantiles
#

quantile(DF,probs=c(0.01,0.05,0.1))
quantile(NB,probs=c(0.01,0.05,0.1))

#
# plot histograms and QQ plots of densities
#

par(mfrow=c(2,2))
hist(DF,main="Simulated DF distribution")
hist(NB,main="Simulated normalized bias")
plot(density(DF),type="l",
	main="Density estimate of simulated DF distribution",
	xlab="DF values",ylab="Density")
plot(density(NB),type="l",
	main="Density estimate of simulated normalized bias",
	xlab="Normalized bias",ylab="Density")
par(mfrow=c(1,1))

#
# Example: p-values and quantiles of DF and NB distributions using
# MacKinnon's programs
#

args(qunitroot)

# compute usual asymptotic critical values for t-stat and normalized bias

qunitroot(c(0.01,0.05,0.10),trend="nc",statistic="t")
qunitroot(c(0.01,0.05,0.10),trend="nc",statistic="n")

# compute critical values appropriate for sample of size 100

qunitroot(c(0.01,0.05,0.10),trend="nc",statistic="t",n.sample=100)
qunitroot(c(0.01,0.05,0.10),trend="nc",statistic="n",n.sample=100)

# compute p-values for specific values

args(punitroot)

punitroot(-1.645,trend="nc",statistic="t")
pnorm(-1.645)						# p-value from N(0,1)

#
# simulate data from trend cases I, II 
#

set.seed(201)
e = rnorm(250)
y1.H0 = cumsum(e)
y1.H1 = arima.sim(model=list(ar=0.75),
		 innov=e)
y2.H0 = 5 + y1.H0
y2.H1 = 5+ y1.H1
y3.H0 = 5 + 0.1*seq(250)+ y1.H0
y3.H1 = 5+ 0.1*seq(250) + y1.H1
		
par(mfrow=c(2,2))
# tsplot(y1.H0,main="Case I: I(1) data")
# tsplot(y1.H1,main="Case I: I(0) data")
tsplot(y2.H0,main="Case I: I(1) data")
tsplot(y2.H1,main="Case I: I(0) data")
tsplot(y3.H0,main="Case II: I(1) data")
tsplot(y3.H1,main="Case II: I(0) data")
par(mfrow=c(1,1))

# note: must add legend to above graph

#
# generate p-values and quantiles for different trend cases
#
qunitroot(c(0.01,0.05,0.10),trend="c",statistic="t",
n.sample=100)
qunitroot(c(0.01,0.05,0.10),trend="c",statistic="n",
n.sample=100)
punitroot(-1.645,trend="c",statistic="t",
n.sample=100)
punitroot(-1.645,trend="c",statistic="n",
n.sample=100)

qunitroot(c(0.01,0.05,0.10),trend="ct",statistic="t",
n.sample=100)
qunitroot(c(0.01,0.05,0.10),trend="ct",statistic="n",
n.sample=100)
punitroot(-1.645,trend="ct",statistic="t",
n.sample=100)
punitroot(-1.645,trend="ct",statistic="n",
n.sample=100)

#
# Example: unit root testing on exchange rates (no drift)
# see Meese and Diebold papers
#

#
# first do ADF tests
#

uscn.spot = lexrates.dat[,"USCNS"]
uscn.spot@title = "Log US/CN spot exchange rate"
par(mfrow=c(2,2))
plot.timeSeries(uscn.spot,
main="Log of US/CN spot exchange rate",
reference.grid=F)
xx = acf(uscn.spot)
plot.timeSeries(diff(uscn.spot),
main="First difference of log US/CN spot exchange rate",
reference.grid=F)
xx = acf(diff(uscn.spot))
par(mfrow=c(1,1))

args(unitroot)

# set max lag = 6

adft.out = unitroot(uscn.spot,trend="c",statistic="t",
method="adf",lags=6)
class(adft.out)
names(adft.out)
adft.out
summary(adft.out)

adft.out = unitroot(uscn.spot,trend="c",lags=2)
summary(adft.out)

adfn.out = unitroot(uscn.spot,trend="c",lags=2,statistic="n")
adfn.out

#
# next do pp tests
#
unitroot(uscn.spot,trend="c",method="pp")
unitroot(uscn.spot,trend="c",method="pp",statistic="n")

#
# Example: unit root testing on log-monthly stock prices
#

lnp = log(singleIndex.dat[,1])
lnp@title = "Log of S&P 500 Index"
par(mfrow=c(2,2))
plot.timeSeries(lnp,main="Log of S&P 500 index",
reference.grid=F)
acf.plot(acf(lnp,plot=F))
plot.timeSeries(diff(lnp),main="First difference of log S&P 500 Index",
reference.grid=F)
acf.plot(acf(diff(lnp),plot=F))
par(mfrow=c(1,1))

adft.out = unitroot(lnp,trend="ct",lags=4)
summary(adft.out)

#
# Example: stationary tests on interest rate differentials (get reference - Crowder?)
#

#
# simulate asymptotic distributions for KPSS distribution
#

wiener2 = function(nobs) {
	e = rnorm(nobs)
# create detrended errors
	e1 = e - mean(e)
	e2 = residuals(OLS(e~seq(1,nobs)))
# compute simulated Brownian Bridges
	y1 = cumsum(e1)
	y2 = cumsum(e2)
	intW2.1 = nobs^(-2) * sum(y1^2)
	intW2.2 = nobs^(-2) * sum(y2^2)
	ans = list(intW2.1=intW2.1,
				  intW2.2=intW2.2)
	ans
}

#
# simulate KPSS distributions:
# This could take some time.
#

nobs = 1000
nsim = 1000
KPSS1 = rep(0,nsim)
KPSS2 = rep(0,nsim)
for (i in 1:nsim) {
	BN.moments = wiener2(nobs)
	KPSS1[i] = BN.moments$intW2.1
	KPSS2[i] = BN.moments$intW2.2
}

#
# compute quantiles of distribution
#
quantile(KPSS1,probs=c(0.90,0.925,0.95,0.975,0.99))
quantile(KPSS2,probs=c(0.90,0.925,0.95,0.975,0.99))

#
# stationary test on exchange rates
#
args(stationaryTest)
kpss.out = stationaryTest(uscn.spot,trend="c")
class(kpss.out)
kpss.out