rm(list=ls()) # Caution: this clears the Environment
library(httr)
library(xlsx)
library(dynlm)
library(broom)
library(caTools)
library(uroot)
library(fUnitRoots)
library(tseries)
library(forecast)
library(lmtest)
library(car)
#library(vars)

url <- "http://www.ekonometria.wne.uw.edu.pl/uploads/Main/bezrob_BAEL.xlsx"
destfile <- "bezrob_BAEL.xlsx"
curl::curl_download(url, destfile)
Dane <- read.xlsx(destfile,sheetIndex = 1)

Dane_ts <- ts(data=Dane, frequency = 4, start = c(1993, 1))
U = Dane_ts[,"UO"]/100
dU = diff(U,4)                                 #zmiana bezrobocia q/q

dGDP =Dane_ts[,"PKB"]/100-1                       #wzrost pkb q/q
dgdp = log(Dane_ts[,"PKB"]/100)                  #log zmiany gdp q/q

#Dane na logarytmach
u = -log(1-U)                  #log (ln(aktywni zawodowo)-ln(pracujący))
du = diff(u,4)                             #zmiana bezrobocia q/

plot(U)
plot(u)
cor(U,u,use = "complete.obs")
lag(u,1)

#interpolacja brakujących danych
library(imputeTS)
U1=na_seadec(U)
u1=na_seadec(u)
plot(ts.union(U1,U))
mydata = ts.union(U1,u1)


#Średnia ruchoma
ma_u = ts(runmean(u,k=4), frequency = 4, start = c(1993, 1))

plot(ma_u)

#dekompozycja szeregu na trend, składnik sezonowy i losowy
ud = decompose(u1)
plot(ud)
ud1=stl(u1,s.window=4)
plot(ud1)

#białego szumu, procesu trendu liniowego, błądzenia losowego, błądzenia losowego z dryfem

lin_trend = ts()

#Funkcje ACF i PACF

acf(na.exclude(u))
pacf(na.exclude(u))

acf(na.exclude(diff(u,4)))
pacf(na.exclude(diff(u,4)))

acf(na.exclude(diff(u,1)))
pacf(na.exclude(diff(u,1)))

#ARIMA model

arima(u,order=c(2,1,2))
arima_model <- arima(u,order=c(2,1,2))

#przekształcenie wyników w ramkę™
arima_model_table = tidy(arima_model,conf.int=TRUE)
#dodanie statystyk z
arima_model_table[,"z"]=arima_model_table$estimate/arima_model_table$std.error
#dodanie p-value
arima_model_table[,"p.value"]=pnorm(-abs(arima_model_table$z))*2
arima_model_table = arima_model_table[,c(1,2,3,6,7,4,5)]

arima_model_table

#Test ADF
adf.test(na.exclude(u),k=8)
adfTest(na.exclude(u), lags = 8, type = "nc")
adfTest(na.exclude(u), lags = 8, type = "c")
adfTest(na.exclude(u), lags = 8, type = "ct")

adf.test(na.exclude(diff(u,4),k=8))

#Test KPSS
kpss.test(na.exclude(u))
kpss.test(na.exclude(u),null="Trend")
kpss.test(na.exclude(diff(u,4)))
kpss.test(na.exclude(diff(u,4)),null="Trend")

hegy.test(u, deterministic = c(1,1,1), lag.method = "fixed", maxlag = 2)

#generowanie prognoz za pomocą całego zbioru danych
auto_model= auto.arima(u)
summary(auto_model)
fcst=forecast(auto_model,h=10)
summary(fcst) #statystyki ilustrujące jakość prognoz
plot(fcst)

#porównanie jakości prognoz (10 ostatnich okresów jest wyłączona z estymacji)
horiz=8
u_training = ts(u1[1:(length(u)-horiz)], frequency = 4, start = c(1993, 1)) #zbiór treningowy
auto_model= auto.arima(u_training)  #SARIMA
fcst_sarima = forecast(auto_model,h=horiz)
accuracy(fcst_sarima, u)
plot(ts(cbind(u,fcst_sarima$mean)[(length(u)-horiz):length(u),],frequency = 4, end = c(2017, 4)),plot.type = c("single"))

