Merge pull request #24 from bsvars/21-truncated-normal-forecasts

21 truncated normal forecasts

R/bvarPANELs-package.R

@@ -155,9 +155,16 @@
 #' unknown values to be predicted, \eqn{\mathbf{Y}_f}, given data, \eqn{\mathbf{Y}}
 #' closely following Karlsson (2013):
 #' \deqn{p(\mathbf{Y}_f | \mathbf{Y})}
-#' The package offers the possibility of conditional predictions given provided 
-#' future projections for some of the variables. The forecasting is performed 
-#' using function \code{\link{forecast.PosteriorBVARPANEL}}.
+#' The package offers the possibility of:
+#' \describe{
+#'   \item{forecasting for models with exogenous variables}{given the provided 
+#'   future values of the exogenous variables.}
+#'   \item{conditional predictions}{given provided future projections for some 
+#'   of the variables.}
+#'   \item{trucated forecasts}{for variables that represents rates from the 
+#'   interval \eqn{[0,100]}.}
+#' }
+#' The forecasting is performed using function \code{\link{forecast.PosteriorBVARPANEL}}.
 #' 
 #' @seealso \code{\link{specify_bvarPANEL}}, \code{\link{estimate.BVARPANEL}}, 
 #' \code{\link{forecast.PosteriorBVARPANEL}}
@@ -206,7 +213,7 @@
 #' plot(predictive, which_c = "POL")
 #' 
 #' # Full estimation and forecasting example with 
-#' #   exogenous variables and conditional forecasts
+#' #   exogenous variables, conditional forecasts, and truncation for rates
 #' ############################################################
 #' data(ilo_dynamic_panel)                                 # load the data
 #' data(ilo_exogenous_variables)                           # load the exogenous variables
@@ -215,7 +222,8 @@
 #' set.seed(123)
 #' specification = specify_bvarPANEL$new(
 #'                   ilo_dynamic_panel,
-#'                   exogenous = ilo_exogenous_variables
+#'                   exogenous = ilo_exogenous_variables,
+#'                   type = c("real", rep("rates", 3))
 #'                 )
 #' burn_in       = estimate(specification, S = 10)         # run the burn-in; use say S = 10000
 #' posterior     = estimate(burn_in, S = 10)               # estimate the model; use say S = 10000

R/forecast.R

@@ -52,6 +52,13 @@
 #' of these exogenous variables in the argument \code{exogenous_forecast} of the
 #' \code{\link{forecast.PosteriorBVARPANEL}} function.
 #' 
+#' \strong{Truncated forecasts for variables of type 'rates'.}
+#' The package provides the option to truncate the forecasts for variables of 
+#' for which the corresponding element of argument \code{type} of the function 
+#' \code{specify_bvarPANEL$new()} is set to \code{"rates"}. The one-period-ahead
+#' predictive normal density for such variables is truncated to values from 
+#' interval  \eqn{[0,100]}. 
+#' 
 #' @method forecast PosteriorBVARPANEL
 #' 
 #' @param posterior posterior estimation outcome - an object of class 
@@ -81,7 +88,7 @@
 #' \emph{Review of Economics and Statistics}, \bold{81}(4), 639-651,
 #' \doi{10.1162/003465399558508}.
 #'
-#' @seealso \code{\link{estimate.PosteriorBVARPANEL}}, 
+#' @seealso \code{\link{specify_bvarPANEL}}, \code{\link{estimate.PosteriorBVARPANEL}}, 
 #' \code{\link{summary.ForecastsPANEL}}, \code{\link{plot.ForecastsPANEL}}
 #'
 #' @author Tomasz Woźniak \email{wozniak.tom@pm.me}
@@ -111,10 +118,13 @@
 #' 
 #' # conditional forecasting 6 years ahead conditioning on 
 #' #  provided future values for the Gross Domestic Product 
-#' #  growth rate
+#' #  and truncated forecasts for the rates
 #' ############################################################
 #' data(ilo_conditional_forecasts)                        # load the conditional forecasts of dgdp
