Classic instrumental variables regression model for Stan

classic_iv.stan

data {
  int N;
  int PX; // dimension of exogenous covariates
  int PN; // dimension of endogenous covariates
  int PZ; // dimension of instruments
  matrix[N, PX] X_exog; // exogenous covariates
  matrix[N, PN] X_endog; // engogenous covariates
  matrix[N, PZ] Z; // instruments
  vector[N] Y_outcome; // outcome variable
  int<lower=0,upper=1> run_estimation; // simulate (0) or estimate (1)
}
transformed data {
  matrix[N, 1 + PN] Y;
  Y[,1] = Y_outcome;
  Y[,2:] = X_endog;
}
parameters {
  vector[PX + PN] gamma1;
  matrix[PX + PZ, PN] gamma2;
  vector[PN + 1] alpha;
  vector<lower = 0>[1 + PN] scale;
  cholesky_factor_corr[1 + PN] L_Omega;
}
transformed parameters {
  matrix[N, 1 + PN] mu; // the conditional means of the process

  mu[:,1] = rep_vector(alpha[1], N) + append_col(X_endog,X_exog)*gamma1;
  mu[:,2:] = rep_matrix(alpha[2:]', N) + append_col(X_exog, Z)*gamma2;

}
model {
  // priors
  to_vector(gamma1) ~ normal(0, 1);
  to_vector(gamma2) ~ normal(0, 1);
  to_vector(alpha) ~ normal(0, 1);
  scale ~ cauchy(0, 2);
  L_Omega ~ lkj_corr_cholesky(4);

  // likelihood
  if(run_estimation ==1){
    for(n in 1:N) {
      Y[n] ~ multi_normal_cholesky(mu[n], diag_pre_multiply(scale, L_Omega));
    }
  }
}
generated quantities {
  matrix[N, 1 + PN] Y_simulated;

  for(n in 1:N) {
    Y_simulated[n, 1:(1+PN)] = multi_normal_cholesky_rng(mu[n]', diag_pre_multiply(scale, L_Omega))';
  }
}