# Stochastic Search Variable Selection

#### Franz X. Mohr

#### 2020-07-23

## Introduction

A general drawback of vector autoregressive (VAR) models is that the number of estimated coefficients increases disproportionately with the number of lags. Therefore, fewer information per parameter is available for the estimation as the number of lags increases. In the Bayesian VAR literature one approach to mitigate this so-called curse of dimensionality is stochastic search variable selection (SSVS) as proposed by George et al. (2008). The basic idea of SSVS is to assign commonly used prior variances to parameters, which should be included in a model, and prior variances close to zero to irrelevant parameters. By that, relevant parameters are estimated in the usual way and posterior draws of irrelevant variables are close to zero so that they have no significant effect on forecasts and impulse responses. This is achieved by adding a hierarchial prior to the model, where the relevance of a variable is assessed in each step of the sampling algorithm.1

This vignette presents code for the estimation of a Bayesian vector autoregressive (BVAR) model with SSVS. It uses dataset E1 from Lütkepohl (2007), which contains data on West German fixed investment, disposable income and consumption expenditures in billions of DM from 1960Q1 to 1982Q4. Following a related example in Lütkepohl (2007, Section 5.2.10) only the first 71 observations of a VAR(4) model are used. The `bvartools` package can be used to load the data and generate the data matrices:

    # devtools::install_github("franzmohr/bvartools")
    library(bvartools)

    # Load and transform data
    data("e1")
    e1 <- diff(log(e1))

    # Generate VAR
    data <- gen_var(e1, p = 4, deterministic = "const")

    # Get data matrices
    y <- data$Y[, 1:71]
    x <- data$Z[, 1:71]

## Estimation

The prior variances of the parameters are set in accordance with the semiautomatic approach described in George et al. (2008). Hence, the prior variance of the \(i\)th parameter is set to \(\tau_{1,i}^2 = (10 \hat{\sigma}_i)^2\) if this parameter should be included in the model and to \(\tau_{0,i}^2 = (0.1 \hat{\sigma}_i)^2\) if it should be excluded. \(\hat{\sigma}_i\) is the standard error associated with the unconstrained least squares estimate of parameter \(i\). For all variables the prior inclusion probabilities are set to 0.5. The prior of the error variance-covariance matrix is uninformative.

    # Reset random number generator for reproducibility
    set.seed(1234567)

    t <- ncol(y) # Number of observations
    k <- nrow(y) # Number of endogenous variables
    m <- k * nrow(x) # Number of estimated coefficients

    # Coefficient priors
    a_mu_prior <- matrix(0, m) # Vector of prior means

    # SSVS priors (semiautomatic approach)
    vs_prior <- ssvs_prior(data, semiautomatic = c(.1, 10))
    tau0 <- vs_prior$tau0
    tau1 <- vs_prior$tau1

    # Prior for inclusion parameter
    prob_prior <- matrix(0.5, m)

    # Prior for variance-covariance matrix
    u_sigma_df_prior <- 0 # Prior degrees of freedom
    u_sigma_scale_prior <- diag(0, k) # Prior covariance matrix
    u_sigma_df_post <- t + u_sigma_df_prior # Posterior degrees of freedom

The initial parameter values are set to zero and their corresponding prior variances are set to \(\tau_1^2\), which implies that all parameters should be estimated relatively freely in the first step of the Gibbs sampler.

    # Initial values
    a <- matrix(0, m)
    a_v_i_prior <- diag(1 / c(tau1)^2, m) # Inverse of the prior covariance matrix

    # Data containers for posterior draws
    iter <- 15000 # Number of total Gibs sampler draws
    burnin <- 5000 # Number of burn-in draws

    store <- iter - burnin
    draws_a <- matrix(NA, m, store)
    draws_lambda <- matrix(NA, m, store)
    draws_sigma <- matrix(NA, k^2, store)

