rm(list=ls()) #clear the workspace
set.seed(2020) #so that everyone gets the same results

#Set the working directory to the folder containing this lecture's files, e.g.:
#In RStudio: Session -> Set Working Directory -> To Source File Location

#note: require() alone does not load the package after installing it - hence library() in the braces
if(!require("tseries")) {install.packages("tseries"); library("tseries")}
if(!require("forecast")) {install.packages("forecast"); library("forecast")}


#####################################
#1. White noise - a stationary process
x <- rnorm(100, mean = 0, sd = 1)
plot(x, type="l", main="White noise", sub="", xlab="", ylab="")
Acf(x, main="ACF", sub="", xlab="", ylab="", ylim=c(-1,1))


#####################################
#2. Random walk - a non-stationary process
y<-seq(from = 0, by=0, length.out = 100) #create an empty object
#use a loop to add the previous value of the series to the current shock
for (i in 1:99){
  y[1]=x[1]
  y[i+1]=y[i]+x[i+1]
}

#alternatively we can use the function
y2<-cumsum(x)

plot(y, type="l")
points(y2, col=2) #sanity check - both methods give an identical result

plot(y, type="l", main="Random walk", sub="", xlab="", ylab="")
Acf(y, main="ACF", sub="", xlab="", ylab="", ylim=c(-1,1))


#####################################
#3. Let's run the tests and check their power
#ADF test
adf.test(x) #We reject H0 - the series is stationary
adf.test(y) #No grounds to reject H0 - the series is NON-stationary

#KPSS test
kpss.test(x) #No grounds to reject H0 - the series is stationary
kpss.test(y) #We reject H0 - the series is NON-stationary


#####################################
#4. What about a highly persistent series (0.9<rho<1)?
#rho - how persistent the process is
#M - number of Monte Carlo simulations

p.value<-function(rho, M, adf=c(), kpss=c()){

for (i in 1:M){
  #set.seed(i)  #uncomment if you want to get the same results every time
  e<-rnorm(50, mean = 0, sd = 1)
  y<-seq(from = 0, by=0, length.out = 50)
  for (j in 1:49){
    y[1]=e[1]
    y[j+1]=rho*y[j]+e[j+1]
  }
  adf<-c(adf,adf.test(y)$p.value)
  kpss<-c(kpss,kpss.test(y)$p.value)
}
  pval<-as.data.frame(cbind(adf,kpss))
  colnames(pval)<-c("adf","kpss")
  return(pval)
}

#How often did the test find stationarity (rejecting H0):
pval<-p.value(rho=0.1, M=120)
pval$adf
pval$kpss

#How often did the test find stationarity (%):
sum(pval$adf<(0.05))/length(pval$adf)*100
sum(pval$kpss>(0.05))/length(pval$adf)*100

#INTERPRETATION NOTE: for ADF (H0: non-stationarity) the percentage above is the actual
#POWER of the test, because H0 is false here (the simulated process IS stationary for rho<1).
#For KPSS (H0: stationarity) the share of "correct" calls in this experiment measures the
#test's ability to NOT REJECT a true H0 (i.e. its actual size), not power in the strict
#sense - we could only talk about KPSS's power by simulating genuinely non-stationary
#data and checking how often KPSS correctly rejects it.

#let's build a summary
adf=c()
kpss=c()

for (i in seq(from=0.9, to=0.99, by=0.01)){
  pval<-p.value(rho=i, M=250)
  adf<-c(adf,sum(pval$adf<(0.05))/length(pval$adf)*100)
  kpss<-c(kpss,sum(pval$kpss>(0.05))/length(pval$adf)*100)
}

wynik<-rbind(adf, kpss)
colnames(wynik)<-seq(from=0.9, to=0.99, by=0.01)
wynik

#let's plot the power curve - this best shows what happens as rho -> 1
plot(seq(from=0.9, to=0.99, by=0.01), wynik["adf",], type="b", col="darkgreen",
     ylim=c(0,100), xlab=expression(rho), ylab="% of stationarity calls",
     main="Power curve: ADF and KPSS")
lines(seq(from=0.9, to=0.99, by=0.01), wynik["kpss",], type="b", col="grey40", lty=2)
abline(h=5, lty=3, col="grey")
legend("topright", legend=c("ADF (power)","KPSS (correct stationarity call)"),
       col=c("darkgreen","grey40"), lty=c(1,2), bty="n")


