#************************************************************ 
# Unit Root Test Simulation
# Jin Man Lee
# Spectral Analysis on ARMA type of Data
#install.packages("urca")
# install.packages("FinTS")
#install.packages("fUnitRoots")
#------------------------------------------------------------


library("urca")
library("tseries")

# 1. Unit Root for Random Walk 
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 
}

nobs = 1000
nsim = 10000
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) 
}

quantile(DF,probs=c(0.01,0.05,0.1))
quantile(NB,probs=c(0.01,0.05,0.1))
par(mfrow=c(2,1))
hist(DF)
hist(NB)


#KPSS Type
quantile(DF,probs=c(0.9,0.95,0.99))
quantile(NB,probs=c(0.9,0.95,0.99))


# 2. Unit Root with Serial Correlation 
wiener2 <- function(nobs) {
  e = rnorm(nobs)
  # create detrended errors
  e1 = e - mean(e)
  e2 = residuals(lm(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
#
nobs = 100
nsim = 100
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
}

#KPSS Type
quantile(KPSS1,probs=c(0.9,0.95,0.99))
quantile(KPSS2,probs=c(0.9,0.95,0.99))
par(mfrow=c(2,1))
hist(KPSS1,breaks=12, col="blue")
hist(KPSS2,breaks=12, col="red")


# DF test Using urca based on MacKinnon’s (1996)
" No constant No Trend"
qunitroot(c(0.01,0.05,0.1), trend="nc", statistic = "t", N =100)
qunitroot(c(0.01,0.05,0.1), trend="nc", statistic = "n", N =100)

" Constant No Trend"
qunitroot(c(0.01,0.05,0.1), trend="c", statistic = "t", N =100)
qunitroot(c(0.01,0.05,0.1), trend="c", statistic = "n", N =100)

" Constant and Trend"
qunitroot(c(0.01,0.05,0.1), trend="ct", statistic = "t", N =100)
qunitroot(c(0.01,0.05,0.1), trend="ct", statistic = "n", N =100)

# UniRoot Tests using urca
y = rnorm(100)
df1 <- ur.df(y, type = "drift",  selectlags ="AIC")
pp1 <- ur.pp(y, type = "Z-tau", model = "constant",lags ="short")
kp1 <- ur.kpss(y, type = "tau",  lags ="short")
summary(df1)
summary(pp1)
summary(kp1)

# UnitRoot Tests using tseries
df2 <- adf.test(y, alternative = "stationary", k = 0)
df2$statistic 
df2$p.value 

pp2 <- pp.test(y, type= "Z(t_alpha)")
pp2$statistic
pp2$p.value

kp2 <- kpss.test(y)
kp2$statistic
kp2$p.value