#porównanie z rozbudowaną wersją wyrównywania wykładniczego
fit_ets = ets(u_training)
fcst_ets = forecast(fit_ets,h=horiz)
plot(fcst_ets)
accuracy(fcst_ets,u)
plot(ts(cbind(u,fcst_ets$mean)[(length(u)-horiz):length(u),],frequency = 4, end = c(2017, 4)),plot.type = c("single"))

#testowanie diagnostyczne model arima(2,1,2)
acf(na.exclude(arima_model$residuals))
pacf(na.exclude(arima_model$residuals))

Box.test(na.exclude(arima_model$residuals),type="Box-Pierce")
Box.test(na.exclude(arima_model$residuals),type="Ljung-Box")

Box.test(na.exclude(arima_model$residuals),lag=4,type="Box-Pierce")
Box.test(na.exclude(arima_model$residuals),lag=4,type="Ljung-Box")

#Okun Law

mydata<-ts(cbind(U,dU,u,du,dGDP,dgdp), frequency = 4, start = c(1993, 1))

#Modele DL, Okun Law
model1 <- dynlm(dU~dGDP,data=mydata)
summary(model1)
bgtest(model1, order=4)

model2 <- dynlm(dU~dGDP+L(dGDP)+L(dGDP,2)+L(dGDP,3)+L(dGDP,4),data=mydata)
summary(model2)
bgtest(model2, order=4)

#Model ADL, Okun Law
model3 <- dynlm(dU~L(dU,1)+L(dU,2)+dGDP, data=mydata)
summary(model3)
bgtest(model3, order=4)

#Test przyczynowości w sensie Grangera
model3 <- dynlm(dU~L(dU,1)+L(dU,2)+L(dGDP) +L(dGDP,2),data=mydata)
summary(model3)
bgtest(model3, order=4)

linearHypothesis(model3, c("L(dGDP) = 0", "L(dGDP, 2) = 1"))

#Model SARIMAX
regr = cbind(dGDP,lag(dGDP,-1))[1:length(dGDP),]
auto.arima(U,xreg=regr)

#VAR
library(vars)
vardata = ts(cbind(U1,dGDP)[5:length(dGDP),],frequency = 4, start = c(1994, 1))
var_model=VAR(vardata,p=4,type = c("const"), season=4)
roots(var_model,modulus=TRUE)
serial.test(var_model)

#Wybór liczby opóźnień
VARselect(vardata,lag.max=8,type = c("const"),season=4)
var_model_auto = VAR(vardata,lag.max = 8,ic = c("SC")) #procedura zautomatyzowana
summary(var_model_auto)
roots(var_model_auto,modulus=TRUE)
serial.test(var_model_auto)

causality(var_model_auto, cause=c("dGDP"))
causality(var_model_auto, cause=c("U1"))

#nieskumulowane funkcje reakcji
irf_unit=irf(var_model_auto,ortho=FALSE)
plot(irf_unit)

#skumulowane funkcje reakcji
irf_unit_cumm=irf(var_model_auto,response="dGDP",ortho=FALSE,cummulative=TRUE)
plot(irf_unit_cumm)

#funkcje reakcji dla szoków zortogonalizowanych
irf_ortho=irf(var_model_auto,ortho=TRUE)
plot(irf_ortho)

#Test Johansena
vecm_model=ca.jo(vardata, ecdet = c("const"), type="trace", K=5, spec="longrun",season=4)
summary(vecm_model)
var_vecm = vec2var(vecm_model, r = 1)
serial.test(var_vecm)

#prognozowanie za pomocą VAR
fcst_vecm = predict(var_vecm, n.ahead = 10)
fanchart(fcst_vecm,plot.type = "single")

#dekompozycja błędu prognozy
fevd_unit=fevd(var_model_auto,ortho=FALSE)
plot(fevd_unit,plot.type = "single")
fevd_ortho=fevd(var_model_auto,ortho=TRUE)
plot(fevd_ortho,plot.type = "single")