@ --------------------------------------------------------------------------- @
@ GAUSS code for structural estimation with bootstrapped standard errors      @
@ GARCH case                                                                  @
@ --------------------------------------------------------------------------- @

new;
cls;
library optmum;

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

boot=100;                             // Number of bootstrap samples - for estimation on actual data set boot=1
l=2;                                // Number of regressors structural equation (2)

seed=100998855;	

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

load microdata[31631,4]=micro.txt;
load macrodata[25,5]=macro.txt;
 
/* Calculate growth rates and number of cross sections */

k=rows(microdata);

count=1;
datamid=zeros(1,5);

i=1;
do until i>(k-1);

z1=microdata[i+1,1];
z2=microdata[i+1,2];
z3=microdata[i+1,3]-microdata[i,3];
z4=microdata[i,3];
z5=microdata[i,4];
z=z1~z2~z3~z4~z5;

if microdata[i+1,1]==microdata[i,1];datamid=datamid|z;endif;
if microdata[i+1,1]>microdata[i,1];count=count+1;endif;

i=i+1;
endo;

kk=rows(datamid);
datamid=datamid[2:kk,.];
kkk=rows(datamid);

/* Calculate number of periods for every cross section */

t=ones(2069,1);
j=1;

i=1; do until i>(kkk-1);

if datamid[i+1,1]>datamid[i,1];j=j+1;endif;
if datamid[i+1,1]==datamid[i,1];t[j,1]=t[j,1]+1;endif;

i=i+1;
endo;

/* Datamid sorted by time */

datamid2=sortc(datamid,2);

/* Calculate number of cross sections for every period */

n=ones(25,1);
j=1;

i=1; do until i>(kkk-1);

if datamid2[i+1,2]>datamid2[i,2];j=j+1;endif;
if datamid2[i+1,2]==datamid2[i,2];n[j,1]=n[j,1]+1;endif;

i=i+1;
endo;

/* Calculate cross section fixed effects */

lambdai=zeros(2069,1);
sto=0;

j=1;
do until j>2069;

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

lambdai[j,1]=meanc(datamid[sta:sto,3]);

j=j+1;
endo;

/* Calculate time fixed effects */

lambdat=zeros(25,1);
sto=0;

j=1;
do until j>25;

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

lambdat[j,1]=meanc(datamid2[sta:sto,3]);

j=j+1;
endo;

/* Calculate transformed growth rates */

datafinal=zeros(kkk,4);
datafinal[.,1:2]=datamid[.,1:2];
year=datamid[1,2];
datafinal[1,3]=datamid[1,3]-lambdat[year-67,1];
datafinal[1,4]=lambdat[year-67,1];

i=1;do until i>(kkk-1);

year=datamid[i+1,2];
datafinal[i+1,3]=datamid[i+1,3]-lambdat[year-67,1];
datafinal[i+1,4]=lambdat[year-67,1];

i=i+1;
endo;

@ ++++++++++++++++++++++++++++++++++++++++++++++++++ Cross section index +++++++++++++++++++++++++++++++++++++++++++++++++ @

micro=datafinal[.,1:3]~datamid[.,4]~datamid[.,5];   

k=rows(micro);
index=zeros(k,1);

j=1;
i=1; do until i>(k-1);

index[i,1]=j;
if micro[i+1,1]>micro[i,1];j=j+1;endif;

i=i+1;
endo;

index[k,1]=2069;
micro=index~micro[.,2:5];

@ ++++++++++++++++++++++++++++++++++++++++++++ Reshuffle consumers and macro data +++++++++++++++++++++++++++++++++++++++++++ @

// Start loop

coeff_structural_1=zeros(boot,l);
coeff_structural_2=zeros(boot,l);

omega_boot=zeros(25,boot);

m=1;do until m>boot;

// Reshuffle micro

temp=zeros(1,5);indexb=0;
shuffle=floor(1+(rndus(2069,1,seed).*2069)); 

j=1; do until j>2069;

temp1=indexcat(micro[.,1],shuffle[j,1]);
nr=rows(temp1);
temp2=micro[temp1[1,1]:temp1[nr,1],.];
indexb2=j*ones(nr,1);

temp=temp|temp2; 
indexb=indexb|indexb2;

j=j+1;endo;

kbtemp=rows(temp);
microb=temp[2:kbtemp,.];
indexb=indexb[2:kbtemp,.];
microb=indexb~microb[.,2:5];
kb=rows(microb);

// Calculate number of periods for every cross section

tb=ones(2069,1);

j=1;
i=1; do until i>(kb-1);

if microb[i+1,1]>microb[i,1];j=j+1;endif;
if microb[i+1,1]==microb[i,1];tb[j,1]=tb[j,1]+1;endif;

