Sorry for slow reply, I’ve been traveling.
If you have no missing values in the observed time series, running the KF (what pymc-statespace does) is indeed wasteful. I just merged this PR that switches to a fast-path for computing the VAR likelihood. It takes out the loop entirely, so you can get the 1000x scale speedups I report in the pr.
Basically if you have no missing data and no measurement error, you can just write the VAR as a big regression problem. If you do this, you will find that the likelihood you need to calculate only depends on sufficient statistics of the data, which is very low-dimensional compared to the full system. Since you have a static A, it’s even cheaper because you can pre-compute a bunch of stuff. The statespace machine is still nice to have for post-estimation, though, because this setup doesn’t give you IRF or forecasting “for free”, but it does give you a fast path to parameter inference. If you use the latest version of pymc-extras, each model can implement a custom make_likelihood method that lets it override the default Kalman Filter recursions. That gives you the best of both worlds, you can do fast inference and have the post-estimation machinery.
The specific details are as follows. You have a VAR(p) that looks like:
\begin{align}
x_t &= A_1 x_{t-1} + A_2 x_{t-2} + \ldots + B z_t + \varepsilon_t \\
\varepsilon_t &\sim N(0, Q)
\end{align}
z_t are the exogenous regression data, in your case:
z_t = \begin{bmatrix}1 \\ t \end{bmatrix}
Again, all of the x_{t-\tau} are observed, so there’s no missing data. There’s also no measurement error, so the likelihood is just going to be based on the distribution of the residuals. Stack all this into a single system:
x_t = \bar{A} \tilde{x_t} + B z_t + \varepsilon_t
Where \bar{A} = \begin{bmatrix} A_1 & A_2 & \ldots & A_p \end{bmatrix} and \tilde{x_t} = \begin{bmatrix} x_{t-1} & x_{t-2} & \ldots & x_{t-p}\end{bmatrix}^T
This equation describes one time-step, but we want to do the whole thing in one shot. So take the input data and do an im2col operation on it to stack it so that each row corresponds to a regression.
The regression targets are going to be the “current” values, starting from p (since you don’t have data for the right-hand side before then)
Y = \begin{bmatrix}x_p^T \\ x_{p+1}^T \\ \vdots \end{bmatrix}
While the regressors are going to be the corresponding shifted vector \tilde{x}_t, again starting at time p:
X = \begin{bmatrix} \tilde{x}^T_p \\ \tilde{x}^T_{p-1} \\ \vdots \end{bmatrix}
You will also stack up all the z_t terms the same way into Z, giving:
Y = X\bar{A}^T + ZB^T + E
For shapes, you have A of shape (k, kp) and B of shape k, 2. Y is shape (T-p+1, k) and X is shape (T-p+1, kp). E of shape (T-p+1, k) is the matrix of all the residuals for each time step, which are assumed to be independent between steps.
For brevity, define N = T - p + 1 as the number of “usable time steps” you VAR(p) implies.
Now we want a likelihood. This could be written as a matrix normal distribution but that sucks so let’s ravel everything with the column-wise \text{vec} operator, so vec(x) = x.T.ravel() It’s linear so:
\text{vec}(Y) = \text{vec}(X \bar{A}^T + Z B^T) + \text{vec}(E)
As mentioned, the rows of E are iid with constant covariance Q, so:
\text{vec}(E) \sim N(0, Q \otimes I_N)
Where \otimes is the kronecker product. This gives a single multivariate normal likelihood for the whole system:
\text{vec}(Y) \sim N(\text{vec}(X \bar{A}^T + Z B^T), Q \otimes I_N)
Or, defining residuals R = Y - X \bar{A}^T - Z B^T:
\text{vec}(R) \sim N(0, Q \otimes I_N)
So the likelihood we’re after is:
\ln p(R) = -\frac{N k}{2} \ln 2 \pi - \frac{1}{2} |Q \otimes I_N| - \frac{1}{2} \text{vec}(R)^T (Q \otimes I_N)^{-1} \text{vec}(R)
The determinant is immediately cheap because of kronecker rules, it’s just \frac{N}{2} \ln |Q|, so just one k \times k calculation (or easier if we have a cholesky factor lying around).
For the quadratic part, some vec and kron tricks. First, you can push the inverse inside the kron, so (Q \otimes I_N)^{-1} = (Q^{-1} \otimes I_N^{-1}) = Q^{-1} \otimes I_N. Then you have to know that (A \otimes B) \text{vec}(C) = \text{vec}(BCA^T), so:
\begin{align}
\text{vec}(R)^T(Q^{-1} \otimes I_N)\text{vec}(R) &= \text{vec}(R)^T\text{vec}(I_N R (Q^{-1})^T) \\
&= \text{vec}(R)^T \text{vec}(R Q^{-1})
\end{align}
The next fun-fact is the “Frobenius inner-product identity”, which states that \text{vec}(A)^T \text{vec}(B) = \text{trace}(A^T B) so we actually have:
\begin{align}
\text{vec}(R)^T \text{vec}(R Q^{-1}) &= \text{trace}(R^T R Q^{-1}) \\
&= \text{trace}(Q^{-1} R^T R)
\end{align}
The last step used the cyclicality of the trace operator. Anyway all this matters because now the matrix multiplication we’re left with is contracting over the big dimension. R has shape (N, k) so R^T R has shape just (k,k), which is great because now we’re solving against a tiny matrix.
So what is R^T R? Expanding back out:
\begin{align}
R^T R &= (Y - X \bar{A}^T - Z B^T)^T(Y - X \bar{A}^T - Z B^T) \\
&= (Y^T - \bar{A} X^T - B Z^T) (Y - X \bar{A}^T - Z B^T) \\
&= Y^T Y - Y^T X \bar{A}^T - Y^T Z B^T - \bar{A} X^T Y +\bar{A} X^T X \bar{A}^T + \bar{A} X^T Z B^T - B Z^T Y + B Z^T X \bar{A}^T + B Z^T Z B^T)
\end{align}
Define S_{YY} = Y^T Y, S_{XY} = X^T Y, and so on as pure functions of input data you only need to compute once and who don’t vary between MCMC draws. Also note that each S_{AB} does not include the large dimension N, only k and 2, so the matrix multiplications you do have to do on each draw stay cheap. So you compute each of the 6 S matrices just once and re-use them. The residuals are then:
R^T R = S_{YY} - S^T_{XY} \bar{A}^T - S_{YZ} B^T - \bar{A} S_{XY} + \bar{A} S_{XX} \bar{A}^T + \bar{A} S_{XZ} B^T - B S_{YZ}^T + B S_{XZ}^T \bar{A}^T + B S_{ZZ} B^T
And the final likelihood is:
\ln p(Y | X, Z, \bar{A}, B, Q) = -\frac{N k}{2} - \frac{N}{2} \ln |Q| - \frac{1}{2} \text{trace}(Q^{-1}R^T R)
This is what the StationaryVAR distribution gives you. You can also hand-code it with a potential, skipping statespace entirely, but again you will miss out on the post-estimation tools.