Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
96,618 characters · 35 sections · 56 citation commands
Forecasting with Bayesian Panel Vector Autoregressions Using the R Package bpvars (Version 2.0)
Modelling dynamic panel data requires application-specific approaches due to the data's distinctive features, including temporal persistence and cross-sectional dependence, and often necessitates handling missing observations. Consequently, the literature is full of models attempting to strike a balance between the efficient extraction of data informational content and handling the large number of parameters that this task would require Canova2007. Bayesian inference offers a wide range of tools for addressing these challenges canova2013. Its successful applications in economics and finance provide evidence at a high level of generality and offer improved forecasting precision jarocinski_responses_2010.
This paper introduces the \pkg{bpvars} package by bpvars for \proglang{R} Rcore for Bayesian Forecasting with Panel Vector Autoregressions. It was developed for the International Labour Organization for forecasting labour market outcomes, including unemployment, employment, and labour force participation rates, for 189 countries. The model's formulation was inspired by existing implementations of Bayesian hierarchical modelling of global data for other United Nations agencies, which include forecasting of carbon dioxide emissions and temperature raftery2017, the world's population raftery_bayesian_2014, or fertility rates fosdick_regional_2014,liu_bayesian_2024. Additionally, Bayesian inference plays an essential role in global modelling of climate change\footnote{See sokolov_probabilistic_2009,director_probabilistic_2021}, demographics\footnote{See gerland2014world,raftery_joint_2014,raftery_probabilistic_2023,yu_probabilistic_2023}, and epidemiology\footnote{See godwin_bayesian_2017,irons_estimating_2021}.
Our new family of panel Vector Autoregressions (VARs) builds on the model proposed by jarocinski_responses_2010. The common feature of our models is a country-specific VAR model with autoregressive and error term covariance parameter matrices that features prior distributions centred around their global counterparts. In other words, under the prior mean, the country-specific data follow a VAR model with global parameters that are also estimated. With that respect, our framework relates to the literature on the so-called exchangeable priors proposed by LindleySmith1972 and random coefficient modelling rendon2013, extending them to all country-specific model parameters in a dynamic panel modelling context. Our models treat all variables as endogenous and contemporaneously related, capturing the persistence of individual variables, their dynamic interactions, and supplementing equations for individual countries with information from other countries' data through the global prior. This is complemented by a multiple-level hierarchical prior structure that allows the model to adjust to various types of data without the necessity of making arbitrary choices regarding the hyper-parameters of the prior distributions. These features lead to a simple model formulation that benefits from the flexibility of hierarchical priors, thereby improving its forecasting performance.
The package includes four panel VAR models, each featuring several variations. The benchmark model exhibits all the features mentioned above. It can be estimated in a version proposed by jarocinski_responses_2010, which introduces a minimal prior structure and the exchangeable prior for the autoregressive parameters only. Two additional models introduce country groupings: one for country-specific parameters and the other for global parameters. All of these models exhibit similar hierarchical priors and feature two alternatives for the global autoregressive parameters prior mean: setting it to a multivariate random walk, as in the Minnesota prior by Doan1984, or to the pooled estimator, as proposed by zellner_hong. Finally, the package implements the estimation of a single-country VAR model, enabling users to perform comparisons with panel counterparts with ease.
Importantly, the package provides Bayesian treatment of missing observations. It considers three types of such observations: the initial conditions, missing observations throughout the sampling period, and those arising from late data releases at the end of the sample resulting in rugged edges. In line with the approach initiated by guttman_bayesian_1983, Bayesian inference treats all of these missing observations as random variables, and estimates them alongside with the parameters, using a likelihood-consistent joint posterior distribution given the observed data. In other words, the missing observations are sampled in each iteration of the estimation algorithm from the distribution coherent with the model equations and its distributional assumptions. Consequently, the forecasts based on a model with missing observations accommodate the uncertainty arising from their existence, facilitating adequate risk accountability.
To authors' best knowledge, the \pkg{bpvars} package is the first \proglang{R} package dedicated specifically to Bayesian panel VAR models.\footnote{The only other similar facility for handling dynamic panel data is provided by the \pkg{BGVAR} package by BGVAR,boeck_bgvar_2022 implementing Bayesian Global VAR modelling proposed by CFH2016. Extensive code is available through \pkg{The BEAR Toolbox} by BEAR and Dieppe2016 for \proglang{MATLAB}, offering a range of Bayesian panel VARs, including that by jarocinski_responses_2010. We acknowledge other packages focused on frequentist estimation of models for dynamic panel data such as the \pkg{pvars} by pvars, \pkg{panelvar} by panelvar,panelvar_pap, and \pkg{plm} by plm,plm_pap.} Its distinctive feature is the combination of hierarchical modelling for panel VARs with missing variables treatment and a wide range of forecasting tools. This functionality is unmatched in the existing software. Additionally, the \pkg{bpvars} package is compatible in terms of the code structure, workflow and object design, and user experience with the \pkg{bsvars} package by bsvars,wozniak2024 and \pkg{bsvarSIGNs} package by WangWozniak2025,bsvarSIGNs. This compatibility provides an opportunity for synergies by simplifying the learning process for the user and expanding their toolset by alternative models. All of these packages implement frontier econometric and numerical techniques written in \proglang{C++} facilitated using the package \pkg{Rcpp} by Rcpp,eddelbuettel_seamless_2013 and eddelbuettel2011rcpp.\footnote{Similarly, linear algebra and random number generation is implemented in \proglang{C++} code using the package \pkg{RcppArmadillo} by RcppArmadillo and eddelbuettel_rcpparmadillo2014 providing headers to the library \pkg{Armadillo} by sanderson2016armadillo. Other similarities include using the progress bar from the package \pkg{RcppProgress} by RcppProgress, and the structure--rich input and output out-letting using the objects from the \proglang{R} package \pkg{R6} by R6.}
This article proceeds as follows. Section (ref) introduces our novel benchmark model and its detailed formulation, whereas Section (ref) focuses on its extensions. Section (ref) scrutinises Bayesian missing observation treatment and Section (ref) presents Bayesian forecasting for dynamic panel data. The \pkg{bpvars} package workflow detailing each of the stages of the forecasting exercise are highlighted in Section (ref), whereas the pseudo-out-of-sample forecasting exercise and the missing observation treatment are outlined in Sections (ref) and (ref), respectvely. Finally, Appendix (ref) presents the detail of the estimation algorithm for the sake of transparency and reproducibility.
Consider a benchmark model that consists of a set of country-specific Vector Autoregressions (VARs) with global prior specification. In this model, labour market outcomes are forecasted with the country-specific VAR parameters. These parameters, however, entertain a panel model feature by sharing the same prior distribution, which we call a global prior distribution. Such formulation of the model can be understood as a country-specific VAR whose parameters follow a global VAR under the prior mean. Below, the country specific model is presented and the global prior is explained.
Let an \(N\)-vector \(\mathbf{y}_{c.t} =
'\) collect the dependent variables for country \(c\) at time \(t\), where the country indicator takes values \(c \in\{ 1,\dots,C\}\) with a total number of \(C\) countries in the sample data, and the time indicator \(t \in\{ 1,\dots,T_c\}\), with the country-specific sample size \(T_c\). These variables follow a multivariate dynamic specification, namely, the Gaussian VAR model Sims1980 given by
where \(\mathbf{A}_{c.l}\) are \(N\times N\) are autoregressive matrices at lag \(l\), \(\mathbf{d}_{c.t}\) is the \(D\)-vector of deterministic terms, and \(\mathbf{A}_{c.d}\) is the \(N\times D\) matrix of the corresponding coefficients, \(\boldsymbol\epsilon_{c.t}\) is the \(N\)-vector of error terms that is normally distributed with covariance matrix \(\boldsymbol\Sigma_c\). Additionally, this dynamic model relies on initial conditions \(\mathbf{y}_{c.0}, \dots, \mathbf{y}_{c.-(p-1)}\) which are considered in Section (ref).
In the model from expressions (ref) and (ref), all the variables are treated as endogenous through the specification of the joint conditional normal distribution with covariance \(\boldsymbol\Sigma_c\) and their dynamics and temporal inter-dependencies are captured by the autoregressive parameters \(\mathbf{A}_{c.l}\). These features decide on the improved forecasting performance of VAR models compared to their univariate or static alternatives.
Rewrite this model in a matrix notation. Define a \(T_c\times N\) matrix \(\mathbf{Y}_c =
'\), a \(T_c\times K\) matrix \(\mathbf{X}_c =
'\), where \(\mathbf{x}_{c.t} =
'\), and \(K=Np+D\), a \(T_c\times N\) matrix \(\mathbf{E}_c =
'\), and a \(K\times N\) matrix \(\mathbf{A}_c =
'\). Then, the model from (ref)--(ref) can be written in an equivalent form as
where \(\mathcal{MN}_{T_c\times N}()\) denotes a matrix-variate normal distribution for a \(T_c\times N\) matrix blr1999,wozniak2016.\footnote{Let operator $\text{vec}()$ stack the columns of $\mathbf{E}_c$ one under another in a $T_cN\times 1$ vector $\text{vec}(\mathbf{E}_c)$. Then the distribution specification in (ref) with the mean matrix $\mathbf{0}_{T_c\times N}$, the row-specific covariance parameter $\boldsymbol\Sigma_c$, and the column-specific covariance parameter $\mathbf{I}_{T_c}$, is equivalent to the multivariate normal distribution $\text{vec}(\mathbf{E}_c)\sim\mathcal{N}_{T_cN}\left(\mathbf{0}_{T_cN\times1}, \boldsymbol\Sigma_c\otimes\mathbf{I}_{T_c}\right)$, where $\otimes$ denotes the Kronecker product of two matrices.} Note that the assumptions in matrix specification (ref) correspond to those from assumption (ref). It is further used to introduce the prior structure for the model.
Bayesian inference and estimation requires the specification of prior distributions of the model's parameters. These distributions are specified based on the feasibility of estimation, interpretability, and to optimise forecasting performance. In what follows, the complete prior specification is presented for this benchmark model and the corresponding estimation procedure is explained in Appendix (ref).
The main objectives motivating the prior distributions selection are:
The prior distribution for the country-specific parameters \(\mathbf{A}_c\) and \(\boldsymbol\Sigma_c\) is specified in a hierarchical manner, which allows for the estimation of the prior hyper-parameters addressing the objectives stated above. Its main characteristic is the interpretation of the prior means for country-specific autoregressive parameters \(\mathbb{E}\left[\mathbf{A}_c\right]=\mathbf{A}\) and error term covariance matrix \(\mathbb{E}\left[\boldsymbol\Sigma_c\right]=\boldsymbol\Sigma\) as global parameters, that is, invariant over countries \(c\). These prior expected values imply a VAR model with global parameters is given by:
This prior specification is implemented by assuming a convenient matrix-variate normal inverse Wishart distribution karlsson2013,wozniak2016 given by
that additionally includes hyper-parameters determining the scale and shape of the prior distribution for \(\mathbf{A}_c\) and \(\boldsymbol\Sigma_c\), namely, a \(K\times K\) column-specific covariance matrix \(\mathbf{V}\) and the shape parameter \(\nu\). Additional advantage of this specification is that it leads to a convenient Gibbs sampler for the estimation of the model.
The global parameters are estimated, which is facilitated by Bayesian hierarchical modelling. Therefore, the global autoregressive matrix, \(\mathbf{A}\), follows a matrix-variate normal distribution with the mean \(K\times N\) matrix \(m\underline{\mathbf{M}}\), \(K\times K\) column-specific covariance \(\mathbf{V}\), and \(N\times N\) row-specific covariance \(s\underline{\mathbf{S}}\) denoted by
where \(\underline{\mathbf{M}}\) and \(\underline{\mathbf{S}}\) are fixed matrices of appropriate sizes and types while scalar hyper-parameters \(m\) and \(s\) are further estimated.
The global error term covariance matrix, \(\mathbf{\Sigma}\), follows a Wishart distribution with \(N\times N\) scale matrix \(s\underline{\mathbf{S}}_\Sigma\) and shape parameter \(\underline{\mu}_\Sigma\)
where the matrix \(\underline{\mathbf{S}}_\Sigma\) and shape parameter \(\underline{\mu}_\Sigma\) are fixed, while the positive scalar \(s\) is estimated.
In the hierarchical Panel VAR model proposed here, all of the hyper-parameters of the prior in (ref) are estimated. This gives the model the advantage of fitting the data closely, while avoiding the necessity of making arbitrary choices regarding the values of these hyper-parameters.
Consequently, a prior distribution is assumed for the column-specific covariance \(\mathbf{V}\) that controls the level of shrinkage of the autoregressive parameters \(\mathbf{A}_c\) and \(\mathbf{A}\) around their respective prior means. It is set to the inverse-Wishart distribution with scale \(w\underline{\mathbf{W}}\) and shape \(\underline{\eta}\)
with the \(K\times K\) scale matrix \(\underline{\mathbf{W}}\) and the shape parameter \(\underline{\eta}\) being fixed and the positive scalar \(w\) estimated. The shape parameter \(\nu\) follows an exponential distribution with mean \(\underline\lambda\) denoted by
The prior specification is complemented by the average global persistence hyper-parameters \(m\), and scaling factors \(w\) and \(s\) following the normal, gamma, and inverted gamma 2 prior distributions respectively
To summarise, the joint prior distribution for the parameters of the model is given by
where the particular distributions are as follows:
These prior distributions depend on fixed hyper-parameters which in our notation are underscored. In what follows, we provide a justification for their default values, which are then utilized in the \pkg{bpvars} package.
The package offers two alternative values of the matrix contributing to the global autoregressive parameters prior mean, \(\underline{\mathbf{M}}\). The first choice is motivated by the interpretability of the Minnesota Prior proposed by Doan1984, in which case this matrix is set to \(
'\). It implies that the prior mean of the own lag for each of the variables is estimated by the value of another hyper-parameter \(m\) pre-multiplying this matrix in (ref). The alternative choice draws on the idea by zellner_hong where this matrix is set to the polled estimator equal to \(\left( \sum_{c=1}^{C}\mathbf{X}_c'\mathbf{X}_c \right)^{-1}\left( \sum_{c=1}^{C}\mathbf{X}_c'\mathbf{Y}_c \right)\). The prior mean specification is complemented by a normal prior distribution for the hyper-parameter \(m\) with mean \(\underline{\mu}_m = 1\) and the variance \(\underline{\sigma}_m^2 = 1\). This distribution centres the prior around the matrix \(\underline{\mathbf{M}}\), that is around random walk process in the Minnesota prior, and the pooled estimator otherwise.
The \(\mathbf{V}\) parameter prior scale matrix is set to \(\mathbf{W} = diag
\), where \(\mathbf{p}\) is a \(p\)-vector of values from 1 to \(p\), implements the feature of the Minnesota prior for global parameters where the shrinkage towards the prior mean becomes exponentially stronger for autoregressive matrices \(\mathbf{A}_i\) with increasing lag order \(i = 1,\dots,p\). The estimated hyper-parameter pre-multiplying this matrix, namely \(w\), features a gamma prior with the scale \(\underline{s}_w = 1\) and shape \(\underline{a}_w = 1\). These values make the prior distribution little informative and lets the data decide on the underlying estimate. The row-specific covariance matrix of \(\mathbf{A}\) is a product of the identity matrix \(\underline{\mathbf{S}} = \mathbf{I}_N\) and the estimated hyper-parameter \(s\), featuring the inverted gamma 2 prior with the scale and shape set to \(\underline{s}_s = 1\) and \(\underline{\nu}_s = 3\) respectively.
Similar choices are made for the prior scale of the global covariance matrix \(\mathbf{\Sigma}\) being the identity matrix \(\underline{\mathbf{S}}_\Sigma = \mathbf{I}_N\) pre-multiplied by the estimated \(s\). The shape parameter of this Wishart distribution is set to \(\underline{\mu}_\Sigma = N + 1\), which ensures finite prior variance of \(\mathbf{\Sigma}\). Similarly, the value of the shape parameter for the prior distribution in (ref) is set to \(\underline{\eta} = N + 1\).
Finally, the exponential prior for the degrees of freedom parameter \(\nu\) is set to \(\underline\lambda = 72\), which assigns 50% of the prior probability to the degrees of freedom parameter being less than 50. This choice makes the prior span the part of the parameter space implying Student-t like distribution for low values of \(\nu\), as well as close approximations of the normal distribution for values of \(\nu>30\).
The benchmark specification presented in Section (ref) has an alternative prior specification proposed by jarocinski_responses_2010. This model facilitates the estimation of the country-specific parameters and the global autoregressive matrix assuming a minimal prior structure. It features equations (ref) and (ref) with the prior distribution for the country-specific parameters given by
In this model, the global autoregressive parameters follow an improper prior distribution:
and \(s\) follows the inverted gamma 2 prior distribution as in (ref) with \(\underline{s}_s=\underline{\nu}_s= 0.001\).
As a variation on the country-specific VAR model, we consider a model with country grouping. Consider a set of \(C\) countries grouped into \(G\) groups, where each group \(g\in\{1,\dots,G\}\) contains \(C_g\) countries. Each group \(g\) has its own parameters, \(\mathbf{A}_{g}\) and \(\boldsymbol\Sigma_{g}\), defining the VAR model with country grouping for country \(c\) from group \(g\) given by
In other words, this model is specified by imposing restrictions on the country-specific parameters such that \(\mathbf{A}_c = \mathbf{A}_g\) and \(\boldsymbol\Sigma_c = \boldsymbol\Sigma_g\). The group allocations can be fixed and determined, for example, by geographical location or economic development, or they can be estimated from the data given a pre-specified number of groups \(G\).
The group-specific parameters follow the hierarchical prior distribution given by
that is complemented by the hierarchical structure as presented in in Section (ref).
Another variation is a model featuring equations (ref) and (ref) in which the country-specific parameters have group-specific global parameters. This model specification is implemented by setting the prior expectations to \(\mathbb{E}\left[\mathbf{A}_c\right]=\mathbf{A}_g\) and \(\mathbb{E}\left[\boldsymbol\Sigma_c\right]=\boldsymbol\Sigma_g\) and to let the country groupings to be fixed and specified by the user or estimated for a fixed group number \(G\). Consequently, the country-specific parameters follow the hierarchical prior distribution given by
with the priors for the group-specific global parameters given by:
Similarly, to other specifications, the remaining prior hierarchy is as described in Section (ref) with the same choices for the values of matrix \(\underline{\mathbf{M}}\).
Finally, the package allows the estimation of VAR models for individual countries. This facility is provided for comparisons with the Hierarchical Panel models. In this model, the country specific parameters are estimated independently for each country \(c\) by setting the model equations as in (ref) and (ref), and assuming the following prior distribution:
where the hyper-parameters \(m\), \(s\), \(w\), and \(\nu\) may be pre-specified, or estimated. In the latter case, the hyper-parameters \(m\) and \(\nu\) follow the distributions specified in (ref) and (ref), respectively, while those for \(s\) and \(w\) are:
An alternative prior specification for these country-specific models is the diffuse prior set as:
ensuring that the posterior mean estimator for the country-specific parameters is equal to the corresponding maximum likelihood estimator karlsson2013.
Dynamic panel data often contain missing observations due to varying starting periods of data collection in various jurisdictions, temporary breaks in data collection, sluggish data announcements, and other reasons. The \pkg{bpvars} package comprehensively implements solutions to two aspects of the problem. Firstly, it takes advantage of all of the provided observations by specifying the likelihood function for all of them and estimating the initial values \(\mathbf{y}_{c.0}, \dots, \mathbf{y}_{c.-(p-1)}\). Secondly, it treats the missing observations within the sample period as additional parameters to be estimated. Consequently, Bayesian inference specifies the joint distribution for the initial conditions, missing and observed data based on the parametric assumptions of the model and estimates the initial conditions and missing data along with the model parameters by using their implied conditional distribution given the observed data. This sampler is integrated into the Gibbs sampler presented in Appendix (ref). All the models in the \pkg{bpvars} package estimate the initial values, whereas the missing observations treatment is automatically applied if the provided data include missing values.
Begin by rewriting the model in a convenient form. Consider a vector collecting the initial conditions and data, observed and missing, \(\mathbf{y}_c\) of dimension \(T_c + p\). It stacks the country-specific vectorised data \(\mathbf{y}_c =
'\). Let the vector of initial conditions and missing observations, \(\mathbf{y}_{c.m}\), be created using a \(T_{c.m}\times T_c+p\) selection matrix \(\mathbf{S}_{c.m}\), such that \(\mathbf{y}_{c.m} = \mathbf{S}_{c.m}\mathbf{y}_c\). Similarly, let the vector of observed data, \(\mathbf{y}_{c.o}\), be created using a \(T_{c.o}\times T_c+p\) selection matrix \(\mathbf{S}_{c.o}\), such that \(\mathbf{y}_{c.o} = \mathbf{S}_{c.o}\mathbf{y}_c\).
The joint distribution of the initial conditions, missing and observed data, \(\mathbf{y}_c\), is based on rewritten model equations (ref) and (ref) that take the form:
where \(\text{bdiag}\) constructs a block-diagonal matrix from the provided blocks. For simplicity of exposition we present the notation for \(p=1\). Vector \(\boldsymbol\mu_c =
'\) stacks the deterministic components, and \(\epsilon_c =
'\) stacks the error terms. The square matrix of order \(N(T_c+p)\), \(\mathbf{H}_{\mathbf{A}_c}\), is constructed based on the autoregressive parameters such that:
The first row of matrix \(\mathbf{H}_{\mathbf{A}_c}\) is specified for the initial condition \(\mathbf{y}_0\) in line with the global prior assumption in equation (ref). Similarly, the first element of \(\boldsymbol\mu_c\) and the top-left-hand-side block of the covariance matrix in (ref) include global parameters. The remaining rows of matrix \(\mathbf{H}_{\mathbf{A}_c}\) rely on country-specific parameters in line with the model equations.
The resulting joint distribution of the initial conditions, missing and observed data, \(\mathbf{y}_c\), is the following multivariate normal distribution:
It is subsequently used to provide a conditional distribution of the missing values given the observed data:
A sampler from the distribution (ref) is used in every iteration of the Gibbs sampler, and the model's parameters are sampled conditionally on the observed data and the current draw of the initial values and missing values.
Bayesian forecasting in this panel VAR model is performed by sampling from its predictive density. In this section, the conditional predictive density -- that is, the common element with frequentist approach -- is derived, a numerical sampling procedure from this distribution is presented, and a range of techniques, including marginal, restricted, and conditional forecasting are discussed. Then the pseudo-out-of-sample forecasting is explained and complemented by the measures of forecast accuracy.
The joint conditional predictive density of the unknown future values to be predicted at the forecast horizon \(h\) from 1 to \(H\) periods ahead, given the forecast origin at period \(T_c\), denoted by \(\mathbf{y}_{c.T_c + H}, \dots, \mathbf{y}_{c.T_c + 1}\), is formed on the basis of the Bayesian Hierarchical Panel Vector Autoregressive model. This joint conditional density is factorised in terms of the one-period-ahead conditional densities
Notably, each of the conditional predictive densities on the RHS of the expression above is a one-period-ahead forecasting density determined by the model in (ref) and (ref) that can be presented as
In what follows, the density in equation (ref) is used to form the Bayesian predictive density in Section (ref) and to facilitate a straightforward numerical algorithm to obtain the random draws from it in Section (ref).
The purpose of Bayesian forecasting is to obtain a sample of random numbers drawn from the predictive density. This density is defined as the joint predictive density of the future unknown values, \(\mathbf{y}_{c.T_c + H}, \dots, \mathbf{y}_{c.T_c + 1}\), given sample data in matrices \(\mathbf{Y}_c\) and \(\mathbf{X}_c\). An important feature is that this density does not involve conditioning on the parameter values, which distinguishes Bayesian from frequentist forecasting. This lack of dependence is obtained by integrating out the parameter values from the predictive density using the posterior distribution
where the forecast origin, \(T_c\), is a period for which the last observation is available. Therefore, Bayesian forecasting accounts for estimation uncertainty by integrating out the parameters \(\mathbf{A}_c\) and \(\boldsymbol\Sigma_c\) from the predictive density karlsson2013. The estimation outcome, namely, the posterior distribution of the parameters given data, \(p\left( \mathbf{A}_c, \boldsymbol\Sigma_c \mid \mathbf{Y}_c, \mathbf{X}_c\right)\), is used in the integration. Finally, the predictive density can easily be formed using the conditional predictive density from equation (ref).
The numerical integration algorithm provides random draws from the predictive density that can be used to provide forecast summaries of interest, such as the mean and median forecast values, and the predictive interval. These summaries are easy-to-communicate characteristics of the predictive density.
The density on the LHS of equation (ref) is usually not of known analytical form and, thus, the forecasting is performed by generating random draws from the predictive density using numerical integration techniques. The predictive density in (ref) together with the conditional density as defined in Section (ref) and the one-period-ahead density in (ref) lead to the following algorithm. It provides a sample from the predictive density for country \(c\) using the sample of draws from the posterior distribution of the country-specific parameters \(\{\mathbf{A}_c^{(s)},\mathbf{\Sigma}_c^{(s)}\}_{s=1}^S\). Let \(T\) denote both the sample size for the estimation of the model and the forecast origin. We are interested in forecasting at horizons \(h\) going from 1 to \(H\). The algorithm is as follows:
Note that the algorithm is iterative in nature, that is, the draws of forecasts at shorter horizons are used to construct the mean of the density to sample those at further horizons. Also, note that the future values of exogenous variables \(\mathbf{d}_{c.T+1},\dots,\mathbf{d}_{c.T+H}\) are treated as given and must be provided.
Institutional standards for forecast reporting often require the presentation of forecasts for a subset of variables of interest despite the modelling and forecasting frameworks including more variables. Bayesian forecasting offers a simple and well-founded way to report marginal forecasts for the subset of variables. To illustrate its workings, consider a situation in which labour market forecasts are reported but forecasted using also Gross Domestic Product. Divide the vector \(\mathbf{y}_{c.t}\) into the GDP, \(gdp_{c.t}\), and labour market variables, \(\mathbf{l}_{c.t}\), such that \(\mathbf{y}_{c.t} =
'\). Then the joint predictive density can be expressed considering the corresponding block of the mean and covariance matrix:
where \(\mathbf{m}_{c.t+1}^{(l)}\) and \(m_{c.t+1}^{(gdp)}\) are appropriately factorised elements of the predictive density mean, and \(\boldsymbol\Sigma_{c}^{(l)}\), \(\boldsymbol\Sigma_{c}^{(l.gdp)}\), and \(\sigma_c^{2(gdp)}\) are those for the covariance matrix \(\boldsymbol\Sigma_c\).
Reporting marginal forecasts is based on two results. Firstly, the marginal predictive density for the labour market outcomes when the joint predictive density is as in equation (ref) is simply given by
The density in (ref) is the marginal predictive density of labour market variables, \(\mathbf{l}_{c.t+1}\). However it still depends on the past GDP values.
For one-period-ahead forecasting the task is finished because the past values of GDP are sample data. Labour market variables forecasts at more distant horizons will, however, depend on GDP forecasts. This is addressed by marginalisation of the predictive density over these GDP forecasts. This marginalisation is formally performed by integration of the joint density of GDP and labour market variables forecasts over the former values. To illustrate this integration and for simplicity of exposition without loosing the generality, consider two-period-ahead forecasts, \(h=2\), performed by a model with one autoregressive lag, \(p=1\). Then the marginal density forecast for the labour variables at horizon one and two is given by:
It is obtained by integrating out one-period-ahead forecast of GDP from the joint predictive density. This joint predictive density is constructed in (ref) by factorising it into the marginal two-period ahead predictive density for \(\mathbf{l}_{c.t+2}\) given by \(p\left( \mathbf{l}_{c.t+2}\mid \mathbf{l}_{c.t+1}, gdp_{c.t+1}, \mathbf{A}_c, \boldsymbol\Sigma_c\right)\) as defined in (ref), and the joint one-period-ahead density \(p\left(\mathbf{l}_{c.t+1}, gdp_{c.t+1} \mid \mathbf{l}_{c.t}, gdp_{c.t}, \mathbf{A}_c, \boldsymbol\Sigma_c\right)\) as defined in (ref).
The simplicity of this solution is that the integral is solved numerically automatically in the sampling procedure for Bayesian forecasting presented in Section (ref). It suffices to report the summary statistics of labour market outcomes forecasts, without the DGP forecasts, to obtain the demanded effect. The interpretation of such forecast is that labour market forecasts are averaged over all \(S\) paths of future GDP simulated from its predictive density.
The International Monetary Fund's GDP projections for all 189 countries are borrowed into our VAR setup. At the same time, including GDP as one of the endogenous variables in the Vector Autoregression is a necessary feature of a reliable forecasting model. In this context, following the argument by Sims1980, we include GDP as if it was an endogenous variable in the Vector Autoregression for modelling labour market outcomes. We proceed this way in order to estimate the model and form the joint predictive density. This enables us to reliably predict labour market outcomes based on the future trajectories of GDP.
Such conditional forecasting is an indispensable feature of Bayesian forecasting Doan1984,waggoner_conditional_1999. The one-period-ahead predictive density in expression (ref) as factorised in equation (ref) is the basis for our implementation of such forecasts. Then, forecasting of the labour market outcomes, \(\mathbf{l}_{c.t+1}\), given the contemporaneous projection of GDP, \(gdp_{c.t}\), is performed using the conditional predictive density:
In conditional forecasting as in equation (ref) GDP is never predicted. Instead, its projections are used to forecast labour market outcomes.
The approach to forecasting described in the current paper implements optimal choices based on the best Bayesian forecasting practices, thus leading to reliable predictions. This means we forecast labour market rates without applying any transformations such as taking the logarithm or differencing. The rates must be within the interval \([0,100]\) and the forecasts must respect this constraint. This is ensured following the approach by waggoner_conditional_1999, which transforms the one-period-ahead predictive densities specified in equations (ref) and (ref) into truncated multivariate normal distributions between \([0,100]\). The numerical implementation is based on the function mvrandn() from the package \pkg{TruncatedNormal} by TruncatedNormal implementing the constrained sampler by Botev. This forecasting method can be implemented without further adjustments to the estimation algorithm.
We implement pseudo-out-of-sample forecasting exercise to check the external validity of the models. This task is facilitated by using the available data to verify the forecast performance. The sample is split into training and forecast evaluation samples with the cut-off point \(T_{f_1} < T_c\), where the observations up to time \(T_{f_1}\) are used for the estimation of the model. Then, we forecast at the application-driven horizons of interest, say one and two periods ahead, \(T_{f_1} + 1\) and \(T_{f_1} + 2\) respectively. The maximum value of the forecasting horizon is denoted by \(H\) and in this example is equal to \(H=2\). The observations at periods \(T_{f_1} + 1\) and \(T_{f_1} + 2\) are used to verify prediction precision according to selected measures. This estimation and forecasting exercise is repeated iteratively using an expanding window of the training sample, that is, the cut-off point \(T_{f_1}\) is moved forward by one period at a time. For instance, in the second iteration of the exercise the cut-off point for the estimation sample is set to \(T_{f_2} = T_{f_1} + 1\). Altogether, the exercise is repeated \(F = T_c - T_{f_1} - H + 1\) times. At each iteration, the models are re-estimated, and forecasts are saved. This procedure uses parallel computations using \pkg{OpenMP} by dagum1998openmp for each of the iterations.
The pseudo-out-of-sample forecasting exercise is complemented by calculating the forecast performance measures. Let the forecasting horizon of interest be denoted by \(h\) and the observations used to verify forecast precision for a country \(c\) by \(\mathbf{y}_{c.T_{f_1}+h}, \dots, \mathbf{y}_{c.T_{f_{F}}+h}\). Consider point forecasting with the predictions given by the mean of the predictive density denoted by \(\bar{\mathbf{y}}_{c.T_{f_i}+h} = S^{-1}\sum_{s=1}^{S}\mathbf{y}_{c.T_{f_i}+h}^{(s_)}\) for \(i = 1,\dots,F\). Then the root-mean-squared-forecast error (\(RMSFE\)) for the \(n^{\text{th}}\) variable for forecasting horizon \(h\) is defined as
The RMSFE can be averaged across all variables using \(RMSFE = \sqrt{N^{-1}\sum_{n=1}^N RMSFE_n^2}\). As an alternative measure of point forecast performance, the mean-absolute-forecast error (\(MAFE\)) for the \(n^{\text{th}}\) variable for horizon \(h\) is defined as
and is aggregated to a joint measure using \(MAFE = N^{-1}\sum_{n=1}^N MAFE_n\). Both of the measures, in their single-variable and averaged versions, can be averaged further across countries by applying similar averaging formulae. For both, \(RMSFE\) and \(MAFE\), the lower the value, the better the forecast performance.
As a performance measure of density forecasts consider the predictive log-score (PLS) of a model defined by geweke_comparing_2010 as
where the predictive density ordinates \(\hat{p}\left(\mathbf{y}_{c.T_{f_i}+h} \mid \mathbf{Y}_{c.T_{f_i}}, \mathbf{X}_{c.T_{f_i}}\right)\), following gelfand1990sampling, are estimated through numerical integration using the sample from the posterior distribution of the parameters
where the conditional densities on the right-hand side are defined in expression (ref). \(PLS\) can be reported for individual variables using their marginal predictive densities or for the whole system as above. It is also subject to aggregation over countries. The higher the \(PLS\) value, the better the forecast performance.
The methods proposed in this paper are implemented in the \proglang{R} Package \pkg{bpvars} by bpvars. The package offers a simple workflow for the model specification, estimation, forecasting, and their summaries and visualisations. This section first presents the basic workflow of the package and then focuses on its particular steps. This is followed by the presentation of the pseudo-out-of-sample forecasting exercise.
Load the package to the \proglang{R} environment by executing the following line:
The dynamic panel data to be used for estimation must be formatted as a list containing in its elements ts matrices with country-specific time series with \(T_c\) rows and \(N\) columns. The package provides appropriately formatted sample data for a system containing gross domestic product, unemployment rate, employment rate, and labour force participation rate in object ilo_dynamic_panel. It is automatically loaded upon package loading. Display several first observations for Australia by running
Specify the model by running a simple function that specifies the model and all the required values, such as the number of lags, data matrices for individual countries, prior distribution hyper-parameters, starting values for the Gibbs sampler, and some characteristics for specific steps of the estimation algorithm. For instance, the code below relies on the default model setup and will specify a model with one lag, and the Minnesota prior with the prior mean for the autoregressive parameters reflecting the unit-root non-stationarity of the variables and save the model specification in object spec that is of class BVARPANEL.
This object is then provided as the first argument of the function estimate() that runs the initial 1000 iterations of the Gibbs sampler described in Appendix (ref) in the burn-in phase:
The code above reads the first argument and applies appropriate algorithms to estimate the model. More precisely, the fact that the first argument, that is the object spec, is of class BVARPANEL triggers the execution of method estimate.BVARPANEL() that reads the starting values and runs the Gibbs sampler. The object burn is of class PosteriorBVARPANEL.
As the first run achieves convergence, the second execution of the function estimate() will provide the final 1000 draws from the posterior distribution. To facilitate this provide object burn as the first argument to the estimate() function:
This execution of the function estimate() is determined by the class of the object burn and triggers the estimation algorithm using method estimate.PosteriorBVARPANEL(). This method reads the last draw from burn, sets it as the starting value, and continues the Gibbs sampler to obtain the final 1000 draws from the target posterior distribution. The object post is also of class \texttt{PosteriorBVARPANEL}.
Investigate the estimates of the autoregressive parameters for the second variable, that is unemployment rate UR, for Australia denoted by the ISO country code AUS by applying the method summary() to the estimation outcome, saving it in object post_summ, and extracting the list element called after the ISO code of the country and then extracting element A standing for the autoregressive matrix:
The displayed characteristics of the posterior distribution of the parameters in this equation include the mean, standard deviation, as well as the 5% and 95% quantiles.
Forecast the labour market outcomes and the GDP three years ahead by running the function forecast(), setting its argument horizon to 3, and saving its output in an object fore:
This function uses the draws from the posterior distribution in object post and applies Bayesian forecasting procedure described in Section (ref). The object fore is of class ForecastBVARPANEL and can be used to summarise and visualise the forecasting results. Apply the summary() function choosing the country-specific forecasts by setting the argument which_c to the ISO code of the country, and display the forecasts for the second variable, that is, the unemployment rate. The table can be interpreted as reporting the characteristics of the marginal predictive density of the unemployment rate.
Visualise the forecasts for Australia in Figure 2 by applying method plot() and setting its argument which_c to the ISO code of the country:
Compute the forecast error variance decomposition by running the function compute_variance_decompositions() and saving its output in object fevd. The forecast horizon is determined by argument horizon set to 3 and the choice of the country is made using the argument which_c. The plot is given in Figure 3.
This basic workflow can be coded using the pipe operator \textbar\textgreater that provides the object obtained by the code preceding it as the first argument of the function following it. We propose to separate the estimation part from the forecasting one. The workflow below first provides the data to the function specifying the model, which is provided to the estimate() function to run the initial 1000 draws of the burn-in stage, which in turn is provided to the estimate() function to obtain the final 1000 iterations of the Gibbs sampler. The single output of this part is saved in object post:
Having estimated the model, forecasting is executed and its outcome can be plotted using:
Similarly, the forecast error variance decompositions can computed and plotted:
In what follows, we focus on particular elements of the workflow described above. They are presented expanding the estimation and forecasting to a more realistic scenario using a model specification with exogenous variables, and with conditional forecasts given the projections of GDP.
The package is installed once prior to working with it from the CRAN repository by running the command:
The correct functioning of the \pkg{bpvars} requires the installation of the \proglang{R} package \pkg{bsvars} and other dependencies, which will proceed automatically during the installation of the \pkg{bpvars} package. Before every use of the package load it to the memory by running the command:
In order to implement the model with dummy variables, the package provides data objects formatted as required by the package functions. The dynamic panel data set contains annual observations from 1991 to 2024 for the logarithm of gross domestic product (\(gdp\)), unemployment rate (\(UR\)), employment rate (\(EPR\)), labour force participation rate (\(LFPR\)) for 189 countries are downloaded from ILO using the \proglang{R} Package \pkg{Rilostat} by Rilostat. The GDP data are taken from the IMF database. Figure 4 plots all time series contained in the data object ilo_dynamic_panel and highlights the series for Australia in each of the plots.
Inspect the dummy variables for year 2008 of the global financial crisis, and years 2020 and 2021 of the COVID pandemic:
The forecasted values of these dummy variables are necessary for forecasting:
Each of these data objects features documentation that can be accessed by the ? operator, for example:
The function specifying the model presented in Section (ref) provides ample possibilities of modifying the model. In order to implement the model with exogenous variables, one autoregressive lag, for the labour market variables and GDP, one needs to execute:
The spec object is of class BVARPANEL and includes the model specification. The arguments of the function specify_bvarPANEL\$new() provide the possibility of implementing the basic setup of the model. The argument data reads appropriately formatted data, the argument p specifies the number of autoregressive lags, and the argument exogenous reads the exogenous variables. Finally, the argument \texttt{stationary} specifies whether the variables are stationary, which subsequently sets the Minnesota prior mean for the global autoregressive matrix \(\mathbf{A}\). A value \texttt{FALSE} sets the prior mean to the random walk process for the particular variable, whereas \texttt{TRUE} sets it to white noise process.
Further investigation or customisation of the model specification can be made by displaying or altering the elements of the spec object. For instance, to display the matrix \(\underline{\mathbf{M}}\) from equation (ref) showing the prior mean of the global autoregressive matrix, execute:
This \(8\times 4\) matrix includes the autoregressive matrix in the top \(4\times 4\), the prior mean for the constant term vector in the \(5^{\text{th}}\) row, and the prior means for the global slopes on the dummies in the last three rows. The ones on the main diagonal mean that the overall global persistence hyper-parameter \(m\) is estimated. This global autoregressive matrix prior mean can be changed to the pooled estimator by zellner_hong by running spec\$set_global2pooled() similarly and the model with the prior by jarocinski_responses_2010 can be set by running spec\$set_to_Jarocinski().
Any element of the model specification can be altered but this should be done with care, making certain that the class of the provided object aligns with the class of the element being altered, and that the provided value is appropriately solicited. This includes all the fixed prior hyper-parameters from Section (ref) that are coded as elements of list in object spec\$prior. For instance, consider the prior mean \(\underline{\lambda}\) of the shape parameter \(\nu\) from equation (ref) whose default value can be accessed by:
Consider setting the median of this exponential prior distribution to thirty, which can be implemented and verified by running:
Finally, the documentation for the model specification can be accessed by running ?specify_bvarPANEL with details on its particular parts, such as the function generating the model specification specify_bvarPANEL\$new(), available following the links in the documentation just displayed.
The specification of the model with country grouping is similar to the one presented in Section (ref) and uses facility specify_bvarGroupPANEL. To specify a model with fixed group allocations the argument group_allocation must be provided as a vector of length \(C\) with integers denoting the groups \(g\). The package provides several such vectors that group countries by their geographical location with various levels of aggregation, see vectors country_grouping_region, country_grouping_subregionbroad, and country_grouping_subregiondetailed, and by their income level implmented in country_grouping_incomegroup. The example below uses the vector \texttt{country_grouping_region}.
Alternatively, the group allocation can be estimated. Such a model requires the specification of the number of groups \(G\) using the argument G as in the example below where it is set to value 2.
The spec objects created in the listings above can be used in the workflows that follow. Consult the package documentation on the model specification using ?specify_bvarGroupPANEL or on the provided country groupings using, for instance, ?country_grouping_region.
The function estimate() runs the Gibbs sampler and is implemented as two methods, estimate.BVARPANEL() and estimate.PosteriorBVARPANEL(). It starts the estimation respectively at the starting values from the model specification object of class BVARPANEL or from the last draw of the previous run provided in an object of class PosteriorBVARPANEL. The Gibbs sampler presented in Appendix (ref) performs iterations sampling random draws from the full conditional posterior distributions of the parameters of the model, namely, the country-specific parameter matrices \(\mathbf{A}_c\) and \(\mathbf{\Sigma}_c\), the global parameter matrices \(\mathbf{A}\) and \(\mathbf{\Sigma}\), the hierarchical prior estimated parameters \(\mathbf{V}\), \(\nu\), \(m\), \(w\), and \(s\), and missing observations if required.
The estimate() function uses argument S to set the number of Gibbs sampler iterations, and argument thin to set the thinning parameter. Thinning is a procedure to reduce the autocorrelation in the posterior draws, rendering the posterior estimates more efficient. Alternatively, this function is used to manage the computer memory used for estimation. It is obtained by returning every thin draw in the final sample. This procedure reduces the number of draws returned by the estimate() function to S/thin. The last argument of this function is \texttt{show_progress} that gives users the choice of whether or not to display the progress bar. Finally, the documentation can be accessed by running \texttt{?estimate.BVARPANEL} or \texttt{?estimate.PosteriorBVARPANEL}.
The model is estimated in two stages, the burn-in and the final one both run for 1000 iterations:
and the posterior estimates can be computed using the summary() function, for instance, executing:
reports the estimates of the fourth row of Poland's error term covariance matrix.
Based on authors' experience of working with the package, the data set ilo_dynamic_panel, and a model with exogenous variables ilo_exogenous_variables and \(p=1\) lag, using at least S\ =\ 1000 iterations is recommended for both burn-in and final estimation. More iterations for both of the estimation stages might be required for data sets with a larger number of variables or a model with more lags.
The function forecast() implements method forecast.PosteriorBVARPANEL() that provides draws from the predictive density as discussed in Section (ref). This forecasting routine performs forecasting for all countries and facilitates conditional forecasting given the future trajectories of some variables as well as forecasting for models with exogenous variables.
The example below presents the full extent of the forecasting capabilities of the forecast() function. Its first argument is the object post of class PosteriorBVARPANEL that contains the posterior draws from the model's posterior distribution. The remaining arguments must be aligned in terms of the forecasting horizon. The code below sets the argument horizon to value three, which requires the lists provided to the last argument, namely exogenous_forecast, to contain matrices with three rows. Appropriately constructed objects provided as arguments are described in Section (ref). The outcome of the forecasting procedure is stored in the object fore of class \texttt{ForecastsPANEL}.
The format of these arguments and the output is explained in the package documentation available by running ?forecast.PosteriorBVARPANEL.
Having obtained the draws from the predictive density saved in object fore, the user can report the forecasts for the variable of interest. The function summary() provides the summary statistics of the forecasts such as the mean, standard deviation, and the 5% and 95% quantiles for Poland. The choice of country is made by setting argument which_c to value "POL".
Finally, the forecasts can be visualised using the plot() function. The example in Figure 5 plots the forecasts for Poland with a 68% forecast interval, with optional plot's colour, main title, and axis label.
The documentation for the summary() and plot() functions can be accessed by running ?summary.ForecastsPANEL and ?plot.ForecastsPANEL.
The results reported in Figure 5 include unconstrained forecasts. The interval forecasts of unemployment rate, \(UR\), include negative values, which makes the forecasts difficult to interpret and report in official communication. The labour market rates are variables specified within the interval from 0 to 100. The package facilitates constrained forecasting of rates ensuring the draws from the predictive density fall within the interval. To apply the constraint, use the argument type in the specify_bvarPANEL function by setting it to c("real",\ "rate",\ "rate",\ "rate"). This setup results in unconstrained forecasts for the first variable, \(gdp\), and treats the remaining three, that is \(UR\), \(EPR\), and \(LFPR\), as rates. Follow by estimating the model.
In the second step, predict and plot the forecasts in Figure 6. Unemployment rate forecasts are now bind by the restrictions. The variable type can be specified for all the models in the package.
The package facilitates conditional forecasting given the future values of some variables. This feature is illustrated by forecasting labour market outcomes three years ahead given the future values of gross domestic product. The package includes an object ilo_conditional_forecasts with these future values appropriately formatted as a list of matrices. Each of these matrices includes the number of rows equal to the forecast horizon and the number of columns equal to the number of dependent variables. Its missing values coded using NA denote the future values to be forecasted, while any provided numerical values are conditioned on in the forecasting. The listing below is an example of future values of Polish gross domestic product used to conditionally forecast labour market outcomes.
To use the future values in forecasting, provide the list object ilo_conditional_forecasts as the value for argument conditional_forecast of the forecast method.
The plot of the conditional forecast also includes the provided future values.
The \pkg{bpvars} package provides a range of functions to perform a pseudo-out-of-sample forecasting exercise and report forecasting performance measures. Begin, by specifying a model using the function specify_bvarPANEL\$new() or specify_bvarGroupPANEL\$new() as in the listing below. Then, specify the parameters of the pseudo-out-of-sample (poos) forecasting exercise as in the second line of the listing. In this example, we chose the model to be estimated using 1000 draws in the burn-in run of the MCMC algorithm to obtain convergence followed by 2000 draws from the posterior and predictive distributions. The horizons at which the forecast performance is evaluated are set to 1. Finally, the argument \texttt{training_sample} set to \texttt{20} means that the first \(T_{f_1} = 20\) annual observations will be used to estimate the model in the first iteration, and a total of \(F = 14\) iterations will be performed and used to evaluate the forecasts and predictive ability of models. Both the specification of the model object \texttt{spec} and the specification of the forecasting exercise \texttt{poos} are provided as arguments to function \texttt{forecast_poos_recursively()}, which runs the pseudo-out-of-sample forecasting exercise and returns the forecasts in object \texttt{fore}. Depending on the setup of the exercise, the execution of this function might take substantial computational resources and time as it estimates the model thirteen times. One full-sample estimation is used to provide adequate starting values for the estimations in the following fourteen iterations of the exercise.
Finally, one line computes the forecast performance measures using the function compute_forecast_performance() and saves them in object fper.
Two points are due at this stage. Firstly, the forecasts from the pseudo-out-of sample exercise are of class "ForecastsPANEL" and can be as such analysed and visualised using the summary() and plot() methods. For instance, the object fore includes forecasts for Poland from the fourth iteration of the forecasting exercise:
Secondly, one can analyse the forecasting performance of this particular model by looking at the computed performance measures. For instance, the line below displays the overall predictive log-scores for the Panel VAR model
whereas the following line displays, root-mean-squared-forecast error for Poland.
We now move on to practical implementation of the forecast performance comparison with other models, namely, Panel VARs with country grouping where the group allocations are driven by geographical location,
or by the income group.
Below, we report the relative predictive log-scores one period ahead for the models with country groupings, in the order of mentioning them above in columns, to the model with country-specific parameters. The negative signs of these values indicate a higher \(PLS\) value for the latter model.
In another example we check which of the models systematically predicts better the outcomes for Poland one year ahead. Therefore, we report the relative root-mean-squared-forecast error of the models with country groupings to the model with country-specific parameters.
These values being greater than one for all outcomes show that the country-specific Panel VAR also predicts better in terms of point forecast for Poland.
Working with macroeconomic dynamic panel data often requires handling missing observations. The data set from object ilo_dynamic_panel used in the previous sections includes cubic data in which the missing observations are replaced by the values imputed by the Statistics Department at the International Labour organisation. For the sake of illustrating the handling of missing observations in the \pkg{bpvars} package, we estimate the model using data in object ilo_dynamic_panel_missing containing just over 34% of missing observations. It is formatted the same way as object ilo_dynamic_panel, but it contains NA for missing observations. The series are plotted in Figure 8.
In what follows, we proceed with the usual specification and estimation of the model. The data object ilo_dynamic_panel_missing is provided as the first argument to the function specifying the model specify_bvarPANEL. This function extracts the information regarding missing data and provides it in an appropriate form to the estimate function. Any model in the package can handle data objects with missing observations.
For illustrative purposes, Figure 9 reports Brazilian unemployment rate that includes missing observations for 1994, 2000, and 2010. In Figure 9, the black thin line represents the posterior mean of the estimated missing observations and the blue thin line those from object ilo_dynamic_panel. Despite some differences, these observations follow quite closely.
Finally, we illustrate the effect of missing observations on forecasts on the example of Poland. For this country, the data begins in 1992. However, its forecasts are affected by the uncertainty due to missing observations also in other countries. This is reflected in slightly wider predictive intervals.
We propose a new set of forecasting models for labour market outcomes that are based on newest modelling practices and inspired by solutions implemented for global forecasting for United Nation Agencies. This framework produces reliable and interpretable forecasts, as well as a range of other tools supporting diversification and accountability of risks in the process of advising and decision-making. These tools include convenient forecasting reporting and visualisation, analysis of forecast error variance decompositions, and pseudo-out-of-sample forecasting exercises for the forecasting performance assessment.
The models are implemented in the \proglang{R} package \pkg{bpvars} that is available on CRAN. The package provides a simple workflow for the model specification, estimation, forecasting, and their summaries and visualisations. The package is designed to be user-friendly and accessible to a wide range of users, including those with limited experience in Bayesian econometrics. The package builds on the best of the two worlds: excellent computational speed thanks to algorithms written in \proglang{C++} and the convenience of data analysis in \proglang{R}. Both of these features greatly simplify the work with the Bayesian dynamic hierarchical system models excelling in forecasting performance.