i=i+1;
endo;

microb2=sortc(microb,2);

// Calculate number of cross sections for every period 

nb=ones(25,1);

j=1;
i=1; do until i>(kb-1);

if microb2[i+1,2]>microb2[i,2];j=j+1;endif;
if microb2[i+1,2]==microb2[i,2];nb[j,1]=nb[j,1]+1;endif;

i=i+1;
endo;

// Estimation micro income process

if boot==1;data=micro;t=t;n=n;k=k;else;
data=microb;t=tb;n=nb;k=kb;endif;

par=4;     
par_in1=2|1|-3|-3;
{xout,fout,gout,cout}=optmum(&li,par_in1);
par_final=trans(xout);
teta=par_final[1,1];
sigma=par_final[2,1];
alfa=par_final[3,1];
beta=par_final[4,1];

// Output

final=zeros(k,8);
final[.,1:2]=data[.,1:2];
final[.,4]=data[.,4];
final[.,3]=data[.,5];

final=fin(par_final);
final2=sortc(final,2);

// Reshuffle macro (lambdat)

shuffle=floor(1+(rndu(25,1).*25));    @adjust!@

lambdatb=zeros(25,1);
j=1; do until j>25;
ind=shuffle[j,1];
lambdatb[j,1]=lambdat[ind,1];
j=j+1;
endo;

// Estimation macro income process

if boot==1;data_agg=lambdat;
else;data_agg=lambdatb;endif;

par_agg=2;     
par_in1_agg=0|-1;
{xout,fout,gout,cout}=optmum(&li_agg,par_in1_agg);
par_final_agg=trans_agg(xout);
kap=par_final_agg[1,1];  
varpi=par_final_agg[2,1];

// Output

final=fin_agg(par_final_agg);
omega_calc=zeros(k,3);

@ ++++++++++++++++++++++++++++++++++++++++++++ Estimate structural +++++++++++++++++++++++++++++++++++++++++++ @

// GMM estimation structural model

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

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

coeff_structural_1[m,.]=par_out1';

S=vcm(par_out1);  

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

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

coeff_structural_2[m,.]=par_out2';

if par_out2[1,1]>50;continue;endif;

@ +++++ Omega +++++ @

omega_calc2=sortc(omega_calc,2);

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

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

omega_boot[j,m]=meanc(omega_calc2[sta:sto,3]);

j=j+1;
endo;

m=m+1;
endo;

if boot>1;

// Standard errors parameters and 5 - 95 error bounds omega

se_structural_1=stdc(coeff_structural_1);
se_structural_2=stdc(coeff_structural_2);

perc5=zeros(25,1);
perc95=zeros(25,1);

j=1;
do until j>25;