#####################################
#5. Does the test's deterministic specification matter for its power?
#tseries' adf.test() ALWAYS fits a regression with a constant and a linear trend
#(this cannot be switched off). Let's check how much power we lose by testing a
#process with no drift or trend (like ours) as if it had a trend.
if(!require("urca")) {install.packages("urca"); library("urca")}

sim_ar1 <- function(rho, n=50){
  e<-rnorm(n); y<-numeric(n); y[1]<-e[1]
  for (t in 2:n) y[t]<-rho*y[t-1]+e[t]
  y
}

rho_test <- 0.1
M2 <- 300
rej_trend<-0; rej_const<-0; rej_none<-0
for (m in 1:M2){
  yy <- sim_ar1(rho_test)
  if (suppressWarnings(adf.test(yy)$p.value) < 0.05) rej_trend<-rej_trend+1

  u_const <- ur.df(yy, type="drift", lags=trunc((50-1)^(1/3)), selectlags="Fixed")
  if (u_const@teststat[1] < u_const@cval[1,"5pct"]) rej_const<-rej_const+1

  u_none <- ur.df(yy, type="none", lags=trunc((50-1)^(1/3)), selectlags="Fixed")
  if (u_none@teststat[1] < u_none@cval[1,"5pct"]) rej_none<-rej_none+1
}
cat(sprintf("ADF power (constant+trend): %.1f%%\n", 100*rej_trend/M2))
cat(sprintf("ADF power (constant only): %.1f%%\n", 100*rej_const/M2))
cat(sprintf("ADF power (no constant/trend): %.1f%%\n", 100*rej_none/M2))
#CONCLUSION: fitting an overly rich deterministic specification (a trend that isn't
#actually in the data) strongly reduces test power - regardless of the persistence rho.


#####################################
#6. Panel data
#We read in an example dataset
data<-read.csv2("W4_data.csv", dec=".")
head(data) #check what the data looks like
if(!require("plm")) {install.packages("plm"); library("plm")}

#Levin-Lin-Chu (2002)
#method 1 - split the data frame
sle<-data.frame(split(data$SLE, data$Country))
purtest(sle,  pmax = 2, exo = "intercept", test="levinlin" )

#method 2 - build a panel data frame object
data<-pdata.frame(data, index=c("Country","Year"), drop.index = TRUE)
purtest(data$PAFP,  pmax = 4, exo = "intercept", test="levinlin" )
purtest(data$SLE,  pmax = 4, exo = "intercept", test="levinlin" )
purtest(data$X1,  pmax = 4, exo = "intercept", test="levinlin" )
purtest(data$FDIYS,  pmax = 4, exo = "intercept", test="levinlin" )
purtest(data$GOVC1,  pmax = 4, exo = "intercept", test="levinlin" )
purtest(data$POT,  pmax = 4, exo = "intercept", test="levinlin" )

#Im-Pesaran-Shin (2003)
purtest(data$PAFP,  pmax = 4, exo = "intercept", test="ips" )
purtest(data$SLE,  pmax = 4, exo = "intercept", test="ips" )
purtest(data$X1,  pmax = 4, exo = "intercept", test="ips" )
purtest(data$FDIYS,  pmax = 4, exo = "intercept", test="ips" )
purtest(data$GOVC1,  pmax = 4, exo = "intercept", test="ips" )
purtest(data$POT,  pmax = 4, exo = "intercept", test="ips" )

#NOTE: for the variable X1 the LLC and IPS tests give CONFLICTING conclusions (check the p-values!).
#This is not an accident - LLC assumes a COMMON rho parameter for all countries (alternative
#hypothesis: ALL countries have a stationary series), while IPS allows DIFFERENT rho_i across
#countries (alternative hypothesis: SOME countries have a stationary series). When the panel is
#heterogeneous, the two tests can lead to different conclusions.


#####################################
# Homework
#1. Repeat the power-curve simulation for n=100 and n=200 instead of n=50. How does the shape of
#   the power curve change as the sample size grows?
#2. For rho=0.95, check the power of the KPSS test against a FALSE H0 - simulate a non-stationary
#   process (a random walk) and compute how often KPSS correctly rejects it. Compare with ADF's
#   power in the same case.
#3. Find another variable in the panel data (or your own dataset) for which the LLC and IPS tests
#   give different conclusions, and try to explain the reason for the difference.