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

#Ustaw katalog roboczy na folder, w którym zapisany jest plik W2_EURPLN.csv, np.:
#setwd("C:/.../Ekonometria praktyczna/W5_6")
#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("strucchange")) {install.packages("strucchange"); library("strucchange")}
if(!require("tseries")) {install.packages("tseries"); library("tseries")}
if(!require("MSwM")) {install.packages("MSwM"); library("MSwM")}
#tsDyn bywa okresowo zdejmowane z CRAN z powodu drobnych uwag w dokumentacji
#(np. sierpień 2026 - przestarzała składnia w przykładowych danych, nic wspólnego z modelami STAR/SETAR)
#jeśli zwykła instalacja zawiedzie, pobieramy naprawioną wersję z r-universe (niezależnie od CRAN)
if(!require("tsDyn")) {
  install.packages("tsDyn")
  if(!require("tsDyn")) {
    install.packages("tsDyn", repos = c("https://matthieustigler.r-universe.dev", "https://cloud.r-project.org"))
    library("tsDyn")
  }
}


#####################################
#1. Przykład sztucznej zmiany strukturalnej

e<-rnorm(100, mean = 0, sd = 0.1)
x<-seq(from = 0, by=1, length.out = 100)
plot(x, type="l") #utworzyliśmy x (zmienną objaśniającą)
y1<- 2 + 0.1*x +e  #tworzymy y, nie ma zmiany strukturalnej
plot(y1, type="l")

#zmiana strukturalna - jednorazowe przesunięcie o stałą (zmiana beta_0)
y2<-y1
y2[60:100]<-5+ 0.1*x[60:100] +e[60:100]
plot(y2, type="l")

#zmiana strukturalna - zmiana współczynnika regresji (zmiana beta_1)
y3<-y1
y3[60:100]<- -6.9 +0.25*x[60:100] +e[60:100] #stała dobrana tak, aby szereg był ciągły w punkcie zmiany
plot(y3, type="l")

m1<-lm(y1~x)
summary(m1)
plot(resid(m1), type="l")

m2<-lm(y2~x)
summary(m2)
plot(resid(m2), type="l") #trend spadkowy "złamany" w punkcie zmiany

m3<-lm(y3~x)
summary(m3)
plot(resid(m3), type="l") #reszty w kształcie litery V


#####################################
#2. Test Chowa
sctest(y1~x, type="Chow", point=60)
sctest(y2~x, type="Chow", point=60)
sctest(y3~x, type="Chow", point=60)

#uogólniony test fluktuacji - NIE wymaga podania punktu zmiany
#(bada, czy oszacowania parametrów są stabilne w całej próbie)
sctest(m3)

#a jeżeli chcemy, aby to dane wskazały MOMENT zmiany - procedura Bai-Perrona:
bp<-breakpoints(y3~x)
summary(bp)
breakdates(bp)      #wykryta data zmiany
confint(bp)         #przedział ufności dla momentu zmiany
plot(bp)            #BIC i RSS dla różnej liczby załamań

plot(y3, type="l")
lines(fitted(bp, breaks=1), col="red")
abline(v=breakpoints(bp)$breakpoints, lty=2, col="blue")


#####################################
#3. Rozwiązanie za pomocą zmiennych binarnych
#tworzymy zmienną binarną (1=występuje zmiana strukturalna)
S<-as.numeric(seq_along(y1) > 59) #odpowiednik pętli: 0 dla obs. 1-59, 1 dla obs. 60-100

#zmiana stałej - wystarczy sama zmienna binarna
m2_2<-lm(y2~x+S)
summary(m2_2)
plot(resid(m2_2), type="l") #reszty wyglądają już poprawnie

#zmiana nachylenia - potrzebna interakcja
m3_2<-lm(y3~x*S) #x*S rozwija się do: x + S + x:S
summary(m3_2)
plot(resid(m3_2), type="l")


#####################################
#4. Modele progowe (TAR/SETAR/STAR)



#zróbmy ćwiczenie z kursem walutowym
eurpln<-read.csv2("W2_EURPLN.csv", dec=".")
eurpln<-ts(eurpln[,2], start = c(2000,1), frequency = 12)
plot(eurpln, type="l") #Zawsze warto narysować wykres po wczytaniu danych (sprawdzenie)

horizon<-12
f<-length(eurpln)-horizon

mod.fx <- star(eurpln[1:f],  d=2, control=list(maxit=3000))
fct<-predict(mod.fx, eurpln[1:f], horizon, type="bootstrap")

eurpln2<-ts(eurpln[180:length(eurpln)], end =c(2020,2), frequency = 12) #skrócony szereg do wykresu
point.fct<-ts(fct$pred, end =c(2020,2), frequency = 12)
#fct$ci zawiera dwie kolumny z granicami przedziału - sprawdźcie, która jest która:
head(fct$ci)
ci1<-ts(fct$ci[,1], end =c(2020,2), frequency = 12)
ci2<-ts(fct$ci[,2], end =c(2020,2), frequency = 12)

plot(eurpln2, type="l", ylim=c(3,5))
lines(point.fct, col="red")
lines(ci1, col=3)
lines(ci2, col=3)


#####################################
#5. Model przełącznikowy (Markov switching)
#Zamiast danych symulowanych używamy tego samego szeregu co przy modelu STAR,
#dzięki czemu można porównać oba podejścia do zmiany reżimu na tych samych danych.

#pracujemy na logarytmicznych stopach zwrotu (szereg stacjonarny)
r_eur <- diff(log(eurpln))*100
plot(r_eur, main="EUR/PLN - miesięczne stopy zwrotu (%)", ylab="%")

#model bazowy: AR(0) ze stałą - interesuje nas przełączanie średniej i wariancji
d_eur <- data.frame(r = as.numeric(r_eur))
mod_lin <- lm(r ~ 1, data = d_eur)

#k = 2 reżimy, p = 1 opóźnienie, sw = które parametry przełączają się między reżimami
#sw ma długość: liczba współczynników + 1 (wariancja); tutaj: stała, AR(1), wariancja
mod_ms <- msmFit(mod_lin, k = 2, p = 1, sw = c(TRUE, TRUE, TRUE),
                 control = list(parallel = FALSE))
summary(mod_ms)

#prawdopodobieństwa wygładzone - w którym reżimie znajdował się kurs w danym miesiącu
plotProb(mod_ms, which = 1)
plotProb(mod_ms, which = 2)

#interpretacja: zwykle jeden reżim ma niską, a drugi wysoką wariancję
#(okresy spokoju vs okresy napięć na rynku walutowym - 2008/09, 2020, 2022)
plotDiag(mod_ms, which = 1)

#PORÓWNANIE PODEJŚĆ:
#STAR/SETAR - reżim zależy DETERMINISTYCZNIE od przekroczenia progu przez obserwowaną zmienną
#Markov switching - reżim jest NIEOBSERWOWANY, wnioskujemy o nim probabilistycznie