-#' specification = specify_bvarPANEL$new(ilo_dynamic_panel)   # specify the model
+#' specification = specify_bvarPANEL$new(
+#'                   ilo_dynamic_panel,
+#'                   type = c("real", rep("rates", 3))
+#'                 )   # specify the model
 #' burn_in       = estimate(specification, 10)            # run the burn-in; use say S = 10000
 #' posterior     = estimate(burn_in, 10)                  # estimate the model; use say S = 10000
 #' # forecast 6 years ahead
@@ -124,7 +134,7 @@
 #' ############################################################
 #' set.seed(123)
 #' ilo_dynamic_panel |>
-#'   specify_bvarPANEL$new() |>
+#'   specify_bvarPANEL$new(type = c("real", rep("rates", 3))) |>
 #'   estimate(S = 10) |> 
 #'   estimate(S = 20) |> 
 #'   forecast(
@@ -192,14 +202,26 @@ forecast.PosteriorBVARPANEL = function(
     )
   }
 
+  type      = posterior$last_draw$data_matrices$type
+  LB        = rep(-Inf, N)
+  UB        = rep(Inf, N)
+  rates_id  = which(type == "rates")
+  if (length(rates_id) > 0) {
+    LB[rates_id] = 0
+    UB[rates_id] = 100
+  }
+
   # perform forecasting
   fff           = .Call(`_bvarPANELs_forecast_bvarPANEL`, 
                         posterior_A_c_cpp, 
                         posterior_Sigma_c_cpp, 
                         X_c, 
                         conditional_forecast, 
                         exogenous_forecast, 
-                        horizon
+                        horizon,
+                        LB,
+                        UB,
+                        TRUE
                        )
 
   forecasts       = list()

R/specify_bvarpanel.R

