***
*** Example 3.7. Contours and netmaps (Figures 3.13 and 3.14)
***
*** The code for the contours and netmaps is given in GAUSS 
*** because these facilities do not exist in STATA. I have 
*** also calculated the kernel estimates in GAUSS, though it 
*** would also have been possible to go to this stage in STATA.
*** 
*** Note: The following provides 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.
***
/* GAUSS program for contouring joint density of log pce for first
two years of CDI data, results shown in Figures 3.13 and 3.14 */ 

/* the data pan8586 of log pce in the two years was
created in STATA and converted to GAUSS format */
library pgraph;
open f1="pan8586";
x=readr(f1,10000);
f1=close(f1);
x=packr(x);
n=rows(x);
v=(x-meanc(x)')'*(x-meanc(x)')/n;
dt=det(v);
v_1=invpd(v);
/* the inverse variance-covariance matrix of the data */
xmin=3;
xmax=8;
ymin=3;
ymax=8; 

grid=99;
stx=(xmax-xmin)/(grid-1);
sty=(ymax-ymin)/(grid-1);
ans=zeros(grid,grid);
ggx=seqa(xmin,stx,grid);
ggy=seqa(ymin,sty,grid); 

/* setting the bandwidth */
h=1.0;
/* doing the density estimation */
ic=1;
do until ic==grid+1;
ic;
jc=1;
do until jc==grid+1;
locate 10,1; ic;;jc;
xx=ggx[ic,1]|ggy[jc,1];
df=(x-xx')/h;
t2=v_1[1,1]*df[.,1]^2+2*v_1[1,2]*df[.,1].*df[.,2]
+v_1[2,2]*df[.,2]^2;
ans[ic,jc]=sumc(selif((1-t2),t2 .lt 1));
jc=jc+1;
endo;
ic=ic+1;
endo;
ans=ans*2/(pi*n*h^2*sqrt(dt));
ans=missrv(ans,0);
/* save for future use */
save ansp1h1=ans; 

/* making Figure 3.13 */
graphset;
begwind;
xtics(3,8,1.0,2);
ytics(3,8,1.0,2);
_pdate="";
_plctrl=-1;
_pstype=4;
_psymsiz=1;
_plev={0,0.01,0.02,0.03,0.05,0.1,0.2,0.4,0.6,0.8,1.0,1.2,1.4};
makewind(10,8,0,0,1);
contour(ggx',ggy,ans');
makewind(10,8,0,0,1);
xy(x[.,1],x[.,2]);
endwind;
gkp=seqa(1,2,50);
graphset;
xtics(3,8,1.0,2);
ytics(3,8,1.0,2);
ztics(0,0.4,0.1,1);
_pdate="";
_pticout=1;
begwind;
window(1,2,0);
surface(ggx[gkp,.]',ggy[gkp,.],ans[gkp,gkp]');
nextwind;
_pview={1,10,-3,4};
surface(ggx[gkp,.]',ggy[gkp,.],ans[gkp,gkp]');
endwind;
end; 

[Return to Program List]