Bayesian Forecasting of Labour Market Indicators using the R package bvarPANELs

by Tomasz Woźniak

Coming up next

\[ \]

modelling and forecasting framework

the R package bvarPANELs

Materials

\[ \]

Lecture Slides as a website

R script for the easy results reproduction

GitHub repo to reproduce the slides and results

bvarPANELs package repo

bvarPANELs package website

modelling and forecasting framework

modelling and forecasting framework

characterisation

• contemporary Bayesian modelling and institutional setup
• incorporates best knowledge and practices
• a balance between model size, flexibility, and its capacity to extract signal from data
• highly computational, application-specific modelling
• inspirations: UN, IPCC, ECB, FED, Christopher Sims

modelling and forecasting framework

modelling features

• Bayesian nonstationary variables handling
• system modelling
• dynamic approach
• global–to–local formulation
• embedded flexibility
• parameter estimation risk accountability

modelling and forecasting framework

forecasting features

• original non-stationary variables
• density forecasting
• conditional forecasting given \(gdp\) projections
• forecasting for models with exogenous variables
• restricted forecasting of rates

the model

Bayesian hierarchical panel VAR

• country-specific vector autoregression
• panel modelling: global model for prior mean
• flexible 3-level hierarchical prior structure
• density forecasting

the model

country-specific vector autoregression

