/*
**	Monte Carlo Program for 2SLS and Hausman Tests 
*/


@ Data generation under Strong Ideal Conditions @

/*
**    True Model: y = beta(1) + beta(2)*x* + e 
**    Observed x = x* + v
**    Instruments: z1, z2
**    x* = 1 + rho*z1 + rho*z2 + u, u iid N(0,1)
*/ 

/*
**  By running this program, we will learn:
**  (i)     effects of the size of measurement errors on OLS
**  (ii)    consistency of 2SLS estimator even under measurement errors
**  (iii)   the weak instruments problem.
*/

seed  = 1;
beta2 = 0.5;
beta1 = 0.5;
vvar  = 1; @ variance of measurement errors @
rho   = 1; @ correlation between x and z @

tt      = 100; @ # of observations @
kk      = 2;  @ # of betas @
iter    = 3000; @ # of sets of different data @

storb  = zeros(iter,1);
storse = zeros(iter,1);
stort  = zeros(iter,1);

stortb  = zeros(iter,1);
stortse = zeros(iter,1);
stortt  = zeros(iter,1);

hausman1 = zeros(iter,1);
hausman2 = zeros(iter,1);


i = 1; do while i <= iter;

@ Generating y  and x with measurement errors @

z1 = rndns(tt,1,seed) ;
z2 = rndns(tt,1,seed) ;
xt = 1 + rho*z1 + rho*z2 + rndns(tt,1,seed);
mea = rndns(tt,1,seed);
y  = beta1 + beta2*xt + 2*rndns(tt,1,seed);
x  = xt + sqrt(vvar)*mea ;

@	OLS using yy and xx  @

yy = y;
xx = ones(tt,1)~x;

b  = invpd(xx'xx)*(xx'yy);
e  = yy - xx*b ;

s2 = (e'e)/(tt-kk);
v  = s2*invpd(xx'xx);

se = sqrt(diag(v));

storb[i,1]  = b[2,1];
storse[i,1] = se[2,1];
stort[i,1] = (b[2,1]-beta2)/se[2,1];

@   2SLS using zz as instruments  @

zz  = ones(tt,1)~z1~z2;

xpx = (xx'zz)*invpd(zz'zz)*(zz'xx) ;
xpy = (xx'zz)*invpd(zz'zz)*(zz'yy) ;
tb  = invpd(xpx)*xpy;
te  = yy - xx*tb ;

ts2 = (te'te)/(tt-kk);
tv  = ts2*invpd(xpx);

tse = sqrt(diag(tv));

stortb[i,1]  = tb[2,1];
stortse[i,1] = tse[2,1];
stortt[i,1] = (tb[2,1]-beta2)/tse[2,1];

@  Hausman Test: Original @

haus = (tb-b)'pinv(s2*invpd(xpx)-s2*invpd(xx'xx))*(tb-b) ;
df   = rank(invpd(xpx)-invpd(xx'xx)) ;
pval = cdfchic(haus,df);
if pval >= 0.05; hausman1[i,1]= 0; else; hausman1[i,1] = 1; endif;

@  Hausman Test: t-test version @

res  = x - zz*invpd(zz'zz)*(zz'x);
xres = xx~res ;
augb = invpd(xres'xres)*(xres'yy);
auge = yy - xres*augb ;
augs2 = (auge'auge)/(tt-cols(xres)) ;
augv  = augs2*invpd(xres'xres);

haus = augb[3,1]/sqrt(augv[3,3]);
df   = tt-4;
pval = 2*cdftc(abs(haus),df);
if pval >= 0.05; hausman2[i,1]= 0; else; hausman2[i,1] = 1; endif;

i = i + 1; endo;

@ Reporting Monte Carlo results @

output file = tslsmonte.out reset;

format /rd 12,3;

"Monte Carlo results";
"-----------";
"Mean of OLS b(2)                   ="  meanc(storb);
"s.e. of OLS b(2)                   ="  stdc(storb);
"mean of estimated s.e. of OLS b(2) ="  meanc(storse) ;
"-----------";
"Mean of 2SLS b(2)                   ="  meanc(stortb);
"s.e. of 2SLS b(2)                   ="  stdc(stortb);
"mean of estimated s.e. of 2SLS b(2) ="  meanc(stortse);
"-----------";
"Rejection Rate of Hausman Test 1    ="  meanc(hausman1);
"Rejection Rate of Hausman Test 2    ="  meanc(hausman2);

output off ;