SSVS can be added to a standard Gibbs sampler algorithm for VAR models in a straightforward manner. The `ssvs` function can be used to obtain a draw of inclusion parameters and its corresponding inverted prior variance matrix. It requires the current draw of parameters, standard errors \(\tau_0\) and \(\tau_1\), and prior inclusion probabilities as arguments. In this example constant terms are excluded from SSVS, which is achieved by specifying `include = 1:36`. Hence, only parameters 1 to 36 are considered by the function and the remaining three parameters have prior variances that correspond to their values in \(\tau_1^2\).

    # Reset random number generator for reproducibility
    set.seed(1234567)

    # Start Gibbs sampler
    for (draw in 1:iter) {
      # Draw variance-covariance matrix
      u <- y - matrix(a, k) %*% x # Obtain residuals
      # Scale posterior
      u_sigma_scale_post <- solve(u_sigma_scale_prior + tcrossprod(u))
      # Draw posterior of inverse sigma
      u_sigma_i <- matrix(rWishart(1, u_sigma_df_post, u_sigma_scale_post)[,, 1], k)
      # Obtain sigma
      u_sigma <- solve(u_sigma_i)

      # Draw conditional mean parameters
      a <- post_normal(y, x, u_sigma_i, a_mu_prior, a_v_i_prior)

      # Draw inclusion parameters and update priors
      temp <- ssvs(a, tau0, tau1, prob_prior, include = 1:36)
      a_v_i_prior <- temp$v_i # Update prior

      # Store draws
      if (draw > burnin) {
        draws_a[, draw - burnin] <- a
        draws_lambda[, draw - burnin] <- temp$lambda
        draws_sigma[, draw - burnin] <- u_sigma
      }
    }

