23 update the forecast method

NAMESPACE

@@ -1,6 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
 S3method(estimate,BSVARSIGN)
+S3method(forecast,PosteriorBSVARSIGN)
 export(specify_bsvarSIGN)
 export(specify_identification_bsvarSIGN)
 export(specify_posterior_bsvarSIGN)

R/forecast.R

@@ -0,0 +1,132 @@
+
+#' @title Forecasting using Structural Vector Autoregression
+#'
+#' @description Samples from the joint predictive density of all of the dependent 
+#' variables for models from packages \pkg{bsvars}, \pkg{bsvarSIGNs} or 
+#' \pkg{bvarPANELs} at forecast horizons from 1 to \code{horizon} specified as 
+#' an argument of the function.
+#' 
+#' @method forecast PosteriorBSVARSIGN
+#' @param posterior posterior estimation outcome - an object of class 
+#' \code{PosteriorBSVARSIGN} obtained by running the \code{estimate} function.
+#' @param horizon a positive integer, specifying the forecasting horizon.
+#' @param exogenous_forecast a matrix of dimension \code{horizon x d} containing 
+#' forecasted values of the exogenous variables. 
+#' @param conditional_forecast a \code{horizon x N} matrix with forecasted values 
+#' for selected variables. It should only contain \code{numeric} or \code{NA} 
+#' values. The entries with \code{NA} values correspond to the values that are 
+#' forecasted conditionally on the realisations provided as \code{numeric} values.
+#' 
+#' @return A list of class \code{Forecasts} containing the
+#' draws from the predictive density and data. The output list includes element:
+#' 
+#' \describe{
+#'  \item{forecasts}{an \code{NxhorizonxS} array with the draws from predictive density}
+#'  \item{Y}{an \eqn{NxT} matrix with the data on dependent variables}
+#' }
+#' 
+#' @author Tomasz Woźniak \email{[email protected]} and Xiaolei Wang \email{[email protected]}
+#' 
+#' @examples
+#' # upload data
+#' data(oil)
+#' 
+#' # specify the model and set seed
+#' set.seed(123)
+#' sign_irf       = array(matrix(c(-1, -1, 1, rep(0, 6)), nrow = 3), dim = c(3, 3, 1))
+#' specification  = specify_bsvarSIGN$new(oil, p = 12, sign_irf = sign_irf)
+#' 
+#' # estimate the model
+#' posterior      = estimate(specification, 20)
+#' 
+#' # sample from predictive density 1 year ahead
+#' predictive     = forecast(posterior, 4)
+#' 
+#' # workflow with the pipe |>
+#' ############################################################
+#' set.seed(123)
+#' oil |>
+#'   specify_bsvarSIGN$new(p = 12, sign_irf = sign_irf) |>
+#'   estimate(S = 20) |> 
+#'   forecast(horizon = 4) -> predictive
+#' 
+#' # conditional forecasting 2 quarters ahead conditioning on 
+#' #  provided future values for the Gross Domestic Product 
+#' ############################################################
+#' cf        = matrix(NA , 2, 3)
+#' cf[,3]    = tail(oil, 1)[3]   # conditional forecasts equal to the last gdp observation
+#' predictive    = forecast(posterior, 2, conditional_forecast = cf)
+#' 
+#' # workflow with the pipe |>
+#' ############################################################
+#' set.seed(123)
+#' oil |>
+#'   specify_bsvarSIGN$new(p = 12, sign_irf = sign_irf) |>
+#'   estimate(S = 10) |> 
+#'   forecast(horizon = 2, conditional_forecast = cf) -> predictive
+#' 
+#' @export
+forecast.PosteriorBSVARSIGN = function(
+    posterior, 
+    horizon = 1, 
+    exogenous_forecast = NULL,
+    conditional_forecast = NULL
+) {
+
+  posterior_Sigma = posterior$posterior$Sigma
+  posterior_A     = posterior$posterior$A
+  T               = nrow(posterior$last_draw$data_matrices$X)
+  X_T             = posterior$last_draw$data_matrices$X[T,]
+  Y               = posterior$last_draw$data_matrices$Y
+
+  N               = nrow(posterior_Sigma)
+  K               = length(X_T)
+  d               = K - N * posterior$last_draw$p - 1
+  S               = dim(posterior_Sigma)[3]
+
+  # prepare forecasting with exogenous variables
+  if (d == 0 ) {
+    exogenous_forecast = matrix(NA, horizon, 1)
+  } else {
+    stopifnot("Forecasted values of exogenous variables are missing." = (d > 0) & !is.null(exogenous_forecast))
+    stopifnot("The matrix of exogenous_forecast does not have a correct number of columns." = ncol(exogenous_forecast) == d)
+    stopifnot("Provide exogenous_forecast for all forecast periods specified by argument horizon." = nrow(exogenous_forecast) == horizon)
+    stopifnot("Argument exogenous has to be a matrix." = is.matrix(exogenous_forecast) & is.numeric(exogenous_forecast))
+    stopifnot("Argument exogenous cannot include missing values." = sum(is.na(exogenous_forecast)) == 0 )
+  }
+
+  # prepare forecasting with conditional forecasts
+  if ( is.null(conditional_forecast) ) {
+    # this will not be used for forecasting, but needs to be provided
+    conditional_forecast = matrix(NA, horizon, N)
+  } else {
+    stopifnot("Argument conditional_forecast must be a matrix with numeric values."
+              = is.matrix(conditional_forecast) & is.numeric(conditional_forecast)
+    )
+    stopifnot("Argument conditional_forecast must have the number of rows equal to 
+              the value of argument horizon."
+              = nrow(conditional_forecast) == horizon
+    )
+    stopifnot("Argument conditional_forecast must have the number of columns 
+              equal to the number of columns in the used data."
+              = ncol(conditional_forecast) == N
+    )
+  }
+
+  # perform forecasting
+  for_y       = .Call(`_bsvarSIGNs_forecast_bsvarSIGNs`, 
+                      posterior_Sigma,
+                      posterior_A,
+                      X_T,
+                      exogenous_forecast,
+                      conditional_forecast,
+                      horizon
+  ) # END .Call
+
+  fore            = list()
+  fore$forecasts  = for_y
+  fore$Y          = Y
+  class(fore)     = "Forecasts"
+
+  return(fore)
+} # END forecast.PosteriorBSVARSIGN

