SSVS-VAR — R cross-check¶

Bayesian variable selection, Part 3 (R)¶

Validates gsnvar_python.ipynb. There is no drop-in SSVS-VAR routine on this R toolchain (the bvartools building blocks would have to be assembled by hand), so — as with the horseshoe and heteroskedasticity notebooks — we implement the same George-Sun-Ni spike-and-slab Gibbs from scratch in base R and check it reproduces the coefficient restriction map: own lags and the Taylor-rule terms kept, cross-terms switched off. The coefficient block is identical to the Python engine (vectorised GLS + semiautomatic spike/slab); for brevity the error covariance is drawn from an inverse-Wishart rather than the selected $\Psi$ factor, which has a negligible effect on the coefficient inclusions.

Data: gsnvar_data.csv — US real GDP growth, CPI inflation, 3-month T-bill; quarterly 1959–2009 (statsmodels macrodata).

In [1]:
.libPaths('C:/Users/user/R/win-library/4.6')
d <- read.csv('gsnvar_data.csv'); vn <- colnames(d); Y <- as.matrix(d); p <- 2; m <- ncol(Y)
T0 <- nrow(Y); Yt <- Y[(p+1):T0, ]                                  # build the VAR design
X  <- cbind(1, Y[p:(T0-1), ], Y[(p-1):(T0-2), ]); Tn <- nrow(Yt); k <- ncol(X)
labs <- c('const', paste0(rep(vn, 2), '.L', rep(1:2, each=m)))
XtX <- crossprod(X); Bols <- solve(XtX, crossprod(X, Yt)); E <- Yt - X %*% Bols
Sig <- crossprod(E)/Tn
se <- sqrt(outer(diag(solve(XtX)), diag(Sig)))                      # k x m OLS std errors
c0 <- 0.1; c1 <- 10; kap0 <- c0*se; kap1 <- c1*se; w <- 0.5         # semiautomatic spike/slab scales
isint <- matrix(FALSE, k, m); isint[1, ] <- TRUE

set.seed(1); ndraw <- 6000; burn <- 3000
B <- Bols; Sinv <- solve(Sig); gam <- matrix(1, k, m); GAM <- matrix(0, k, m); nk <- 0
for (it in 1:(ndraw+burn)) {
  d2 <- ifelse(gam==1, kap1^2, kap0^2); d2[isint] <- 100            # prior variances (spike/slab)
  Prec <- kronecker(Sinv, XtX) + diag(1/as.vector(d2))
  rhs  <- as.vector(crossprod(X, Yt) %*% Sinv)                      # vec(X'Y Sigma^{-1})
  L <- chol(Prec); mu <- backsolve(L, forwardsolve(t(L), rhs))
  B <- matrix(mu + backsolve(L, rnorm(k*m)), k, m)                  # (1) beta | Sigma, gamma
  l1 <- log(w)-log(kap1)-0.5*B^2/kap1^2; l0 <- log(1-w)-log(kap0)-0.5*B^2/kap0^2
  gam <- matrix(as.numeric(runif(k*m) < 1/(1+exp(l0-l1))), k, m); gam[isint] <- 1   # (2) gamma | beta
  E <- Yt - X %*% B
  Sinv <- rWishart(1, Tn+m+1, solve(diag(m) + crossprod(E)))[, , 1]  # (3) Sigma^{-1} | B  (inverse-Wishart)
  if (it > burn) { GAM <- GAM + gam; nk <- nk + 1 }
}
incl <- GAM/nk; dimnames(incl) <- list(labs, vn)
cat('coefficient inclusion probabilities (rows = regressor, cols = equation):\n'); print(round(incl, 2))
cat('\ncross-check -> Python GSN: own lags + Taylor-rule terms kept, ~half of lags dropped -> same restriction map.\n')

options(repr.plot.width=6.5, repr.plot.height=5)
zz <- t(incl[k:1, ])                                                # heatmap: equations x regressors
image(1:m, 1:k, zz, col=hcl.colors(20,'RdYlBu',rev=TRUE), zlim=c(0,1), axes=FALSE, xlab='', ylab='',
      main='SSVS-VAR restriction map (base R): inclusion probabilities')
axis(1, 1:m, vn); axis(2, 1:k, rev(labs), las=1)
for (i in 1:m) for (j in 1:k) text(i, j, sprintf('%.2f', zz[i, j]), cex=0.8)

coefficient inclusion probabilities (rows = regressor, cols = equation):
         dlgdp infl  ffr
const     1.00 1.00 1.00
dlgdp.L1  0.98 0.11 0.30
infl.L1   0.10 1.00 0.10
ffr.L1    0.09 0.45 1.00
dlgdp.L2  0.62 0.22 0.65
infl.L2   0.71 0.99 0.71
ffr.L2    0.16 0.29 0.10

cross-check -> Python GSN: own lags + Taylor-rule terms kept, ~half of lags dropped -> same restriction map.
No description has been provided for this image

Same restriction map¶

  • Base R reproduces the George-Sun-Ni coefficient selection. The from-scratch Gibbs finds the same structure as the Python engine and PyMC: each variable's own lags are retained, the interest-rate equation keeps its Taylor-rule loadings on lagged output and inflation, and roughly half of the lag coefficients are pushed into the spike. Two independent implementations agreeing on which restrictions the data support confirms the selection is a property of the model and data, not the code. (This base-R check exercises the coefficient layer only; the $\Psi$ covariance-selection layer lives in the Python notebook.)

References¶

  • George, E. I., Sun, D. & Ni, S. (2008). Bayesian stochastic search for VAR model restrictions. J. Econometrics 142, 553–580.
  • Koop, G. & Korobilis, D. (2010). Bayesian multivariate time series methods for empirical macroeconomics. Found. Trends Econometrics 3, 267–358.