Skip to content

Instantly share code, notes, and snippets.

@khakieconomics
Created October 13, 2017 19:33
Show Gist options
  • Save khakieconomics/e6e34469b5ac47bbf0ed1f5b22a16ecb to your computer and use it in GitHub Desktop.
Save khakieconomics/e6e34469b5ac47bbf0ed1f5b22a16ecb to your computer and use it in GitHub Desktop.
Classic instrumental variables regression model for 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))';
}
}
@cjgb
Copy link

cjgb commented Sep 7, 2020

I am trying to make sense of this code and there is something I do not understand: shouldn't line 27 be something like

mu[:,1] = rep_vector(alpha[1], N) + append_col(mu[:,2:],X_exog)*gamma1;

You want your estimation of the exogenous variables (and not the exogenous variables themselves) into the estimation of the main effect, don't you?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment