rm(list=ls()) #czyszczenie obszaru roboczego
set.seed(2020) #aby wszyscy otrzymali takie same wyniki

#Ustaw katalog roboczy na folder, w którym zapisane są pliki tego wykładu, np.:
#W RStudio: Session -> Set Working Directory -> To Source File Location

#uwaga: require() samo nie ładuje pakietu po instalacji - stąd library() w nawiasie klamrowym
if(!require("tseries")) {install.packages("tseries"); library("tseries")}
if(!require("forecast")) {install.packages("forecast"); library("forecast")}


#####################################
#1. Biały szum - proces stacjonarny
x <- rnorm(100, mean = 0, sd = 1)
plot(x, type="l", main="Biały szum", sub="", xlab="", ylab="")
Acf(x, main="ACF", sub="", xlab="", ylab="", ylim=c(-1,1))


#####################################
#2. Błądzenie losowe - proces niestacjonarny
y<-seq(from = 0, by=0, length.out = 100) #tworzę pusty obiekt
#za pomocą pętli sumuję poprzednią wartość szeregu i bieżącą wartość składnika losowego
for (i in 1:99){
  y[1]=x[1]
  y[i+1]=y[i]+x[i+1]
}

#alternatywnie możemy wykorzystać funkcję
y2<-cumsum(x)

plot(y, type="l")
points(y2, col=2) #sanity check - obie metody dają identyczny wynik

plot(y, type="l", main="Błądzenie losowe", sub="", xlab="", ylab="")
Acf(y, main="ACF", sub="", xlab="", ylab="", ylim=c(-1,1))


#####################################
#3. Wykonajmy testy i sprawdźmy ich moc
#Test ADF
adf.test(x) #Odrzucamy H0 - szereg jest stacjonarny
adf.test(y) #Brak podstaw do odrzucenia H0 - szereg jest NIEstacjonarny

#Test KPSS
kpss.test(x) #Brak podstaw do odrzucenia H0 - szereg jest stacjonarny
kpss.test(y) #Odrzucamy H0 - szereg jest NIEstacjonarny


#####################################
#4. A co z szeregiem o wysokiej persystencji (0.9<rho<1)?
#rho - jaka jest persystencja
#M - liczba symulacji Monte Carlo

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

for (i in 1:M){
  #set.seed(i)  #do odznaczenia, jeśli chcemy za każdym razem dostawać te same wyniki
  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)
}

#Jak często test wykazał stacjonarność (odrzucenie H0):
pval<-p.value(rho=0.1, M=120)
pval$adf
pval$kpss

#Jak często test wykazał stacjonarność (%):
sum(pval$adf<(0.05))/length(pval$adf)*100
sum(pval$kpss>(0.05))/length(pval$adf)*100

#UWAGA interpretacyjna: dla ADF (H0: niestacjonarność) powyższy odsetek to rzeczywista
#MOC testu, bo H0 jest tu fałszywe (symulowany proces JEST stacjonarny dla rho<1).
#Dla KPSS (H0: stacjonarność) odsetek "poprawnych" wskazań mierzy w tym eksperymencie
#zdolność testu do NIEODRZUCENIA prawdziwej H0 (czyli jego rzeczywisty rozmiar/size),
#a nie moc w ścisłym sensie - o mocy KPSS moglibyśmy mówić dopiero symulując
#dane rzeczywiście niestacjonarne i sprawdzając, jak często KPSS je poprawnie odrzuca.

#zróbmy podsumowanie
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

#narysujmy krzywą mocy - to najlepiej pokazuje, co się dzieje przy 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="% wskazań stacjonarności",
     main="Krzywa mocy: ADF i 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 (moc)","KPSS (poprawne wskazanie stacjonarności)"),
       col=c("darkgreen","grey40"), lty=c(1,2), bty="n")


#####################################
#5. Czy specyfikacja deterministyczna testu ma znaczenie dla jego mocy?
#adf.test() z pakietu tseries ZAWSZE dopasowuje regresję ze stałą i trendem liniowym
#(nie da się tego wyłączyć). Sprawdźmy, jak dużo mocy tracimy, testując proces
#bez dryfu i trendu (taki jak nasz) tak, jakby miał 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("Moc ADF (stała+trend): %.1f%%\n", 100*rej_trend/M2))
cat(sprintf("Moc ADF (tylko stała): %.1f%%\n", 100*rej_const/M2))
cat(sprintf("Moc ADF (bez stałej/trendu): %.1f%%\n", 100*rej_none/M2))
#WNIOSEK: dopasowanie zbyt bogatej specyfikacji deterministycznej (trend, którego
#w danych nie ma) mocno obniża moc testu - niezależnie od persystencji rho.


#####################################
#6. Dane panelowe
#Wczytujemy przykładowy zbiór danych
data<-read.csv2("W4_data.csv", dec=".")
head(data) #sprawdzanie, jak wyglądają dane
if(!require("plm")) {install.packages("plm"); library("plm")}

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

#2 sposób - tworzymy obiekt panelowy dataframe
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" )

#UWAGA: dla zmiennej X1 testy LLC i IPS dają SPRZECZNE wnioski (sprawdźcie p-value!).
#To nieprzypadkowe - LLC zakłada WSPÓLNY parametr rho dla wszystkich krajów (hipoteza
#alternatywna: WSZYSTKIE kraje mają szereg stacjonarny), a IPS pozwala na RÓŻNE rho_i
#w różnych krajach (hipoteza alternatywna: CZĘŚĆ krajów ma szereg stacjonarny). Gdy
#panel jest niejednorodny, oba testy mogą prowadzić do różnych wniosków.


#####################################
# Praca domowa
#1. Powtórz symulację krzywej mocy dla n=100 i n=200 zamiast n=50. Jak zmienia się
#   kształt krzywej mocy wraz ze wzrostem próby?
#2. Dla rho=0.95 sprawdź moc testu KPSS wobec FAŁSZYWEJ H0 - zasymuluj proces
#   niestacjonarny (błądzenie losowe) i policz, jak często KPSS go poprawnie odrzuca.
#   Porównaj z mocą ADF dla tego samego przypadku.