The following is a simplistic multilevel model to estimate effects of locations and products on transaction distributions. Results are good (using parameters to simulate and test actual transactions), but it takes a LONG time to run. Hoping for suggestions on how to improve approach and performance. Thanks:

```
data {
int N; // number of observations
vector[N] y; // observations
int J; // number of groups- Prod
int PGrp[N]; // group indexes
int K; // number of grps - Location
int LGrp[N]; // Location level grp indexes
real p_Mu;
real p_sigma;
real p_alpha;
vector[J] p_PGrpMu;
vector<lower=0>[J] p_sigma_PGrp;
vector[J] p_alpha_PGrp;
vector[K] p_LGrpMu;
vector<lower=0>[K] p_sigma_LGrp;
vector[K] p_alpha_LGrp;
real muSigma;
real sigmaSigma;
real alphaSigma;
real pGrpMuSigma;
real pGrpSigmaSigma;
real pGrpAlphaSigma;
real lGrpMuSigma;
real lGrpSigmaSigma;
real lGrpAlphaSigma;
}
parameters {
real mu;
real<lower=0> sigma;
real alpha;
vector[J] PGrpMu;
vector<lower=0>[J] sigma_PGrp;
vector[J] alpha_PGrp;
vector[K] LGrpMu;
vector<lower=0>[K] sigma_LGrp;
vector[K] alpha_LGrp;
}
model {
target += cauchy_lpdf(mu | p_Mu, muSigma);
target += cauchy_lpdf(PGrpMu | p_PGrpMu, pGrpMuSigma);
target += cauchy_lpdf(LGrpMu | p_LGrpMu, lGrpMuSigma);
target += cauchy_lpdf(sigma | p_sigma, sigmaSigma);
target += cauchy_lpdf(sigma_PGrp | p_sigma_PGrp, pGrpSigmaSigma);
target += cauchy_lpdf(sigma_LGrp | p_sigma_LGrp, lGrpSigmaSigma);
target += cauchy_lpdf(alpha | p_alpha, alphaSigma);
target += cauchy_lpdf(alpha_PGrp | p_alpha_PGrp, pGrpAlphaSigma);
target += cauchy_lpdf(alpha_LGrp | p_alpha_LGrp, lGrpAlphaSigma);
target += skew_normal_lpdf(y | mu+PGrpMu[PGrp]+LGrpMu[LGrp], sigma+sigma_PGrp[PGrp]+sigma_LGrp[LGrp], alpha+alpha_PGrp[PGrp]+alpha_LGrp[LGrp] );
}
```