***
*** Chapter 5
***
*** I provide the code for calculating the system of demand 
*** equations, including the own and cross-price elasticities, 
*** for completing the system, and for calculating the 
*** symmetry-constrained estimates. The code here is for the 
*** Maharashtran case; except for minor details-the number and 
*** names of the goods, and the definition of the other 
*** variables-the Pakistani code is the same. There are four 
*** separate programs: the first, allindia.do, is for estimating 
*** the demand ssystem. Appended to it is a program mkmats.do, that 
*** calculates the commutation and selection matrices required for 
*** the symmetry-constrained estimates, as well as procedures for 
*** making the "vec" of a matrix, and for reversing the operation. 
*** The code bootall.do bootstraps the procedure in order to obtain 
*** measures of sampling variability. Finally, the program policy.do 
*** calculates the efficiency and equity components of the cost-benefit 
*** ratios for price reform and was used to give the results in Tables 
*** 5.12 and 5.13.
*** 
*** Note: The following provides the STATA code used to produce 
*** results used in "The Analysis of Household Surveys: 
*** A Microeconometric Approach to Development Policy", by Angus Deaton.
*** This book, published for the World Bank by The Johns Hopkins University
*** Press and scheduled for release in 1997.
***
program allindia.do
* program for doing eveything in one shot
* requires a data set indinr on log unit values
* budget shares, lnx, and lnn, as well as subrounds
* and regions and household characteristics
*
version 4.0
#delimit ;
drop _all;
set matsize 400;
set maxvar 1724 width 1832;
set more 1;
capture log close;
log using allindia, replace;
*The log unit values begin with luv, e.g. luvwhe;
*The budget shares begin with w, e.g. wric;
*lnexp is log of outlay, lnhhs log of household size;
*regional and subround (seasonal) variables;
*Demographics and other variables to taste;
*Here they are rm1-rm5 and rf1-rf4, age-sex ratios;
*Caste, religion, and labor types;
use indinr;
*These are the commodity identifiers: used as three letter prefixes;
global goods "ric whe jow oce pul dai oil mea veg fru sug";
*Number of goods in the system;
global ngds=11; 

matrix define sig=J($ngds,1,0);
matrix define ome=J($ngds,1,0);
matrix define lam=J($ngds,1,0);
matrix define wbar=J($ngds,1,0);
matrix define b1=J($ngds,1,0);
matrix define b0=J($ngds,1,0); 

