#####################
# Econometrics SzD  #
# user functions    #
#####################



#####################################################################
# Function that plots IRF                                           #
#####################################################################

MyPlot <- function(tab,main="Impulse response function",xlab="horizon",ylab="",hline=0){
  require(ggplot2)
  ggplot(data=tab, aes(x=tab[,1])) + 
    geom_line(aes(y=tab[,2]), colour="blue", size=2) + 
    theme_light()+
    theme(text = element_text(size=12))+
    labs(title=main, y=ylab, x=xlab, caption="")+
    geom_hline(yintercept=hline)
}

#####################################################################
# Function returns realizations of y                                 #
# Output:                                                           #                       
# act   - realization for y  [matrix of dimension (T-T1)xh]         #
# Inputs:                                                           #
# y     - data (size Tx1)                                           #                      
# horiz - horizon of forecasts                                      # 
# T1    - period for first forecast (forecast evalutaion from T1+1) #
#####################################################################

Actfct <- function(y){
  
  # definition of matrices with forecasts
  act             <- matrix(NA,T-T1,horiz)
  for (n in 0:(T-T1-1)){
    h <- min(horiz,T-T1-n)
    act[n+1,1:h] <- y[(T1+n+1):((T1+n+h))]
  }
  return(act) 
}

#####################################################################
# Funtion returns recursive forecasts for y from RW                 #
# Output:                                                           #                       
# fct   - forecasts for y  [matrix of dimension (T-T1)xh]           #
# Inputs:                                                           #
# y     - data (size Tx1)                                           #                      
# T1    - period for first forecast (forecast evalutaion from T1+1  #
#####################################################################

RWfct <- function(y){
  T <- as.numeric(length(y))
  
  # definition of matrices with forecasts
  fct             <- matrix(NA,T-T1,horiz)
  fctSD           <- matrix(NA,T-T1,horiz)
  
  for (n in 0:(T-T1-1)){
    fct[n+1,]    <- rep(1,horiz)*y[(T1+n)]
    fctSD[n+1, ] <- sqrt(1:horiz)*sd(y[1:(T1+n)])
  }
  return(list(fct=fct,fctSD = fctSD))
}



#####################################################################
# Funtion returns recursive forecasts for y from ARMA(p,q) model    #
# Output:                                                           #                       
# fct   - forecasts for y  [matrix of dimension (T-T1)xh]           #
# fctSD   - forecasts std dev. for y                                #
# Inputs:                                                           #                       
# y     - data (size Tx1)                                           #                      
# p,q   - ARMA model lags                                           #
#####################################################################

ARMAfct <- function(y, p, q){
  require(forecast)
  y      <- ts(y)
  fct    <- matrix(NA,T-T1,horiz); rownames(fct)=rownames(act)
  fctSD <- matrix(NA,T-T1,horiz);  rownames(fctSD)=rownames(act)
  for (n in 0:(T-T1-1)){
    z            <- ts(y[1:(T1+n)])   # recursive sample
    temp          <- Arima(z, order=c(p,0,q))
    # point forecast
    fct[n+1, ]    <- forecast(object=temp, h=horiz)$mean
    # standard deviation
    IRF          <- c(1,ARMAtoMA(ar = coef(temp)[1:p],ma = coef(temp)[(p+1):(p+q)],horiz-1))
    fctSD[n+1, ] <- sqrt(temp$sigma2*cumsum(IRF^2))           
  }
  return(list(fct=fct,fctSD = fctSD))
}


#####################################################################
# Funtion returns recursive forecasts for y from VAR(p) model       #
# Output:                                                           #                       
# fct   - forecasts for y  [matrix of dimension (T-T1)xh]           #
# fctSD   - forecasts std dev. for y                                #
# Inputs:                                                           #
# y      - vector of observable (size TxN)                          #                      
# varN   - variable number                                          #
# p      - lags of VAR model                                        #
#####################################################################

VARfct <-function(y, p, varN){
  y <- ts(y)
  
  # definition of matrices with forecasts
  fct    <- matrix(NA,T-T1,horiz); rownames(fct)=rownames(act)
  fctSD <- matrix(NA,T-T1,horiz);  rownames(fctSD)=rownames(act)
  for (n in 0:(T-T1-1)){
    z       <- y[1:(T1+n),]   # recursive sample
    temp    <- VAR(z, p)
    # point forecast
    tempF   <- predict(temp, n.ahead=horiz)
    fct[n+1,1:horiz] <- tempF$fcst[[varN]][,"fcst"]
    # forecast standard deviation
    fctSD[n+1,1:horiz] <- tempF$fcst[[varN]][,"CI"]/qnorm(0.975)
  }
  return(list(fct=fct,fctSD = fctSD))
}


#####################################################################
# Calculates DM test (two-sided)                                    #
# Example:                                                          #                       
# stat <- DMstat(err1, err2)                                        #
# DMstat                                                            # 
#####################################################################

DMstat <- function(err1, err2, type="two.sided", power=2){
  require(forecast)
  h       <- ncol(err1)
  DMstat  <- rep(NA,h)
  for (n in 1:h){
    e1 <- na.omit(err1[,n])
    e2 <- na.omit(err2[,n])
    temp       <- dm.test(e1, e2, alternative=type, h=n, power=power)
    DMstat[n]  <- temp$p.value
  }
  return(DMstat)
}


#######################################
# Getting matrices y and x for VAR(P) #
#######################################

VARmat <- function(z,P){
  K     <- ncol(z)  
  temp  <- z
  names <- colnames(z)
  for(p in 1:P){
    temp  <- cbind(temp,stats::lag(z,-p))
    names <- c(names,paste0(colnames(z),p))
  }
  temp <- na.omit(temp)
  y           <- temp[,1:K]
  colnames(y) <- names[1:K]
  x           <- temp[,(K+1):ncol(temp)]
  x           <- cbind(1,x)
  colnames(x) <- c("const",names[(K+1):ncol(temp)])
  return(list(y=y, x=x))
}                            

#################################################################
# IRF for BVAR with N=2 and P=2, given B and Sigma draw         #
#################################################################

singleIRF <- function(B,Sig,H){
  IRF   <- array(0,dim=c(K,K,H))
  Pstar <- t(chol(Sig))
  A     <- t(B[2:L,1:K])
  Atemp <- cbind(diag(1,K*(P-1)),matrix(0,K*(P-1),K))
  Acan  <- rbind(A,Atemp)
  temp  <- rbind(Pstar,diag(0,K*(P-1)))
  for(k in 1:H){
    IRF[,,k] <- temp[1:K,1:K]
    temp     <- Acan%*%temp
  }
  return(IRF)
}


#################################################################
# single draw from predictive density for BVAR with N=2 and P=2 #
#################################################################

singleFCT <- function(B,Sig,H,y,shock=FALSE){
  require(MASS)
  Fct  <- matrix(0,H,K)
  y  <- as.matrix(y)
  yT <- c(1,vec(t(tail(y,P)[P:1,])))
  for(k in 1:H){
    temp    = t(B)%*%yT
    if (shock){
      eps = mvrnorm(n=1, mu = rep(0,K), Sigma = Sig)
      temp=temp+eps
    }
    Fct[k,] = t(temp) 
    yT = c(1,temp,yT[2:(L-K)])
  }
  return(Fct)
}