### BREAKS IN LEVEL-TREND REGRESSION
rm(list=ls())
graphics.off()
setwd(dirname(parent.frame(2)$ofile)) #Set working directory to folder with R file
setwd('../Bases')

library(fracdiff)
library(strucchange)
library(urca)
library("readxl")
library(memochange)
library(tseries)
library(LongMemoryTS)


##### DATA
database <- read_excel("manycountries.xlsx",sheet = "Hoja1")
series <- ts(database, frequency = 1, start=c(1270))

############################
# COUNTRY
############################

ind = 9    #Country selection
bdw = 0.80 #Bandwidth for GPH y ELW

country <- series[,ind]
TSdata_raw <- ts(country,frequency=1,start=1300)
TSdata = na.remove(TSdata_raw)
TSdata = log(TSdata)

###############################
# Breaks with Bai-Perron
###############################
tt <- 1:length(TSdata) 
TSdata_brk <- breakpoints(TSdata ~ 1+tt, h = 0.17)
TSdata_brk

nbreaks = length(TSdata_brk$breakpoints)  #Numero de breaks
###################Plot
par(mar=c(2.5,2.5,1.5,1)) 
plot(TSdata,ylim=c(min(TSdata),max(TSdata)),ylab="",xlab="",lwd = 3,main = "Real GDP per capita in the Netherlands, 1348-2019") 
lines(fitted(TSdata_brk, breaks = length(TSdata_brk$breakpoints)), col = "red",lty = 1, lwd = 2)
lines(confint(TSdata_brk, breaks = length(TSdata_brk$breakpoints)),lwd=2,col="blue")

# par(mar=c(2.5,2.5,1.5,1)) 
# plot(TSdata,ylim=c(min(TSdata),max(TSdata)),ylab="",xlab="",lwd = 3,main = names(database)[ind]) #names(database)[ind]
# # grid(nx = NULL, ny = NULL, col = "lightgray", lty = "dotted",  lwd = par("lwd"), equilogs = TRUE)
# abline(v=breakdates(TSdata_brk), col="gray70", lty=3, lwd=2, untf=TRUE)
# lines(fitted(TSdata_brk, breaks = length(TSdata_brk$breakpoints)), col = "red",lty = 1, lwd = 2)

#####################
# Detrending
#####################
Reg<-lm(TSdata ~ 1+tt*breakfactor(TSdata_brk, breaks = nbreaks))
detre = Reg$residuals

detre = TSdata-fitted(TSdata_brk, breaks = length(TSdata_brk$breakpoints))
# detre_brk <- breakpoints(detre ~ 1, h = 0.1)


brkpts = c(1,TSdata_brk$breakpoints,length(TSdata))  #Inicio y final de regimenes
#####################
# Long memory estimations without trend
#####################

est_longmemo = data.frame(gph_est=rep(0,nbreaks+2),gph_std=rep(0,nbreaks+2),elw_est=rep(0,nbreaks+2),elw_std=rep(0,nbreaks+2)) #,hp_est=rep(0,nbreaks+2),hp_std=rep(0,nbreaks+2)
####################
#Memory of the complete series
####################
est_longmemo$gph_est[1] = fdGPH(detre,bdw)[1]
est_longmemo$gph_std[1] = fdGPH(detre,bdw)[2]
est_longmemo$elw_est[1] = ELW(detre,m=floor(1+length(TSdata) ^bdw))[1]
est_longmemo$elw_std[1] = ELW(detre,m=floor(1+length(TSdata) ^bdw))[2]
# est_longmemo$hp_est[1] = Hou.Perron(detre,m=floor(1+length(TSdata) ^bdw))[1]
# est_longmemo$hp_std[1] = Hou.Perron(detre,m=floor(1+length(TSdata) ^bdw))[2]

####################
# Memory in regimes
#################
for(j in 1:(nbreaks+1)){
  T = brkpts[j+1]-brkpts[j]+1
  est_longmemo$gph_est[j+1] = fdGPH(detre[brkpts[j]:brkpts[j+1]],bdw)[1]
  est_longmemo$gph_std[j+1] = fdGPH(detre[brkpts[j]:brkpts[j+1]],bdw)[2]
  est_longmemo$elw_est[j+1] = ELW(detre[brkpts[j]:brkpts[j+1]],m=floor(1+T^bdw))[1]
  est_longmemo$elw_std[j+1] = ELW(detre[brkpts[j]:brkpts[j+1]],m=floor(1+T^bdw))[2]
  # est_longmemo$hp_est[j+1] = Hou.Perron(detre[brkpts[j]:brkpts[j+1]],m=floor(1+T^bdw))[1]
  # est_longmemo$hp_std[j+1] = Hou.Perron(detre[brkpts[j]:brkpts[j+1]],m=floor(1+T^bdw))[2]
}

est_longmemo


############################
# Presistence change test of Martins & Rodrigues
############################
test_stat = data.frame(p90 = rep(0,nbreaks), p95 = rep(0,nbreaks), p99 = rep(0,nbreaks), stat = rep(0,nbreaks))

for(j in 1:nbreaks){
  test_stat[j,] = MR_test(detre[brkpts[j]:brkpts[j+2]],trend = "none")[3,1:4] #serial="TRUE",simu=1,M=1000
}
test_stat

TSdata_brk
est_longmemo
test_stat