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

#Set the working directory to the folder where W2_EURPLN.csv is saved, e.g.:
#setwd("C:/.../Applied econometrics/W2_3")
#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("strucchange")) {install.packages("strucchange"); library("strucchange")}
if(!require("tseries")) {install.packages("tseries"); library("tseries")}
if(!require("MSwM")) {install.packages("MSwM"); library("MSwM")}
#tsDyn is periodically archived from CRAN over minor documentation notes
#(e.g. August 2026 - deprecated syntax in the example datasets, unrelated to the STAR/SETAR models)
#if the normal install fails, we fetch the fixed version from r-universe (independent of 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. Example of an artificial structural change

e<-rnorm(100, mean = 0, sd = 0.1)
x<-seq(from = 0, by=1, length.out = 100)
plot(x, type="l") #we created x (the explanatory variable)
y1<- 2 + 0.1*x +e  #we create y, with no structural change
plot(y1, type="l")

#structural change - a one-off shift by a constant (a change in beta_0)
y2<-y1
y2[60:100]<-5+ 0.1*x[60:100] +e[60:100]
plot(y2, type="l")

#structural change - a change in the regression coefficient (a change in beta_1)
y3<-y1
y3[60:100]<- -6.9 +0.25*x[60:100] +e[60:100] #constant chosen so the series is continuous at the break point
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") #downward trend "broken" at the point of change

m3<-lm(y3~x)
summary(m3)
plot(resid(m3), type="l") #V-shaped residuals


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

#generalized fluctuation test - does NOT require specifying the break point
#(tests whether the parameter estimates are stable over the whole sample)
sctest(m3)

#and if we want the data itself to reveal the MOMENT of change - the Bai-Perron procedure:
bp<-breakpoints(y3~x)
summary(bp)
breakdates(bp)      #detected date of the change
confint(bp)         #confidence interval for the moment of change
plot(bp)            #BIC and RSS for different numbers of breaks

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


#####################################
#3. Solving it with binary variables
#we create a binary variable (1 = structural change is present)
S<-as.numeric(seq_along(y1) > 59) #equivalent to a loop: 0 for obs. 1-59, 1 for obs. 60-100

#a change in the constant - the binary variable alone is enough
m2_2<-lm(y2~x+S)
summary(m2_2)
plot(resid(m2_2), type="l") #the residuals now look correct

#a change in the slope - an interaction is needed
m3_2<-lm(y3~x*S) #x*S expands to: x + S + x:S
summary(m3_2)
plot(resid(m3_2), type="l")


#####################################
#4. Threshold models (TAR/SETAR/STAR)

#let's do an exercise with an exchange rate
eurpln<-read.csv2("W2_EURPLN.csv", dec=".")
eurpln<-ts(eurpln[,2], start = c(2000,1), frequency = 12)
plot(eurpln, type="l") #it's always worth plotting the data right after loading it (a sanity check)

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) #shortened series for the plot
point.fct<-ts(fct$pred, end =c(2020,2), frequency = 12)
#fct$se contains two columns with the interval bounds - check which one is which:
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. A switching model (Markov switching)
#Instead of simulated data we use the same series as for the STAR model,
#so that the two approaches to regime change can be compared on the same data.

#we work on log returns (a stationary series)
r_eur <- diff(log(eurpln))*100
plot(r_eur, main="EUR/PLN - monthly returns (%)", ylab="%")

#baseline model: AR(0) with a constant - we are interested in switching of the mean and variance
d_eur <- data.frame(r = as.numeric(r_eur))
mod_lin <- lm(r ~ 1, data = d_eur)

#k = 2 regimes, p = 1 lag, sw = which parameters switch between regimes
#sw has length: number of coefficients + 1 (variance); here: constant, AR(1), variance
mod_ms <- msmFit(mod_lin, k = 2, p = 1, sw = c(TRUE, TRUE, TRUE),
                 control = list(parallel = FALSE))
summary(mod_ms)

#smoothed probabilities - which regime the exchange rate was in during a given month
plotProb(mod_ms, which = 1)
plotProb(mod_ms, which = 2)

#interpretation: usually one regime has low variance and the other has high variance
#(calm periods vs. periods of tension in the FX market - 2008/09, 2020, 2022)
plotDiag(mod_ms, which = 1)

#COMPARING THE TWO APPROACHES:
#STAR/SETAR - the regime depends DETERMINISTICALLY on an observed variable crossing a threshold
#Markov switching - the regime is UNOBSERVED, we infer it probabilistically