##############################
# Econometrics SzD           #
# Theme 8 - Ex-post Foreacst #
##############################
rm(list=ls())
source("libs.R")
source("utils.R")

# a. Loading data
####################

# filename <- "dataJIE.csv"
# filename <- "http://web.sgh.waw.pl/~mrubas/AdvEcon/data/dataJIE.csv"
# data     <- read.csv(filename, header=TRUE, sep=";", dec=",", row.names="daty")

filename <- "dataJIE.RData"
load(filename)

x     <- cbind(data["Y_US"],data["Y_UK"])   
dates <- as.yearqtr(as.Date(rownames(data),format="%d.%m.%Y"))

z  <- zoo(x, dates)
y  <- ts(x,start=c(1975,1), freq=4)

# Start of forefcasting competition

# b. settings
T          <- dim(y)[1] 	                            # number of observations
T1         <- 80                                      # 20 years to estimate the first model
horiz      <- 12                                      # forecast horizon
varN       <- 2                                       # which variable 

# c. Calculating realizations, point forecasts and forecast standatrd deviation

act           <- Actfct(y[,varN])
rownames(act) <- as.character(index(z)[T1:(T-1)])
# View(act)

# d. Calculating point forecasts and ex-ante forecast error 

fctRW    <- RWfct(y[,varN])            # here it makes little sense, why?
fctARMA  <- ARMAfct(y[,varN], 2, 0)
fctVAR   <- VARfct(y, 2, varN)

# relative standard deviation
# View(fctVAR$fctSD / fctARMA$fctSD)

# e. Sequantial forecast graph

# Choose model "RW", "ARMA", VAR"
fct   = fctVAR$fct

datesF = c(index(z),last(index(z)) + (1:horiz)/4)
plot_data    = data.frame(t=index(z[,varN]),y=coredata(z[,varN]))
temp_plot   <- MyPlot(plot_data,main="Sequential forecast",xlab="",hline=mean(z[,varN]))
for (n in 0:(T-T1-1)) {
  fct_data = data.frame(yf = fct[n+1,], tf = datesF[(T1+n+1):(T1+n+horiz)])
  temp_plot <- temp_plot + geom_line(data=fct_data, aes(y=yf,x=tf), colour="red", size=1)
}
temp_plot + 
  xlim(1990,2018) +
  ylim(-3,3)


# f. Effectiveness test graph

h     = 12                          # Horizon
fct   = fctVAR$fct                  # Choose model "RW", "ARIMA", "VAR"

a  <- act[1:(T-T1-h+1),h] 
f  <- fct[1:(T-T1-h+1),h] 

par(mfrow=c(1,1), cex = 0.7, bty="l")
a1 <- min(cbind(a,f))-0.01; a2 <- max(cbind(a,f))+0.01; 
plot(f,a, main="level", xlab="forecast", ylab="realization", pch=19, col="red",  xlim=c(a1,a2), ylim = c(a1,a2))
abline(0,1,lwd=2)

summary(lm(a ~ f))

# out of sample forecast evaluation
Hsel = c(1,2,3,6,9,12)

# g. Point forecasts evaluation 
errRW     <- act - fctRW$fct  
errARMA   <- act - fctARMA$fct  
errVAR    <- act - fctVAR$fct 

# mean error
ME <- cbind(colMeans(errRW, na.rm = T) , 
            colMeans(errARMA, na.rm = T), 
            colMeans(errVAR, na.rm = T))
ME   <- ME[Hsel,]
colnames(ME) <- c("RW", "ARMA","VAR")
rownames(ME) <- Hsel
kable(ME,digits=2)

# root mean squared error
RMSE <- cbind(  sqrt(colMeans(errRW^2, na.rm = T)),
                sqrt(colMeans(errARMA^2, na.rm = T)),
                sqrt(colMeans(errVAR^2, na.rm = T)))
RMSE   <- RMSE[Hsel,]
colnames(RMSE) <- c("RW", "ARMA","VAR")
rownames(RMSE) <- Hsel
kable(RMSE,digits=2)

# We prefer to look at relaive RMSE
RMSErel       = RMSE
RMSErel[,2:3] = RMSE[,2:3]/RMSE[,1] 
kable(RMSErel,digits=2)

