R/foref.R

In gamalamboy/stresstest: Implement the ETI Stress Testing Simulation

#' Update the design matrix with a forecasted data point.
#' 
#' Given a model structure and a \code{OneStep} object to use as
#' new data value, create a new design matrix \eqn{Z_{t}} from
#' \eqn{ Z_{t-1}}.
#' 
#' @details this function operates in three (3) parts:\enumerate{
#' \item shift old values of the variables by 1 using the lag operator.
#' \item reinsert shifted values in the proper positions of \eqn{Z_{t}}.
#' \item replace lag 1 values with values from the \code{OneStep} object.}
#' This produces a design matrix \eqn{Z_{t}} for the forecasted year.
#' 
#' @param dats_withexo the dataset used to create the forecast. It consists of
#' an enlarged \code{zoo} object with NA's in the forecast rows (except for
#' exogenous variables). The dataset must be generated by \code{DataInput()}.
#' @param period the number of years to forecast. As of yet, the function
#' does not support backcasting.
#' @param ztprop1 the model structure (a \code{Ztprops} object) to be used
#' in determining the proper positions of each value.
#' @param gibbs_mean the \emph{Mean} component of the \code{Gibbs} simulation.
#' @param gibbs_var the \emph{Var} component of the \code{Gibbs} simulation.
#' 
#' @return an array with dimensions (<number_of_forecasted_years>, 
#' <number_of_equations>, <number_of_iterations_in_Gibbs>).
#' @export
ForecastGibbs <- function( dats_withexo, period, ztprop1, gibbs_mean, gibbs_var)
{
    freqs <- attr(ztprop1, "nterms")
    gme <- as.vector( gibbs_mean, mode = "numeric")
    gva <- as.matrix( gibbs_var)
    coln <- colnames(dats_withexo)
    maxlag <- max(ztprop1$lags)
    fl <- length(freqs); fs <- sum(freqs); pl <- length(period)
    dexo0 <- zoo::coredata(dats_withexo)
    attr(dexo0, "yrs") <- stats::time(dats_withexo)

# ofore_index is the index matching the columns of the data set with the
# design matrix rows.

    vv <- match(names(freqs), coln)
    ofore_index <- sub2inds( maxlag + 1, maxlag + 1, vv) - 1
    
    datindex <- sub2inds( maxlag+1, -ztprop1$lags + maxlag + 1,
    match( ztprop1$vars_indep,  coln)) - 1
    
    Zindex <- sub2inds( fl, rep( 1:fl, freqs), 1:fs ) - 1
    
    basep <- which( stats::time( dats_withexo) == min(period) - 1)
    
    Zmat0 <- Zarray( datindex, dexo0[ seq( to = basep, length.out = 
    maxlag + 1),, drop = FALSE], Zindex, matrix(0, fl, fs))
    
    res1 <- parallel::mclapply( 1:ncol(gibbs_mean), function(x)
    {   
        res1 <- OneStep( dexo0, Zmat0, pl, maxlag, gibbs_mean[,x], 
                         gibbs_var[,,x], ofore_index, datindex, Zindex)
    })
    
    res1 <- array( unlist(res1), dim = c( dim( res1[[1]]), length(res1)))
    dimnames(res1) <- list( period, names(freqs), NULL)
    
    res1
}