*Average budget shares;
cap program drop mkwbar;
program def mkwbar;
local ig=1;
while "`1'" ~= ""{;
qui summ w`1';
matrix wbar[`ig',1]=_result(3);
local ig=`ig'+1;
mac shift};
end; 

mkwbar $goods; 

*First stage regressions: within village;
cap program drop st1reg;
program def st1reg;
local ig=1;
while "`1'" ~= ""{;
*Cluster fixed effect regression;
areg luv`1' rm* rf* lnexp lhhs
schcaste hindu buddhist labtyp*, absorb(cluster);
*Measurement error variance;
matrix ome[`ig',1]=$S_E_sse/$S_E_tdf; 

*Quality elasticity;
matrix b1[`ig',1]=_coef[lnexp];
*These residuals still have cluster effects in;
predict ruv`1', resid;
*Purged y's for next stage;
gen y1`1'=luv`1'-_coef[lnexp]*lnexp-_coef[lhhs]*lhhs-_coef[rm1]*
rm1-_coef[rm2]*rm2-_coef[rm3]*rm3-_coef[rm4]*rm4-_coef[rm5]*
rm5-_coef[rf1]*rf1-_coef[rf2]*rf2-_coef[rf3]*rf3-_coef[rf4]*rf4
-_coef[schcaste]*schcaste-_coef[hindu]*hindu-_coef[buddhist]*
buddhist-_coef[labtyp1]*labtyp1-_coef[labtyp2]*labtyp2
-_coef[labtyp3]*labtyp3-_coef[labtyp4]*labtyp4;
drop luv`1';
*Repeat for budget shares;
areg w`1' rm* rf* lnexp lhhs
schcaste hindu buddhist labtyp*, absorb(cluster);
predict rw`1', resid;
matrix sig[`ig',1]=$S_E_sse/$S_E_tdf;
matrix b0[`ig',1]=_coef[lnexp];
gen y0`1'=w`1'-_coef[lnexp]*lnexp-_coef[lhhs]*lhhs-_coef[rm1]*
rm1-_coef[rm2]*rm2-_coef[rm3]*rm3-_coef[rm4]*rm4-_coef[rm5]*
rm5-_coef[rf1]*rf1-_coef[rf2]*rf2-_coef[rf3]*rf3-_coef[rf4]*rf4
-_coef[schcaste]*schcaste-_coef[hindu]*hindu-_coef[buddhist]*
buddhist-_coef[labtyp1]*labtyp1-_coef[labtyp2]*labtyp2
-_coef[labtyp3]*labtyp3-_coef[labtyp4]*labtyp4;
*This next regression is necessary to get covariance of residuals;
qui areg ruv`1' rw`1' lnexp lhhs rm* rf*
schcaste hindu buddhist labtyp*, absorb(cluster);
matrix lam[`ig',1]=_coef[rw`1']*sig[`ig',1];
drop w`1' rw`1' ruv`1';
local ig=`ig'+1;
mac shift};
end; 

st1reg $goods; 

matrix list sig;
matrix list ome;
matrix list lam;
matrix list b0;
matrix list b1; 

drop lnexp rm* rf*;
drop labtyp* schcaste hindu buddhist lhhs;
*Saving so far as a protection;
save tempa, replace; 

drop _all;
use tempa; 

*Averaging by cluster;
*Counting numbers of obs in each cluster for n and n+;
cap program drop clustit;
program def clustit;
local ig=1;
while "`1'" ~= ""{;
egen y0c`ig'=mean(y0`1'), by(cluster);
egen n0c`ig'=count(y0`1'), by(cluster);
egen y1c`ig'=mean(y1`1'), by(cluster);
egen n1c`ig'=count(y1`1'), by(cluster);
drop y0`1' y1`1';
local ig=`ig'+1;
mac shift };
end; 

clustit $goods;
sort cluster;
*keeping one obs per cluster; 

*NB subround and region are constant within cluster;
qui by cluster: keep if _n==1; 

*Saving cluster level information;
*Use this for shortcut bootstrapping;
save tempclus, replace; 

*Removing province and quarter effects;
qui tab region, gen(regiond);
qui tab subround, gen(quard);
drop regiond6 quard4;
cap program drop purge;
program define purge;
local ig=1;
while `ig' <= $ngds {;
qui regress y0c`ig' quard* regiond*;
predict tm, resid;
replace y0c`ig'=tm;
drop tm;
qui regress y1c`ig' quard* regiond*;
predict tm, resid;
replace y1c`ig'=tm;
drop tm;
local ig=`ig'+1;
};
end; 

purge;
drop regiond* quard*; 

matrix define n0=J($ngds,1,0);
matrix define n1=J($ngds,1,0); 

*Averaging (harmonically) numbers of obs over clusters;
cap program drop mkns;
program define mkns;
local ig=1;
while `ig' <= $ngds {;
replace n0c`ig'=1/n0c`ig';
replace n1c`ig'=1/n1c`ig';
qui summ n0c`ig';
matrix n0[`ig',1]=(_result(3))^(-1);
qui summ n1c`ig';
matrix n1[`ig',1]=(_result(3))^(-1);
drop n0c`ig' n1c`ig';
local ig=`ig'+1;
};
end; 

mkns; 

*Making the intercluster variance and covariance matrices;
*This is done in pairs because of the missing values; 

matrix s=J($ngds,$ngds,0);
matrix r=J($ngds,$ngds,0);
cap program drop mkcov;
program def mkcov;
local ir=1;
while `ir' <= $ngds {;
local ic=1;
while `ic' <= $ngds {;
qui corr y1c`ir' y1c`ic', cov;
matrix s[`ir',`ic']=_result(4);
qui corr y1c`ir' y0c`ic', cov;
matrix r[`ir',`ic']=_result(4);
local ic=`ic'+1;
};
local ir=`ir'+1; 

};
end; 

mkcov; 

*We don't need the data any more;
drop _all; 

matrix list s;
matrix list r; 

*Making OLS estimates;
matrix bols=syminv(s);
matrix bols=bols*r;
display("Second-stage OLS estimates: B-matrix");
matrix list bols;
display("Column 1 is coefficients from 1st regression, etc");

*Corrections for measurement error;
cap program drop fixmat;
program def fixmat;
matrix def sf=s;
matrix def rf=r;
local ig=1;
while `ig' <= $ngds {;
matrix sf[`ig',`ig']=sf[`ig',`ig']-ome[`ig',1]/n1[`ig',1];
matrix rf[`ig',`ig']=rf[`ig',`ig']-lam[`ig',1]/n0[`ig',1];
local ig=`ig'+1;
};
end; 

fixmat; 

matrix invs=syminv(sf);
matrix bhat=invs*rf;
*Estimated B matrix without restrictions;
matrix list bhat; 

*Housekeeping matrices, including elasticities; 

cap program drop mormat;
program def mormat;
matrix def xi=J($ngds,1,0);
matrix def el=J($ngds,1,0);
local ig=1;
while `ig' <= $ngds {;
matrix xi[`ig',1]=b1[`ig',1]/(b0[`ig',1]+
(1-b1[`ig',1]*wbar[`ig',1]));
matrix el[`ig',1]=1-b1[`ig',1]+b0[`ig',1]/wbar[`ig',1];
local ig=`ig'+1;
};
end; 

mormat; 

global ng1=$ngds+1;
matrix iden=I($ngds);
matrix iden1=I($ng1);
matrix itm=J($ngds,1,1);
matrix itm1=J($ng1,1,1);
matrix dxi=diag(xi);
matrix dwbar=diag(wbar);
matrix idwbar=syminv(dwbar);
display("Average budget shares");
matrix tm=wbar';
matrix list tm;
display("Expenditure elasticities");
matrix tm=el'; 

matrix list tm;
display("Quality elasticities");
matrix tm=b1';
matrix list tm; 

*This all has to go in a program to use it again later;
*Basically uses the b matrix to form price elasticity matrix;
cap program drop mkels;
program define mkels;
matrix cmx=bhat';
matrix cmx=dxi*cmx;
matrix cmx1=dxi*dwbar;
matrix cmx=iden-cmx;
matrix cmx=cmx+cmx1;
matrix psi=inv(cmx);
matrix theta=bhat'*psi;
display("Theta matrix");
matrix list theta;
matrix ep=bhat';
matrix ep=idwbar*ep;
matrix ep=ep-iden;
matrix ep=ep*psi;
display("Matrix of price elasticities");
matrix list ep;
end; 

mkels; 

**Completing the system by filling out the matrices;

cap program drop complet;
program define complet;
*First extending theta;
matrix atm=theta*itm;
matrix atm=-1*atm;
matrix atm=atm-b0;
matrix xtheta=theta,atm;
matrix atm=xtheta';
matrix atm=atm*itm;
matrix atm=atm';
matrix xtheta=xtheta\atm;
*Extending the diagonal matrices;
matrix wlast=wbar'*itm;
matrix won=(1);
matrix wlast=won-wlast;
matrix xwbar=wbar\wlast;
matrix dxwbar=diag(xwbar);
matrix idxwbar=syminv(dxwbar);
matrix b1last=(0.25);
matrix xb1=b1\b1last;
matrix b0last=b0'*itm;
matrix b0last=-1*b0last;
matrix xb0=b0\b0last;
matrix xe=itm1-xb1;
matrix tm=idxwbar*xb0;
matrix xe=xe+tm;
matrix tm=xe';
display("extended outlay elasticities");
matrix list tm;
matrix xxi=itm1-xb1;
matrix xxi=dxwbar*xxi;
matrix xxi=xxi+xb0;
matrix tm=diag(xb1);
matrix tm=syminv(tm);
matrix xxi=tm*xxi;
matrix dxxi=diag(xxi);
*Extending psi;
matrix xpsi=dxxi*xtheta; 

matrix xpsi=xpsi+iden1;
matrix atm=dxxi*dxwbar;
matrix atm=atm+iden1;
matrix atm=syminv(atm);
matrix xpsi=atm*xpsi;
matrix ixpsi=inv(xpsi);
*Extending bhat & elasticity matrix;
matrix xbhatp=xtheta*ixpsi;
matrix xep=idxwbar*xbhatp;
matrix xep=xep-iden1;
matrix xep=xep*xpsi;
display("extended matrix of elasticities");
matrix list xep;
end; 

complet; 

**Calculating symmetry restricted estimators;
**These are only approximately valid & assume no quality effects;
run mkmats;
vecmx bhat vbhat;
** R matrix for restrictions;
lmx $ngds llx;
commx $ngds k;
global ng2=$ngds*$ngds;
matrix bigi=I($ng2);
matrix k=bigi-k;
matrix r=llx*k;
matrix drop k;
matrix drop bigi;
matrix drop llx;
** r vector for restrictions, called rh;
matrix rh=b0#wbar;
matrix rh=r*rh;
matrix rh=-1*rh;
**doing the constrained estimation;
matrix iss=iden#invs;
matrix rp=r';
matrix iss=iss*rp;
matrix inn=r*iss;
matrix inn=syminv(inn);
matrix inn=iss*inn;
matrix dis=r*vbhat;
matrix dis=rh-dis;
matrix dis=inn*dis;
matrix vbtild=vbhat+dis;
unvecmx vbtild btild; 

**the following matrix should be symmetric;
matrix atm=b0';
matrix atm=wbar*atm;
matrix atm=btild+atm;
matrix list atm; 

**going back to get elasticities and complete sytem;
matrix bhat=btild;
mkels;
complet; 

log close;

program mkmats.do
**calculates two matrices, the commutation matrix and
**the lower diagonal selection matrix that are needed
**in the main calculations; these are valid only for
**square matrices
**also a routine for taking the vec of a matrix
**and a matching unvec routine 

**for calculating the commutation matrix k
**the matrix is defined by K*vec(A)=vec(A') 

cap program drop commx
program define commx
local n2=`1'^2
matrix `2'=J(`n2',`n2',0)
local i=1
local ik=0
while `i' <= `1'{
local j=1
local ij=`i'
while `j' <= `1'{
local ir=`j'+`ik'
matrix `2'[`ir',`ij']=1
local ij=`ij'+`1'
local j=`j'+1
}
local i=`i'+1
local ik=`ik'+`1'
}
end 

**for vecing a matrix, i.e. stacking it into a column vector 

cap program drop vecmx
program def vecmx
local n=rowsof(`1')
local n2=`n'^2
matrix def `2'=J(`n2',1,0)
local j=1
while `j' <= `n' {
local i=1
while `i' <= `n' {
local vcel=(`j'-1)*`n'+`i'
matrix `2'[`vcel',1]=`1'[`i',`j']
local i=`i'+1
}
local j=`j'+1
}
end 

*program for calculating the matrix that extracts
*from vec(A) the lower left triangle of the matrix A 

cap program drop lmx
program define lmx
local ng2=`1'^2
local nr=0.5*`1'*(`1'-1)
matrix def `2'=J(`nr',`ng2',0)
local ia=2
local ij=1
while `ij' <= `nr'{
local ik=0
local klim=`1'-`ia'
while `ik' <= `klim' {
local ip=`ia'+(`ia'-2)*`1'+`ik'
matrix `2'[`ij',`ip']=1
local ij=`ij'+1
local ik=`ik'+1
}
local ia=`ia'+1
}
end 

**program for unvecing the vec of a square matrix; 

cap program drop unvecmx
program def unvecmx 

local n2=rowsof(`1')
local n=sqrt(`n2')
matrix def `2'=J(`n',`n',0)
local i=1
while `i' <= `n' {
local j=1
while `j' <= `n' {
local vcel=(`j'-1)*`n'+`i'
matrix `2'[`i',`j']=`1'[`vcel',1]
local j=`j'+1
}
local i=`i'+1
}
end 

