/*Calculate Marginal Effects and their standard deviations
for bivariate probit with sample selection.
By Mark Schreiner.
Last Update: 2/6/96
*/
/* The order of regressors and parameters is as follows:
First regressors unique to the first (selected sample) equation.
The intercept is included as a unique variable so that even if nothing
else is unique there will still be one explanatory unique to eq. 1.
Second, regressors in first equation that are also in second.
Third, regressors unique to the second equation (include intercept).
Fourth, regressors in the second equation that are also in the first.
Notation:
kx=Number of regressors in a group.
px=Position of sub-vector within a larger vector.
bx=Parameters for group x.
xx=Regressors for group x.
x1=Regressors/parameters unique to equation 1.
x2=Resgressors/parameters unique to equation 2.
x1_=Parameters on common regressors in equation 1.
x2_=Parameters on common regressors in equation 2.
x_=Regressors that appear in both equations.
n=Number of observations.
k=Number of parameters, including rho.
i=index over observations.
xi=index over regressors.
bi=index over parameters.
x=regressor matrix.
b=parameter vector.
covb=estimated covaraiance matrix of parameters.
*/
/*User must modify the hard-code below*/
k1 =1;
k1_=6;
k2 =2;
k2_=k1_;
k =k1+k1_+k2+k2_+1;
k_ =k1_;
kx =k1+k2+k_;
stop;
end;
/*px indicates positions of sub-vectors within vectors*/
p1 =1;
p1_=p1+k1;
p2 =p1+k1+k_;
p2_=p1+k1+k_+k2;
/*Set up vector of estimated coefficients*/
/*These are from Feb 6 1996, 5:28 pm LIMDEP output*/
load b[] = c:\graham\gambia\b;
/*q. Name Type k p */
/*Eq. 1 Constant Unique 1 1 */
/*Eq. 1 a26_25 Common 6 2 */
/*Eq. 1 a36_plus Common 6 2 */
/*Eq. 1 Noread Common 6 2 */
/*Eq. 1 Koranic Common 6 2 */
/*Eq. 1 Female Common 6 2 */
/*Eq. 1 HHSize Common 6 2 */
/*Eq. 2 Constant Unique 2 8 */
/*Eq. 2 Informal Unique 2 8 */
/*Eq. 2 a26_25 Common 6 10 */
/*Eq. 2 a36_plus Common 6 10 */
/*Eq. 2 Noread Common 6 10 */
/*Eq. 2 Koranic Common 6 10 */
/*Eq. 2 Female Common 6 10 */
/*Eq. 2 HHSize Common 6 10 */
/*Eq. 2 Rho 16 */
/*Set up beta subvectors*/
print b;
rho=b(k);
b1 =b[p1 :p1 +k1-1,1];
b1_=b[p1_:p1_+k_-1,1];
b2 =b[p2 :p2 +k2-1,1];
b2_=b[p2_:p2_+k_-1,1];
/*Read in and set up data matrix in form x1 x* x2*/
load xin[] = c:\graham\gambia\data;
n=rows(xin);
x=ones(n,kx)*0;
/*Adjust hard-code below so that x has
first unique regressors for equation 1,
second common regressors for both equations, and
third unique regressors for equation 2. */
x[.,p1]=ones(n,1); /*Constant eq 1 unique*/
x[.,p2]=ones(n,1); /*Constant eq 2 unqiue*/
x[.,p2+1]=xin[.,12]; /*Informal eq 2 unique*/
x[.,p1_] =xin[.,4]; /*a26_35 common*/
x[.,p1_+1]=xin[.,5]; /*a36_plus common*/
x[.,p1_+2]=xin[.,6]; /*Noread common*/
x[.,p1_+3]=xin[.,7]; /*Koranic common*/
x[.,p1_+4]=xin[.,15]; /*Female common*/
x[.,p1_+5]=xin[.,17]; /*HHSize common*/
/*Read in the covariance matrix of the estimated parameters*/
load covb[] = c:\graham\gambia\covb;
/*dummy[x,1]=1 if the xth regressor is a dummy variable*/
/*When a regressor that is a dummy appears alone in a derivitive,
I must make sure it equals one or the derivitive vanishes*/
dummy=ones(k,1)*0;
/*Adjust hard-code below*/
dummy[2:6,1]=ones(5,1);
dummy[10:14,1]=ones(5,1);
/*Initialize accumulators of marginal effects and standard errors*/
effects1=ones(1,k)*0;
effects2=effects1;
sd1 =ones(1,k)*0;
sd2 =sd1;
/*Loop over observations and set up intermediate values*/
i=0;
do while i<=n;
/*Set up regressor subvectors*/
x1 =x[i,p1 :p1 +k1-1];
x1_=x[i,p1_:p1_+k_-1];
x2 =x[i,p2 :p2 +k2-1];
w1=x1*b1+x1_*b1_;
w2=x2*b2+x2_*b2_;
delta=1/(sqrt(1-rho^2));
v1=delta*(w2-rho*w1);
v2=delta*(w1-rho*w2);
a=(1/2)*(1/PI)*delta;
bb=(1/2)*(delta^2)*(-1)*(w1^2+w2^2-2*rho*w1*w1);
phiw1=pdfn(w1);
phiw2=pdfn(w2);
phiv1=pdfn(v1);
phiv2=pdfn(v2);
phi2=a*exp(bb);
bigphiv1=cdfn(v1);
bigphiv2=cdfn(v2);
/*Now loop over the regressors*/
xi=1;
/*Addition to sum of marginal effects of standard errors for each obs*/
dm1=ones(1,k)*0;
dm2=ones(1,k)*0;
dsd1=ones(1,k)*0;
dsd2=ones(1,k)*0;
/*Process marginal effects for the first equation first*/
do while xi<=kx;
dj=ones(k,kx)*0;
if (xi>=p1) and (xi<pl_) then do;
dm1[1,xi]=b1[xi-p1+1,1]*phiw1*bigphiv1;
/*Change in x1*/
bi=1;
xj=xi-p1+1;
do while bi<=k;
if (bi>=p1) and (bi<p1_) then do;
/*Change in marginal effect of x1 with respect to b1*/
bj=bi-p1+1;
if dummy[bi,1]=1 then do;
xx=1;
else;
xx=x1[1,bj];
endif;
dj[bi,xi]=phiw1*(bigphiv1*(1-w1*b1[xj,1]*xx)
+b1[xj,1]*xx*phiv1);
endif;
if (bi>=p1_) and (bi<p2) then do;
/*Change in effect of x1 with respect to b1_*/
bj=bi-p1_+1;
if dummy[bi,1]=1 then do;
xx=1;
else;
xx=x_[1,bj];
endif;
dj[bi,xi]=b1[xj,1]*xx*phiw1*(phiv1-w1*bigphiv1);
endif;
if (bi>=p2) and (bi<p2_) then do;
/*Change in effect of x1 with respect to b2*/
bj=bi-p2+1;
if dummy[bi,1]=1 then do;
xx=1;
else;
xx=x2[1,bj];
endif;
dj[bi,xi]=b1[xj,1]*xx*phiw1*phiv1;
endif;
if (bi>=p2_) and (bi<k) then do;
/*Change in effect of x1 with respect to b2_*/
bj=bi-p2_+1;
if dummy[bi,1]=1 then do;
xx=1;
else;
xx=x_[1,bj];
endif;
dj[bi,xi]=b1[xj,1]*xx*phiw1*phiv1;
endif;
if (bi=k) then do;
/*Change in effect of x1 with respect to rho*/
dj[bi,xi]=b1[xj,1]*phiw1*phiv1*delta*
(delta*rho*v1-w1);
endif;
bi=bi+1;
endo;
endif;
if (xi>=p1_) and (xi<p2) then do;
dm1[1,xi]=b2_[xi-p1_+1,1]*b1_[xi-p1_+1,1]*phi2;
/*Change in x_*/
bi=1;
xj=xi-p1_+1;
do while bi<=k;
if (bi>=p1) and (bi<p1_) then do;
/*Change in marginal effect of x_ with respect to b1*/
bj=bi-p1+1;
if dummy[bi,1]=1 then do;
xx=1;
else;
xx=x1[1,bj];
endif;
dj[bi,xi]=b1_[xj,1]*b2_[xj,1]*xx*2*(rho*w2-w1)*phi2;
endif;
if (bi>=p1_) and (bi<p2) then do;
/*Change in effect of x_ with respect to b1_*/
bj=bi-p1_+1;
if dummy[bi,1]=1 then do;
xx=1;
else;
xx=x_[1,bj];
endif;
dj[bi,xi]=b2_[xj,1]*phi2*(1+b1_[xj,1]*xx*2*(rho*w2-w1));
endif;
if (bi>=p2) and (bi<p2_) then do;
/*Change in effect of x_ with respect to b2*/
bj=bi-p2+1;
if dummy[bi,1]=1 then do;
xx=1;
else;
xx=x2[1,bj];
endif;
dj[bi,xi]=b1_[xj,1]*b2_[xj,1]*xx*2*phi2*(rho*w1-w2);
endif;
if (bi>=p2_) and (bi<k) then do;
/*Change in effect of x_ with respect to b2_*/
bj=bi-p2_+1;
if dummy[bi,1]=1 then do;
xx=1;
else;
xx=x_[1,bj];
endif;
dj[bi,xi]=b1_[xj,1]*phi2*(1+b2_[xj,1]*xx*2*(rho*w1-w2));
endif;
if (bi=k) then do;
/*Change in effect of x_ with respect to rho*/
dj[bi,xi]=b1_[xj,1]*b2_[xj,1]*phi2*(delta^2)*
(rho+w1*w2+2*(delta^2)*rho*bb);
endif;
bi=bi+1;
endo;
endif;
if (xi>=p2) and (xi<k) then do;
dm1[1,xi]=b2[xi-p2+1,1]*phiw2*bigphiv2;
/*Change in x2*/
bi=1;
xj=xi-p2+1;
do while bi<=k;
if (bi>=p1) and (bi<p1_) then do;
/*Change in marginal effect of x2 with respect to b1*/
bj=bi-p1+1;
if dummy[bi,1]=1 then do;
xx=1;
else;
xx=x1[1,bj];
endif;
dj[bi,xi]=b2[xj,1]*xx*phiw2*phiv2;
endif;
if (bi>=p1_) and (bi<p2) then do;
/*Change in effect of x2 with respect to b1_*/
bj=bi-p1_+1;
if dummy[bi,1]=1 then do;
xx=1;
else;
xx=x_[1,bj];
endif;
dj[bi,xi]=b2_[xj,1]*xx*phiw2*phiv2;
endif;
if (bi>=p2) and (bi<p2_) then do;
/*Change in effect of x2 with respect to b2*/
bj=bi-p2+1;
if dummy[bi,1]=1 then do;
xx=1;
else;
xx=x2[1,bj];
endif;
dj[bi,xi]=phiw2*(bigphiv2*(1-b2[xj,1]*xx*w2)+
b2[xj,1]*xx*phiv2);
endif;
if (bi>=p2_) and (bi<k) then do;
/*Change in effect of x2 with respect to b2_*/
bj=bi-p2_+1;
if dummy[bi,1]=1 then do;
xx=1;
else;
xx=x_[1,bj];
endif;
dj[bi,xi]=b2_[xj,1]*xx*phiw2*(phiv2-w2*bigphiv2);
endif;
if (bi=k) then do;
/*Change in effect of x2 with respect to rho*/
dj[bi,xi]=b2_[xj,1]*phiw2*phiv2*delta)*
(delta*rho*v2-w2);
endif;
bi=bi+1;
endo;
endif;
xi=xi+1;
endo;
/*Accumulate marginal effects*/
effects1=effects1+dm1;
/*Calculate the variance matrix of marginal effects*/
vareff=dj'*covb*dj;
sd1=sd1+sqrt(trace(vareff));
/*Now calculate the effects for the second equation*/
xi=p1_;
do while xi<=kx;
dj=ones(k,kx)*0;
if (xi>=p2) and (xi<k) then do;
dm2[1,xi]=b2[xi-p2+1,1]*phiw2;
/*Change in x2*/
bi=p2;
xj=xi-p2+1;
do while bi<=k;
if (bi>=p2) and (bi<p2_) then do;
/*Change in effect of x2 with respect to b2*/
bj=bi-p2+1;
if dummy[bi,1]=1 then do;
xx=1;
else;
xx=x2[1,bj];
endif;
dj[bi,xi]=phiw2*(1-w2*b2[xj,1]*xx);
endif;
if (bi>=p2_) and (bi<k) then do;
/*Change in effect of x2 with respect to b2_*/
bj=bi-p2_+1;
if dummy[bi,1]=1 then do;
xx=1;
else;
xx=x2[1,bj];
endif;
dj[bi,xi]=(-1)*b2*xx*w2*phiw2;
endif;
bi=bi+1;
endo;
endif;
if (xi>=p1_) and (xi<p2) then do;
dm2[1,xi]=b2_[xi-p1_+1,1]*phiw2;
/*Change in x_*/
bi=p2;
xj=xi-p1_+1;
do while bi<=k;
if (bi>=p2) and (bi<p2_) then do;
/*Change in effect of x_ with respect to b2*/
bj=bi-p2+1;
if dummy[bi,1]=1 then do;
xx=1;
else;
xx=x2[1,bj];
endif;
dj[bi,xi]=(-1)*b2_[xj,1]*xx*w2*phiw2;
endif;
if (bi>=p2_) and (bi<k) then do;
/*Change in effect of x_ with respect to b2_*/
bj=bi-p2_+1;
if dummy[bi,1]=1 then do;
xx=1;
else;
xx=x2[1,bj];
endif;
dj[bi,xi]=phiw2*(1-w2*b2_[xj,1]*xx);
endif;
bi=bi+1;
endo;
endif;
xi=xi+1;
endo;
/*Accumulate marginal effects*/
effects2=effects2+dm2;
/*Calculate the variance matrix of marginal effects*/
vareff=dj'*covb*dj;
sd2=sd2+sqrt(trace(vareff));
i=i+1;
endo;
/*Display the results*/
print 'Equation One';
print 'Average marginal effect Average standard error';
i=1;
do while i<=kx;
print i ' ' effects1[1,i]/n ' ' sd1[1,i]/n;
i=i+1;
endo;
print 'Equation Two';
print 'Average marginal effect Average standard error';
i=1;
do while i<=kx;
print i ' ' effects2[1,i]/n ' ' sd2[1,i]/n;
i=i+1;
endo;