The output of a Gibbs sampler with SSVS can be further analysed in the usual way. With the `bvartools` package the posterior draws can be colleced in a `bvar` object and the `summary` function provides summary statistics:

    bvar_est <- bvar(y = y, x = x, A = draws_a[1:36,],
                     C = draws_a[37:39, ], Sigma = draws_sigma)

    summary(bvar_est)
    #>
    #> Model:
    #>
    #> y ~ invest.1 + income.1 + cons.1 + invest.2 + income.2 + cons.2 + invest.3 + income.3 + cons.3 + invest.4 + income.4 + cons.4 + const
    #>
    #> Variable: invest
    #>
    #>              Mean      SD  Naive SD Time-series SD      2.5%       50%   97.5%
    #> invest.1 -0.09876 0.13697 0.0013697      0.0072845 -0.416464 -0.013457 0.01825
    #> income.1  0.04485 0.20903 0.0020903      0.0050995 -0.122589  0.009710 0.74880
    #> cons.1    0.10960 0.34404 0.0034404      0.0113191 -0.148747  0.018986 1.22862
    #> invest.2 -0.01384 0.05605 0.0005605      0.0016562 -0.203125 -0.002278 0.02820
    #> income.2  0.00883 0.16594 0.0016594      0.0023227 -0.256993  0.001545 0.39849
    #> cons.2    0.02677 0.20634 0.0020634      0.0037796 -0.177573  0.006213 0.63160
    #> invest.3  0.04347 0.09345 0.0009345      0.0034445 -0.023887  0.006273 0.31679
    #> income.3 -0.01103 0.14980 0.0014980      0.0020603 -0.352795 -0.002995 0.16288
    #> cons.3   -0.06246 0.25504 0.0025504      0.0060852 -0.903669 -0.013772 0.15090
    #> invest.4  0.27507 0.15086 0.0015086      0.0086043 -0.008008  0.296641 0.53061
    #> income.4 -0.05563 0.22012 0.0022012      0.0051174 -0.804125 -0.010786 0.11592
    #> cons.4   -0.03362 0.21575 0.0021575      0.0035806 -0.696050 -0.006548 0.18149
    #> const     0.01286 0.01265 0.0001265      0.0002661 -0.013699  0.012809 0.03932
    #>
    #> Variable: income
    #>
    #>                Mean       SD  Naive SD Time-series SD      2.5%        50%
    #> invest.1  1.102e-02 0.023340 2.334e-04      9.765e-04 -0.005557  1.626e-03
    #> income.1 -3.207e-02 0.093648 9.365e-04      4.740e-03 -0.334427 -3.963e-03
    #> cons.1    1.471e-01 0.200524 2.005e-03      1.380e-02 -0.027539  2.186e-02
    #> invest.2  2.319e-03 0.012462 1.246e-04      3.318e-04 -0.008327  3.713e-04
    #> income.2  4.523e-03 0.048441 4.844e-04      8.959e-04 -0.055257  1.473e-03
    #> cons.2   -1.959e-04 0.039987 3.999e-04      5.228e-04 -0.048945 -3.272e-04
    #> invest.3 -4.794e-05 0.008515 8.515e-05      9.442e-05 -0.014517 -1.062e-05
    #> income.3  1.933e-02 0.064395 6.439e-04      1.960e-03 -0.030042  4.399e-03
    #> cons.3    8.511e-03 0.052209 5.221e-04      1.179e-03 -0.040069  2.029e-03
    #> invest.4  1.632e-03 0.010241 1.024e-04      1.774e-04 -0.007799  3.558e-04
    #> income.4 -9.724e-03 0.043825 4.383e-04      1.080e-03 -0.151364 -2.703e-03
    #> cons.4    1.145e-04 0.038461 3.846e-04      4.782e-04 -0.058874 -1.345e-04
    #> const     1.714e-02 0.004044 4.044e-05      1.508e-04  0.008481  1.760e-02
    #>            97.5%
    #> invest.1 0.08080
    #> income.1 0.03235
    #> cons.1   0.60885
    #> invest.2 0.04358
    #> income.2 0.13748
    #> cons.2   0.04947
    #> invest.3 0.01230
    #> income.3 0.23939
    #> cons.3   0.16225
    #> invest.4 0.03414
    #> income.4 0.02893
    #> cons.4   0.06297
    #> const    0.02413
    #>
    #> Variable: cons
    #>
    #>                Mean       SD  Naive SD Time-series SD      2.5%        50%
    #> invest.1 -0.0022297 0.010621 1.062e-04      0.0002869 -0.037301 -0.0004459
    #> income.1  0.1540046 0.139680 1.397e-03      0.0103872 -0.016656  0.1654882
    #> cons.1   -0.2569915 0.196688 1.967e-03      0.0155472 -0.592388 -0.2902832
    #> invest.2  0.0059985 0.015759 1.576e-04      0.0006758 -0.005214  0.0010501
    #> income.2  0.3012156 0.107148 1.071e-03      0.0061845  0.005264  0.3110944
    #> cons.2    0.0137197 0.065816 6.582e-04      0.0040666 -0.040179  0.0022285
    #> invest.3  0.0003970 0.006939 6.939e-05      0.0001015 -0.007895  0.0001237
    #> income.3  0.0098598 0.043903 4.390e-04      0.0012604 -0.027608  0.0025765
    #> cons.3    0.0251849 0.066412 6.641e-04      0.0026261 -0.027473  0.0058177
    #> invest.4 -0.0041075 0.012570 1.257e-04      0.0003991 -0.046010 -0.0008125
    #> income.4  0.0249975 0.059880 5.988e-04      0.0020598 -0.021079  0.0057097
    #> cons.4   -0.0007571 0.031803 3.180e-04      0.0005640 -0.050184 -0.0001114
    #> const     0.0137934 0.003673 3.673e-05      0.0001372  0.006585  0.0138125
    #>             97.5%
    #> invest.1 0.006613
    #> income.1 0.408646
    #> cons.1   0.017079
    #> invest.2 0.056528
    #> income.2 0.479963
    #> cons.2   0.237534
    #> invest.3 0.013710
    #> income.3 0.154372
    #> cons.3   0.242239
    #> invest.4 0.005168
    #> income.4 0.213962
    #> cons.4   0.038692
    #> const    0.020933
    #>
    #> Variance-covariance matrix:
    #>
    #>                    Mean        SD  Naive SD Time-series SD       2.5%       50%
    #> invest_invest 2.175e-03 3.966e-04 3.966e-06      5.973e-06  1.529e-03 2.128e-03
    #> invest_income 4.938e-05 7.491e-05 7.491e-07      1.038e-06 -9.466e-05 4.773e-05
    #> invest_cons   1.406e-04 6.269e-05 6.269e-07      8.748e-07  2.907e-05 1.365e-04
    #> income_invest 4.938e-05 7.491e-05 7.491e-07      1.038e-06 -9.466e-05 4.773e-05
    #> income_income 1.505e-04 2.717e-05 2.717e-07      2.893e-07  1.067e-04 1.471e-04
    #> income_cons   6.644e-05 1.816e-05 1.816e-07      2.313e-07  3.604e-05 6.473e-05
    #> cons_invest   1.406e-04 6.269e-05 6.269e-07      8.748e-07  2.907e-05 1.365e-04
    #> cons_income   6.644e-05 1.816e-05 1.816e-07      2.313e-07  3.604e-05 6.473e-05
    #> cons_cons     9.773e-05 1.833e-05 1.833e-07      2.958e-07  6.825e-05 9.562e-05
    #>                   97.5%
    #> invest_invest 0.0030697
    #> invest_income 0.0002051
    #> invest_cons   0.0002767
    #> income_invest 0.0002051
    #> income_income 0.0002129
    #> income_cons   0.0001075
    #> cons_invest   0.0002767
    #> cons_income   0.0001075
    #> cons_cons     0.0001396

