Merge pull request #9 from bsvars/conditional-forecast

Conditional forecast

R/forecast.R

@@ -11,6 +11,11 @@
 #' @param horizon a positive integer, specifying the forecasting horizon.
 #' @param exogenous_forecast not used here ATM; included for compatibility with 
 #' generic \code{forecast}.
+#' @param conditional_forecast a list of length \code{C} containing 
+#' \code{horizon x N} matrices with forecasted values for selected variables. 
+#' These matrices 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{PanelForecasts} containing the
 #' draws from the predictive density and data. The output list includes element:
@@ -40,18 +45,71 @@
 #'   estimate(S = 20) |> 
 #'   forecast(horizon = 2) -> predictive
 #' 
+#' # conditional forecasting 6 years ahead conditioning on 
+#' #  provided future values for the Gross Domestic Product 
+#' #  growth rate
+#' ############################################################
+#' #' data(ilo_conditional_forecast)                        # load the conditional forecasts of dgdp
+#' predictive    = forecast(posterior, 6, conditional_forecast = ilo_conditional_forecast)
+#' 
+#' # workflow with the pipe |>
+#' ############################################################
+#' set.seed(123)
+#' ilo_cubic_panel |>
+#'   specify_bvarPANEL$new() |>
+#'   estimate(S = 10) |> 
+#'   estimate(S = 20) |> 
+#'   forecast(
+#'     horizon = 6, 
+#'     conditional_forecast = ilo_conditional_forecast
+#'   ) -> predictive
+#' 
 #' @export
