#delimit ;
set more off;
*log close;
clear;
set memory 20m;
use nhis_mcod, replace;               * load data;
log using aproblemset6.log, replace;  * write results to log;

********** question 1;
xi i.educ i.income i.race;
dprobit diedin5 _I* age married bmi;


********** question 2;
probit diedin5 _I* age married bmi;
prchange, x(_Ieduc_2=1 _Ieduc_3=0 _Iincome_2=0 _Iincome_3=1 _Iincome_4=0
_Iincome_5=0 _Iincome_6=0 _Irace_2=0 _Irace_3=0 _Irace_4=0 age1=30 married=1 bmi=28);

prchange, x(_Ieduc_2=1 _Ieduc_3=0 _Iincome_2=0 _Iincome_3=1 _Iincome_4=0
_Iincome_5=0 _Iincome_6=0 _Irace_2=0 _Irace_3=0 _Irace_4=0 age1=60 married=1 bmi=28);


********** question 3;
gen underweight=bmi<20;
gen obese=bmi>=30;
dprobit diedin5 _I* age married underweight obese;   * unrestricted model;
test obese underweight;
dprobit diedin5 _I* age married;                     * restricted model;


********** question 4;
logistic diedin5 _I* age married underweight obese;   * unrestricted model;
sum obese;

* this subroutine generates the marginal effects for a probit model;
* the notation follows greene so beta is (k x 1) and xbar is (k x 1);
* so beta1`xbar is a scalar;

logit diedin5 married age1 bmi;
mfx compute;
matrix betat=e(b);    * get beta from probit (1 x k);
matrix beta=betat';
matrix covp=e(V);     * get v/c matric from probit (k x k);


* get means of x -- call it xbar (k x 1);
* must be the same order as in the probit statement;
matrix accum zz = married age1 bmi, mean(xbart);

matrix xbar=xbart';                 * transpose beta; 
matrix xbeta=beta'*xbar;            * get xbeta (scalar);
matrix prob=exp(xbeta[1,1])/(1+exp(xbeta[1,1]));   
matrix pdf=prob*(1-prob);
matrix k=rowsof(beta);              * get number of covariates;
matrix Ik=I(k[1,1]);                * construct I(k);
matrix scale=(1-exp(xbeta[1,1]))/(1+exp(xbeta[1,1]));
matrix G=Ik+scale*beta*xbar';             * construct G;
matrix v_c=(pdf*pdf)*G*covp*G';     * get v-c matrix of marginal effects;
matrix me= beta*pdf;                * get marginal effects;
matrix se_me1=cholesky(diag(vecdiag(v_c)));  * get square root of main diag;
matrix se_me=vecdiag(se_me1)';      *take diagonal values;
matrix z_score=vecdiag(diag(me)*inv(diag(se_me)))';  * get z score;
matrix results=me,se_me,z_score;    *  construct results matrix;
matrix colnames results=marg_eff std_err z_score;     * define column names;
matrix list results;                * list results;


* this is an example of a marginal effect for a dichotomous outcome;
* in this case, set the 1st variable worka as 1 or 0;
matrix x1=xbar;
matrix x1[1,1]=1;
matrix x0=xbar;
matrix x0[1,1]=0;
matrix xbeta1=beta'*x1;
matrix xbeta0=beta'*x0;
matrix prob1=exp(xbeta1[1,1])/(1+exp(xbeta1[1,1]));
matrix prob0=exp(xbeta0[1,1])/(1+exp(xbeta0[1,1]));
matrix me_1=prob1-prob0;
matrix pdf1=prob1/(1+exp(xbeta1[1,1]));
matrix pdf0=prob0/(1+exp(xbeta0[1,1]));
matrix G1=pdf1*x1 - pdf0*x0;
matrix v_c1=G1'*covp*G1;
matrix se_me_1=sqrt(v_c1[1,1]);
* marginal effect of workplace bans;
matrix list me_1;
* standard error of workplace a;
matrix list se_me_1;