# Structural Bayesian Vector Autoregression

**URL:** <https://discourse.pymc.io/t/structural-bayesian-vector-autoregression/16603>\
**Category:** General\
**Tags:** modeling\
**Created:** [February 28, 2025, 11:20am UTC](https://discourse.pymc.io/t/structural-bayesian-vector-autoregression/16603 "2025-02-28T11:20:15Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![aeturrell](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/aeturrell/32/3489_2.png) [@aeturrell](https://discourse.pymc.io/u/aeturrell)\
**Post date:** [February 28, 2025, 11:20am UTC](https://discourse.pymc.io/t/structural-bayesian-vector-autoregression/16603/1 "2025-02-28T11:20:15Z")

</div>

Really enjoyed the posts by, eg, @ricardoV94 on BVARs.

I’m interested in creating a simple structural bayesian vector autoregression using PyMC. These are like VAR(p) models but have a number of potential restrictions imposed:

- sign restrictions
- long run restrictions
- theoretically informed or narrative restrictions

In an economic context, that last set of restrictions might include rules that:

- Only shocks to output can shift output contemporaneously.
- Only shocks to output and inflation can shift inflation contemporaneously.
- All shocks can affect the interest rate contemporaneously.

A good example of that combines all of this functionality can be found in the R packages [bsvarSIGNs](https://bsvars.org/bsvarSIGNs/). There’s also a neat worked through example [here](https://kevinkotze.github.io/ts-8-tut/).

The closest example I’ve been able to find in the PyMC literature is this [BVAR page](https://www.pymc.io/projects/examples/en/latest/time_series/bayesian_var_model.html) but it doesn’t have any ‘structural’ restrictions in it. Statsmodels supports SVARS too; the page is [here](https://www.statsmodels.org/dev/generated/statsmodels.tsa.vector_ar.svar_model.SVAR.html).

Is it possible to implement SBVARs in PyMC? And how would one do it?

---

<div class="post-metadata">

**Author:** ![jessegrabowski](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jessegrabowski/32/5010_2.png) [@jessegrabowski](https://discourse.pymc.io/u/jessegrabowski)\
**Post date:** [February 28, 2025, 11:32am UTC](https://discourse.pymc.io/t/structural-bayesian-vector-autoregression/16603/2 "2025-02-28T11:32:32Z")

</div>

You can use the `BayesianVARMA` model in `pymc_extras.statespace`. Example notebooks [here](https://github.com/pymc-devs/pymc-extras/blob/main/notebooks/VARMAX%20Example.ipynb) and [here](https://github.com/jessegrabowski/statespace-presentation/blob/main/02.%20VAR%20Models%20and%20Bayesian%20IRFs.ipynb), and I also talk about them in [this video](https://youtu.be/G9VWXZdbtKQ?feature=shared&t=1916).

You have full control over the coefficient matrix and covariance matrix when doing inference (so you can impose zeros in either), and over the covariance matrix used when doing impulse response function. There’s an option to do cholesky decomposition on the covariance matrix out of the box, but it needs a bit of work. PRs welcome 🙂

---

<div class="post-metadata">

**Author:** ![aeturrell](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/aeturrell/32/3489_2.png) [@aeturrell](https://discourse.pymc.io/u/aeturrell)\
**Post date:** [February 28, 2025, 12:01pm UTC](https://discourse.pymc.io/t/structural-bayesian-vector-autoregression/16603/3 "2025-02-28T12:01:21Z")

</div>

Super helpful, thank you (and for the very prompt reply!)

So, if I’ve understood correctly, to implement the equivalent of the `a.mat` and `b.mat` from [this example](https://kevinkotze.github.io/ts-8-tut/):

```plaintext
# a.mat
## [,1] [,2] [,3]
## [1,] NA 0 0
## [2,] NA NA 0
## [3,] NA NA NA
# b.mat
## [,1] [,2] [,3]
## [1,] NA 0 0
## [2,] 0 NA 0
## [3,] 0 0 NA

```

The lines in your first example one would have to change would be:

```python
    state_chol, _, _ = pm.LKJCholeskyCov(
        "state_chol", eta=1, n=bvar_mod.k_posdef, sd_dist=pm.Exponential.dist(lam=1)
    )

    ar_params = pm.Normal("ar_params", mu=0, sigma=1, dims=ar_dims)
    state_cov = pm.Deterministic("state_cov", state_chol @ state_chol.T, dims=state_cov_dims)

```

Am I right in saying that a.mat is ar\_params and b.mat is state\_cov in your example? And apologies if this is a stupid question, but how do you modify the pm.Normal and pm.Deterministic lines to implement matrix structures with zeros in the right places?

---

<div class="post-metadata">

**Author:** ![jessegrabowski](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jessegrabowski/32/5010_2.png) [@jessegrabowski](https://discourse.pymc.io/u/jessegrabowski)\
**Post date:** [February 28, 2025, 2:03pm UTC](https://discourse.pymc.io/t/structural-bayesian-vector-autoregression/16603/4 "2025-02-28T14:03:19Z")

</div>

The statespace package is always working with a formula like this:

\begin{align} y\_t &=Z x\_t + d + \eta\_t, & \eta\_t \sim N(0, H) \\ x\_t &= T x\_{t-1} + c + R \varepsilon\_t, & \varepsilon \sim N(0, Q) \end{align}

Looking at this formula from the statsmodel docs:

Ay\_t = A\_1 y\_{t-1} + \ldots + A\_p y\_{t-p} + B \varepsilon\_t, \quad \varepsilon\_t \sim N(0, I\_k)

I write I\_k because we assume y\_t has shape k \times 1 (there are k endogenous variables) and we’re working with a VAR(p, 0) model. I assume the errors are spherical in statsmodels (that is anything we know about the covariance is pre-factored out into B), but I actually don’t know – that’s an important detail to check.

Anyway, to get from one to the other, you’re going to have to do some rewriting:

y\_t = A^{-1} A\_1 y\_{t-1} + \dots + A^{-1} A\_p y\_{t-p} + A^{-1} B \varepsilon\_t

Call A^{-1} A\_p = \Gamma\_p, so the system is y\_t = \sum\_{s=1}^p \Gamma\_p y\_{t-p} + A^{-1} B \varepsilon\_t

In the statespace framework, you work with everything “stacked” into a single state vector. So define:

x\_t = \begin{bmatrix} y\_t \\ y\_{t-1} \\ \vdots \\ y\_{t-p} \end{bmatrix}

Next you need a map from \{\Gamma\_1, \dots, \Gamma\_p\} to T. This is done with vectorization. ~~Define \text{rvec}(X) as a “row vectorize” operator that stacks the rows of a matrix into a single row vector from top to bottom (the usual \text{vec} operator works on columns, so I’m try to avoid confusion)~~

`<BEGIN EDIT>`

The formula I had before for the T matrix was wrong. I should have taken my advice and looked at the digression. Here’s an example of T for a VAR(2, 0). The indices are in the form (equation, coefficient, lag order):

\begin{bmatrix} a\_{x,x,1} & a\_{x,y,1} & a\_{x,z,1} & a\_{x,x,2} & a\_{x,y,2} & a\_{x,z,2} \\ a\_{y,x,1} & a\_{y,y,1} & a\_{y,z,1} & a\_{y,x,2} & a\_{y,y,2} & a\_{y,z,2} \\ a\_{z,x,1} & a\_{z,y,1} & a\_{z,z,1} & a\_{z,x,2} & a\_{z,y,2} & a\_{z,z,2} \\ 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \end{bmatrix}

You see there there are 4 square blocks. On top there’s one for the lag 1 coefficients, with equations on the rows and coefficients on the columns, then the same for lag 2. Then there’s a k \times k identity matrix, and a k \times k zero matrix. So in block form, we have:

T = \begin{bmatrix} \Gamma\_1 & \Gamma\_2 \\ I\_k & 0\_k \end{bmatrix}

Continuing the pattern, for VAR(3) one has:

T = \begin{bmatrix} \Gamma\_1 & \Gamma\_2 & \Gamma\_3 \\ I\_k & 0\_k & 0\_k \\ 0\_k & I\_k & 0\_k \end{bmatrix}

And so on.

`<END EDIT>`

There’s a section in the linked notebook called " Digression: A note on representation" that would be useful to look at now, to see how this would look in a VAR(1,0) case with 3 endogenous variables.

Now we just need R and Q. \varepsilon\_t isn’t going to change. So the job of R is to make it conform in shape to the new state vector x\_t. Nothing in x\_t gets shocks except the first k entries (the ones associated with t-1, so R is going to have a lot of zeros. One choice is to define R = \begin{bmatrix} A^{-1} B \\ 0 \\ \vdots \\ 0 \end{bmatrix}, which implies that Q = I\_k.

Alternatively, you could decide to define R = \begin{bmatrix} I\_k \\ 0 \\ \vdots \\ 0 \end{bmatrix}, in which case you’re now forced to pick Q = A^{-1} B B^T A^{-T} (that follows by analolgy from multiplying a scalar random variable a \sim N(0, 1) by a scalar b gives a \cdot b \sim N(0, b^2)

With that we’re basically done. Define Z = \begin{bmatrix} I\_k & 0 & \dots & 0 \end{bmatrix} (it just picks out the top k elements of the x\_t vector), and we have everything in statespace form.

From here, you have two choices:

Option 1 is to just use the pre-defined `BayesianVARMA` model directly. This will force you to go with the covariance definition of R and Q, because it internally sets R to the one with the I\_k on top.  
a. Going with the example you gave, you need 6 coefficients in the A matrix, Sample those, then arrange them into a lower triangular matrix.  
b. Next, sample your VAR coefficient matrices. Then construct the \Gamma matrices, and finally ravel, concatenate, and reshape into the form expected by the model  
c. Finally, do the same for the covariance, computing A^{-1} B B^T A^{-T}

This model would look something like this (priors are dumb on purpose):

```py
from pymc_extras.statespace import BayesianVARMAX
import pymc as pm
import pytensor.tensor as pt
import numpy as np

mod = BayesianVARMAX(order=(2, 0), 
                     endog_names=['x', 'y', 'z'],
                     stationary_initialization=True)

k = mod.k_posdef
A_mask = np.array([[1, 0, 0],
                   [1, 1, 0],
                   [1, 1, 1]])
B_mask = np.eye(3)

n_nonzero_A = A_mask.ravel().sum()
n_nonzero_B = B_mask.ravel().sum()

with pm.Model(coords=mod.coords) as m:
    A_coefs = pm.Normal('A_coefs', 0, 1, shape=(n_nonzero_A, ))
    B_coefs = pm.Normal('B_coefs', 0, 1, shape=(n_nonzero_B, ))
    
    # .set is necessary in pytensor to allocate a value inside a matrix. Here, we use A_mask 
    # and B_mask to do fancy indexing 
    A = pm.Deterministic('A',
                         pt.zeros((k, k))[A_mask.astype(bool)].set(A_coefs), 
                         dims=['observed_state', 'observed_state_aux'])
    B = pm.Deterministic('B', 
                         pt.zeros((k, k))[B_mask.astype(bool)].set(B_coefs),
                         dims=['observed_state', 'observed_state_aux'])
    
    # According to the build report, ar_lags needs to be shape ['observed_state', 'ar_lags', 'observed_state_aux']. 
    # But we need compute A_inv @ A_i for each lag, so first we use ar_lag as a "batch dim",
    # then we'll fix it later.
    lag_coefs = pm.Normal('lag_matrices', 0, 1, dims=['ar_lag', 'observed_state', 'observed_state_aux'])

    # Pytensor will automatically vectorize solve across the left-most dimension. We have to tell it
    # that b has 2 dims, otherwise it will solve A against column vectors with 2 batch dims, which
    # isn't what we want.
    Gammas = pt.linalg.solve(A, lag_coefs, b_ndim=2, assume_a='gen')
    
    # Now that we have Gammas, we can use swapaxes to put the dims into the right order.      
    # Here, the name of the variable is extremely important -- this is what will actually go into the T matrix.
    ar_params = pm.Deterministic('ar_params', pt.swapaxes(Gammas, 0, 1), dims=['observed_state', 'ar_lag', 'observed_state_aux'])
    
    # Similarly, we want to build to covariance matrix, then give it the right name in a deterministic
    # First, get A^-1 @ B
    A_inv_B = pt.linalg.solve(A, B, assume_a='gen')

    # Then "square" it and put it in a deterministic with the right name.
    state_cov = pm.Deterministic('state_cov', A_inv_B @ A_inv_B.T, dims=['shock', 'shock_aux'])
    
    # Using all nan data will draw unconditional trajectories based on the prior (the kalman filter, having nothing to learn from, will just do nothing)
    mod.build_statespace_graph(np.full((100, 3), np.nan),
                               mode='JAX')
    
# This model is a mess -- the priors are really nasty and you'll have a hard time sampling numerically stable trajectories
# JAX has somewhat more stable linalg routines, so we want to use it to do forward sampling too.
with freeze_dims_and_data(m):
    prior = pm.sample_prior_predictive(compile_kwargs={'mode':'JAX'})

```

Option 2 is to make a custom statespace model to handle all of these reshapes and solves for you, and ideally open a PR to include it for others to use in the future! Explanation [here](https://github.com/pymc-devs/pymc-extras/blob/main/notebooks/Making%20a%20Custom%20Statespace%20Model.ipynb), and some conversations [here](https://discourse.pymc.io/t/does-the-state-space-model-support-time-dependent-transition-matrices/16412/23).

---

<div class="post-metadata">

**Author:** ![aeturrell](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/aeturrell/32/3489_2.png) [@aeturrell](https://discourse.pymc.io/u/aeturrell)\
**Post date:** [February 28, 2025, 2:19pm UTC](https://discourse.pymc.io/t/structural-bayesian-vector-autoregression/16603/5 "2025-02-28T14:19:13Z")

</div>

Amazing, thank you so much for this. Think I would have been quite stuck without your example. I will have a go with it. Yeah, it would be amazing to have this is all wrapped up nicely, at least into a custom state space model but ideally into a package that handled other typical SVAR type actions (eg computing and plotting impulse response functions). Really appreciate your help!

---

<div class="post-metadata">

**Author:** ![jessegrabowski](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jessegrabowski/32/5010_2.png) [@jessegrabowski](https://discourse.pymc.io/u/jessegrabowski)\
**Post date:** [February 28, 2025, 2:20pm UTC](https://discourse.pymc.io/t/structural-bayesian-vector-autoregression/16603/6 "2025-02-28T14:20:41Z")

</div>

IRFs are handled by the `BayesianVARMA` class, see the example notebooks for details.