I’m trying to use reduce_sum with a multilevel VAR, and I’m having trouble trying to figure out how to wrangle my data structure to work with reduce_sum (or if it’s even possible).

My units are regions, nested within countries. The model partially pools at the country level (country-level random intercepts and slopes). The data are structured as a 11x5x120 dimension array, where each 11x5 matrix is a region.

I’m running into trouble on how to define the slice, given that it needs to slice an array and each region’s time series must be kept separate (for the lags to make sense). (This is my first time using reduce_sum, so I may also be misunderstanding something about how this works).

This is the code I have so far:

```
functions{
real partial_sum(real[] slice_y,
int start, int end,
matrix beta,
matrix alpha,
matrix Omega,
vector tau,
int t,
int country,
int R,
int K
)
for(r in 1:R){
matrix[t, K] y_temp = to_matrix(slice_y[:, :, r]);
for(tt in 1:t){
vector[K] lagvars; //lagged variables
//linear predictor
if(tt > 1){
lagvars = to_vector(y_temp[tt-1,:]);
real mu[K]; //linear predictor
mu = alpha[country[r]] + beta[country[r]]*lagvars;
return multi_normal_lpdf(y[t]|mu,quad_form_diag(Omega[country[r]], tau[country[r]]))
}
}
}
}
data{
int N; //# observations
int K; //# dimensions of Y
int C; //# of countries
int R; //# of regions
int t; //# of time periods in the panel
int<lower = 1, upper=C> country[R]; //country id for each region
int<lower = 1, upper=t> time[N]; //time period id for each obs
int<lower = 1, upper=R> region[N]; //region id for each obs
real Y[t, K, R]; //the outcome array - each variable's time series stacked by region
}
parameters{
//individual level
vector<lower = 0>[K] tau[C]; //scale for residuals
matrix[K, K] z_beta[C]; //untransformed betas
vector[K] z_alpha[C]; //untransformed intercepts
//hierarchical parameters
corr_matrix[K] Omega[C]; //country level correlation matrix
vector[K] tau_loc; //mean for variance scaling factor
vector<lower=0>[K] tau_scale; //scale for tau
matrix[K, K] bhat_location; //mean for prior on beta
matrix<lower = 0>[K, K] bhat_scale; //scale for prior on beta
vector[K] ahat_location; //means for prior on intercepts
vector<lower = 0>[K] ahat_scale; //variance for intercept prior
}
transformed parameters{
matrix[K, K] beta[C]; //VAR(1) coefficients, country specific
vector[K] alpha[C]; //country specific intercepts
for(c in 1:C){
//recentering random effects
alpha[c] = ahat_location + ahat_scale .*z_alpha[c];
beta[c] = bhat_location + bhat_scale*z_beta[c];
}
}
model{
//hyperpriors
tau_loc ~ cauchy(0,1);
tau_scale ~ cauchy(0,1);
ahat_location ~ normal(0,1);
ahat_scale ~ cauchy(0, 1);
to_vector(bhat_location) ~ normal(0, 0.5);
to_vector(bhat_scale) ~ cauchy(0, 0.5);
//hierarchical priors
for(c in 1:C){
//non-centered paramaterization to avoid convergence issues
z_alpha[c] ~ normal(0, 1);
to_vector(z_beta[c]) ~ normal(0, 1);
tau[c] ~ normal(tau_loc, tau_scale);
Omega[c] ~ lkj_corr(5);
}
//reduce sum likelihood
target += reduce_sum(partial_sum, Y, grainsize, beta, alpha,
Omega, tau, t, country, R, K)
}
```

Where I run into trouble is the `matrix[t, K] y_temp = to_matrix(slice_y[:, :, r])`

line. In the original model, this was designed to subset a particular region out of the array. I believe this is where the start and end integers need to go, but I’m not sure how to incorporate that while still respecting the multi-level nature of the data.

I’ve looked at this post, but the nature and structure of their model is different enough from mine (and the code is complex enough) that I haven’t been able to find a clear way to alter it to work with my model.

How can I get this model running with sum_reduce, if possible?

Thanks in advance for any help!