/*
**	Monte Carlo Program for Multicollinearity
*/


@ Data generation under Strong Ideal Conditions @

/*    Model: y = beta(1) + beta(2)*x2 + beta(3)*x3 + e    */ 

seed   = 1;
beta3  = 0.5;
beta2  = 0.5;
beta1  = 0.5;

tt        = 100; @ # of observations @
kk      = 3;  @ # of betas @
iter    = 5000; @ # of sets of different data @
rho   = 0.5; @ control correlation between x(2) and x(3) @

x2   = rndns(tt,1,seed);
x3  = rho*x2 + (1-rho)*rndns(tt,1,seed) ;

storb    = zeros(iter,1);
storse   = zeros(iter,1);
stort_tr = zeros(iter,1);
stort_fa = zeros(iter,1);
rejt_tr  = zeros(iter,1);
rejt_fa  = zeros(iter,1);

i = 1; do while i <= iter;

@ Generating y  @

y  = beta1 + beta2*x2 + beta3*x3 + rndns(tt,1,seed);

@	OLS using yy and xx  @

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

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_tr[i,1] = (b[2,1]-beta2)/se[2,1]; @ testing H_o: b(2) = 0.5 @
stort_fa[i,1] = (b[2,1])/se[2,1]; @ testing H_o: b(2) = 0 @
if cdftc(abs(stort_tr[i,1]),tt-kk) <= 0.025; rejt_tr[i,1] = 1; endif;
if cdftc(abs(stort_fa[i,1]),tt-kk) <= 0.025; rejt_fa[i,1] = 1; endif;

i = i + 1; endo;

@ Reporting Monte Carlo results @

output file = mulmonte.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) ;
"% of rejection of true H at 5%       =" meanc(rejt_tr)*100;
"% of rejection of incorrect H at 5%  =" meanc(rejt_fa)*100;

library pgraph;
graphset;
{a1,a2,a3}=hist(storb,50);
output off ;