rm(list=ls()) # Caution: this clears the Environment

library(forecast)
library(tseries)
library(fUnitRoots)
library(dynlm)
library(aTSA)
set.seed(100)

#Tworzenie trendu liniowego i zmiennych sezonowych
lin_trend = ts(1:100,frequency = 4, start = c(1990, 1))
lin_trend
q = ts(seasonaldummy(lin_trend),frequency = 4, start = c(1990, 1))
q

#opóźnianie zmiennych
lag(q,-1)               #zwróc uwagę na parametr lag(q) to wyprzedzenie! 

#różnicowanie zmiennych
diff(lin_trend,1)

#generowanie białego szumu
e=ts(rnorm(100),frequency = 4, start = c(1990, 1))
plot(e)

#model trendu liniowego
lin_trend_model = 50+lin_trend+10*e
plot(lin_trend_model)

#generowanie błądzenia przypadkowego
random_walk=ts(cumsum(rnorm(100)),frequency = 4, start = c(1990, 1))
plot(random_walk)

#generowanie błądzenia przypadkowego z dryfem
random_walk_drift=ts(cumsum(rnorm(100,mean=0.1)),frequency = 4, start = c(1990, 1))
plot(random_walk_drift)

#model z trendem stochastycznym

stoch_trend_model = 10 + random_walk_drift + e
plot(stoch_trend_model)

#dodawanie wahań sezonowych

lin_trend_season_model = lin_trend_model -25*q[,"Q1"] - 15*q[,"Q2"] - 4*q[,"Q3"]
stoch_trend_season_model = stoch_trend_model -4*q[,"Q1"] - 3*q[,"Q2"] - 1*q[,"Q3"]

plot(lin_trend_season_model)
plot(stoch_trend_season_model)

#dekompozycja na trend, składnik sezonowy i losowy
d_lint_model = decompose(lin_trend_season_model)
d_stocht_model = decompose(stoch_trend_season_model)
plot(d_lint_model)
plot(d_stocht_model)

plot(ts.union(lin_trend,d_lint_model$trend))
plot(ts.union(random_walk_drift,d_stocht_model$trend))

acf(e)
pacf(e)

acf(random_walk_drift)
pacf(random_walk_drift)

#acf i pacf dla modelu trendu liniowego
acf(lin_trend_model)
pacf(lin_trend_model)

#acf i pacf dla zmiennej (model z trendem deterministycznym) z usuniętym trendem
trend_model = tslm(lin_trend_model~lin_trend)
summary(trend_model)
acf(trend_model$residuals)
pacf(trend_model$residuals)

#acf i pacf dla modelu trendu stochastycznego
acf(stoch_trend_model)
pacf(stoch_trend_model)

#acf i pacf dla zmiennych z usuniętym trendem linowym
strend_model = tslm(stoch_trend_model~lin_trend)
summary(strend_model)
acf(strend_model$residuals)
pacf(strend_model$residuals)

#Testowanie istnienia pierwiastka jednostkowego

#biały szum
adf.test(e)
adf.test(e,k=1)
#błądzenie losowe
adf.test(random_walk)

#ADF model trendu liniowego
adf.test(lin_trend_model)
adf.test(lin_trend_model,k=1)

#ADF model trendu stochastycznego
adf.test(stoch_trend_model)
library(fUnitRoots)
adfTest(stoch_trend_model, lags = 2, type = "nc")
adfTest(stoch_trend_model, lags = 2, type = "c")
adfTest(stoch_trend_model, lags = 2, type = "ct")
adfTest(lin_trend_model, lags = 2, type = "nc")
adfTest(lin_trend_model, lags = 2, type = "c")
adfTest(lin_trend_model, lags = 2, type = "ct")

#ADF dla różnicowanie szeregu
adf.test(diff(stoch_trend_model,1))



#Test KPSS 

kpss.test(lin_trend_model)                 #H0: stacjonarność
kpss.test(lin_trend_model,null="Trend")    #H0: trendostacjonarność

#KPSS dla modelu trendu stochastycznego
kpss.test(stoch_trend_model)               #H0: stacjonarność
kpss.test(stoch_trend_model,"Trend")       #H0: trendostacjonarność

kpss.test(diff(stoch_trend_model,1))       #H0: trendostacjonarność

#spektrum dla modelu z sezonowością roczną i tygodniową
periodic=ts(1/7*((1:3650)%%7)+1/365*((1:3650)%%365)+2*rnorm(3650))
plot(periodic)
delta=1/365
mspect=spectrum(periodic,log="no")
spectx <- mspect$freq/delta
specty <- 2*mspect$spec
plot(spectx, specty,type="l")
plot(spectx[1:30], specty[1:30],type="l")

#Kointegracja

x=random_walk_drift
e1=ts(rnorm(100),frequency = 4, start = c(1990, 1))
y = 5+0.9*x+e1
plot(ts.union(y,x),plot.type = c("single"))
resid = y-0.9*x
plot(resid)

#ECM
T = 100
rm(x)
x=ts(cumsum(rnorm(T,mean=0.1)),frequency = 4, start = c(1990, 1)) #random walk with drift
rm(e1)
e1 = ts(rnorm(T),frequency = 4, start = c(1990, 1))
rm(y)
y = vector()
y[1:T] = 0
for(t in 3:T) {
  y[t] = y[t-1] - 0.50*(y[t-1]-0.9*x[t-1])+0.1*lag(y[t-1]-y[t-2],-1)+0.1*(x[t-1]-x[t-2])+e1[t]
}
y = ts(y,frequency = 4, start = c(1990, 1))

plot(ts.union(y,x),plot.type = c("single"))

#Dwustopniowa procedura Engla-Grangera
#1 stopień
adf.test(y)
adf.test(diff(y))
adf.test(x)
adf.test(diff(x))
ce_model = lm(y~x)
summary(ce_model)
ce = ts(ce_model$residuals,frequency = 4, start = c(1990, 1))
#2 stopień
coint.test(y,x)
ecm = dynlm(d(y) ~  L(ce) + L(d(y,1)) + L(d(x,1)))
summary(ecm)