The package bvartools
implements functions for Bayesian inference of linear vector autoregressive (VAR) models. It separates a typical BVAR analysis workflow into multiple steps:
In each step researchers are provided with the opportunitiy to fine-tune a model according to their specific requirements or to use the default framework for commonly used models and priors. Since version 0.1.0 the package comes with posterior simulation functions that do not require to implement any further simulation algorithms. For Bayesian inference of stationary VAR models the package covers
For Bayesian inference of cointegrated VAR models the package implements the algorithm of Koop, León-González and Strachan (2010) [KLS] -- which places identification restrictions on the cointegration space -- in the following variants
For Bayesian inference of dynamic factor models the package implements the althorithm used in the textbook of Chan, Koop, Poirer and Tobias (2019).
Similar packages worth checking out are
install.packages("bvartools")
# install.packages("devtools") devtools::install_github("franzmohr/bvartools")
This example covers the estimation of a simple Bayesian VAR (BVAR) model. For further examples on time varying parameter (TVP), stochastic volatility (SV), and vector error correction (VEC) models as well as shrinkage methods like stochastic search variable selection (SSVS) or Bayesian variable selection (BVS) see the vignettes of the package and r-econometrics.com.
To illustrate the estimation process the dataset E1 from Lütkepohl (2006) is used. It contains data on West German fixed investment, disposable income and consumption expenditures in billions of DM from 1960Q1 to 1982Q4. Like in the textbook only the first 73 observations of the log-differenced series are used.
library(bvartools) # Load data data("e1") e1 <- diff(log(e1)) * 100 # Reduce number of oberservations e1 <- window(e1, end = c(1978, 4)) # Plot the series plot(e1)
The gen_var
function produces an object, which contains information on the specification of the VAR model that should be estimated. The following code specifies a VAR(2) model with an intercept term. The number of iterations and burn-in draws is already specified at this stage.
model <- gen_var(e1, p = 2, deterministic = "const", iterations = 5000, burnin = 1000)
Note that the function is also capable of generating more than one model. For example, specifying p = 0:2
would result in three models.
Function add_priors
produces priors for the specified model(s) in object model
and augments the object accordingly.
model_with_priors <- add_priors(model, coef = list(v_i = 0, v_i_det = 0), sigma = list(df = 1, scale = .0001))
If researchers want to fine-tune individual prior specifications, this can be done by directly accessing the respective elements in object model_with_priors
.
The output of add_priors
can be used as the input for user-written algorithms for posterior simulation. However, bvartools
also comes with built-in posterior simulation functions, which can be directly applied to the output of the prior specification step by using function draw_posterior
:
bvar_est <- draw_posterior(model_with_priors)
The following code sets up a simple Gibbs sampler algorithm.
# Reset random number generator for reproducibility set.seed(1234567) iterations <- 10000 # Number of saved iterations of the Gibbs sampler burnin <- 5000 # Number of burn-in draws draws <- iterations + burnin # Total number of MCMC draws y <- t(model_with_priors$data$Y) x <- t(model_with_priors$data$Z) tt <- ncol(y) # Number of observations k <- nrow(y) # Number of endogenous variables m <- k * nrow(x) # Number of estimated coefficients # Set (uninformative) priors a_mu_prior <- model_with_priors$priors$coefficients$mu # Vector of prior parameter means a_v_i_prior <- model_with_priors$priors$coefficients$v_i # Inverse of the prior covariance matrix u_sigma_df_prior <- model_with_priors$priors$sigma$df # Prior degrees of freedom u_sigma_scale_prior <- model_with_priors$priors$sigma$scale # Prior covariance matrix u_sigma_df_post <- tt + u_sigma_df_prior # Posterior degrees of freedom # Initial values u_sigma_i <- diag(1 / .00001, k) # Data containers for posterior draws draws_a <- matrix(NA, m, iterations) draws_sigma <- matrix(NA, k^2, iterations) # Start Gibbs sampler for (draw in 1:draws) { # Draw conditional mean parameters a <- post_normal(y, x, u_sigma_i, a_mu_prior, a_v_i_prior) # Draw variance-covariance matrix u <- y - matrix(a, k) %*% x # Obtain residuals u_sigma_scale_post <- solve(u_sigma_scale_prior + tcrossprod(u)) u_sigma_i <- matrix(rWishart(1, u_sigma_df_post, u_sigma_scale_post)[,, 1], k) u_sigma <- solve(u_sigma_i) # Invert Sigma_i to obtain Sigma # Store draws if (draw > burnin) { draws_a[, draw - burnin] <- a draws_sigma[, draw - burnin] <- u_sigma } }
bvar
objectsFunction bvar
can be used to collect relevant output of the Gibbs sampler in a standardised object, which can be used by further applications such as predict
to obtain forecasts or irf
for impulse respons analysis.
bvar_est <- bvar(y = model_with_priors$data$Y, x = model_with_priors$data$Z, A = draws_a[1:18,], C = draws_a[19:21, ], Sigma = draws_sigma)
Summary statistics can be obained in the usual manner:
summary(bvar_est)
The means of the posterior draws are very close to the results of the frequentist estimatior in Lütkepohl (2006).
Posterior draws can be visually inspected by using the plot
function. By default, it produces a series of histograms of all estimated coefficients.
plot(bvar_est)
Alternatively, the trace plot of the post-burnin draws can be draws by adding the argument type = "trace"
:
plot(bvar_est, type = "trace")
Summary statistics can be obtained in the usual way using the summary
method.
summary(bvar_est)
The MCMC series in object est_bvar
can be thinned using
bvar_est <- thin(bvar_est, thin = 10)
Forecasts can be obtained with the function predict
. If the model contains deterministic terms, new values have to be provided in the argument new_D
, which must be of the same length as the argument n.ahead
.
bvar_pred <- predict(bvar_est, n.ahead = 5, new_D = rep(1, 5)) plot(bvar_pred)
IR <- irf(bvar_est, impulse = "income", response = "cons", n.ahead = 8) plot(IR, main = "Forecast Error Impulse Response", xlab = "Period", ylab = "Response")
OIR <- irf(bvar_est, impulse = "income", response = "cons", n.ahead = 8, type = "oir") plot(OIR, main = "Orthogonalised Impulse Response", xlab = "Period", ylab = "Response")
GIR <- irf(bvar_est, impulse = "income", response = "cons", n.ahead = 8, type = "gir") plot(GIR, main = "Generalised Impulse Response", xlab = "Period", ylab = "Response")
bvar_fevd <- fevd(bvar_est, response = "cons") plot(bvar_fevd, main = "FEVD of consumption")
Eddelbuettel, D., & Sanderson C. (2014). RcppArmadillo: Accelerating R with high-performance C++ linear algebra. Computational Statistics and Data Analysis, 71, 1054-1063. https://doi.org/10.1016/j.csda.2013.02.005
George, E. I., Sun, D., & Ni, S. (2008). Bayesian stochastic search for VAR model restrictions. Journal of Econometrics, 142(1), 553-580. https://doi.org/10.1016/j.jeconom.2007.08.017
Kim, S., Shephard, N., & Chib, S. (1998). Stochastic volatility: Likelihood inference and comparison with ARCH models. Review of Economic Studies 65(3), 361-396.
Koop, G., León-González, R., & Strachan R. W. (2010). Efficient posterior simulation for cointegrated models with priors on the cointegration space. Econometric Reviews, 29(2), 224-242. https://doi.org/10.1080/07474930903382208
Koop, G., León-González, R., & Strachan R. W. (2011). Bayesian inference in a time varying cointegration model. Journal of Econometrics, 165(2), 210-220. https://doi.org/10.1016/j.jeconom.2011.07.007
Korobilis, D. (2013). VAR forecasting using Bayesian variable selection. Journal of Applied Econometrics, 28(2), 204-230. https://doi.org/10.1002/jae.1271
Lütkepohl, H. (2006). New introduction to multiple time series analysis (2nd ed.). Berlin: Springer.
Pesaran, H. H., & Shin, Y. (1998). Generalized impulse response analysis in linear multivariate models. Economics Letters, 58, 17-29. https://doi.org/10.1016/S0165-1765(97)00214-0
Sanderson, C., & Curtin, R. (2016). Armadillo: a template-based C++ library for linear algebra. Journal of Open Source Software, 1(2), 26. https://doi.org/10.21105/joss.00026
[^cpp]: RcppArmadillo
is the Rcpp
bridge to the open source 'Armadillo' library of Sanderson and Curtin (2016).
[^tvpvec]: In contrast to Koop et al. (2011) version 0.2.1 assumes a fixed value for the autocorrelation coefficient of the time varying cointegration space. A step for drawing this coefficient will be introduced in a future release.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.