/*SIMULTANEUS EQUATIONS AS IN CLASS */

new;
 output file = weak.out reset ;
library optmum;

N = 150;

let alpha1= .5; alpha2=.3 ; alpha3=0.3 ; alpha4=.2 ; alpha5=.01 ;
sigma1 = 1; sigma2=2 ;

/*generate some regressors*/

x1=2+ seqa(1,1,N)/N.*rndn(N,1) ;
x2=3+0.5*x1+sigma1*rndn(N,1);

Xvec=x1~x2;

GAM_inv=zeros(2,2) ;

detGAM=1-alpha1*alpha3;

GAM_inv[1,1]=1/detGAM;
GAM_inv[1,2]=alpha3/detGAM;
GAM_inv[2,1]=alpha1/detGAM;
GAM_inv[2,2]=1/detGAM;



Bmat=zeros(2,2) ;

Bmat[1,1]=alpha2 ; Bmat[1,2]=alpha4 ; Bmat[2,1]=0 ;
Bmat[2,2]=alpha5 ;

/* The model is

y1=alpha1*y2+alpha2*x1+u1
y2=alpha3*y1+alpha4*x1+alpha5*x2+u2 */

Nsim =10 ; "Number of simulations:" Nsim ;

/* create a matrix to hold the estimated values */

results = zeros(Nsim,3) ;

rsqrsim=zeros(Nsim,1);

Fsim=zeros(Nsim,1) ;

 estim=0 ; do while estim < Nsim ;
estim=estim+1 ;

/* generate correlated error terms */

u1=sigma1*rndn(N,1) ;

u2=sigma2*rndn(N,1)+0.2*u1 ;

/*generate linear relations */

ymat=xvec*Bmat*GAM_inv+(u1~u2);
y1vec=ymat[.,1];
y2vec=ymat[.,2];

/*2SLS for first equation*/

z1=ones(N,1)~x1~x2 ;

z2=ones(N,1)~x1 ; @instrument only (plus constant)@

/*first stage*/

y2hat=z1*inv(z1'z1)*z1'y2vec ;

y2hat2=z2*inv(z2'z2)*z2'y2vec ;

 ybarvec=meanc(y2vec)*ones(N,1);

 res1=y2vec-y2hat ; res2=y2vec-y2hat2 ;

/*calculate F for weak instrument */

sighat1=res1'res1/(N-3) ; vcov1=sighat1*inv(z1'z1);

t_inst=y2hat[2]/sqrt(vcov1[2,2]) ;

Fsim[estim,.]=t_inst^2 ;


Rsqr1=1-res1'res1/((y2vec-ybarvec)'(y2vec-ybarvec)) ;

Rsqr2=1-res2'res2/((y2vec-ybarvec)'(y2vec-ybarvec)) ;

PartR2=Rsqr1-Rsqr2 ; @Partial R-square for the instrument@

Rsqrsim[estim,.]=PartR2 ;

Xhat=ones(N,1)~y2hat~x1;


/*Second stage*/

alp_hat=inv(Xhat'Xhat)*XHat'y1vec ;

results[estim,.]=alp_hat' ;

endo ; " " ;

/*results of last estimation*/

/*"estimated parameters:" alp_hat'; */

"true param are:" 0 " " alpha1 " " alpha2 ;


/* averages and std dev of estimates */

meanres=meanc(results); stdres=stdc(results) ;

"mean estimates" meanres' ;
 "std of estimates" stdres' ;
"Partial R-squares for instrument" meanc(Rsqrsim) ;

"F test for instrument:" meanc(Fsim);


output off ;