/*SIMULTANEUS EQUATIONS AS IN CLASS */

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

N = 150;

let alpha1= .5; alpha2=.3 ; alpha3=.3 ; alpha4=.2 ; alpha5=1 ;
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]=alpha3 ;
Bmat[1,2]=alpha4 ;
Bmat[2,1]=0 ;
Bmat[2,2]=alpha5 ;



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

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

results = zeros(Nsim,3) ;

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 ;
RHS1=ones(N,1)~y2vec~x1 ;

/*First Stage*/

PUT IN THE FIRST STAGE HERE

/*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 */ "mean estimates"
meanc(results); "std of estimates" stdc(results) ;

" " ;


output off ;