// stan_code/var_1_pois.stan // a poisson var model for count data generated by a VAR poisson process data{ int forecast_len; // length of the forecast to generate int N; // length of the time series int K; // number of series // array notation is [rows, columns] int Y[K,N]; // output matrix } parameters{ vector[K] tau; // LJK parameterize the covariance matrix cholesky_factor_corr[K] L_Omega; // covariance matrix in cholesky form matrix[K,K] beta; vector[K] alpha; // latent poisson parameters vector[K] lambda[N]; } transformed parameters{ matrix[K,K] L_Sigma; L_Sigma = diag_pre_multiply(tau, L_Omega); } model{ vector[K] mu[N-1]; for (n in 2:N){ mu[n-1] = (alpha + beta*lambda[n-1]); } tau ~ cauchy(0,10); lambda[1] ~ cauchy(0,8); lambda[2:N] ~ multi_normal_cholesky(mu, L_Sigma); L_Omega ~ lkj_corr_cholesky(1.0); // prior on Omega to_vector(beta) ~ normal(0,0.8); alpha ~ cauchy(0,10); // intercepts can be all over // try this with map_rect, but maybe it can be better vectorized too for(k in 1:K) Y[k] ~ poisson_log(lambda[1:N,k]); } generated quantities{ // stan doesn't have a 64 bit rng, so let's just forecast lambdas // and then we'll generate y in R/python downstream. vector[K] lambda_new[forecast_len]; cov_matrix[K] Sigma; { // just so we output the covariance matrix Sigma = multiply_lower_tri_self_transpose(L_Sigma); lambda_new[1] = multi_normal_cholesky_rng(alpha + beta*lambda[N], L_Sigma); for (n in 2:forecast_len){ lambda_new[n] = multi_normal_cholesky_rng(alpha + beta*lambda_new[n-1], L_Sigma); } } }