@@ -305,16 +305,24 @@ specify_panel_data_matrices = R6::R6Class(
     #' regressors, \eqn{\mathbf{X}_c}. 
     X     = list(),
 
+    #' @field type an \code{N} character vector with elements set to "rates" or "real"
+    #' determining the truncation of the predictive density to \code{[0, 100]} and
+    #' \code{(-Inf, Inf)} (no truncation) for each of the variables.
+    type  = character(),
+
     #' @description
     #' Create new data matrices DataMatricesBVARPANEL
     #' @param data a list containing \code{(T_c+p)xN} matrices with country-specific
     #' time series data.
     #' @param p a positive integer providing model's autoregressive lag order.
     #' @param exogenous a list containing \code{(T_c+p)xd} matrices with 
-    #' country-specificof exogenous variables. This matrix should not include a 
+    #' country-specific of exogenous variables. This matrix should not include a 
     #' constant term.
+    #' @param type an \code{N} character vector with elements set to "rates" or "real"
+    #' determining the truncation of the predictive density to \code{[0, 100]} and
+    #' \code{(-Inf, Inf)} (no truncation) for each of the variables.
     #' @return New data matrices DataMatricesBVARPANEL
-    initialize = function(data, p = 1L, exogenous = NULL) {
+    initialize = function(data, p = 1L, exogenous = NULL, type = rep("real", ncol(data[[1]]))) {
       if (missing(data)) {
         stop("Argument data has to be specified")
       } else {
@@ -324,6 +332,11 @@ specify_panel_data_matrices = R6::R6Class(
       }
       stopifnot("Argument p must be a positive integer number." = p > 0 & p %% 1 == 0)
 
+      stopifnot("Argument type must include elements 'rates' or 'real'." 
+                = all(unique(type) %in% c("rates", "real")) )
+      stopifnot("Argument type must be of length corresponding to the numbers of variables."
+                = length(type) == ncol(data[[1]]) )
+
       C             = length(data)
       if (is.null(exogenous)) {
         d = 0
@@ -358,6 +371,7 @@ specify_panel_data_matrices = R6::R6Class(
         self$X[[c]]   = X
       } # END c loop
       names(self$Y)   = names(self$X) = names(data)
+      self$type       = type
     }, # END initialize
 
     #' @description
@@ -371,7 +385,8 @@ specify_panel_data_matrices = R6::R6Class(
     get_data_matrices = function() {
       list(
         Y = self$Y,
-        X = self$X
+        X = self$X,
+        type = self$type
       )
     } # END get_data_matrices
   ) # END public
@@ -428,15 +443,21 @@ specify_bvarPANEL = R6::R6Class(
     #' @param stationary an \code{N} logical vector - its element set to 
     #' \code{FALSE} sets the prior mean for the autoregressive parameters of the 
     #' \code{N}th equation to the white noise process, otherwise to random walk.
+    #' @param type an \code{N} character vector with elements set to "rates" or "real"
+    #' determining the truncation of the predictive density to \code{[0, 100]} and
+    #' \code{(-Inf, Inf)} (no truncation) for each of the variables.
     #' @return A new complete specification for the Bayesian Panel VAR model BVARPANEL.
     initialize = function(
     data,
     p = 1L,
     exogenous = NULL,
-    stationary = rep(FALSE, ncol(data[[1]]))
+    stationary = rep(FALSE, ncol(data[[1]])),
+    type = rep("real", ncol(data[[1]]))
     ) {
-      stopifnot("Argument data has to contain matrices with the same number of columns." = length(unique(simplify2array(lapply(data, ncol)))) == 1)
-      stopifnot("Argument p has to be a positive integer." = ((p %% 1) == 0 & p > 0))
+      stopifnot("Argument data has to contain matrices with the same number of columns." 
+                = length(unique(simplify2array(lapply(data, ncol)))) == 1)
+      stopifnot("Argument p has to be a positive integer." 
+                = ((p %% 1) == 0 & p > 0))
 
       self$p    = p
       C         = length(data)
@@ -445,8 +466,8 @@ specify_bvarPANEL = R6::R6Class(
       if (!is.null(exogenous)) {
         d           = ncol(exogenous[[1]])
       }
-      
-      self$data_matrices   = specify_panel_data_matrices$new(data, self$p, exogenous)
+
+      self$data_matrices   = specify_panel_data_matrices$new(data, self$p, exogenous, type)
       self$prior           = specify_prior_bvarPANEL$new(C, N, self$p, d, stationary)
       self$starting_values = specify_starting_values_bvarPANEL$new(C, N, self$p, d)
       self$adaptiveMH      = c(0.44, 0.6)

inst/include/bvarPANELs_RcppExports.h

@@ -25,17 +25,59 @@ namespace bvarPANELs {
         }
     }
 
-    inline Rcpp::List forecast_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, Rcpp::List& exog_forecasts, const int horizon) {
-        typedef SEXP(*Ptr_forecast_bvarPANEL)(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP);
+    inline arma::vec mvnrnd_truncated(arma::vec mu, arma::mat Sigma, arma::vec LB, arma::vec UB) {
+        typedef SEXP(*Ptr_mvnrnd_truncated)(SEXP,SEXP,SEXP,SEXP);
+        static Ptr_mvnrnd_truncated p_mvnrnd_truncated = NULL;
+        if (p_mvnrnd_truncated == NULL) {
+            validateSignature("arma::vec(*mvnrnd_truncated)(arma::vec,arma::mat,arma::vec,arma::vec)");
+            p_mvnrnd_truncated = (Ptr_mvnrnd_truncated)R_GetCCallable("bvarPANELs", "_bvarPANELs_mvnrnd_truncated");
+        }
+        RObject rcpp_result_gen;
+        {
+            RNGScope RCPP_rngScope_gen;
+            rcpp_result_gen = p_mvnrnd_truncated(Shield<SEXP>(Rcpp::wrap(mu)), Shield<SEXP>(Rcpp::wrap(Sigma)), Shield<SEXP>(Rcpp::wrap(LB)), Shield<SEXP>(Rcpp::wrap(UB)));
+        }
+        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 arma::vec mvnrnd_cond_truncated(arma::vec x, arma::vec mu, arma::mat Sigma, arma::vec LB, arma::vec UB) {
+        typedef SEXP(*Ptr_mvnrnd_cond_truncated)(SEXP,SEXP,SEXP,SEXP,SEXP);
+        static Ptr_mvnrnd_cond_truncated p_mvnrnd_cond_truncated = NULL;
+        if (p_mvnrnd_cond_truncated == NULL) {
+            validateSignature("arma::vec(*mvnrnd_cond_truncated)(arma::vec,arma::vec,arma::mat,arma::vec,arma::vec)");
+            p_mvnrnd_cond_truncated = (Ptr_mvnrnd_cond_truncated)R_GetCCallable("bvarPANELs", "_bvarPANELs_mvnrnd_cond_truncated");
+        }
+        RObject rcpp_result_gen;
+        {
+            RNGScope RCPP_rngScope_gen;
+            rcpp_result_gen = p_mvnrnd_cond_truncated(Shield<SEXP>(Rcpp::wrap(x)), Shield<SEXP>(Rcpp::wrap(mu)), Shield<SEXP>(Rcpp::wrap(Sigma)), Shield<SEXP>(Rcpp::wrap(LB)), Shield<SEXP>(Rcpp::wrap(UB)));
+        }
+        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_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, Rcpp::List& exog_forecasts, const int horizon, arma::vec LB, arma::vec UB, const bool show_progress) {
+        typedef SEXP(*Ptr_forecast_bvarPANEL)(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP);
         static Ptr_forecast_bvarPANEL p_forecast_bvarPANEL = NULL;
         if (p_forecast_bvarPANEL == NULL) {
-            validateSignature("Rcpp::List(*forecast_bvarPANEL)(arma::field<arma::cube>&,arma::field<arma::cube>&,Rcpp::List&,Rcpp::List&,Rcpp::List&,const int)");
+            validateSignature("Rcpp::List(*forecast_bvarPANEL)(arma::field<arma::cube>&,arma::field<arma::cube>&,Rcpp::List&,Rcpp::List&,Rcpp::List&,const int,arma::vec,arma::vec,const bool)");
             p_forecast_bvarPANEL = (Ptr_forecast_bvarPANEL)R_GetCCallable("bvarPANELs", "_bvarPANELs_forecast_bvarPANEL");
         }
         RObject rcpp_result_gen;
         {
             RNGScope RCPP_rngScope_gen;
-            rcpp_result_gen = p_forecast_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(exog_forecasts)), Shield<SEXP>(Rcpp::wrap(horizon)));
+            rcpp_result_gen = p_forecast_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(exog_forecasts)), Shield<SEXP>(Rcpp::wrap(horizon)), Shield<SEXP>(Rcpp::wrap(LB)), Shield<SEXP>(Rcpp::wrap(UB)), Shield<SEXP>(Rcpp::wrap(show_progress)));
         }
         if (rcpp_result_gen.inherits("interrupted-error"))
             throw Rcpp::internal::InterruptedException();

inst/tinytest/test_forecast.R

@@ -136,3 +136,32 @@ expect_error(
   info = "exogenous forecast: provided exogenous forecasts contain missing values."
 )
 
+
+
+# truncated forecasts
+
+expect_identical(
+  class(specify_bvarPANEL$new(ilo_dynamic_panel, type = c("real",rep("rates",3))))[1], 
+  "BVARPANEL",
+  info = "truncated forecast: good specification of argument type."
+)
+
+expect_error(
+  specify_bvarPANEL$new(ilo_dynamic_panel, type = rep("rates",3)),
+  pattern = "length",
+  info = "truncated forecast: wrong specification of argument type."
+)
+
+set.seed(1)
+suppressMessages(
+  specification_no1 <- specify_bvarPANEL$new(ilo_dynamic_panel, type = c("real",rep("rates",3)))
+)
+run_no1             <- estimate(specification_no1, 3, 1, show_progress = FALSE)
+suppressMessages(
+  ff                  <- forecast(run_no1, 6)
+)
+
+expect_true(
+  all(ff$POL$forecasts[2,,] >= 0),
+  info = "truncated forecast: unemployment rates forecasts are non-negative."
+)