Hello Stan users,

I’m trying to write Stan code that is both a) a categorical logit and b) is a ‘structural model’ (economics terminology). I can find useful tutorials/code on either the former (e.g. https://github.com/ksvanhorn/ART-Forum-2017-Stan-Tutorial) or the latter (e.g. https://github.com/CCS-Lab/hBayesDM/blob/master/exec/ra_prospect.stan) but not the combination, and have got stuck in adapting the code for my purpose.

My data comes from 300 people choosing one of 11 options, in two rounds. It is a risky choice, with a simplified utility function below. (I’m aware ‘normal’ categorical logits need K-1 for identification, I don’t think that is true here.)

Previously, I calculated the utility of each option and the code took ~3 hours to run. It appears vectorising is a sensible option to improve speed, but I’ve come up against the problem that you can’t use exponents on vectors. [error: No matches for: pow(vector, real)]

The answer appears to be looping over choices in the model block, but I’m getting confused between the different indices, as R respondents make a choice between C alternatives in S rounds. Furthermore, this seems to defeat the point of vectorising. If I *should* loop over choices to use the exponential, how do I do this? If I shouldn’t, or a better set up would be optimal, what combination of vector/array/indexing options seems best? Parameters will either be the same for everyone, or have some distribution over respondents (but constant across scenarios and alternatives).

Current code below

```
data {
int<lower=2> C; // choices/alternatives per scenario=11
int<lower=2> R; // respondents=300ish
int<lower=2> S; // scenarios=2
int<lower=1> N; // observations=R*S
int<lower=1,upper=C> Y[N]; // number of chosen alternative
vector[C] gain[N];
vector[C] loss[N];
vector[C] EV[N];
}
parameters {
real<lower=0> lambda[N];
real mean_lambda;
real<lower=0> sigma_lambda;
real<lower=1> r;
}
model {
mean_lambda ~ cauchy(5, 1);
sigma_lambda ~ cauchy(0, 1);
lambda ~ normal(mean_lambda, sigma_lambda);
r ~ cauchy(1,1);
for (i in 1:N){
Y[i] ~ categorical_logit(EV[i] + gain[i] - lambda[i]*pow(loss[i],r));
}}
```