-forecast.PosteriorBVARPANEL = function(posterior, horizon = 1, exogenous_forecast = NULL) {
-
+forecast.PosteriorBVARPANEL = function(
+    posterior, 
+    horizon = 1, 
+    exogenous_forecast = NULL,
+    conditional_forecast = NULL
+) {
 
   posterior_A_c_cpp     = posterior$posterior$A_c_cpp
   posterior_Sigma_c_cpp = posterior$posterior$Sigma_c_cpp
   X_c             = posterior$last_draw$data_matrices$X
   Y_c             = posterior$last_draw$data_matrices$Y
+  N               = dim(Y_c[[1]])[2]
 
-  fore            = .Call(`_bvarPANELs_forecast_bvarPANEL`, posterior_A_c_cpp, posterior_Sigma_c_cpp, X_c, horizon)
+  do_conditional_forecasting = !is.null(conditional_forecast)
+
+  if (!do_conditional_forecasting) {
+    # perform forecasting
+    fore          = .Call(`_bvarPANELs_forecast_bvarPANEL`, 
+                          posterior_A_c_cpp, posterior_Sigma_c_cpp, X_c, horizon)
+  } else {
+    stopifnot("Argument conditional_forecast must be a list with the same countries 
+              as in the provided data." 
+              = is.list(conditional_forecast) & length(conditional_forecast) == length(Y_c)
+    )
+    stopifnot("Argument conditional_forecast must be a list with the same countries 
+              as in the provided data."
+              = all(names(Y_c) == names(conditional_forecast))
+    )
+    stopifnot("Argument conditional_forecast must be a list with matrices with numeric values."
+              = all(sapply(conditional_forecast, function(x) is.matrix(x) & is.numeric(x)))
+    )
+    stopifnot("All the matrices provided in argument conditional_forecast must have 
+              the same number of rows equal to the value of argument horizon."
+              = unique(sapply(conditional_forecast, function(x) nrow(x) )) == horizon
+    )
+    stopifnot("All the matrices provided in argument conditional_forecast must have 
+              the same number of columns equal to the number of columns in the used data."
+              = unique(sapply(conditional_forecast, function(x) ncol(x) )) == N
+    )
+
+    # perform conditional forecasting
+    fore          = .Call(`_bvarPANELs_forecast_conditional_bvarPANEL`, 
+                          posterior_A_c_cpp, posterior_Sigma_c_cpp, X_c, conditional_forecast, horizon)
+  }
 
-  N               = dim(Y_c[[1]])[2]
   S               = dim(posterior_A_c_cpp)[1]
   C               = length(Y_c)
   forecasts      = array(NA, c(horizon, N, S, C))

R/ilo_conditional_forecast.R

@@ -0,0 +1,29 @@
+
+#' @title Data containing conditional projections for growth rate of GDP (dgdp)
+#' for 189 United Nations countries from 2024 to 2029
+#'
+#' @description For each of the countries a time series of 6 observations on 
+#' GDP growth rates (sgdp) #' formatted so they is provided to generate 
+#' conditional forecasts of labour market outcomes given the provided projected 
+#' paths of output. Last data update was implemented on 2024-05-11.
+#'
+#' @usage data(ilo_conditional_forecast)
+#' 
+#' @format A list of 189 \code{ts} objects with time series of 6 observations 
+#' on 4 variables:
+#' \describe{
+#'   \item{UR}{unemployment rate - contains missing values}
+#'   \item{EPR}{annual employment rate - contains missing values}
+#'   \item{LFPR}{annual labour force participation rate - contains missing values}
+#'   \item{dgdp}{annual growth rate of gross domestic product - contains projected 
+#'   values}
+#' }
+#' 
+#' @source
+#' International Labour Organization. (2020). ILO modelled estimates database,
+#' ILOSTAT [database]. Available from \url{https://ilostat.ilo.org/data/}.
+#' 
+#' @examples 
+#' data(ilo_conditional_forecast)   # upload the data
+#' 
+"ilo_conditional_forecast"

inst/include/bvarPANELs_RcppExports.h

@@ -46,6 +46,48 @@ namespace bvarPANELs {
         return Rcpp::as<Rcpp::List >(rcpp_result_gen);
     }
 
+    inline arma::vec mvnrnd_cond(arma::vec x, arma::vec mu, arma::mat Sigma) {
+        typedef SEXP(*Ptr_mvnrnd_cond)(SEXP,SEXP,SEXP);
+        static Ptr_mvnrnd_cond p_mvnrnd_cond = NULL;
+        if (p_mvnrnd_cond == NULL) {
+            validateSignature("arma::vec(*mvnrnd_cond)(arma::vec,arma::vec,arma::mat)");
+            p_mvnrnd_cond = (Ptr_mvnrnd_cond)R_GetCCallable("bvarPANELs", "_bvarPANELs_mvnrnd_cond");
+        }
+        RObject rcpp_result_gen;
+        {
+            RNGScope RCPP_rngScope_gen;
+            rcpp_result_gen = p_mvnrnd_cond(Shield<SEXP>(Rcpp::wrap(x)), Shield<SEXP>(Rcpp::wrap(mu)), Shield<SEXP>(Rcpp::wrap(Sigma)));
+        }
+        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::vec >(rcpp_result_gen);
+    }
+
+    inline Rcpp::List forecast_conditional_bvarPANEL(arma::field<arma::cube>& posterior_A_c_cpp, arma::field<arma::cube>& posterior_Sigma_c_cpp, Rcpp::List& X_c, Rcpp::List& cond_forecasts, const int horizon) {
+        typedef SEXP(*Ptr_forecast_conditional_bvarPANEL)(SEXP,SEXP,SEXP,SEXP,SEXP);
+        static Ptr_forecast_conditional_bvarPANEL p_forecast_conditional_bvarPANEL = NULL;
+        if (p_forecast_conditional_bvarPANEL == NULL) {
+            validateSignature("Rcpp::List(*forecast_conditional_bvarPANEL)(arma::field<arma::cube>&,arma::field<arma::cube>&,Rcpp::List&,Rcpp::List&,const int)");
+            p_forecast_conditional_bvarPANEL = (Ptr_forecast_conditional_bvarPANEL)R_GetCCallable("bvarPANELs", "_bvarPANELs_forecast_conditional_bvarPANEL");
+        }
+        RObject rcpp_result_gen;
+        {
+            RNGScope RCPP_rngScope_gen;
+            rcpp_result_gen = p_forecast_conditional_bvarPANEL(Shield<SEXP>(Rcpp::wrap(posterior_A_c_cpp)), Shield<SEXP>(Rcpp::wrap(posterior_Sigma_c_cpp)), Shield<SEXP>(Rcpp::wrap(X_c)), Shield<SEXP>(Rcpp::wrap(cond_forecasts)), 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<Rcpp::List >(rcpp_result_gen);
+    }
+
 }
 
 #endif // RCPP_bvarPANELs_RCPPEXPORTS_H_GEN_

inst/tinytest/test_forecast.R

@@ -38,3 +38,46 @@ expect_error(
   forecast(run_no1, horizon = 1.5),
   info = "forecast: specify horizon as integer."
 )
+
+
+# conditional forecasting
+data(ilo_conditional_forecast)
+
+set.seed(1)
+suppressMessages(
+  specification_no1 <- specify_bvarPANEL$new(ilo_cubic_panel)
+)
+run_no1             <- estimate(specification_no1, 3, 1, show_progress = FALSE)
+ff                  <- forecast(run_no1, 6, conditional_forecast = ilo_conditional_forecast)
+
+set.seed(1)
+suppressMessages(
+  ff2              <- ilo_cubic_panel |>
+    specify_bvarPANEL$new() |>
+    estimate(S = 3, thin = 1, show_progress = FALSE) |>
+    forecast(horizon = 6, conditional_forecast = ilo_conditional_forecast)
+)
+
+
+expect_identical(
+  ff$forecasts[1,1,1,1], ff2$forecasts[1,1,1,1],
+  info = "conditional forecast: forecast identical for normal and pipe workflow."
+)
+
+expect_true(
+  is.numeric(ff$forecasts) & is.array(ff$forecasts),
+  info = "conditional forecast: returns numeric array."
+)
+
+expect_error(
+  forecast(run_no1, horizon = 4, conditional_forecast = ilo_conditional_forecast),
+  pattern = "horizon",
+  info = "conditional forecast: provided forecasts different from horizon."
+)
+
+expect_error(
+  forecast(run_no1, horizon = 6, conditional_forecast = ilo_conditional_forecast[-1]),
+  info = "conditional forecast: uneven number of countries in forecasts and data."
+)
+
+

inst/varia/ilo_conditional_forecast.R

@@ -0,0 +1,26 @@
+
+# Create ILO conditional forecasts from provided files
+library(dplyr)
+
+# this file contains a cubic panel dynamic dataset
+all_cv    <- read.csv("inst/varia/ilo_cubic_panel.csv")
+colnames(all_cv) = c("year", "iso3code", "country", "UR", "EPR", "LFPR", "dgdp")
+
+# all variables all countries
+data_cv <- all_cv %>% 
+  filter(year >= 2024) %>% 
+  select(year, iso3code, UR, EPR, LFPR, dgdp)
+
+# Create a list with the country data
+countries = unique(data_cv$iso3code)
+countries = countries[order(countries)]
+ilo_conditional_forecast = list()
+for (i in 1:length(countries)) {
+  ilo_conditional_forecast[[i]] <- data_cv %>% 
+    filter(iso3code == countries[i]) %>% 
+    select(UR, EPR, LFPR, dgdp) %>% 
+    ts(start = 2024, frequency = 1)
+  names(ilo_conditional_forecast)[i] <- countries[i]
+}
+
+save(ilo_conditional_forecast, file = "data/ilo_conditional_forecast.rda")

inst/varia/ilo_cubic_panel.R

@@ -1,6 +1,5 @@
 
-# A first look at the data
-
+# Create ILO dataset from provided files
 library(dplyr)
 
 # this file contains a cubic panel dynamic dataset