Skip to content

Commit

Permalink
Merge pull request #28 from bsvars/23-update-the-forecast-method
Browse files Browse the repository at this point in the history
23 update the forecast method
  • Loading branch information
donotdespair authored Jul 15, 2024
2 parents c78c165 + 713025d commit 93075c0
Show file tree
Hide file tree
Showing 8 changed files with 409 additions and 0 deletions.
1 change: 1 addition & 0 deletions NAMESPACE
Original file line number Diff line number Diff line change
@@ -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)
Expand Down
132 changes: 132 additions & 0 deletions R/forecast.R
Original file line number Diff line number Diff line change
@@ -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
21 changes: 21 additions & 0 deletions inst/include/bsvarSIGNs_RcppExports.h
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
41 changes: 41 additions & 0 deletions inst/tinytest/test_forecast.R
Original file line number Diff line number Diff line change
@@ -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."
)

84 changes: 84 additions & 0 deletions man/forecast.PosteriorBSVARSIGN.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

42 changes: 42 additions & 0 deletions src/RcppExports.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand Down Expand Up @@ -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&)");
Expand All @@ -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);
Expand All @@ -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},
Expand Down
Loading

0 comments on commit 93075c0

Please sign in to comment.