#' @export
#' @title Conditional forecasts
#' @param obj Estimated model
#' @param forecastHorizon forecast horizon
#' @param cfconds An \eqn{h\times n} matrix containing the forecast conditions with \eqn{h} being the forecast horizon and \eqn{n} the number of variables in the VAR. Unconstrained variables should be NA.
#' @param id_obj Identification
#' @param interval forecast bands
#' @param ... not used
#' @return returns an S3 object of the class fcbvar
#' @details
#' Conditional forecasts are forecasts conditional on given values for a subset of variables and are obtained by pre-determining the path of certain variables. Waggoner and Zha (1999) show that the distribution of future shocks is normal with
#' \deqn{\eta\sim N(\bar{\eta},\bar\Gamma)}
#' where
#' \deqn{\bar{\eta}=R'(RR')^{-1}r}
#' and
#' \deqn{\Gamma=I-R'(RR')^{-1}R}
#' with \eqn{\eta} the \eqn{s\times1}-vector of structural shocks and \eqn{r} is the vector of differences between predicted and conditional values.
#' Instead of drawing directly from the above distribution we use a singular value decomposition of \eqn{R} as proposed by
#'
#' @rdname cforecast
#' @references Waggoner, Daniel F. and Tao Zha, Conditional Forecasts in Dynamic Multivariate Models, The Review of Economics and Statistics, Vol. 81, No. 4 (Nov 1999), pp. 639-651
#' @references Jarocinski, M., Conditional forecasts and uncertainty about forecast revisions in vector autoregressions, Economics Letters, Vol. 25, No. 3 (2010), pp. 257-259
cforecast <- function(obj,forecastHorizon,id_obj,cfconds,interval = c(0.05,0.95),...) UseMethod("cforecast")
#' @export
#' @rdname cforecast
cforecast.bvar <- function(obj,forecastHorizon,id_obj,cfconds,interval = c(0.05,0.95),...){
# Step one: unconditional forecasts
# Preliminary Calculations
nVariables <- dim(obj$data_info$data)[2] # Number of variables
nLength <- dim(obj$data_info$data)[1] # Length of time series
nLags <- obj$general_info$nolags # Number of lags
nForecasts <- dim(obj$mcmc_draws$Alpha)[3] # Number of forecasts, depends on the number sampled posteriors
mForecastStorage <- array(NA,dim=c(forecastHorizon,nVariables,nForecasts)) # Matrix to storage forecasts
cdForecast <- array(NA,dim=c(forecastHorizon,nVariables))
fc_tmp <- array(NA,dim=c(forecastHorizon+nLags,nVariables)) # Matrix to temporary store forecasts
# If the data is of class ts, get start date and frequency of data
if(stats::is.ts(obj$data_info$data)){
tsStart <- stats::start(obj$data_info$data)
tsFrequency <- stats::frequency(obj$data_info$data)
}
tsColnames <- colnames(obj$data_info$data)
# Check if mydata is a time series object
# If it is a time series object get frequency, etc.
if(stats::is.ts(obj$data_info$data)){
nFreq <- frequency(obj$data_info$data)
nFirstDate <- min(stats::time(obj$data_info$data))
nFirstYear <- floor(nFirstDate)
nFirstMonthQuarter <- (nFirstDate-nFirstYear) * nFreq
nLastDate <- max(stats::time(obj$data_info$data))
}
# Do the forecasts
for(ii in 1:nForecasts){
fc_tmp <- array(NA,dim=c(forecastHorizon+nLags,nVariables)) # Matrix to temporary store forecasts
fc_tmp[1:nLags,] <- (obj$data_info$data[nLength:(nLength-nLags+1),])
# Step 1: Unconditional forecast for one path
for(jj in 1:forecastHorizon){
nstart <- jj
nEnd <- jj + nLags - 1
y <- fc_tmp[nstart:nEnd,]
randDraw <- rnorm(nVariables) %*% t(chol(obj$mcmc_draws$Sigma[,,ii]))
if(obj$general_info$intercept == TRUE){
y <- c(1,t(y))
tempForecast <- y %*% obj$mcmc_draws$Alpha[,,ii] + randDraw
}
else{
y <- t(y)
tempForecast <- y %*% obj$mcmc_draws$Alpha[,,ii] + randDraw
}
# Store the unconditional
fc_tmp[jj + nLags,] <- tempForecast
}
# Step 2: Calculate Impulse-Response Function
ortirf <- compirf(Alpha = obj$mcmc_draws$Alpha[,,ii], Sigma = obj$mcmc_draws$Sigma[,,ii], id_obj = id_obj,
nolags = nLags, intercept = obj$general_info$intercept, nhor = forecastHorizon)
# Step 3: create R matrix and r vector
R <- vector()
r <- vector()
for(jj in 1:forecastHorizon){
for(kk in 1:nVariables){
if(!is.na(cfconds[jj,kk])){
# r-vector condition - unconditional forecast
r_tmp <- matrix(cfconds[jj,kk] - fc_tmp[jj + nLags,kk],nrow=1)
r <- rbind(r,r_tmp)
# R-Matrix
R_tmp <- matrix(0,nrow=1,ncol = nVariables * forecastHorizon)
for(ll in 1:jj){
R_tmp[1,((ll - 1) * nVariables + 1):(ll * nVariables)] <- ortirf[kk,,(jj - ll + 1)]
}
R <- rbind(R,R_tmp)
}
}
}
# Step 4: Draw shocks from the Waggoner-Zha distribution
Rdim <- dim(R)
Q <- Rdim[1]
K <- Rdim[2]
Rsvd <- svd(R,nu=nrow(R),nv=ncol(R))
U <- Rsvd$u
V <- Rsvd$v
S <- Rsvd$d
P <- diag(S)
V1 <- V[,1:Q]
V2 <- V[,-c(1:Q)]
eta <- matrix(pracma::mrdivide(V1,P) %*% t (U) %*% r + V2 %*% rnorm(n = K - Q),ncol=forecastHorizon)
# Step 5: Obtain the conditional forecasts
for(jj in 1:forecastHorizon){
temp <- 0
for(kk in 1:jj){
temp <- temp + ortirf[,,jj - kk + 1] %*% eta[,kk]
}
cdForecast[jj,] <- fc_tmp[nLags + jj,] + t(temp)
}
# Step 6: Save forecasts
mForecastStorage[,,ii] <- cdForecast
}
# Summarize results
forecastMean <- array(0,dim=c(forecastHorizon,nVariables))
forecastUpper <- array(0,dim=c(forecastHorizon,nVariables))
forecastLower <- array(0,dim=c(forecastHorizon,nVariables))
for(ii in 1:nVariables){
for(jj in 1:forecastHorizon){
forecastMean[jj,ii] <- mean(mForecastStorage[jj,ii,])
forecastLower[jj,ii] <- quantile(mForecastStorage[jj,ii,],probs=max(interval))
forecastUpper[jj,ii] <- quantile(mForecastStorage[jj,ii,],probs=min(interval))
}
}
OriginalPath <- array(NA,dim=c(nLength + forecastHorizon, nVariables))
forecastFinalMean <- OriginalPath
forecastFinalUpper <- OriginalPath
forecastFinalLower <- OriginalPath
OriginalPath[1:nLength,] <-obj$data_info$data
forecastFinalMean[(nLength + 1):(nLength + forecastHorizon),] <- forecastMean
forecastFinalUpper[(nLength + 1):(nLength + forecastHorizon),] <- forecastUpper
forecastFinalLower[(nLength + 1):(nLength + forecastHorizon),] <- forecastLower
if(is.ts(obj$data_info$data)){
forecastFinalMean <- ts(forecastFinalMean,start=tsStart,frequency=tsFrequency)
forecastFinalUpper <- ts(forecastFinalUpper,start=tsStart,frequency=tsFrequency)
forecastFinalLower <- ts(forecastFinalLower,start=tsStart,frequency=tsFrequency)
OriginalPath <- ts(OriginalPath,start=tsStart,frequency=tsFrequency)
}
colnames(forecastFinalMean) <- colnames(obj$data_info$data)
colnames(forecastFinalUpper) <- colnames(obj$data_info$data)
colnames(forecastFinalLower) <- colnames(obj$data_info$data)
colnames(OriginalPath) <- colnames(obj$data_info$data)
retList <- structure(list(forecast = forecastFinalMean, Upper = forecastFinalUpper, Lower = forecastFinalLower,
Original = OriginalPath),class="fcbvar")
return(retList)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.