/*SIMULTANEUS EQUATIONS AS IN CLASS */
new;
 output file = simu.out reset ;
library optmum;

N = 100;

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);

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


bootrep=200; "Number of bootstrap replications:" bootrep ;

Nsim =1 ; "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) ;

/*generate linear relation */


y1vec=xmat*(alpha1|alpha2|alpha3)+u1 ;

alphhat=inv(xmat'xmat)*xmat'y1vec ;

resid=y1vec-xmat*alphhat ;

sigmhat=resid'resid/(N-3) ;

stderr=sqrt(diag(sigmhat*inv(xmat'xmat))) ;

 stdboot=bootstr(bootrep,xmat,alphhat,resid) ;

"stderr" stderr'; "" ;

"stdboot" stdboot' ;

/*results[estim,.]=alp_hat' ;*/

endo ; " " ;

/*

/* averages and std dev of estimates */ "mean estimates"
meanc(results); "std of estimates" stdc(results) ;

" " ;

*/


proc bootstr(bootrep,xmat,alphhat,resid);
  local uni,pick,bootres,i;


bootres=zeros(bootrep,rows(alphhat)) ;

i=0;

do while i<bootrep ;

i=i+1 ;

/* Draw errors from empirical distribution */
uni=rndu(N,1);/*uniform on [0,1] */

uni=ones(N,1)+N*uni;/*uniform on [1,N+1] */

pick=trunc(uni) ;/*equal probability on 1,...,N */




y1vec=xmat*alphhat+resid[pick,1] ;

/*"pick" pick';*/

bootres[i,.]=(inv(xmat'xmat)*xmat'y1vec)' ;

endo;

retp(stdc(bootres)) ;

endp; output off ;