Hi,
I have code for the Steady State BVAR(p) model by Mattias Villani (Villani, 2009), which is one of the main macroeconomic forecasting models used by the Swedish central bank (Ankargren, Unosson, and Yang, 2020) . Unfortunately, it takes quite a long time to estimate in Stan (compared to a Gibbs sampler) and I am wondering if it is possible to make it faster. Any ideas are very welcome! Thanks :-)
The model is
where y_t is a k-dimensional vector of time series at time t, and x_t is a q-dimensional vector of deterministic variables at time t (for example a constant and/or a dummy), and u_t \sim N_k(0,\Sigma_u) with independence between time periods. Here A_\ell for \ell=1,\dots,p is (k \times k), and \Psi is (k \times q).
Now
is the âsteady stateâ. With this parametrization of a VAR model, we can specify our priors beliefs about the long run forecasts for the variables in the system, since they converge to the steady state (unconditional mean). As Villani points out: âPrior information on the steady state is typically available, and quite often in strong form. The forecasts of inflation undertaken at central banks operating under an explicit inflation target is an apparent example.â (Villani, 2009).
We can stack the A matrices in the (kp \times k) matrix \beta
We can then rewrite the model as a nonlinear regression (Karlsson, 2013)
where w_t'=(y_{t-1}',\dots,y_{t-p}') is a kp-dimensional vector of lagged endogenous variables, and q_t'=(x_{t-1}',\dots,x_{t-p}') is a qp-dimensional vector of lagged exogenous (deterministic) variables, I_p is the (p \times p) identity matrix and \otimes denotes the Kronecker product.
Following Villani (2009), starting with \beta the usual âMinnesotaâ prior is used
Iâm leaving out the specifics here, but \Omega_\beta is a diagonal matrix. And then for the steady state coefficients, we use
where \Omega_\Psi also is a diagonal matrix. At last for \Sigma_u Jeffreys prior is used
References
Ankargren, S., Unosson, M., and Yang, Y. (2020). A Flexible Mixed-Frequency Vector Autoregression with a Steady-State Prior. Journal of Time Series Econometrics. 12(2), Article 20180034.
Karlsson, S. (2013). Forecasting with Bayesian Vector Autoregression. In: Elliott, G. and Timmerman, A. (eds) Handbook of Economic Forecasting. Elsevier B.V. Vol 2, Part B., pp. 791-897.
Villani, M. (2009). Steady-state priors for vector autoregressions.Journal of Applied Econometrics. 24(4), pp. 630-650.
The stan code is
functions {
matrix kron(matrix A, matrix B) { //kronecker product
matrix[rows(A) * rows(B), cols(A) * cols(B)] C;
int m;
int n;
int p;
int q;
m = rows(A);
n = cols(A);
p = rows(B);
q = cols(B);
for (i in 1:m) {
for (j in 1:n) {
int row_start;
int row_end;
int col_start;
int col_end;
row_start = (i - 1) * p + 1;
row_end = (i - 1) * p + p;
col_start = (j - 1) * q + 1;
col_end = (j - 1) * q + q;
C[row_start:row_end, col_start:col_end] = A[i, j] * B;
}
}
return C;
}
}
data {
int<lower=0> N; //number of observations
int<lower=0> k; //number of variables
int<lower=0> p; //lag order
int<lower=0> q; //number of exogenous variables
matrix[N, k] Y; //endogenous variables (y_t)
matrix[N, q] X; //exogenous variables (x_t)
matrix[N, k*p] W; //lagged endogenous variables
matrix[N, q*p] Q; //lagged exogenous variables
vector[k*p*k] theta_beta; //vec_beta prior mean
matrix[k*p*k, k*p*k] Omega_beta; //vec_beta prior covariance matrix
vector[k*q] theta_Psi; //vec_Psi prior mean
matrix[k*q, k*q] Omega_Psi; //vec_Psi prior covariance matrix
int<lower=0> H; // Forecast horizon
matrix[H, q] X_pred; //future exogenous variables
}
transformed data {
matrix[p, p] I_p = diag_matrix(rep_vector(1, p)); // Identity matrix
}
parameters {
matrix[k*p, k] beta; //beta' = (A_1,...,A_p)
matrix[k, q] Psi; //Psi * x_t = steady state
cov_matrix[k] Sigma_u;
}
model {
for(t in 1:N){
vector[k] u_t = (Y[t] - (X[t]*Psi' + (W[t]-Q[t]*(kron(I_p,Psi')))*beta))';
u_t ~ multi_normal(rep_vector(0,k), Sigma_u);
}
to_vector(beta) ~ multi_normal(theta_beta, Omega_beta);
to_vector(Psi) ~ multi_normal(theta_Psi, Omega_Psi);
target += -0.5 * (k + 1) * log_determinant(Sigma_u);
}
generated quantities {
matrix[k, k] A[p];
for (i in 1:p) {
A[i] = (beta[((i - 1) * k + 1):(i * k), :])'; //extract A_1, ..., A_p
}
matrix[H, k] Y_pred;
for (h in 1:H) {
vector[k] u_t = multi_normal_rng(rep_vector(0, k), Sigma_u);
vector[k] yhat_t = (X_pred[h]*Psi')';
if (h > 1) {
for (i in 1:min(h-1, p)) {
yhat_t += to_vector((Y_pred[h-i] - X_pred[h-i]*Psi') * A[i]');
}
}
if (h <= p) {
for (i in h:p) {
yhat_t += to_vector((Y[N + h - i] - X[N + h - i]*Psi') * A[i]');
}
}
Y_pred[h] = (yhat_t + u_t)'; //forecasts
}
}