/*
**	Monte Carlo Program 
**  for OLS and GLS
**  under Autocorrelation: AR(1)
*/

seed   = 1;
tt     = 100; @ # of observations @
kk     = 2;  @ # of betas @
iter   = 1000; @ # 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  = .9  ; @ e(t) = rho*e(t-1) + v(t) @
vvar = 4;
bandw = 5; 

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

storpwgb   = zeros(iter,1);
storpwgse  = zeros(iter,1);
storpwgsec = zeros(iter,1);
storpwgt   = zeros(iter,1);

storcogb   = zeros(iter,1);
storcogse  = zeros(iter,1);
storcogsec = zeros(iter,1);
storcogt   = zeros(iter,1);

i = 1; do while i <= iter;

@ Generating y with autocorrelation @

ep= sqrt(vvar/(1-rho^2))*rndns(tt,1,seed);
j = 2; do while j <= tt;
	ep[j,1] = rho*ep[j-1,1]+ep[j,1] ;
j = j + 1 ; endo ;

y = x*tb + 2*ep ;

@	OLS using y and x  @

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

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

proc nw(uv,l) ;
local j,omeo,n,w,ome  ;

    n    = rows(uv)  ;
    omeo = uv'uv/n   ;

    j = 1; do while j <= l  ;

        w    = 1 - j/(l+1)       ;
        ome  = uv[j+1:n,.]'uv[1:n-j,.]/n ;
        omeo = omeo + w*(ome+ome')       ;

    j = j + 1 ; endo ;

retp(omeo) ;
endp       ;

vc = invpd(x'x)*(tt*nw(x.*e,bandw))*invpd(x'x);

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

storb[i,1]   = b[2,1];
storse[i,1]  = se[2,1];
storsec[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; stort[i,1]= 0; else; stort[i,1] = 1; endif;

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


@	GLS using y and x  @

gxx = y[1:rows(y)-1,1]~x[2:rows(x),.]~x[1:rows(x)-1,2];
gyy = y[2:rows(y),1];

durb  = invpd(gxx'gxx)*(gxx'gyy);
drho  = durb[1,1];

gyy = (sqrt(1-drho^2)*y[1,.])|(y[2:rows(y),.]-drho*y[1:rows(y)-1,.]);
gxx = (sqrt(1-drho^2)*x[1,.])|(x[2:rows(y),.]-drho*x[1:rows(y)-1,.]);

pwglsb = invpd(gxx'gxx)*(gxx'gyy);
pwe  = gyy - gxx*pwglsb;
pws2 = (pwe'pwe)/(tt-kk);
pwv  = pws2*invpd(gxx'gxx);

pwse  = sqrt(diag(pwv))  ;

storpwgb[i,1]   = pwglsb[2,1];
storpwgse[i,1]  = pwse[2,1];

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


gyy = gyy[2:rows(gyy),.];
gxx = gxx[2:rows(gxx),.];

coglsb = invpd(gxx'gxx)*(gxx'gyy);
coe  = gyy - gxx*coglsb;
cos2 = (coe'coe)/(tt-kk-1);
cov  = cos2*invpd(gxx'gxx);

cose  = sqrt(diag(cov))  ;

storcogb[i,1]   = coglsb[2,1];
storcogse[i,1]  = cose[2,1];

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

i = i + 1; endo;

@ Reporting Monte Carlo results @

output file = automonte.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 fo NW-corrected s.e. of OLS b(2) ="  meanc(storsec) ;
"----------";
"Mean of P-W feasible GLS b(2)                     ="  meanc(storpwgb);
"s.e. of P-W feasible GLS b(2)                     ="  stdc(storpwgb);
"mean of P-W estimated s.e. of feasible GLS b(2)   ="  meanc(storpwgse) ;
"----------";
"Mean of C-O feasible GLS b(2)                     ="  meanc(storcogb);
"s.e. of C-O feasible GLS b(2)                     ="  stdc(storcogb);
"mean of C-O estimated s.e. of feasible GLS b(2)   ="  meanc(storcogse) ;
"----------";
"Rejection Rate of the t-test with usual OLS           ="  meanc(stort);
"Rejection Rate of the t-test with NW-corrected OLS    ="  meanc(stortc);
"Rejection Rate of the t-test with PW feasible GLS     ="  meanc(storpwgt);
"Rejection Rate of the t-test with Co feasible GLS     ="  meanc(storcogt);


output off ;