program bootall.do
**for bootstrapping Indian demand estimates 

version 4.0
#delimit ;
capture log close;
set more 1;
drop _all;
do allindia;
log using bootall, replace;
drop _all;
set maxvar 300;
vecmx xep vxep;
set obs 1;
gen reps=0;
global nels=$ng1*$ng1;
global nmc=1000;
cap program drop vtodat;
program define vtodat;
local ic=1;
while `ic' <= $nels {;
gen e`ic'=vxep[`ic',1];
local ic=`ic'+1;
};
end;
vtodat;
save bootall, replace;
drop _all; 

cap program drop purge;
program define purge;
local ig=1;
while `ig' <= $ngds {;
qui regress y0c`ig' quard* regiond*;
predict tm, resid;
replace y0c`ig'=tm;
drop tm;
qui regress y1c`ig' quard* regiond*;
predict tm, resid;
replace y1c`ig'=tm;
drop tm;
local ig=`ig'+1;
};
end; 

cap program drop mkns;
program define mkns;
local ig=1;
while `ig' <= $ngds {;
replace n0c`ig'=1/n0c`ig';
replace n1c`ig'=1/n1c`ig';
qui summ n0c`ig';
matrix n0[`ig',1]=(_result(3))^(-1);
qui summ n1c`ig'; 

matrix n1[`ig',1]=(_result(3))^(-1);
drop n0c`ig' n1c`ig';
local ig=`ig'+1;
};
end; 

cap program drop mkcov;
program def mkcov;
local ir=1;
while `ir' <= $ngds {;
local ic=1;
while `ic' <= $ngds {;
qui corr y1c`ir' y1c`ic', cov;
matrix s[`ir',`ic']=_result(4);
qui corr y1c`ir' y0c`ic', cov;
matrix r[`ir',`ic']=_result(4);
local ic=`ic'+1;
};
local ir=`ir'+1;
};
end; 

cap program drop fixmat;
program def fixmat;
matrix def sf=s;
matrix def rf=r;
local ig=1;
while `ig' <= $ngds {;
matrix sf[`ig',`ig']=sf[`ig',`ig']-ome[`ig',1]/n1[`ig',1];
matrix rf[`ig',`ig']=rf[`ig',`ig']-lam[`ig',1]/n0[`ig',1];
local ig=`ig'+1;
};
end; 