\[\begin{align} &\\ \mathbf{y}_{c.t} &= \begin{bmatrix} gdp_{c.t} & UR_{c.t} & EPR_{c.t} & LFPR_{c.t} \end{bmatrix}'\\[3ex] \mathbf{y}_{c.t} &= \mathbf{A}_{c.1} \mathbf{y}_{c.t-1} + \mathbf{A}_{d.c}\mathbf{x}_{c.t} + \boldsymbol\epsilon_{c.t}\\[1ex] \boldsymbol\epsilon_{c.t}\mid \mathbf{y}_{c.t-1} & \sim N_4\left(\mathbf{0}_4, \boldsymbol\Sigma_c\right)\\[2ex] \end{align}\]

• subscript \(c\) is for country, and \(t\) is for time

the model

global model for the prior mean

\[\begin{align} &\\ E_\pi\left[\mathbf{A}_{c}\right] &= \mathbf{A}, \qquad \mathbf{A}_{c} = \begin{bmatrix} \mathbf{A}_{c.1} & \mathbf{A}_{d.c} \end{bmatrix}'\\[1ex] E_\pi\left[\boldsymbol\Sigma_c\right] &= \boldsymbol\Sigma\\[3ex] \mathbf{y}_{c.t} &= \mathbf{A}_{1} \mathbf{y}_{c.t-1} + \mathbf{A}_{d}\mathbf{x}_{c.t} + \boldsymbol\epsilon_{c.t}\\[1ex] \boldsymbol\epsilon_{c.t}\mid \mathbf{y}_{c.t-1} & \sim N_4\left(\mathbf{0}_4, \boldsymbol\Sigma\right) \end{align}\]

the model

flexible 3-level hierarchical prior structure

• estimate country-specific parameters: \(\mathbf{A}_c\) and \(\mathbf{\Sigma}_c\)
• estimate global parameters: \(\mathbf{A}\) and \(\mathbf{\Sigma}\)
• estimate other prior means and shrinkage levels

advantages

• flexible modelling for various types of data
• improved forecasting performance
• robustness to different prior specifications
• convenient estimation using the Gibbs sampler

the model

the local-global prior

\[\begin{align} &\\ \mathbf{A}_c, \boldsymbol\Sigma_c | \mathbf{A}, \mathbf{V}, \mathbf{\Sigma}, \nu &\sim MNIW_{K\times N}\left(\mathbf{A}, \mathbf{V}, (N - \nu - 1)\mathbf{\Sigma}, \nu\right)\\[2ex] \mathbf{A}', \mathbf{V} \mid m, w, s &\sim MNIW_{N\times K}\left(m\underline{\mathbf{M}}', w\underline{\mathbf{W}}, s\underline{\mathbf{S}}, \underline{\eta}\right)\label{eq:pgA}\\[2ex] \mathbf{\Sigma}\mid s &\sim W_{N}\left(s\underline{\mathbf{S}}_\Sigma,\underline{\mu}_\Sigma\right)\\[3ex] \end{align}\]

• \(MNIW\) is the matrix normal-inverse Wishart distribution (see Woźniak (2016))
• \(W\) is the Wishart distribution

the model

hierarchical prior

\[\begin{align} &\\ \nu &\sim\exp\left(\underline\lambda\right)\\[1ex] m &\sim N\left(\underline{\mu}_m, \underline{\sigma}_m^2\right)\\[1ex] w &\sim G\left(\underline{s}_w, \underline{a}_w\right)\\[1ex] s &\sim IG2\left(\underline{s}_s, \underline{\nu}_s\right)\\[3ex] \end{align}\]

• \(IG2\) is the inverted gamma 2 distribution

forecasting

one-period-ahead predictive density

\[\begin{align} &\\ {\color{lig}p\left(\mathbf{y}_{c.t+1}\mid \mathbf{y}_{c.t},\mathbf{A}_{c},\boldsymbol\Sigma_c\right)} & = N_4\left(\mathbf{A}_{c.1} \mathbf{y}_{c.t} + \mathbf{A}_{d.c}\mathbf{x}_{c.t+1}, \boldsymbol\Sigma_c\right)\\[5ex] \end{align}\]

• is implied by the model
• is the same as frequentist predictive density

forecasting

predictive density

\[\begin{align} p\left(\mathbf{y}_{c.t+h},\dots,\mathbf{y}_{c.t+1}\mid \mathbf{Y}_{c.t}\right) &= \int p\left(\mathbf{y}_{c.t+h},\dots,\mathbf{y}_{c.t+1},\mathbf{A}_{c},\boldsymbol\Sigma_c\mid \mathbf{Y}_{c.t}\right)d\left(\mathbf{A}_{c},\boldsymbol\Sigma_c\right)\\[1ex] &= \int{\color{lig} p\left(\mathbf{y}_{c.t+h}\mid \mathbf{y}_{c.t+h-1},\mathbf{A}_{c},\boldsymbol\Sigma_c\right)}\\[1ex] &\qquad\times\dots\\[1ex] &\qquad\times {\color{lig}p\left(\mathbf{y}_{c.t+1}\mid \mathbf{y}_{c.t},\mathbf{A}_{c},\boldsymbol\Sigma_c\right)}\\[1ex] &\qquad\times p\left(\mathbf{A}_{c},\boldsymbol\Sigma_c\mid \mathbf{Y}_{c.t}\right)d\left(\mathbf{A}_{c},\boldsymbol\Sigma_c\right) \end{align}\]

• conditional density structure determines the sampling algorithm

the R package bvarPANELs

the R package bvarPANELs

features

• precise estimation and forecasting
• simple workflows in R
• excellent computational speed
  • frontier econometric and numerical techniques
  • algorithms written in C++
• extensive documentation
• up-to standards: ready for publication
• install the package from the GitHub repo

devtools::install_github("bsvars/bvarPANELs")

the R package bvarPANELs

the R package bvarPANELs

the R package bvarPANELs

load data

library(bvarPANELs)             # load the package
data(ilo_dynamic_panel)         # load the data
ilo_dynamic_panel$COL           # show the data for Colombia

Time Series:
Start = 1991 
End = 2023 
Frequency = 1 
          gdp        UR      EPR     LFPR
1991 25.53464 10.120000 59.96793 66.71999
1992 25.57429  9.440000 60.61868 66.93758
1993 25.62675  7.800000 62.04586 67.29486
1994 25.68326  8.250000 62.11328 67.69840
1995 25.73398  8.720000 62.11020 68.04361
1996 25.75433 11.810000 60.03138 68.07051
1997 25.78805 12.139999 59.96012 68.24507
1998 25.79374 15.000000 57.90753 68.12650
1999 25.75079 20.059999 53.95368 67.49272
2000 25.77962 20.520000 53.75205 67.62965
2001 25.79626 15.040001 57.46694 67.64000
2002 25.82099 14.481000 56.95565 66.60000
2003 25.85942 13.220999 58.52376 67.44000
2004 25.91138 13.717001 57.54213 66.69000
2005 25.95854 11.061999 58.11209 65.34000
2006 26.02355 11.091076 57.45149 64.61836
2007 26.08876 11.204000 56.73176 63.89000
2008 26.12106 11.273000 56.92724 64.16000
2009 26.13240 12.066000 58.87181 66.95000
2010 26.17636 11.153000 59.59857 67.08000
2011 26.24353 10.288000 60.56457 67.51000
2012 26.28191  9.959000 61.29091 68.07000
2013 26.33198  9.246000 61.05022 67.27000
2014 26.37599  8.799000 61.09555 66.99000
2015 26.40512  8.572000 61.27505 67.02000
2016 26.42578  8.922000 60.49401 66.42000
2017 26.43928  9.086000 60.01233 66.01000
2018 26.46460  9.360000 59.13354 65.24000
2019 26.49597 10.280000 57.70790 64.32000
2020 26.42140 15.983000 53.28358 63.42000
2021 26.52396 13.897999 55.38942 64.33000
2022 26.59432 10.547000 56.49851 63.16000
2023 26.60042  9.565347 57.65713 63.75558

the R package bvarPANELs

specify and estimate the model

spec = specify_bvarPANEL$new(                           # specify the model
  ilo_dynamic_panel,                                    # data
  exogenous = ilo_exogenous_variables,                  # exogenous variables
  stationary = c(FALSE, FALSE, FALSE, TRUE),            # stationarity (determines prior mean)
  type = c("real", "rate", "rate", "rate")              # variable types
)

burn = estimate(spec, S = 10000, show_progress = FALSE) # run the burn-in
post = estimate(burn, S = 10000)                        # estimate the model

**************************************************|
bvarPANELs: Forecasting with Bayesian Hierarchical|
            Panel Vector Autoregressions          |
**************************************************|
 Progress of the MCMC simulation for 10000 draws
    Every draw is saved via MCMC thinning
 Press Esc to interrupt the computations
**************************************************|

the R package bvarPANELs

forecast labour market outcomes

fore = forecast(                                     # forecast the model 
  post,                                              # estimation output
  horizon = 6,                                       # forecast horizon
  exogenous_forecast = ilo_exogenous_forecasts,      # forecasts for exogenous variables
  conditional_forecast = ilo_conditional_forecasts   # gdp projections
)

**************************************************|
bvarPANELs: Forecasting with Bayesian Hierarchical|
            Panel Vector Autoregressions          |
**************************************************|
 Progress of sampling 10000 draws from
    the predictive density for 189 countries
    Press Esc to interrupt the computations
**************************************************|

the R package bvarPANELs

forecast labour market outcomes

plot(fore, "COL", main = "Forecasts for Colombia")

the R package bvarPANELs

forecast labour market outcomes

summary(fore, "COL")$variable2

 **************************************************|
 bsvars: Bayesian Structural Vector Autoregressions|
 **************************************************|
   Posterior summary of forecasts                  |
 **************************************************|

       mean       sd 5% quantile 95% quantile
1 10.294857 1.476505    7.872032     12.70846
2  9.841204 1.955355    6.707715     13.01481
3  9.270488 2.261769    5.611791     12.99134
4  8.846041 2.524810    4.745539     13.01309
5  8.492730 2.729845    4.030065     12.93011
6  8.222295 2.887869    3.420855     12.95301

the R package bvarPANELs

forecast error variance decomposition

post |>                                              # estimation output
  compute_variance_decompositions(horizon = 6) |>    # compute variance decompositions
  plot(which_c = "COL")                              # plot variance decompositions