halloween_vassavage.stan

functions {
  matrix make_pairwise_logit_mat(vector theta) {
    matrix[rows(theta), rows(theta)] out;
    
    for(r in 1:rows(out)) {
      for(c in 1:r) {
        if(r==c) {
          out[r, c] = 0;
        } else {
          out[r, c] = exp(theta[c])/(exp(theta[c]) + exp(theta[r]));
          out[c, r] = 1.0 - out[r, c];
        }
      }
    }
    return(out);
  }
  vector get_probs(vector delta, matrix x, matrix eta, matrix L) {
    matrix[rows(eta), rows(delta)] utils;
    matrix[rows(eta), rows(delta)] probs;
    vector[rows(delta)] shares;
    real scaler;
    scaler = pow(cols(probs), -1);
    
    for(i in 1:rows(eta)) {
      matrix[rows(delta), rows(delta)] tmpmat;
      utils[i] = delta' + eta[i] * L * x';
      tmpmat = make_pairwise_logit_mat(utils[i]');
      
      for(j in 1:cols(tmpmat)) {
       probs[i, j] = sum(tmpmat[:,j])/(rows(delta) - 1.0);
      }
    }
    for(i in 1:cols(probs)) {
      shares[i] = mean(probs[:,i]);
    }
    // makes correction for multinomial shares
    return(shares);
  }
  vector get_shares(vector delta, matrix x, matrix eta, matrix L) {
    matrix[rows(eta), rows(delta)] utils;
    matrix[rows(eta), rows(delta)] probs;
    vector[rows(delta)] shares;
    real scaler;
    scaler = pow(cols(probs), -1);
    
    for(i in 1:rows(eta)) {
      utils[i] = delta' + eta[i] * L * x';
      probs[i] = softmax(utils[i]')';
    }
    for(i in 1:cols(probs)) {
      shares[i] = mean(probs[:,i]);
    }
    // makes correction for multinomial shares
    return(shares);
  }
  
}
data {
  int N;
  int Npair;
  int P;
  int NS;
  matrix[NS, P] eta;
  matrix[N, P] X;
  vector[N] price;
  int sales[N];
  int run_estimation;
}
transformed data {
  matrix[N, P-1] X2;
  X2 = X[:,2:];
}
parameters {
  vector[P] beta;
  corr_matrix[P] Omega;
  vector<lower = 0>[P] tau; 
  vector[N] xi;
  real alpha;
  vector[P-1] beta_price;
  real<lower = 0> lambda;
  real<lower = 0> sigma_price;
}
transformed parameters {
  cholesky_factor_cov[P] L;
  L = cholesky_decompose(quad_form_diag(Omega, tau));
}
model {
  vector[N] theta;
  
  // priors
  beta ~ normal(0, 1);
  Omega ~ lkj_corr(3);
  tau ~ normal(0, 1);
  xi ~ normal(0, 1);
  sigma_price ~ normal(0,1);
  lambda ~ normal(0, 1);
  alpha ~ normal(0, 1);
  
  // likelihood
  theta = get_probs(X*beta + xi, X, eta, L);
  
  if(run_estimation==1) {
    price ~ normal(alpha + X2*beta_price+ lambda * xi, sigma_price);
    sales ~ binomial(Npair, theta);
  }
  
}
generated quantities {
  vector[N] shares;
  vector[N] winprob;
  winprob = get_probs(X*beta + xi, X, eta, L);
  shares = get_shares(X*beta + xi, X, eta, L);
}