/*
**	Monte Carlo Program for OLS, OLS with White Correction
**  and GLS under Heteroskedasticity
*/

@ Data generating process @

seed   = 1;    @ seed number @
tt     = 500;  @ # of observations @
kk     = 2;    @ # of betas @
iter   = 4000; @ # of sets of different data @
x  = ones(tt,1)~rndns(tt,kk-1,seed)^2; @ Nonstochastic Regressor @
tb = {1,2}; @ y = x(1)*1 + x(2)*2 + e @
rho =  0.1  ; @ degree of hetyeroskedasticity @

storob   = zeros(iter,1);
storose  = zeros(iter,1);
storosec = zeros(iter,1);
storot   = zeros(iter,1);
storotc  = zeros(iter,1);
storgb   = zeros(iter,1);
storgse  = zeros(iter,1);
storgt   = zeros(iter,1);
storfb   = zeros(iter,1);
storfse  = zeros(iter,1);
storft   = zeros(iter,1);

i = 1; do while i <= iter;

@ Generating y with heteroskedasticity @
y  = x*tb + sqrt(rho*x[.,2]+(1-rho)*1).*rndns(tt,1,seed);

@	OLS using y and x  @

b  = invpd(x'x)*(x'y);
e  = y - x*b ;
ee = e^2;

s2 = (e'e)/(tt-kk);
v  = s2*invpd(x'x);
vc = (x.*e)'(x.*e);
vc = invpd(x'x)*vc*invpd(x'x);

se  = sqrt(diag(v))  ;
sec = sqrt(diag(vc)) ;

storob[i,1]   = b[2,1];
storose[i,1]  = se[2,1];
storosec[i,1] = sec[2,1];

olst = (b[2,1]-tb[2,1])/se[2,1];
df   = tt-2;
pval = 2*cdftc(abs(olst),df);
if pval >= 0.05; storot[i,1]= 0; else; storot[i,1] = 1; endif;

olstc = (b[2,1]-tb[2,1])/sec[2,1];
df   = tt-2;
pval = 2*cdftc(abs(olstc),df);
if pval >= 0.05; storotc[i,1]= 0; else; storotc[i,1] = 1; endif;

@	infeasible GLS using y and x  @

ys = y./sqrt(rho*x[.,2]+(1-rho)*1);
xs = x./sqrt(rho*x[.,2]+(1-rho)*1);

b  = invpd(xs'xs)*(xs'ys);
e  = ys - xs*b ;

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

se  = sqrt(diag(v))  ;

storgb[i,1]   = b[2,1];
storgse[i,1]  = se[2,1];

glst = (b[2,1]-tb[2,1])/se[2,1];
df   = tt-2;
pval = 2*cdftc(abs(glst),df);
if pval >= 0.05; storgt[i,1]= 0; else; storgt[i,1] = 1; endif;

@	feasible GLS using y and x  @

fh = x*invpd(x'x)*(x'ee);

ys = y./sqrt(abs(fh));
xs = x./sqrt(abs(fh));

b  = invpd(xs'xs)*(xs'ys);
e  = ys - xs*b ;

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

se  = sqrt(diag(v))  ;

storfb[i,1]   = b[2,1];
storfse[i,1]  = se[2,1];

glst = (b[2,1]-tb[2,1])/se[2,1];
df   = tt-2;
pval = 2*cdftc(abs(glst),df);
if pval >= 0.05; storft[i,1]= 0; else; storft[i,1] = 1; endif;

i = i + 1; endo;

@ Reporting Monte Carlo results @

output file = hetmonte.out reset;

format /rd 12,3;

"Monte Carlo results";
"-----------";
"Mean of OLS b(2)                              ="  meanc(storob);
"s.e. of OLS b(2)                              ="  stdc(storob);
"mean of estimated s.e. of OLS b(2)            ="  meanc(storose) ;
"mean fo White-corrected s.e. of OLS b(2)      ="  meanc(storosec) ;
"----------";
"Mean of infeasible GLS b(2)                   ="  meanc(storgb);
"s.e. of infeasible GLS b(2)                   ="  stdc(storgb);
"mean of estimated s.e. of infeasible GLS b(2) ="  meanc(storgse) ;
"----------";
"Mean of feasible GLS b(2)                     ="  meanc(storfb);
"s.e. of feasible GLS b(2)                     ="  stdc(storfb);
"mean of estimated s.e. of feasible GLS b(2)   ="  meanc(storfse) ;
"----------";
"Rejection Rate of the t-test with usual OLS           ="  meanc(storot);
"Rejection Rate of the t-test with White-corrected OLS ="  meanc(storotc);
"Rejection Rate of the t-test with infeasible GLS      ="  meanc(storgt);
"Rejection Rate of the t-test with feasible GLS        ="  meanc(storft);

output off ;