It is also possible to obtain the posterior inclusion probabilites of each variable by calculating the means of their posterior draws. As can be seen in the output below, only few variables appear to be relevant in the VAR(4) model, because most inclusion probabilities are relatively low. The inclusion probabilities of the constant terms are 100 percent, because they were excluded from SSVS.

    lambda <- rowMeans(draws_lambda) # Obtain means for every row
    lambda <- matrix(lambda, k) # Transform mean vector into a matrix
    lambda <- round(lambda, 2) # Round values
    dimnames(lambda) <- list(dimnames(y)[[1]], dimnames(x)[[1]]) # Rename matrix dimensions

    t(lambda) # Print
    #>          invest income cons
    #> invest.1   0.41   0.24 0.12
    #> income.1   0.11   0.18 0.63
    #> cons.1     0.17   0.43 0.72
    #> invest.2   0.13   0.11 0.17
    #> income.2   0.10   0.10 0.95
    #> cons.2     0.09   0.06 0.12
    #> invest.3   0.25   0.08 0.08
    #> income.3   0.08   0.12 0.10
    #> cons.3     0.12   0.09 0.16
    #> invest.4   0.86   0.10 0.14
    #> income.4   0.12   0.09 0.18
    #> cons.4     0.11   0.08 0.07
    #> const      1.00   1.00 1.00

Given these values, the researcher could proceed in the usual way and obtain forecasts and impulse responses based on the output of the Gibbs sampler. The advantage of this approach is that it does not only take into account parameter uncertainty, but also model uncertainty. This can be illustrated by the histogram of the posterior draws of the 6th coefficient, which describes the relationship between the first lag of income and the current value of consumption.

    hist(draws_a[6,], main = "Consumption ~ First lag of income", xlab = "Value of posterior draw")

A non-negligible mass of some 23 percent, i.e. 1 - 0.67, of the parameter draws is concentrated around zero. This is the result of SSVS, where posterior draws are close to zero if a constant is assessed to be irrelevant during an iteration of the Gibbs sampler and, therefore, \(\tau_{0,6}^2\) is used as its prior variance. On the other hand, about 67 percent of the draws are dispersed around a positive value, where SSVS suggests to include the variable in the model and the larger value \(\tau_{1,6}^2\) is used as prior variance. Model uncertainty is then described by the two peaks and parameter uncertainty by the dispersion of the posterior draws around them.

However, if the researcher prefers not want to work with a model, where the relevance of a variable can change from one step of the sampling algorithm to another, a different approach would be to work only with a highly probable model. This can be done with a further simulation, where very tight priors are used for irrelevant variables and relatively uninformative priors for relevant parameters. In this example, coefficients with a posterior inclusion probability of above 40 percent are considered to be relevant.2 The prior variance is set to 0.00001 for irrelevant and to 9 for relevant variables. No additional SSVS step is required. Everything else remains unchanged.

    # Select variables that should be included
    include_var <- c(lambda >= .4)

    # Update prior variances
    diag(a_v_i_prior)[!include_var] <- 100000 # Very tight prior close to zero
    diag(a_v_i_prior)[include_var] <- 1 / 9 # Relatively uninformative prior

    # Data containers for posterior draws
    draws_a <- matrix(NA, m, store)
    draws_sigma <- matrix(NA, k^2, store)

    # Start Gibbs sampler
    for (draw in 1:iter) {
      # 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
      }
    }