cap program drop mkels;
program define mkels;
matrix cmx=bhat';
matrix cmx=dxi*cmx;
matrix cmx1=dxi*dwbar;
matrix cmx=iden-cmx;
matrix cmx=cmx+cmx1;
matrix psi=inv(cmx);
matrix theta=bhat'*psi;
display("Theta matrix");
matrix list theta;
matrix ep=bhat';
matrix ep=idwbar*ep;
matrix ep=ep-iden;
matrix ep=ep*psi;
end; 

cap program drop complet;
program define complet;
*First extending theta;
matrix atm=theta*itm;
matrix atm=-1*atm;
matrix atm=atm-b0;
matrix xtheta=theta,atm;
matrix atm=xtheta';
matrix atm=atm*itm;
matrix atm=atm';
matrix xtheta=xtheta\atm;
*Extending the diagonal matrices;
matrix wlast=wbar'*itm;
matrix won=(1);
matrix wlast=won-wlast;
matrix xwbar=wbar\wlast;
matrix dxwbar=diag(xwbar);
matrix idxwbar=syminv(dxwbar); 

matrix b1last=(0.25);
matrix xb1=b1\b1last;
matrix b0last=b0'*itm;
matrix b0last=-1*b0last;
matrix xb0=b0\b0last;
matrix xe=itm1-xb1;
matrix tm=idxwbar*xb0;
matrix xe=xe+tm;
matrix tm=xe';
matrix xxi=itm1-xb1;
matrix xxi=dxwbar*xxi;
matrix xxi=xxi+xb0;
matrix tm=diag(xb1);
matrix tm=syminv(tm);
matrix xxi=tm*xxi;
matrix dxxi=diag(xxi);
*Extending psi;
matrix xpsi=dxxi*xtheta;
matrix xpsi=xpsi+iden1;
matrix atm=dxxi*dxwbar;
matrix atm=atm+iden1;
matrix atm=syminv(atm);
matrix xpsi=atm*xpsi;
matrix ixpsi=inv(xpsi);
*Extending bhat & elasticity matrix;
matrix xbhatp=xtheta*ixpsi;
matrix xep=idxwbar*xbhatp;
matrix xep=xep-iden1;
matrix xep=xep*xpsi;
end; 

