Normal distribution with multiple variance components

I’m trying to model a Generalised Least Squares panel data model with a diagonal variance-covariance matrix, i.e. y_i|X_i,\beta,\Sigma \sim N(X_i\beta,\Sigma),\ \Sigma=diag(\sigma^2_1,\dots,\sigma^2_T). I’ve been able to model a Panel Data model with a common variance, i.e. y_i|X_i,\beta,\sigma^2_\varepsilon\sim N(X_i\beta,\sigma^2_\varepsilon I_T), but I fail with multiple components. I have, however, been able to write the code for an arbitrary \Sigma, but that code failed due to memory limitations (N=2376) and I know that the matrix is diagonal anyway. Could anyone give a hint on how to make this code work for multiple variance components?

data {
int<lower=0> N; \\ Number of cross-sectional units
int<lower=0> K; \\ Number of regressors
real bval[N]; \\ y-value
real btreat[N]; \\ Regressor 1
real bcountrd[N]; \\ Regressor 2
}

parameters {
real beta[K];
real sigma[N];
}

transformed parameters {
real mu[N];
for (i in 1:N) {
mu[i] = beta[1]+beta[2]*btreat[i]+beta[3]*bcountrd[i];
}
}

model {
// Likelihood
bval ~ normal(mu,sigma);
}


I think you’ve coded the model you intended. If you move mu to the model block it’ll use less (there’s less output) but you still have O(N) sigmas:

model {
real mu[N];
for (i in 1:N) {
mu[i] = beta[1]+beta[2]*btreat[i]+beta[3]*bcountrd[i];
}

bval ~ normal(mu,sigma);
}


I think it’s really unlikely this model fits well though honestly. There’s no structure on the sigmas. Do you have multiple samples of each y?

I’d recommend converting those real arrays to vectors and using the vectorized one-liner:

vector[N] mu = beta[1] + beta[2] * btreat + beta[3] * bcountrd;


It’s also more efficient because it reduces the autodiff expression graph.

But as @bbbales2 points out, you can’t let sigma vary by observation without
structure.