The means of the posterior draws are similar to the OLS estimates in Lütkepohl (2007, Section 5.2.10):

    bvar_est <- bvar(y = y, x = x, A = draws_a[1:36,],
                     C = draws_a[37:39, ], Sigma = draws_sigma)

    summary(bvar_est)
    #>
    #> Model:
    #>
    #> y ~ invest.1 + income.1 + cons.1 + invest.2 + income.2 + cons.2 + invest.3 + income.3 + cons.3 + invest.4 + income.4 + cons.4 + const
    #>
    #> Variable: invest
    #>
    #>                Mean       SD  Naive SD Time-series SD      2.5%        50%
    #> invest.1 -2.192e-01 0.109698 1.097e-03      1.115e-03 -0.434669 -2.193e-01
    #> income.1  2.767e-06 0.003208 3.208e-05      3.208e-05 -0.006284 -7.620e-06
    #> cons.1   -3.666e-05 0.003142 3.142e-05      3.142e-05 -0.006293 -9.744e-06
    #> invest.2 -1.348e-04 0.003190 3.190e-05      3.190e-05 -0.006378 -1.505e-04
    #> income.2  5.051e-06 0.003111 3.111e-05      3.111e-05 -0.006065  8.853e-07
    #> cons.2    8.996e-07 0.003182 3.182e-05      3.182e-05 -0.006176 -6.224e-06
    #> invest.3  1.129e-04 0.003181 3.181e-05      3.181e-05 -0.006208  7.075e-05
    #> income.3 -4.310e-05 0.003156 3.156e-05      3.241e-05 -0.006316 -1.907e-05
    #> cons.3    2.501e-05 0.003151 3.151e-05      3.080e-05 -0.006144  2.319e-05
    #> invest.4  3.281e-01 0.108540 1.085e-03      1.085e-03  0.115676  3.279e-01
    #> income.4 -1.727e-05 0.003134 3.134e-05      3.134e-05 -0.006242  5.463e-06
    #> cons.4    5.743e-05 0.003158 3.158e-05      3.098e-05 -0.006085  9.053e-05
    #> const     1.504e-02 0.006084 6.084e-05      6.084e-05  0.003001  1.506e-02
    #>              97.5%
    #> invest.1 -0.007027
    #> income.1  0.006319
    #> cons.1    0.006088
    #> invest.2  0.006088
    #> income.2  0.005990
    #> cons.2    0.006259
    #> invest.3  0.006319
    #> income.3  0.006125
    #> cons.3    0.006178
    #> invest.4  0.538804
    #> income.4  0.006147
    #> cons.4    0.006297
    #> const     0.026850
    #>
    #> Variable: income
    #>
    #>                Mean       SD  Naive SD Time-series SD      2.5%        50%
    #> invest.1  5.500e-04 0.003134 3.134e-05      3.188e-05 -0.005639  5.425e-04
    #> income.1 -4.481e-05 0.003141 3.141e-05      3.167e-05 -0.006252 -2.498e-05
    #> cons.1    2.385e-01 0.139357 1.394e-03      1.394e-03 -0.029708  2.379e-01
    #> invest.2  9.404e-05 0.003097 3.097e-05      3.097e-05 -0.006001  1.572e-04
    #> income.2  3.503e-05 0.003191 3.191e-05      3.391e-05 -0.006151  2.757e-05
    #> cons.2    4.037e-06 0.003157 3.157e-05      3.157e-05 -0.006223 -1.008e-05
    #> invest.3 -9.493e-05 0.003141 3.141e-05      3.141e-05 -0.006293 -7.739e-05
    #> income.3  1.203e-04 0.003160 3.160e-05      3.160e-05 -0.006103  1.173e-04
    #> cons.3   -1.831e-08 0.003142 3.142e-05      3.142e-05 -0.006292  1.169e-05
    #> invest.4  2.402e-04 0.003171 3.171e-05      3.120e-05 -0.005988  2.522e-04
    #> income.4 -5.227e-05 0.003161 3.161e-05      3.179e-05 -0.006333 -1.634e-05
    #> cons.4   -3.993e-05 0.003161 3.161e-05      3.161e-05 -0.006243 -6.883e-05
    #> const     1.542e-02 0.003127 3.127e-05      3.127e-05  0.009222  1.546e-02
    #>             97.5%
    #> invest.1 0.006706
    #> income.1 0.005999
    #> cons.1   0.511573
    #> invest.2 0.006087
    #> income.2 0.006272
    #> cons.2   0.006201
    #> invest.3 0.006096
    #> income.3 0.006222
    #> cons.3   0.006192
    #> invest.4 0.006386
    #> income.4 0.006066
    #> cons.4   0.006134
    #> const    0.021508
    #>
    #> Variable: cons
    #>
    #>                Mean       SD  Naive SD Time-series SD      2.5%        50%
    #> invest.1 -5.229e-04 0.003136 3.136e-05      3.136e-05 -0.006761 -5.321e-04
    #> income.1  2.618e-01 0.085712 8.571e-04      8.793e-04  0.096939  2.611e-01
    #> cons.1   -3.340e-01 0.117457 1.175e-03      1.175e-03 -0.563685 -3.346e-01
    #> invest.2  6.485e-04 0.003115 3.115e-05      3.115e-05 -0.005524  6.700e-04
    #> income.2  3.289e-01 0.076622 7.662e-04      7.831e-04  0.177269  3.292e-01
    #> cons.2   -1.558e-05 0.003176 3.176e-05      3.176e-05 -0.006196  2.898e-06
    #> invest.3  8.561e-05 0.003140 3.140e-05      3.140e-05 -0.006131  1.171e-04
    #> income.3  1.143e-04 0.003196 3.196e-05      3.196e-05 -0.006128  1.196e-04
    #> cons.3    1.699e-04 0.003182 3.182e-05      3.217e-05 -0.005986  1.373e-04
    #> invest.4 -5.370e-04 0.003128 3.128e-05      3.128e-05 -0.006716 -4.964e-04
    #> income.4  1.511e-04 0.003164 3.164e-05      3.234e-05 -0.005999  1.255e-04
    #> cons.4    3.764e-05 0.003128 3.128e-05      3.128e-05 -0.006195  6.791e-05
    #> const     1.408e-02 0.002915 2.915e-05      2.986e-05  0.008236  1.411e-02
    #>              97.5%
    #> invest.1  0.005684
    #> income.1  0.430478
    #> cons.1   -0.100868
    #> invest.2  0.006694
    #> income.2  0.478705
    #> cons.2    0.006196
    #> invest.3  0.006268
    #> income.3  0.006317
    #> cons.3    0.006458
    #> invest.4  0.005552
    #> income.4  0.006441
    #> cons.4    0.006131
    #> const     0.019682
    #>
    #> Variance-covariance matrix:
    #>
    #>                    Mean        SD  Naive SD Time-series SD       2.5%       50%
    #> invest_invest 2.092e-03 3.727e-04 3.727e-06      3.869e-06  1.488e-03 2.052e-03
    #> invest_income 4.582e-05 6.984e-05 6.984e-07      7.258e-07 -8.867e-05 4.469e-05
    #> invest_cons   1.311e-04 5.830e-05 5.830e-07      6.082e-07  2.711e-05 1.273e-04
    #> income_invest 4.582e-05 6.984e-05 6.984e-07      7.258e-07 -8.867e-05 4.469e-05
    #> income_income 1.477e-04 2.642e-05 2.642e-07      2.739e-07  1.041e-04 1.445e-04
    #> income_cons   6.597e-05 1.701e-05 1.701e-07      1.780e-07  3.726e-05 6.427e-05
    #> cons_invest   1.311e-04 5.830e-05 5.830e-07      6.082e-07  2.711e-05 1.273e-04
    #> cons_income   6.597e-05 1.701e-05 1.701e-07      1.780e-07  3.726e-05 6.427e-05
    #> cons_cons     9.455e-05 1.708e-05 1.708e-07      1.811e-07  6.704e-05 9.238e-05
    #>                   97.5%
    #> invest_invest 0.0029335
    #> invest_income 0.0001903
    #> invest_cons   0.0002560
    #> income_invest 0.0001903
    #> income_income 0.0002076
    #> income_cons   0.0001043
    #> cons_invest   0.0002560
    #> cons_income   0.0001043
    #> cons_cons     0.0001333

