@ --------------------------------------------------------------------------- @
@ GAUSS code for reduced form 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
shuffle_res=0;                                  // Reshuffle residuals consumption equation (=1)

l=4;                                            // Number of regressors reduced form equation (2/3/4)

//coeff_con={-0.025,0.428};                     // Value parameters reduced form consumption equation (order:cte,var,lambdat,rate)- only necessary if shuffle_res=1
//coeff_con={0.017,0.091,0.087};
coeff_con={-0.022,0.390,0.048,0.054};

//coeff_con={-0.029,0.465};
//coeff_con={0.013,0.171,0.261};
//coeff_con={-0.020,0.341,0.096,0.192};

iv=0;                                           // OLS (iv=0) or IV (iv=1)
lag_ols=1;
lag_iv=3;

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];   

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:3];

@ +++++++++++++++++++++++++++++++++++++++ Reshuffle consumers ++++++++++++++++++++++++++++++++++++++++ @

// Start loop

coeff_reduced=zeros(boot,l);

m=1;do until m>boot;

// Reshuffle micro

temp=zeros(1,3);indexb=0;
shuffle=floor(1+(rndus(2069,1,seed).*2069));  //random uniform with seed

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:3];
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;

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

// Estimation micro income process

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

outp=zeros(k,3);
outp[.,1:2]=data[.,1:2];

outp=fin(par_final);
outp2=sortc(outp,2);

// Generated regressor h bar

varb=zeros(25,1);

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

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

varb[j,1]=meanc(outp2[sta:sto,3]);

j=j+1;
endo;

@ ++++++++++++++++++++++++++++++++++++++++++ Estimate reduced form ++++++++++++++++++++++++++++++++++ @

rate=macrodata[1:25,4];
cgrow=macrodata[1:25,2];

// 1st lag income

lambdatlag=zeros(25,1);
j=2;do until j>25;
lambdatlag[j,1]=lambdat[j-1,1];
j=j+1;endo;

// Reshuffle residuals reduced form

if shuffle_res==1;

if l==2;
res=cgrow-coeff_con[1,1]-coeff_con[2,1]*varb;	
endif;
if l==3;
res=cgrow-coeff_con[1,1]-coeff_con[2,1]*lambdatlag-coeff_con[3,1]*rate;	
endif;
if l==4;
res=cgrow-coeff_con[1,1]-coeff_con[2,1]*varb-coeff_con[3,1]*lambdatlag-coeff_con[4,1]*rate;	
endif;

resb=zeros(25,1);	
		
if iv==0;
	
shuffle=(zeros(lag_ols,1))|(floor(1+(rndus((25-lag_ols),1,seed).*(25-lag_ols))));  //random uniform with seed

j=(1+lag_ols); 
do until j>25;
ind=shuffle[j,1];
resb[j,1]=res[(ind+lag_ols),1];
j=j+1;
endo;

endif;

if iv==1;
	
shuffle=(zeros(lag_iv,1))|(floor(1+(rndus((25-lag_iv),1,seed).*(25-lag_iv))));    //random uniform with seed

j=(1+lag_iv); 
do until j>25;
ind=shuffle[j,1];
resb[j,1]=res[(ind+lag_iv),1];
j=j+1;
endo;

endif;

if l==2;
cgrowb=coeff_con[1,1]+coeff_con[2,1]*varb+resb;	
endif;
if l==3;
cgrowb=coeff_con[1,1]+coeff_con[2,1]*lambdatlag+coeff_con[3,1]*rate+resb;	
endif;
if l==4;
cgrowb=coeff_con[1,1]+coeff_con[2,1]*varb+coeff_con[3,1]*lambdatlag+coeff_con[4,1]*rate+resb;	
endif;

cgrow=cgrowb;	

endif;

// Specification OLS

if iv==0;
	
T_eff=25-lag_ols;

if l==2;
x=ones(T_eff,1)~varb[(1+lag_ols):25,1];
endif;
if l==3;
x=ones(T_eff,1)~lambdatlag[(1+lag_ols):25,1]~rate[(1+lag_ols):25,1];
endif;
if l==4;
x=ones(T_eff,1)~varb[(1+lag_ols):25,1]~lambdatlag[(1+lag_ols):25,1]~rate[(1+lag_ols):25,1];
endif;
y=cgrow[(1+lag_ols):25,1];

coeff_ols=inv(x'*x)*x'*y;
coeff_reduced[m,.]=coeff_ols';

endif;

// Specification IV

if iv==1;
	
// Instrument lags (-2,-3)

lambdatlag2=zeros(25,1);
j=3;do until j>25;
lambdatlag2[j,1]=lambdat[j-2,1];
j=j+1;endo;

lambdatlag3=zeros(25,1);
j=4;do until j>25;
lambdatlag3[j,1]=lambdat[j-3,1];
j=j+1;endo;

varblag2=zeros(25,1);
j=3;do until j>25;
varblag2[j,1]=varb[j-2,1];
j=j+1;endo;

varblag3=zeros(25,1);
j=4;do until j>25;
varblag3[j,1]=varb[j-3,1];
j=j+1;endo;

ratelag2=zeros(25,1);
j=3;do until j>25;
ratelag2[j,1]=rate[j-2,1];
j=j+1;endo;

ratelag3=zeros(25,1);
j=4;do until j>25;
ratelag3[j,1]=rate[j-3,1];
j=j+1;endo;

cgrowlag2=zeros(25,1);
j=3;do until j>25;
cgrowlag2[j,1]=cgrow[j-2,1];
j=j+1;endo;

cgrowlag3=zeros(25,1);
j=4;do until j>25;
cgrowlag3[j,1]=cgrow[j-3,1];
j=j+1;endo;
	
T_eff=25-lag_iv;

if l==2;
z=ones(T_eff,1)~varblag2[(1+lag_iv):25,1]~varblag3[(1+lag_iv):25,1]~cgrowlag2[(1+lag_iv):25,1]~cgrowlag3[(1+lag_iv):25,1];
x=ones(T_eff,1)~varb[(1+lag_iv):25,1];
endif;
if l==3;
z=ones(T_eff,1)~lambdatlag2[(1+lag_iv):25,1]~lambdatlag3[(1+lag_iv):25,1]
~ratelag2[(1+lag_iv):25,1]~ratelag3[(1+lag_iv):25,1]~cgrowlag2[(1+lag_iv):25,1]~cgrowlag3[(1+lag_iv):25,1];
x=ones(T_eff,1)~lambdatlag[(1+lag_iv):25,1]~rate[(1+lag_iv):25,1];
endif;
if l==4;
z=ones(T_eff,1)~varblag2[(1+lag_iv):25,1]~varblag3[(1+lag_iv):25,1]~lambdatlag2[(1+lag_iv):25,1]~lambdatlag3[(1+lag_iv):25,1]
~ratelag2[(1+lag_iv):25,1]~ratelag3[(1+lag_iv):25,1]~cgrowlag2[(1+lag_iv):25,1]~cgrowlag3[(1+lag_iv):25,1];
x=ones(T_eff,1)~varb[(1+lag_iv):25,1]~lambdatlag[(1+lag_iv):25,1]~rate[(1+lag_iv):25,1];
endif;
y=cgrow[(1+lag_iv):25,1];

coeff_iv=inv(x'*z*inv(z'*z)*z'*x)*x'*z*inv(z'*z)*z'*y;
coeff_reduced[m,.]=coeff_iv';
	
endif;

m=m+1;
endo;

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

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

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 GARCH +++++++++++++++++++++++++++++++++++++ @

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;

outp[sta:sto,3]=variance;

j=j+1;
endo;

retp(outp);
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;


