Normal distribution with multiple variance components

Modeling
,
1

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);   
}

Thanks in advance.

2

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?

3

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.