# Diebold-Mariano -- is the ratio significantly different than 1

DMtest <- cbind(DMstat(errARMA, errRW, type="two.sided", power=2), 
                DMstat(errVAR, errRW, type="two.sided", power=2))
DMtest   <- DMtest[Hsel,]
colnames(DMtest) <- c("ARMA_DM","VAR_DM")
rownames(DMtest) <- Hsel

kable(cbind(RMSErel,DMtest),digits=2)

# DM test, step-by-step
h   = 12
e1h = errRW[,h]
e2h = errARMA[,h]
dh  = na.omit(e1h^2 - e2h^2)
plot.ts(dh)
abline(h=0,col=2)
abline(h=mean(dh),col=3)

dhat = mean(dh)
s    = sd(dh)
Th   = length(dh)

sig2lr = lrvar(dh, type = "Newey-West") # from sandwich package
sig2   = var(dh)/Th

DM = dhat/sqrt(sig2lr)
DM
pnorm(-abs(DM))


# h. Density forecasts evaluation

# ilustration
Tf = 1
h  = 4

x = seq(-5,5,0.01)
temp <- cbind(x,
              dnorm(x,fctRW$fct[Tf,h],fctRW$fctSD[Tf,h]),
              dnorm(x,fctARMA$fct[Tf,h],fctARMA$fctSD[Tf,h]),
              dnorm(x,fctVAR$fct[Tf,h],fctVAR$fctSD[Tf,h]))
colnames(temp) <- c("x","RW","ARMA","VAR")

# temp[,2:4] = log(temp[,2:4]) # uncomment to change from predictive scores to Log PS 

ggplot(data = data.frame(temp), aes(x = x)) + 
  geom_line(aes(y = RW,colour="RW")) + 
  geom_line(aes(y = ARMA,colour="ARMA")) +
  geom_line(aes(y = VAR,colour="VAR"))+
  labs(title="Predictive Score", y="", x="", caption="")+
  theme_light()+
  scale_colour_manual("", 
                        breaks = c("RW","ARMA","VAR"),
                        values = c("red", "blue", "green"))+
  geom_vline(xintercept=act[Tf,h])
  
# which model is best for this single observation?

# log-redictive scores (this is equivalent to errors in point forecasts)
lpsRW   <- log(dnorm(act,fctRW$fct,fctRW$fctSD))
lpsARMA <- log(dnorm(act,fctARMA$fct,fctARMA$fctSD))
lpsVAR  <- log(dnorm(act,fctVAR$fct,fctVAR$fctSD))

LPS <- cbind( colMeans(lpsRW, na.rm = T),
              colMeans(lpsARMA, na.rm = T),
              colMeans(lpsVAR, na.rm = T))
LPS <- LPS[Hsel,]
colnames(LPS) <- c("RW", "ARMA","VAR")
rownames(LPS) <- Hsel
kable(LPS,digits=2)

# We prefer to look at LPS difference (% gain)
LPSrel       = LPS
LPSrel[,2:3] = LPS[,2:3]-LPS[,1] 
kable(LPSrel,digits=2)

# Geweke-Amisano test -- is the difference significantly different from 0 

GAtest <- cbind(DMstat(lpsARMA, lpsRW, type="two.sided", power=1), 
                DMstat(lpsVAR,  lpsRW, type="two.sided", power=1))
GAtest   <- GAtest[Hsel,]
colnames(GAtest) <- c("ARMA_GA","VAR_GA")
rownames(GAtest) <- Hsel

kable(cbind(LPSrel,GAtest),digits=2)



# PIT - Probability Integral Transform
pitRW   <- pnorm(act,fctRW$fct,fctRW$fctSD)
pitARMA <- pnorm(act,fctARMA$fct,fctARMA$fctSD)
pitVAR  <- pnorm(act,fctVAR$fct,fctVAR$fctSD)

# PIT graph
h    <- 1    # select horizon
pitM <- pitVAR
pit  <- na.omit(pitM[,h])
hist(pit, breaks = seq(0,1,0.1), prob = TRUE, col="grey")
abline(h=1, lwd=2)