## Evaluation

Posterior draws can be thinned with function `thin`:

    bvar_est <- thin(bvar_est, thin = 5)

### Forecasts

The `bvar` object can also be used to produce forecasts with `predict` or impulse respones with the `irf` function.

Forecasts with credible bands can be obtained with the function `predict`. If the model contains deterministic terms, new values can be provided in the argument `new_D`. If no values are provided, the function sets them to zero. The number of rows of `new_D` must be the same as the argument `n.ahead`.

    bvar_pred <- predict(bvar_est, n.ahead = 10, new_D = rep(1, 10))

    plot(bvar_pred)

### Impulse response analysis

    OIR <- irf(bvar_est, impulse = "income", response = "cons", n.ahead = 8, type = "oir")

    plot(OIR, main = "Orthogonalised Impulse Response", xlab = "Period", ylab = "Response")

## References

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

Koop, G., & Korobilis, D. (2010). Bayesian multivariate time series methods for empirical macroeconomics. Foundations and trends in econometrics, 3(4), 267-358. https://dx.doi.org/10.1561/0800000013

Lütkepohl, H. (2007). New introduction to multiple time series analysis (2nd ed.). Berlin: Springer.

1. See Koop and Korobilis (2010) for an introduction to Bayesian VAR modelling and SSVS.↩

2. This threshold value is usually set to 50 percent. 40 percent is chosen, because it yields similar results as the restricted model in Lütkepohl (2007, Section 5.2.10).↩