@ ------------------------------------------------------------ @
@ GAUSS code for structural estimation                         @                                              
@ ------------------------------------------------------------ @

/* First: run stage 1 panel IID/ARCH/GARCH to put gmm, t, n in memory */
/* Second: run stage 1 AGGREGATE to put gmm in memory */

@ +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ @

cls;
library optmum pgraph;
load macrodata[25,5]=macro.txt;

par=2; 
df=25-par;
output_file=zeros(kkk,8);   @ output file obtained from step 2 @

@ ++++++++++++++++++++++++++++++++++ Step 1 ++++++++++++++++++++++++++++++++++++++++++++++ @

par_in1=1|-1;
W=eye(25);
gg=zeros(25,1);
{xout,fout,gout,cout}=optmum(&ftn,par_in1);
par_out1=zeros(par,1);
par_out1[1,1]=exp(xout[1,1]);
par_out1[2,1]=exp(xout[2,1]);

gra1=gradp(&mom,par_out1);

S=vcm(par_out1);  

sqn=sqrt(n)~sqrt(n);
vacov1=inv((sqn.*gra1)'*W*(sqn.*gra1))*((sqn.*gra1)'*W*S*W*(sqn.*gra1))*inv((sqn.*gra1)'*W*(sqn.*gra1));
ase1=sqrt(diag(vacov1));

@ ++++++++++++++++++++++++++++++++++++ Step 2 ++++++++++++++++++++++++++++++++++++++++++++ @

par_in2=1|-1;
W=inv(S); 
gg=zeros(25,1);
{xout,fout,gout,cout}=optmum(&ftn,par_in2);
par_out2=zeros(par,1);
par_out2[1,1]=exp(xout[1,1]);
par_out2[2,1]=exp(xout[2,1]);

gra2=gradp(&mom,par_out2);
sqn=sqrt(n)~sqrt(n);
vacov2=inv((sqn.*gra2)'*W*(sqn.*gra2));
ase2=sqrt(diag(vacov2));

test=(sqrt(n).*gg)'*W*(sqrt(n).*gg);
pval=cdfchic(test,df);

@ +++++++++++++++++++++++++++++++++++++ Precautionary component ++++++++++++++++++++++++++ @

output_file2=sortc(output_file,2);

omega=zeros(25,1);
omega_A=zeros(25,1);
omega_B=zeros(25,1);
omega_B2=zeros(25,1);

sto=0;
j=1;
do until j>25;

sta=sto+1;
sto=sto+n[j,1];

omega[j,1]=meanc(output_file2[sta:sto,5]);
omega_A[j,1]=meanc(output_file2[sta:sto,6]);
omega_B[j,1]=meanc(output_file2[sta:sto,7]);
omega_B2[j,1]=meanc(output_file2[sta:sto,8]);

j=j+1;
endo;

@ ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ @

cls;
print "  ";
print "************************************ Estimation output ***************************************";
print "  ";
print "******************************************* Step 1 *******************************************";
print "First step GMM estimates rho and delta (and asymptotic standard errors)";
print par_out1~ase1;
print "   ";
print "******************************************** Step 2 *******************************************";
print "Second step GMM estimates rho and delta (and asymptotic standard errors)";
print par_out2~ase2;
print "  ";
print "Second step test statistic overidentifying restrictions (left) and p-value (middle)";
print test~pval;
print "  ";
print "***********************************************************************************************";

@ ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ @

proc ftn(b);
local rho,delta;
local sta,sto,i,j,year,ini_cons,ini_inc,ratio;
local consdata,consdata2,average,irate;
local g,function;

rho=exp(b[1,1]);
delta=exp(b[2,1]);

consdata=zeros(kkk,4);
consdata[.,1:2]=gmm[.,1:2];
output_file[.,1:2]=gmm[.,1:2];

sto=0;
j=1;
do until j>2069;

sta=sto+1;
sto=sto+t[j,1];
year=gmm[sta,2];

ini_cons=exp(macrodata[year-67,3]);        

i=sta;
do until i>sto;

year=gmm[i,2];
irate=macrodata[year-67,4];
ini_inc=exp(gmm[i,4]);
ratio=ini_inc/ini_cons;

consdata[i,3]=((1/rho)*(irate-delta))+(rho/2)*((((1-teta)^2)*gmm[i,6])+gmm[i,8])*(ratio^2)+((1-teta)*gmm[i,5]+gmm[i,7])*ratio; 
consdata[i,4]=consdata[i,3]+ln(ini_cons);
ini_cons=exp(consdata[i,4]);

output_file[i,3]=ratio;
output_file[i,4]=ratio^2;
output_file[i,5]=(rho/2)*((((1-teta)^2)*gmm[i,6])+gmm[i,8])*(ratio^2);
output_file[i,6]=(rho/2)*((((1-teta)^2)*0.11)+gmm[i,8])*(ratio^2);
output_file[i,7]=(rho/2)*((((1-teta)^2)*gmm[i,6])+gmm[i,8])*(0.75^2);
output_file[i,8]=(rho/2)*((((1-teta)^2)*gmm[i,6])+gmm[i,8])*0.64;

i=i+1;
endo;

j=j+1;
endo;

consdata2=sortc(consdata,2); 
average=zeros(25,1); 
function=0;
g=zeros(25,1); 

sto=0;
j=1;
do until j>25;

sta=sto+1;
sto=sto+n[j,1];

average[j,1]=meanc(consdata2[sta:sto,3]); 
g[j,1]=average[j,1]-macrodata[j,2]; 

j=j+1;
endo;

function=g'*W*g;
gg=g;

retp(function);
endp;

@ +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ @

proc mom(prm);
local rho,delta;
local sta,sto,i,j,year,ini_cons,ini_inc,ratio;
local consdata,consdata2,average,irate;
local g;

rho=prm[1,1];
delta=prm[2,1];

consdata=zeros(kkk,4);
consdata[.,1:2]=gmm[.,1:2];

sto=0;
j=1;
do until j>2069;

sta=sto+1;
sto=sto+t[j,1];
year=gmm[sta,2];
ini_cons=exp(macrodata[year-67,3]);

i=sta;
do until i>sto;

year=gmm[i,2];
irate=macrodata[year-67,4];
ini_inc=exp(gmm[i,4]);
ratio=ini_inc/ini_cons;

consdata[i,3]=((1/rho)*(irate-delta))+(rho/2)*((((1-teta)^2)*gmm[i,6])+gmm[i,8])*(ratio^2)+((1-teta)*gmm[i,5]+gmm[i,7])*ratio;
consdata[i,4]=consdata[i,3]+ln(ini_cons);
ini_cons=exp(consdata[i,4]);

i=i+1;
endo;

j=j+1;
endo;

consdata2=sortc(consdata,2); 
average=zeros(25,1);
g=zeros(25,1); 

sto=0;
j=1;
do until j>25;

sta=sto+1;
sto=sto+n[j,1];

average[j,1]=meanc(consdata2[sta:sto,3]); 
g[j,1]=average[j,1]-macrodata[j,2]; 

j=j+1;
endo;

retp(g);
endp;

@ ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ @

proc vcm(prm);
local rho,delta;
local sta,sto,i,ii,j,year,ini_cons,ini_inc,ratio;
local consdata,consdata2,average,irate;
local S,A,AA,start,k,kk,indS;

rho=prm[1,1];
delta=prm[2,1];

consdata=zeros(kkk,4);
consdata[.,1:2]=gmm[.,1:2];

sto=0;
j=1;
do until j>2069;

sta=sto+1;
sto=sto+t[j,1];
year=gmm[sta,2];

ini_cons=exp(macrodata[year-67,3]);

i=sta;
do until i>sto;

year=gmm[i,2];
irate=macrodata[year-67,4];
ini_inc=exp(gmm[i,4]);
ratio=ini_inc/ini_cons;

consdata[i,3]=((1/rho)*(irate-delta))+(rho/2)*((((1-teta)^2)*gmm[i,6])+gmm[i,8])*(ratio^2)+((1-teta)*gmm[i,5]+gmm[i,7])*ratio; 
consdata[i,4]=consdata[i,3]+ln(ini_cons);
ini_cons=exp(consdata[i,4]);

i=i+1;
endo;

j=j+1;
endo;

S=zeros(25,25);
indS=zeros(25,25);

sto=0;
j=1;
do until j>2069;

sta=sto+1;
sto=sto+t[j,1];

A=zeros(t[j,1],1);
year=consdata[sta,2];
start=year-67;

i=sta;
do until i>sto;

k=i-sta+1;
year=consdata[i,2];
A[k,1]=consdata[i,3]-macrodata[year-67,2]; 
   
i=i+1;
endo;

AA=A*A';

i=sta;
do until i>sto;

ii=sta;
do until ii>sto;

k=i-sta+1;
kk=ii-sta+1;

S[start+(k-1),start+(kk-1)]=S[start+(k-1),start+(kk-1)]+AA[k,kk];
indS[start+(k-1),start+(kk-1)]=indS[start+(k-1),start+(kk-1)]+1;

ii=ii+1;
endo;

i=i+1;
endo;

j=j+1;
endo;

S=S./(indS);

retp(S);
endp;