cap program drop bootindi;
program define bootindi;
local expno=1;
while `expno' <= $nmc {;
display("Experiment number `expno'");
quietly {;
use tempclus;
bsample _N;
qui tab region, gen(regiond);
qui tab subround, gen(quard);
drop regiond6 quard4;
purge;
drop regiond* quard*;
matrix define n0=J($ngds,1,0);
matrix define n1=J($ngds,1,0);
*Averaging (harmonically) numbers of obs over clusters;
mkns;
*Making the intercluster variance and covariance matrices;
*This is done in pairs because of the missing values;
matrix s=J($ngds,$ngds,0);
matrix r=J($ngds,$ngds,0);
mkcov;
*We don't need the data any more;
drop _all;
*Making OLS estimates;
matrix bols=syminv(s);
matrix bols=bols*r;
*Corrections for measurement error;
fixmat;
matrix invs=syminv(sf);
matrix bhat=invs*rf;
global ng1=$ngds+1;
matrix iden=I($ngds);
matrix iden1=I($ng1);
matrix itm=J($ngds,1,1);
matrix itm1=J($ng1,1,1);
matrix dxi=diag(xi);
matrix dwbar=diag(wbar); 

matrix idwbar=syminv(dwbar);
mkels;
**Completing the system by filling out the matrices;
complet;
**Calculating symmetry restricted estimators;
vecmx bhat vbhat;
** R matrix for restrictions;
lmx $ngds llx;
commx $ngds k;
global ng2=$ngds*$ngds;
matrix bigi=I($ng2);
matrix k=bigi-k;
matrix r=llx*k;
matrix drop k;
matrix drop bigi;
matrix drop llx;
** r vector for restrictions, called rh;
matrix rh=b0#wbar;
matrix rh=r*rh;
matrix rh=-1*rh;
**doing the constrained estimation;
matrix iss=iden#invs;
matrix rp=r';
matrix iss=iss*rp;
matrix inn=r*iss;
matrix inn=syminv(inn);
matrix inn=iss*inn;
matrix dis=r*vbhat;
matrix dis=rh-dis;
matrix dis=inn*dis;
matrix vbtild=vbhat+dis;
unvecmx vbtild btild;
**going back to get elasticities and complete sytem;
matrix bhat=btild;
mkels;
complet;
vecmx xep vxep;
set obs 1;
gen reps=`expno';
vtodat;
append using bootall;
save bootall, replace;
drop _all;
local expno=`expno'+1;
};};
end; 