omega_sort=sortc(omega_boot[j,.]',1);
perc5[j,1]=omega_sort[0.05*boot,1];
perc95[j,1]=omega_sort[0.95*boot,1];

j=j+1;
endo;

endif;

@ ---------------------------------------------------------------------------------------------------------------------------------------------------- @

@ ++++++++++++++++++++++++++++++++++ Procedures Micro Income ++++++++++++++++++++++++++++++++++++++++++++++ @
@ ++++++++++++++++++++++++++++++++ Likelihood +++++++++++++++++++++++++++++++++++++ @

proc li(prmtr);
local par_in,teta,vareps,arch,garch,sum,sumli,sto,sta;
local i,j,Y,vt,vt2,vart;
 
par_in=trans(prmtr);
teta=par_in[1,1];
vareps=par_in[2,1];
arch=par_in[3,1];
garch=par_in[4,1];

sum=0;sto=0;

j=1;
do until j>2069;

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

Y=data[sta:sto,3];
vt=0;
vt2=vareps;
vart=vareps;

i=1; do until i>t[j,1];

vart=(1-arch-garch)*vareps+arch*vt2+garch*vart;
vt=Y[i,1]+teta*vt;
vt2=vt^2;

if vart<0;break;endif;

sum=sum+ln(vart)+vt2*inv(vart);

i=i+1; endo;

j=j+1;
endo;

sumli=-(-0.5*(kkk*ln(2*pi))-0.5*sum);

retp(sumli);
endp;

@ ++++++++++++++++++++++++++++++++ Final +++++++++++++++++++++++++++++++++++++ @

proc fin(prmtr);
local par_in,teta,vareps,arch,garch,sto,sta;
local i,j,Y,vt,vt2,vart,shock,variance;
 
par_in=prmtr;
teta=par_in[1,1];
vareps=par_in[2,1];
arch=par_in[3,1];
garch=par_in[4,1];

sto=0;

j=1;
do until j>2069;

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

Y=data[sta:sto,3];
vt=0;
vt2=vareps;
vart=vareps;
shock=zeros(t[j,1],1);
variance=zeros(t[j,1],1);

i=1; do until i>t[j,1];

vart=(1-arch-garch)*vareps+arch*vt2+garch*vart;
variance[i,1]=vart;
vt=Y[i,1]+teta*vt;
shock[i,1]=vt;
vt2=vt^2;

//if vart<0;break;endif;
// not necessary: vart not negative at final parameter values 

i=i+1; endo;

final[sta:sto,5]=shock;
final[sta:sto,6]=variance;

j=j+1;
endo;

retp(final);
endp;

@ ++++++++++++++++++++++++++++++++ Transformation +++++++++++++++++++++++++++++++++++++ @

proc trans(x0); 
local x1;

x1=x0;

x1[1,1]=1/(1+exp(-x0[1,1]));  
x1[2,1]=exp(x0[2,1]);
x1[3,1]=1/(1+exp(-x0[3,1]));
x1[4,1]=(1/(1+exp(-x0[4,1])))-(1/(1+exp(-x0[3,1])));

retp(x1);
endp;

@ ++++++++++++++++++++++++++++++++++ Procedures Macro Income ++++++++++++++++++++++++++++++++++++++++++++++ @
@ ++++++++++++++++++++++++++++++++ Likelihood +++++++++++++++++++++++++++++++++++++ @

proc li_agg(prmtr);
local par_in,kap,varpi,sum,sumli;
local i,Y,vt,vt2;
 
par_in=trans_agg(prmtr);
kap=par_in[1,1];
varpi=par_in[2,1];

sum=0;
Y=data_agg;
vt=0;
vt2=varpi;

i=1;do until i>25;

vt=Y[i,1]-kap;
vt2=vt^2;

sum=sum+ln(varpi)+vt2*inv(varpi);

i=i+1;endo;

sumli=-(-0.5*(25*ln(2*pi))-0.5*sum);

retp(sumli);
endp;

@ ++++++++++++++++++++++++++++++++ Data +++++++++++++++++++++++++++++++++++++ @

proc fin_agg(prmtr);
local par_in,kap,varpi;
local i,Y,vt,vt2,shock_agg,variance_agg;
 
par_in=prmtr;
kap=par_in[1,1];
varpi=par_in[2,1];

Y=data_agg;
vt=0;
vt2=varpi;
shock_agg=zeros(25,1);
variance_agg=zeros(25,1);

i=1;do until i>25;

variance_agg[i,1]=varpi;
vt=Y[i,1]-kap;
shock_agg[i,1]=vt;
vt2=vt^2;

i=i+1;endo;

i=1;do until i>k;

year=final[i,2];
final[i,7]=shock_agg[year-67,1];
final[i,8]=variance_agg[year-67,1];

i=i+1;endo;

retp(final);
endp;

@ ++++++++++++++++++++++++++++++++ Transformation +++++++++++++++++++++++++++++++++++++ @

proc trans_agg(x0); 
local x1;

x1=x0;
x1[2,1]=exp(x0[2,1]);

retp(x1);
endp;

@ ++++++++++++++++++++++++++++ Procedures GMM - Structural estimation +++++++++++++++++++ @
@ ++++++++++++++++++++++++++++++ Criterion function +++++++++++++++++++++++++++++++++++++ @

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(k,4);
consdata[.,1:2]=final[.,1:2];
omega_calc[.,1:2]=final[.,1:2];

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

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

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

i=sta;
do until i>sto;

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

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

omega_calc[i,3]=(rho/2)*((((1-teta)^2)*final[i,6])+final[i,8])*(ratio^2);

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;

@ ++++++++++++++++++++++++++++++++++++++ VCM moment conditions ++++++++++++++++++++++++++++++++++++++++++++++++ @

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,l,ll,indS;

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

consdata=zeros(k,4);
consdata[.,1:2]=final[.,1:2];

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

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

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

i=sta;
do until i>sto;

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

consdata[i,3]=((1/rho)*(irate-delta))+(rho/2)*((((1-teta)^2)*final[i,6])+final[i,8])*(ratio^2)+((1-teta)*final[i,5]+final[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;

l=i-sta+1;
year=consdata[i,2];
A[l,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;

l=i-sta+1;
ll=ii-sta+1;

S[start+(l-1),start+(ll-1)]=S[start+(l-1),start+(ll-1)]+AA[l,ll];
indS[start+(l-1),start+(ll-1)]=indS[start+(l-1),start+(ll-1)]+1;

ii=ii+1;
endo;

i=i+1;
endo;

j=j+1;
endo;

S=S./(indS);

retp(S);
endp;