inst/include/bsvarSIGNs_RcppExports.h

@@ -109,6 +109,27 @@ namespace bsvarSIGNs {
         return Rcpp::as<Rcpp::List >(rcpp_result_gen);
     }
 
+    inline arma::cube forecast_bsvarSIGNs(arma::cube& posterior_Sigma, arma::cube& posterior_A, arma::vec& X_T, arma::mat& exogenous_forecast, arma::mat& cond_forecast, const int& horizon) {
+        typedef SEXP(*Ptr_forecast_bsvarSIGNs)(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP);
+        static Ptr_forecast_bsvarSIGNs p_forecast_bsvarSIGNs = NULL;
+        if (p_forecast_bsvarSIGNs == NULL) {
+            validateSignature("arma::cube(*forecast_bsvarSIGNs)(arma::cube&,arma::cube&,arma::vec&,arma::mat&,arma::mat&,const int&)");
+            p_forecast_bsvarSIGNs = (Ptr_forecast_bsvarSIGNs)R_GetCCallable("bsvarSIGNs", "_bsvarSIGNs_forecast_bsvarSIGNs");
+        }
+        RObject rcpp_result_gen;
+        {
+            RNGScope RCPP_rngScope_gen;
+            rcpp_result_gen = p_forecast_bsvarSIGNs(Shield<SEXP>(Rcpp::wrap(posterior_Sigma)), Shield<SEXP>(Rcpp::wrap(posterior_A)), Shield<SEXP>(Rcpp::wrap(X_T)), Shield<SEXP>(Rcpp::wrap(exogenous_forecast)), Shield<SEXP>(Rcpp::wrap(cond_forecast)), Shield<SEXP>(Rcpp::wrap(horizon)));
+        }
+        if (rcpp_result_gen.inherits("interrupted-error"))
+            throw Rcpp::internal::InterruptedException();
+        if (Rcpp::internal::isLongjumpSentinel(rcpp_result_gen))
+            throw Rcpp::LongjumpException(rcpp_result_gen);
+        if (rcpp_result_gen.inherits("try-error"))
+            throw Rcpp::exception(Rcpp::as<std::string>(rcpp_result_gen).c_str());
+        return Rcpp::as<arma::cube >(rcpp_result_gen);
+    }
+
     inline bool match_sign_narrative(const arma::mat& Epsilon, const arma::mat& sign_narrative, const arma::cube& irf) {
         typedef SEXP(*Ptr_match_sign_narrative)(SEXP,SEXP,SEXP);
         static Ptr_match_sign_narrative p_match_sign_narrative = NULL;

inst/tinytest/test_forecast.R

@@ -0,0 +1,41 @@
+
+data(oil)
+
+# for bsvar
+set.seed(1)
+suppressMessages(
+  specification_no1 <- specify_bsvarSIGN$new(oil, p = 1)
+)
+run_no1             <- estimate(specification_no1, 3, 1, show_progress = FALSE)
+ff                  <- forecast(run_no1, horizon = 2)
+
+set.seed(1)
+suppressMessages(
+  ff2              <- oil |>
+    specify_bsvarSIGN$new(p = 1) |>
+    estimate(S = 3, thin = 1, show_progress = FALSE) |>
+    forecast(horizon = 2)
+)
+
+
+expect_identical(
+  ff$forecasts[1,1,1], ff2$forecasts[1,1,1],
+  info = "forecast: forecast identical for normal and pipe workflow."
+)
+
+expect_true(
+  is.numeric(ff$forecasts) & is.array(ff$forecasts),
+  info = "forecast: returns numeric array."
+)
+
+
+expect_error(
+  specify_bsvar$new(us_fiscal_lsuw) |> forecast(horizon = 3),
+  info = "forecast: wrong input provided."
+)
+
+expect_error(
+  forecast(run_no1, horizon = 1.5),
+  info = "forecast: specify horizon as integer."
+)
+

src/RcppExports.cpp

@@ -173,6 +173,45 @@ RcppExport SEXP _bsvarSIGNs_bsvar_sign_cpp(SEXP SSEXP, SEXP lagsSEXP, SEXP YSEXP
     UNPROTECT(1);
     return rcpp_result_gen;
 }
+// forecast_bsvarSIGNs
+arma::cube forecast_bsvarSIGNs(arma::cube& posterior_Sigma, arma::cube& posterior_A, arma::vec& X_T, arma::mat& exogenous_forecast, arma::mat& cond_forecast, const int& horizon);
+static SEXP _bsvarSIGNs_forecast_bsvarSIGNs_try(SEXP posterior_SigmaSEXP, SEXP posterior_ASEXP, SEXP X_TSEXP, SEXP exogenous_forecastSEXP, SEXP cond_forecastSEXP, SEXP horizonSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::traits::input_parameter< arma::cube& >::type posterior_Sigma(posterior_SigmaSEXP);
+    Rcpp::traits::input_parameter< arma::cube& >::type posterior_A(posterior_ASEXP);
+    Rcpp::traits::input_parameter< arma::vec& >::type X_T(X_TSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type exogenous_forecast(exogenous_forecastSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type cond_forecast(cond_forecastSEXP);
+    Rcpp::traits::input_parameter< const int& >::type horizon(horizonSEXP);
+    rcpp_result_gen = Rcpp::wrap(forecast_bsvarSIGNs(posterior_Sigma, posterior_A, X_T, exogenous_forecast, cond_forecast, horizon));
+    return rcpp_result_gen;
+END_RCPP_RETURN_ERROR
+}
+RcppExport SEXP _bsvarSIGNs_forecast_bsvarSIGNs(SEXP posterior_SigmaSEXP, SEXP posterior_ASEXP, SEXP X_TSEXP, SEXP exogenous_forecastSEXP, SEXP cond_forecastSEXP, SEXP horizonSEXP) {
+    SEXP rcpp_result_gen;
+    {
+        Rcpp::RNGScope rcpp_rngScope_gen;
+        rcpp_result_gen = PROTECT(_bsvarSIGNs_forecast_bsvarSIGNs_try(posterior_SigmaSEXP, posterior_ASEXP, X_TSEXP, exogenous_forecastSEXP, cond_forecastSEXP, horizonSEXP));
+    }
+    Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error");
+    if (rcpp_isInterrupt_gen) {
+        UNPROTECT(1);
+        Rf_onintr();
+    }
+    bool rcpp_isLongjump_gen = Rcpp::internal::isLongjumpSentinel(rcpp_result_gen);
+    if (rcpp_isLongjump_gen) {
+        Rcpp::internal::resumeJump(rcpp_result_gen);
+    }
+    Rboolean rcpp_isError_gen = Rf_inherits(rcpp_result_gen, "try-error");
+    if (rcpp_isError_gen) {
+        SEXP rcpp_msgSEXP_gen = Rf_asChar(rcpp_result_gen);
+        UNPROTECT(1);
+        Rf_error("%s", CHAR(rcpp_msgSEXP_gen));
+    }
+    UNPROTECT(1);
+    return rcpp_result_gen;
+}
 // match_sign_narrative
 bool match_sign_narrative(const arma::mat& Epsilon, const arma::mat& sign_narrative, const arma::cube& irf);
 static SEXP _bsvarSIGNs_match_sign_narrative_try(SEXP EpsilonSEXP, SEXP sign_narrativeSEXP, SEXP irfSEXP) {
@@ -869,6 +908,7 @@ static int _bsvarSIGNs_RcppExport_validate(const char* sig) {
         signatures.insert("arma::cube(*ir1_cpp)(const arma::mat&,const arma::mat&,int,const int&)");
         signatures.insert("arma::field<arma::cube>(*bsvarSIGNs_ir)(arma::cube&,arma::cube&,const int,const int)");
         signatures.insert("Rcpp::List(*bsvar_sign_cpp)(const int&,const int&,const arma::mat&,const arma::mat&,const arma::field<arma::mat>&,const arma::cube&,const arma::mat&,const arma::mat&,const arma::field<arma::mat>&,const Rcpp::List&,const Rcpp::List&,const bool,const int,const int&)");
+        signatures.insert("arma::cube(*forecast_bsvarSIGNs)(arma::cube&,arma::cube&,arma::vec&,arma::mat&,arma::mat&,const int&)");
         signatures.insert("bool(*match_sign_narrative)(const arma::mat&,const arma::mat&,const arma::cube&)");
         signatures.insert("double(*weight_narrative)(const int&,arma::mat,const arma::cube&)");
         signatures.insert("arma::field<arma::mat>(*ZIRF)(const arma::field<arma::mat>&,const arma::mat&)");
@@ -893,6 +933,7 @@ RcppExport SEXP _bsvarSIGNs_RcppExport_registerCCallable() {
     R_RegisterCCallable("bsvarSIGNs", "_bsvarSIGNs_ir1_cpp", (DL_FUNC)_bsvarSIGNs_ir1_cpp_try);
     R_RegisterCCallable("bsvarSIGNs", "_bsvarSIGNs_bsvarSIGNs_ir", (DL_FUNC)_bsvarSIGNs_bsvarSIGNs_ir_try);
     R_RegisterCCallable("bsvarSIGNs", "_bsvarSIGNs_bsvar_sign_cpp", (DL_FUNC)_bsvarSIGNs_bsvar_sign_cpp_try);
+    R_RegisterCCallable("bsvarSIGNs", "_bsvarSIGNs_forecast_bsvarSIGNs", (DL_FUNC)_bsvarSIGNs_forecast_bsvarSIGNs_try);
     R_RegisterCCallable("bsvarSIGNs", "_bsvarSIGNs_match_sign_narrative", (DL_FUNC)_bsvarSIGNs_match_sign_narrative_try);
     R_RegisterCCallable("bsvarSIGNs", "_bsvarSIGNs_weight_narrative", (DL_FUNC)_bsvarSIGNs_weight_narrative_try);
     R_RegisterCCallable("bsvarSIGNs", "_bsvarSIGNs_ZIRF", (DL_FUNC)_bsvarSIGNs_ZIRF_try);
@@ -916,6 +957,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_bsvarSIGNs_ir1_cpp", (DL_FUNC) &_bsvarSIGNs_ir1_cpp, 4},
     {"_bsvarSIGNs_bsvarSIGNs_ir", (DL_FUNC) &_bsvarSIGNs_bsvarSIGNs_ir, 4},
     {"_bsvarSIGNs_bsvar_sign_cpp", (DL_FUNC) &_bsvarSIGNs_bsvar_sign_cpp, 14},
+    {"_bsvarSIGNs_forecast_bsvarSIGNs", (DL_FUNC) &_bsvarSIGNs_forecast_bsvarSIGNs, 6},
     {"_bsvarSIGNs_match_sign_narrative", (DL_FUNC) &_bsvarSIGNs_match_sign_narrative, 3},
     {"_bsvarSIGNs_weight_narrative", (DL_FUNC) &_bsvarSIGNs_weight_narrative, 3},
     {"_bsvarSIGNs_ZIRF", (DL_FUNC) &_bsvarSIGNs_ZIRF, 2},