*SET SEED 1234321.
* For testing Ho: H*beta = delta.
* Note: "Contruct c:\Hmatrix.sav" .
*Note: "Construct c:\delta.sav".
*This does 1000 bootstraps.
*MACRO CALL - bootstr2 y  x1 to xp.

*NOTE: Standardize the independents.

bootstr2 mpg_1 zeng  zhors zweig zaccel.

set printback =none.
set messages = none.
set mprint = off.

DEFINE bootstr2(!positional !TOKENS(1) /!positional !cmdend).
set mxloop 99999.
matrix.
* Get Y.
Get y /variables = !1/file =* /names = namy /missing = omit.

*Get X.
get xindep /variables = !2/file = * /names = namx /missing = omit.
*print x.
compute nr = nrow(xindep).
compute jx =make(nr,1,1).
compute x = {jx,xindep}.
compute ncx = ncol(x).

compute ny = nrow(y).
compute nr1 = {nr,ny}.
* if nr <> ny then missing values are messing you up.

do if (nr <> ny).
print nr1 /Title = "Missing Value Problem".
else.
print nr.
compute mx1 = 0.
compute my = 0.

Compute betahat = ginv(t(x)*x)*T(x)*y.


print betahat.


compute dfnum = trace(ginv(t(x)*x)*(T(x)*x)).
*print dfnum.
compute dfdenom = nr - dfnum.

print dfdenom.
compute jnr = make(nr,1,1).
compute mse = T(y-x*betahat)*(y-x*betahat)/dfdenom.

print mse.

* Get hypothesis matrix   h*beta = delta..
get h /variables = all /file = 'C:\hmatrix.sav'.
compute df_h = nrow(h).
get delta /variables = all /file= 'C:\delta.sav'.
print df_h.	
print h.
print delta.
Compute mstreat = t(h*betahat-delta)*GINV(h*ginv(T(x)*x)*T(h))*(h*betahat-delta)/df_h.

compute nrx = nrow(t(x)*x).
do if nrx >1.
compute iden = mdiag(make(nrx,1,1)).
else.
compute iden =1.
end if.

print mstreat.

compute FCal = mstreat/mse.
print Fcal.


* ______________Now start boot strap.
compute xtx_i = ginv(t(x)*x).

* sample the data.


compute cnt = 0.
* do big loop.
compute bign = 1000.
compute estimb = make(bign,ncx+1,0).
compute b1star = make(bign,1,1).
compute jbig = make(bign,1,1).
compute fstar = make(bign,2,0).
compute bb = make(bign,1,0).
loop i = 1 to bign.

compure rstar = make(nr,1,0).

*get bootstrap sample of x.

compute xstar= make( nr,ncx,1).
compute ystar1= make( nr,1,0).


loop j = 1 to nr.
compute id1 = trunc(nr*uniform(1,1)+1).
compute ystar1(j,1) = y(id1,1).
do if ncx >1.
loop kl = 2 to ncx.
compute xstar(j,kl) = x(id1,kl).
end loop.
end if.
end loop.
*print xstar.

save xstar/outfile = 'c:\xstar.sav'.
compute ystar = ystar1.
compute x1 = xstar.



compute betastar =ginv(T(x1)*x1)*t(x1)*ystar.
*print betastar.

compute msestr = T(ystar - x1*betastar)*(ystar-x1*betastar)/dfdenom.
*print msestr.
compute xtx_i = ginv(T(x1)*x1).

compute mst_str = t(h*betastar-delta)*ginv(h*xtx_i*t(h))*(h*betastar-delta)/df_h.
compute mx1 = 0.

compute fstar(i,1)  =T(h*betastar-h*betahat)*ginv(h*xtx_i*t(h))*(h*betastar-h*betahat)/(df_h*msestr).

compute fstar(i,2) = fcal.

do if fstar(i,1) > fcal.
compute cnt = cnt+1.
end if.

loop j1 = 1 to ncx.
compute estimb(i,j1) = betastar(j1,1).
compute estimb(i,ncx+1) = msestr.
end loop.

end loop.

compute sig = cnt/bign.
print sig /Format "F7.4" /Title "Significance Level".
compute nams = {"Intercpt",namx, "mse"}.
save estimb/outfile = 'c:\estimb.sav' /names = nams.
save fstar/outfile = 'c:\fstar.sav'.

end if.


end matrix.


GET
  FILE='C:\estimb.sav'.

GRAPH
  /HISTOGRAM(NORMAL)=mse .
GRAPH
  /HISTOGRAM(NORMAL)=intercpt .

!do !q !in (!2).

GRAPH
  /HISTOGRAM(NORMAL)=!q .
!doend.



!enddefine.