#####################################################################################
# Function that calculates the value of interest rate at maturity tau from NS model #
# beta - the value of parameters   (3x1 vetor)                                      #
# tau  - matuity                   (vector)                                         #
#####################################################################################

NSrate <-function(beta, tau){
  lambda = 0.0609
  
  if(is.matrix(beta)){
    R = beta[,1]+beta[,2]*((1-exp(-lambda*tau))/(lambda*tau)-exp(-lambda*tau))+beta[,3]*((1-exp(-lambda*tau))/(lambda*tau))  
  } else {
    R = beta[1] +beta[2]*((1-exp(-lambda*tau))/(lambda*tau) -exp(-lambda*tau))+beta[3]*((1-exp(-lambda*tau))/(lambda*tau))  
  }
  return(R)
}

#####################################################################
# Funtion returns recursive forecasts for y from baseline NS model  #
# Output:                                                           #                       
# fct   - forecasts for y  [matrix of dimension (T-T1)xh]           #
# Inputs:                                                           #
# beta  - matrix of NS factors [Tx3]                                #
# daty  - vector of dates [Tx1]                                     #
# horiz - horizon of forecasts (months)                             # 
# T1    - period for first forecast (forecast evalutaion from T1+1) #
# Mat   - maturity (in month)                                       #  
#####################################################################

NSfct <- function(beta, daty, horiz, T1, Mat){
  T <- nrow(beta)
  
  # definition of matrices with forecasts
  fct             <- matrix(NA,T-T1,horiz)
  rownames(fct)   <- as.character(daty[T1:(T-1)])
  
  for (n in 0:(T-T1-1)){
    beta0        = beta[T1+n,]   #NS parameters used for forecast
    for (h in 1:horiz){
      fct[n+1,h]    = (NSrate(beta0,h+Mat)*(h+Mat) - NSrate(beta0,h)*h)/Mat
    }
  }
  
  return(fct) 
}

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

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


#####################################################################
# Funtion returns recursive forecasts for y from ARMA(p,q) model    #
# Output:                                                           #                       
# fct   - forecasts for y  [matrix of dimension (T-T1)xh]           #
# Inputs:                                                           #                       
# y     - data (size Tx1)                                           #                      
# daty  - vector of dates                                           #
# horiz - horizon of forecasts (months)                             # 
# T1    - period for first forecast (forecast evalutaion from T1+1  #
# p,q   - lags of ARMA model                                        #
#####################################################################

ARMAfct <- function(y, daty, horiz,  T1, p, q){
  require(forecast)
	    T <- as.numeric(length(y))
      
      # definition of matrices with forecasts
      fct             <- matrix(NA,T-T1,horiz)
      rownames(fct)   <- as.character(daty[T1:(T-1)])
      

	for (n in 0:(T-T1-1)){
	    z            <- y[1:(T1+n)]   # recursive sample
      junk1        <- Arima(z, order=c(p,0,q))
      junk2        <- forecast::forecast(object=junk1, h=horiz)$mean
      fct[n+1,1:horiz] <- junk2
	                      }
      #detach(package:forecast)
      return(fct) 
}

#####################################################################
# Funtion returns recursive forecasts for y from DL-ARMA(p,q) model #
# Output:                                                           #                       
# fct   - forecasts for y  [matrix of dimension (T-T1)xh]           #
# Inputs:                                                           #
# beta  - vector of DL factors (size Tx3)                           #                      
# daty  - vector of dates                                           #
# horiz - horizon of forecasts                                      #
# T1    - period for first forecast (forecast evalutaion from T1+1  #
# Mat   - maturity of interest rate                                 #
# p,q   - lags of ARMA model                                        #
#####################################################################

DLARMAfct <-function(beta, daty, horiz, T1, Mat, p=2, q=0){
  require(forecast)
  T <- nrow(beta)
  
  # definition of matrices with forecasts
  fct             <- matrix(NA,T-T1,horiz)
  rownames(fct)   <- as.character(daty[T1:(T-1)])
  
  
  for (n in 0:(T-T1-1)){
    z            <- beta[1:(T1+n),]   # recursive sample
    junk         <- arima(z[,1], order=c(p,0,q))
    LF           <- forecast::forecast(object=junk, h=horiz)$mean
    junk         <- arima(z[,2], order=c(p,0,q))
    MF           <- forecast::forecast(object=junk, h=horiz)$mean
    junk         <- arima(z[,3], order=c(p,0,q))
    SF           <- forecast::forecast(object=junk, h=horiz)$mean
    
    betaF        <- cbind(LF,MF,SF)
    fct[n+1,1:horiz] <- NSrate(betaF,Mat)
  }
  #detach(package:forecast)
  return(fct) 
}

#####################################################################
# Funtion returns recursive forecasts for y from DL-VAR(p) model    #
# Output:                                                           #                       
# fct   - forecasts for y  [matrix of dimension (T-T1)xh]           #
# Inputs:                                                           #
# beta  - vector of DL factors (size Tx3)                           #                      
# daty  - vector of dates                                           #
# horiz - horizon of forecasts                                      #
# T1    - period for first forecast (forecast evalutaion from T1+1  #
# Mat   - maturity of interest rate                                 #
# p     - lags of ARMA model                                        #
#####################################################################

DLVARfct <-function(beta, daty, horiz, T1, Mat, p){
  require(vars)
  T <- nrow(beta)
  
  # definition of matrices with forecasts
  fct             <- matrix(NA,T-T1,horiz)
  rownames(fct)   <- as.character(daty[T1:(T-1)])
  
  
  for (n in 0:(T-T1-1)){
    z        <- beta[1:(T1+n),]   # recursive sample
    junk     <- VAR(z, p)	
    fctBeta  <- predict(junk, n.ahead=horiz)
    LF       <- fctBeta$fcst$LT[,"fcst"]
    MF       <- fctBeta$fcst$MT[,"fcst"]
    SF       <- fctBeta$fcst$ST[,"fcst"]
    
    betaF        <- cbind(LF,MF,SF)
    fct[n+1,1:horiz] <- NSrate(betaF,Mat)
  }
  #detach(package:forecast)
  return(fct) 
}