bootindi;
use bootall;
display("Monte Carlo results");
summ; 

log close;

program policy.do
** Does the policy analysis tables starting by running allindia
** and then processes to produce the two tables Table 5.12 and 13
** 

version 4.0
#delimit ;
set more 1;
drop _all; /*
do allindia; */
capture log close;
log using policy, replace; 

drop _all; 

global xgoods "ric whe jow oce pul dai oil mea veg fru sug oth";
matrix ap=(1.5\1.25\1.0\1.0\1.0\1.0\0.67\1.0\1.0\1.0\1.0\1.0); 

cap program drop mkws;
program def mkws;
use indinr;
keep lnexp lhhs w*;
gen xh=exp(lnexp);
gen nh=exp(lhhs);
gen woth=1-wwhe-wric-wjow-woce-wpul-wdai-wmea-woil-wsug-wveg-wfru;
matrix wtild=J(12,1,0);
matrix weps=J(12,4,0);
*Plutocratic and soc rep budget shares;
local ig=1;
while "`1'" ~= ""{;
qui summ w`1' [aweight=xh];
matrix wtild[`ig',1]=_result(3);
matrix weps[`ig',1]=_result(3);
gen tem=w`1'*(xh/nh)^(-0.5);
qui summ tem [aweight=xh];
matrix weps[`ig',2]=_result(3);
drop tem;
gen tem=w`1'*(xh/nh)^(-1);
qui summ tem [aweight=xh];
matrix weps[`ig',3]=_result(3);
drop tem;
gen tem=w`1'*(xh/nh)^(-2);
qui summ tem [aweight=xh];
matrix weps[`ig',4]=_result(3);
drop tem;
local ig=`ig'+1;
mac shift};
*Scaling the matrices;
matrix itm1t=itm1';
matrix wsum=itm1t*weps;
local ic=1;
while `ic'<=$ng1 {;
local id=1;
while `id' <=4 {;
local stmp=1/wsum[1,`id'];
matrix weps[`ic',`id']=`stmp'*weps[`ic',`id'];
local id=`id'+1;
};
local ic=`ic'+1;
};
end;

mkws $xgoods; 

cap program drop equi;
program def equi;
local ic=1;
while `ic' <= 4 {;
matrix wtmp=weps[.,`ic'];
matrix tmp=mxw*wtmp;
matrix tmp1=dtot*tmp;
display("equity effect");
matrix list tmp;
display("Cost-benefit ratio");
matrix list tmp1;
local ic=`ic'+1;
end; 

cap program drop calpol;
program def calpol;
matrix tauotau=itm1-ap;
display("Column one");
matrix list tauotau; 

matrix dth=vecdiag(xtheta);
matrix dth=dth';
matrix mxw=diag(wtild);
matrix mxw=syminv(mxw);
matrix midt=mxw*dth;
matrix col2=midt-itm1;
display("second column");
matrix list col2;
matrix dtau=diag(tauotau);
matrix col3=dtau*col2;
display("third column");
matrix list col3;
matrix col4=xtheta';;
matrix col4=col4*tauotau;
matrix col4=mxw*col4;
matrix tmp=dtau*midt;
matrix col4=col4-tmp;
display("fourth column");
matrix list col4;
matrix tot=itm1+col3;
matrix tot=tot+col4;
display("totals");
matrix list tot;
display("Going on to next Table");
matrix dtot=diag(tot);
matrix dtot=syminv(dtot);
equi;
end;

calpol;

log close; 

[Return to Program List]