Probably a @bgoodri or @anon75146577 question (motivated by previous conversations on the topic with both)
We have some time-series process
x[t]' ~ multi_normal(B * x[t-1]', Sigma);
This process will be stationary – ie x[t] -> 0 as t gets large if the eigenvalues of B are inside the unit circle. This is an attractive property where we don’t want to forecast that the world explodes. Dan mentioned to me a while ago that there was a neat reparameterization that guarantees this property of B
B = matrix_exp(-A);
where A is positive definite. I know I can construct some positive definite matrices from their decompositions (for example a covariance matrix). But I’m wondering if there are more general decompositions of positive definite matrices that would allow me to build them up from parts like we would a covariance matrix?
Yeah, we know how to do VAR(1) in general, but no one has actually sat down to do it yet in Stan. I wrote it down a few years ago, but it is on a different laptop than the one in my lap. It would probably be easier to explain in person anyway (i.e. explaining it is hard), but basically you compose B from its “scaled-SVD” as \mathbf{B} = \rho\mathbf{U}\mathbf{D}\mathbf{V}^\top where \rho \in \left[-1,1\right] is the largest eigenvalue of \mathbf{B}, \mathbf{U} and \mathbf{V} are rotation matrices, and D is a diagonal matrix such that the largest eigenvalue (by magnitude) of \mathbf{U}\mathbf{D}\mathbf{V}^\top is 1.
The hard part is the diagonal of \mathbf{D}. If we followed the usual convention of restricting the diagonal elements of \mathbf{D} to be positive, then in Stan we would have a nasty problem that columns of \mathbf{U} and \mathbf{V} could be multiplied by -1 and still yield an admissible \mathbf{B}. Since \mathbf{U} and \mathbf{V} are asserted to be rotation matrices — i.e. their determinants are each 1 rather than -1 — we have to allow the diagonal elements of \mathbf{D} be negative.
Even still, you have to do iteration in order to obtain a \mathbf{D} such that the largest eigenvalue (by magnitude) of \mathbf{U}\mathbf{D}\mathbf{V}^\top is 1. That is doable if you have a starting value for \mathbf{D} such that the largest implied eigenvalue is less than 1 in magnitude. If the starting vector for the diagonal of \mathbf{D} is a (sorted) unit-vector, then it satisfies the strict inequality on the largest eigenvalue and then you can iterate until the largest eigenvalue is 1 before down-scaling the whole thing by \rho.
Then, to have a jointly uniform density over the admissible values that \mathbf{B} could take, you have to work through the Jacobian determinant of the SVD transformation, which is not as difficult as it sounds. Basically, it only involves the squares of the diagonal elements of \rho \mathbf{D} so our unconventional convention of allowing them to be negative does not change anything with the Jacobian. Also, you have to be able to construct a \mathbf{U} and \mathbf{V} that are uniformly distributed over the space of rotation matrices, but there is another thread somewhere that uses a Cayley transform to do that (with mixed results because the tails are heavy in the unconstrained space, but conditional on the data I think it would be okay).
It’s fairly complex and requires a pretty deep understanding of what stationary really means as the dimensions and lags increase. The arxiv page has the stan code!
You can parameterize the model coefficients via the reflection coefficients, there is a recursion (the levinson recursion or step up recursion, but there may be other names) closely related to the Cholesky factorization that will guarantee the resulting model is stable. I have been experimenting with this as well if it helps (here: https://github.com/RJTK/stanlearn/blob/master/stanlearn/time_series/stan_models/ar_functions.stan) but only so far for the AR case – but the recursion generalizes to the VAR case by including the forwards and backwards recursion. I can also draw attention to the now dated work of Whittle on this problem since it’s not referenced in the work above:
@article{whittle1963fitting,
title={On the fitting of multivariate autoregressions, and the approximate canonical factorization of a spectral density matrix},
author={Whittle, Peter},
journal={Biometrika},
volume={50},
number={1-2},
pages={129--134},
year={1963},
publisher={Oxford University Press}
}
but the Sarah Heaps paper pointed to above looks like exactly what you’re after, and will be useful for myself as well :)