Merge pull request #13 from bsvars/12-exogenous-variables

exogenous variables

NEWS.md

@@ -4,3 +4,7 @@
 2. The package includes the Gibbs sampler for a Bayesian Hierarchical Panel Vector Autoregression model. [#2](https://github.com/bsvars/bvarPANELs/issues/2)
 3. The package includes labour market dynamic panel data. [#3](https://github.com/bsvars/bvarPANELs/issues/3)
 4. The package includes the forecast method. [#5](https://github.com/bsvars/bvarPANELs/issues/5)
+5. The package includes possibility of conditional forecasting given projected paths of some of the variables.
+[#8](https://github.com/bsvars/bvarPANELs/issues/8)
+6. The package includes estimation and forecasting with exogenous variables. [#12](https://github.com/bsvars/bvarPANELs/issues/12)
+

R/bvarPANELs-package.R

@@ -40,19 +40,41 @@
 #' @keywords package models ts
 #' #' @examples
 #' @examples
+#' # Basic estimation and forecasting example
+#' ############################################################
 #' data(ilo_cubic_panel)                                   # load the data
 #' set.seed(123)
 #' specification = specify_bvarPANEL$new(ilo_cubic_panel)  # specify the model
-#' burn_in       = estimate(specification, 20)             # run the burn-in
-#' posterior     = estimate(burn_in, 20)                   # estimate the model
-#' predictive    = forecast(posterior, 2)                  # forecast 2 years ahead
+#' burn_in       = estimate(specification, 10)             # run the burn-in
+#' posterior     = estimate(burn_in, 10)                   # estimate the model
+#' predictive    = forecast(posterior, 2)                  # forecast the future       
 #' 
 #' # workflow with the pipe |>
-#' ############################################################
 #' set.seed(123)
 #' ilo_cubic_panel |>
 #'   specify_bvarPANEL$new() |>
 #'   estimate(S = 20) |> 
 #'   estimate(S = 20) |> 
 #'   forecast(horizon = 2) -> predictive
+#'   
+#' # Full estimation and forecasting example with 
+#' #   exogenous variables and conditional forecasts
+#' ############################################################
+#' data(ilo_cubic_panel)                                   # load the data
+#' data(ilo_exogenous_variables)                           # load the exogenous variables
+#' data(ilo_exogenous_forecasts)                           # load the exogenous forecasts
+#' data(ilo_conditional_forecasts)                         # load the conditional forecasts
+#' set.seed(123)
+#' specification = specify_bvarPANEL$new(
+#'                   ilo_cubic_panel,
+#'                   exogenous = ilo_exogenous_variables
+#'                 )
+#' burn_in       = estimate(specification, 10)             # run the burn-in
+#' posterior     = estimate(burn_in, 10)                   # estimate the model
+#' predictive    = forecast(
+#'                   posterior, 
+#'                   horizon = 6,
+#'                   exogenous_forecast = ilo_exogenous_forecasts,
+#'                   conditional_forecast = ilo_conditional_forecasts
+#'                 )
 "_PACKAGE"

R/estimate.bvarPANEL.R

@@ -38,16 +38,18 @@
 #' 
 #' @examples
 #' data(ilo_cubic_panel)                                   # load the data
+#' data(ilo_exogenous_variables)                           # load the exogenous variables
 #' set.seed(123)
-#' specification = specify_bvarPANEL$new(ilo_cubic_panel)  # specify the model
+#' # specify the model
+#' specification = specify_bvarPANEL$new(ilo_cubic_panel, exogenous = ilo_exogenous_variables)
 #' burn_in       = estimate(specification, 10)             # run the burn-in
 #' posterior     = estimate(burn_in, 10)                   # estimate the model
 #' 
 #' @export
 estimate.BVARPANEL <- function(
     specification, 
     S, 
-    thin = 10L, 
+    thin = 1L, 
     show_progress = TRUE
 ) {
 

R/forecast.R

@@ -30,11 +30,15 @@
 #' 
 #' @examples
 #' data(ilo_cubic_panel)                                   # load the data
+#' data(ilo_exogenous_variables)                           # load the exogenous variables
+#' data(ilo_exogenous_forecasts)                           # load the exogenous forecast
 #' set.seed(123)
-#' specification = specify_bvarPANEL$new(ilo_cubic_panel)  # specify the model
+#' # specify the model
+#' specification = specify_bvarPANEL$new(ilo_cubic_panel, exogenous = ilo_exogenous_variables)
 #' burn_in       = estimate(specification, 10)             # run the burn-in
 #' posterior     = estimate(burn_in, 10)                   # estimate the model
-#' predictive    = forecast(posterior, 2)                  # forecast 2 years ahead
+#' # forecast 6 years ahead
+#' predictive    = forecast(posterior, 6, exogenous_forecast = ilo_exogenous_forecasts)
 #' 
 #' # workflow with the pipe |>
 #' ############################################################
@@ -49,8 +53,12 @@
 #' #  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)
+#' data(ilo_conditional_forecasts)                        # load the conditional forecasts of dgdp
+#' specification = specify_bvarPANEL$new(ilo_cubic_panel)    # specify the model
+#' burn_in       = estimate(specification, 10)               # run the burn-in
+#' posterior     = estimate(burn_in, 10)                     # estimate the model
+#' # forecast 6 years ahead
+#' predictive    = forecast(posterior, 6, conditional_forecast = ilo_conditional_forecasts)
 #' 
 #' # workflow with the pipe |>
 #' ############################################################
@@ -61,7 +69,7 @@
 #'   estimate(S = 20) |> 
 #'   forecast(
 #'     horizon = 6, 
-#'     conditional_forecast = ilo_conditional_forecast
+#'     conditional_forecast = ilo_conditional_forecasts
 #'   ) -> predictive
 #' 
 #' @export
@@ -77,13 +85,30 @@ forecast.PosteriorBVARPANEL = function(
   X_c             = posterior$last_draw$data_matrices$X
   Y_c             = posterior$last_draw$data_matrices$Y
   N               = dim(Y_c[[1]])[2]
+  K               = dim(X_c[[1]])[2]
+  C               = length(Y_c)
 
-  do_conditional_forecasting = !is.null(conditional_forecast)
+  d               = K - N * posterior$last_draw$p - 1
+  if (d == 0 ) {
+    # this will not be used for forecasting, but needs to be provided
+    exogenous_forecast = list()
+    for (c in 1:C) exogenous_forecast[[c]] = 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." 
+              = unique(simplify2array(lapply(exogenous_forecast, function(x){ncol(x)}))) == d)
+    stopifnot("Provide exogenous_forecast for all forecast periods specified by argument horizon." 
+              = unique(simplify2array(lapply(exogenous_forecast, function(x){nrow(x)}))) == horizon)
+    stopifnot("Argument exogenous has to be a matrix." 
+              = all(simplify2array(lapply(exogenous_forecast, function(x){is.matrix(x) & is.numeric(x)}))))
+    stopifnot("Argument exogenous cannot include missing values." 
+              = unique(simplify2array(lapply(exogenous_forecast, function(x){any(is.na(x))}))) == FALSE)
+  }
 
-  if (!do_conditional_forecasting) {
-    # perform forecasting
-    fore          = .Call(`_bvarPANELs_forecast_bvarPANEL`, 
-                          posterior_A_c_cpp, posterior_Sigma_c_cpp, X_c, horizon)
+  if ( is.null(conditional_forecast) ) {
+    conditional_forecast = list()
+    for (c in 1:C) conditional_forecast[[c]] = matrix(NA, horizon, N)
   } else {
     stopifnot("Argument conditional_forecast must be a list with the same countries 
               as in the provided data." 
@@ -104,12 +129,18 @@ forecast.PosteriorBVARPANEL = function(
               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)
   }
 
+  # perform forecasting
+  fore          = .Call(`_bvarPANELs_forecast_bvarPANEL`, 
+                        posterior_A_c_cpp, 
+                        posterior_Sigma_c_cpp, 
+                        X_c, 
+                        conditional_forecast, 
+                        exogenous_forecast, 
+                        horizon
+                       )
+
   S               = dim(posterior_A_c_cpp)[1]
   C               = length(Y_c)
   forecasts      = array(NA, c(horizon, N, S, C))

R/ilo_conditional_forecast.R → R/ilo_conditional_forecasts.R

@@ -5,9 +5,10 @@
 #' @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.
+#' paths of output. 
+#' Last data update was implemented on 2024-05-11.
 #'
-#' @usage data(ilo_conditional_forecast)
+#' @usage data(ilo_conditional_forecasts)
 #' 
 #' @format A list of 189 \code{ts} objects with time series of 6 observations 
 #' on 4 variables:
@@ -24,6 +25,6 @@
 #' ILOSTAT [database]. Available from \url{https://ilostat.ilo.org/data/}.
 #' 
 #' @examples 
-#' data(ilo_conditional_forecast)   # upload the data
+#' data(ilo_conditional_forecasts)   # upload the data
 #' 
-"ilo_conditional_forecast"
+"ilo_conditional_forecasts"

R/ilo_exogenous_forecasts.R

@@ -0,0 +1,25 @@
+
+#' @title Data containing future observations for 189 United Nations countries 
+#' from 2024 to 2029 to be used to forecast with models with 
+#' \code{ilo_exogenous_variables}
+#'
+#' @description For each of the countries a time series of 6 observations on 
+#' On the dummies is provided. These future values are all equal to zero. They 
+#' provide benchmark for the objects to be used when \code{exogenous_variables} 
+#' are used. 
+#' Last data update was implemented on 2024-05-11.
+#'
+#' @usage data(ilo_exogenous_forecasts)
+#' 
+#' @format A list of 189 \code{ts} objects with time series of 6 observations 
+#' on 3 variables:
+#' \describe{
+#'   \item{2008}{the aftermath of the Global Financial Crisis}
+#'   \item{2020}{the COVID pandemic}
+#'   \item{2021}{the aftermath of the COVID pandemic}
+#' }
+#' 
+#' @examples 
+#' data(ilo_exogenous_forecasts)   # upload the data
+#' 
+"ilo_exogenous_forecasts"

R/ilo_exogenous_variables.R

@@ -0,0 +1,23 @@
+
+#' @title A 3-variable annual system for of dummy observations for 2008, 2020, 
+#' and 2021 to be used in the estimation of the Panel VAR model for 
+#' 189 United Nations countries from 1991 to 2023
+#'
+#' @description For each of the countries a time series of 33 observations on 3 
+#' dummy variables for the years 2008, 2020, and 2021 is provided. 
+#' Last data update was implemented on 2024-06-29.
+#'
+#' @usage data(ilo_exogenous_variables)
+#' 
+#' @format A list of 189 \code{ts} objects with time series of 33 observations 
+#' on 3 variables:
+#' \describe{
+#'   \item{2008}{the aftermath of the Global Financial Crisis}
+#'   \item{2020}{the COVID pandemic}
+#'   \item{2021}{the aftermath of the COVID pandemic}
+#' }
+#' 
+#' @examples 
+#' data(ilo_exogenous_variables)   # upload the data
+#' 
+"ilo_exogenous_variables"

R/specify_bvarpanel.R

@@ -76,20 +76,22 @@ specify_prior_bvarPANEL = R6::R6Class(
     #' @param C a positive integer - the number of countries in the data.
     #' @param N a positive integer - the number of dependent variables in the model.
     #' @param p a positive integer - the autoregressive lag order of the SVAR model.
+    #' @param d a positive integer - the number of \code{exogenous} variables in the model.
     #' @return A new prior specification PriorBVARPANEL.
     #' @examples 
     #' # a prior for 2-country, 3-variable example with one lag and stationary data
     #' prior = specify_prior_bvarPANEL$new(C = 2, N = 3, p = 1)
     #' prior$M
     #' 
-    initialize = function(C, N, p){
+    initialize = function(C, N, p, d = 0){
       stopifnot("Argument C must be a positive integer number." = C > 0 & C %% 1 == 0)
       stopifnot("Argument N must be a positive integer number." = N > 0 & N %% 1 == 0)
       stopifnot("Argument p must be a positive integer number." = p > 0 & p %% 1 == 0)
+      stopifnot("Argument d must be a non-negative integer number." = d >= 0 & d %% 1 == 0)
 
-      K                 = N * p + 1
+      K                 = N * p + 1 + d
       self$M            = t(cbind(diag(N), matrix(0, N, K - N)))
-      self$W            = diag(c(kronecker((1:p)^2, rep(1, N) ), rep(10, 1)))
+      self$W            = diag(c(kronecker((1:p)^2, rep(1, N) ), rep(10, 1 + d)))
       self$S_inv        = diag(N)
       self$S_Sigma_inv  = diag(N)
       self$eta          = N + 1
@@ -194,17 +196,19 @@ specify_starting_values_bvarPANEL = R6::R6Class(
     #' @param C a positive integer - the number of countries in the data.
     #' @param N a positive integer - the number of dependent variables in the model.
     #' @param p a positive integer - the autoregressive lag order of the SVAR model.
+    #' @param d a positive integer - the number of \code{exogenous} variables in the model.
     #' @return Starting values StartingValuesBVARPANEL
     #' @examples 
     #' # starting values for Bayesian Panel VAR 2-country model with 4 lags for a 3-variable system.
     #' sv = specify_starting_values_bvarPANEL$new(C = 2, N = 3, p = 4)
     #' 
-    initialize = function(C, N, p){
+    initialize = function(C, N, p, d = 0){
       stopifnot("Argument C must be a positive integer number." = C > 0 & C %% 1 == 0)
       stopifnot("Argument N must be a positive integer number." = N > 0 & N %% 1 == 0)
       stopifnot("Argument p must be a positive integer number." = p > 0 & p %% 1 == 0)
+      stopifnot("Argument d must be a non-negative integer number." = d >= 0 & d %% 1 == 0)
 
-      K               = N * p + 1
+      K               = N * p + 1 + d
       self$A_c        = array(stats::rnorm(C * K * N, sd = 0.001), c(K, N, C))
       self$Sigma_c    = stats::rWishart(C, N + 1, diag(N))
       self$A          = matrix(stats::rnorm(K * N, sd = 0.001), K, N) + diag(K)[,1:N]
@@ -302,8 +306,11 @@ specify_panel_data_matrices = R6::R6Class(
     #' @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 
+    #' constant term.
     #' @return New data matrices DataMatricesBVARPANEL
-    initialize = function(data, p = 1L) {
+    initialize = function(data, p = 1L, exogenous = NULL) {
       if (missing(data)) {
         stop("Argument data has to be specified")
       } else {
@@ -314,6 +321,23 @@ specify_panel_data_matrices = R6::R6Class(
       stopifnot("Argument p must be a positive integer number." = p > 0 & p %% 1 == 0)
 
       C             = length(data)
+      if (is.null(exogenous)) {
+        d = 0
+      } else {
+        stopifnot("Argument exogenous has to be a list of matrices." = is.list(exogenous) & all(simplify2array(lapply(exogenous, function(x){is.matrix(x) & is.numeric(x)}))))
+        stopifnot("Argument exogenous has to contain matrices with the same number of rows as argument data." = unique(simplify2array(lapply(exogenous, function(x){ncol(x)}))) >= 1 & unique(simplify2array(lapply(data, function(x){nrow(x)}))) == unique(simplify2array(lapply(exogenous, function(x){nrow(x)}))))
+        stopifnot("Argument exogenous cannot include missing values." = unique(simplify2array(lapply(exogenous, function(x){any(is.na(x))}))) == FALSE )
+        d = ncol(exogenous[[1]])
+        Td = nrow(exogenous[[1]])
+        test_exogenous = 0
+        for (c in 1:C) {
+          for (i in 1:d) {
+            test_exogenous = test_exogenous + prod(exogenous[[c]][,i] - mean(exogenous[[c]][,i]) == rep(0,Td))
+          }
+        }
+        stopifnot("Argument exogenous cannot include a constant term." = test_exogenous == 0 )
+      }
+
       for (c in 1:C) {
         TT            = nrow(data[[c]])
         T_c           = TT - p
@@ -324,6 +348,9 @@ specify_panel_data_matrices = R6::R6Class(
           X           = cbind(X, data[[c]][(p + 1):TT - i,])
         }
         X             = cbind(X, rep(1, T_c))
+        if (!is.null(exogenous)) {
+          X           = cbind(X, exogenous[[c]][(p + 1):TT,])
+        }
         self$X[[c]]   = X
       } # END c loop
       names(self$Y)   = names(self$X) = names(data)
@@ -393,21 +420,27 @@ specify_bvarPANEL = R6::R6Class(
     #' @param data a list with \code{C} elements of \code{(T_c+p)xN} matrices 
     #' with time series data.
     #' @param p a positive integer providing model's autoregressive lag order.
+    #' @param exogenous a \code{(T+p)xd} matrix of exogenous variables. 
     #' @return A new complete specification for the Bayesian Panel VAR model BVARPANEL.
     initialize = function(
     data,
-    p = 1L
+    p = 1L,
+    exogenous = NULL
     ) {
       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)
       N         = unique(simplify2array(lapply(data, ncol)))
+      d             = 0
+      if (!is.null(exogenous)) {
+        d           = ncol(exogenous[[1]])
+      }
 
-      self$data_matrices   = specify_panel_data_matrices$new(data, self$p)
-      self$prior           = specify_prior_bvarPANEL$new(C, N, self$p)
-      self$starting_values = specify_starting_values_bvarPANEL$new(C, N, self$p)
+      self$data_matrices   = specify_panel_data_matrices$new(data, self$p, exogenous)
+      self$prior           = specify_prior_bvarPANEL$new(C, N, self$p, d)
+      self$starting_values = specify_starting_values_bvarPANEL$new(C, N, self$p, d)
       self$adaptiveMH      = c(0.6, 0.4, 10, 0.1)
     }, # END initialize