SET SEED 1221999999.
* For testing Ho: H*beta = delta.
* Note: "Contruct c:\Hmatrix.sav" .
*Note: "Construct c:\delta.sav".
*This does 1000 bootstraps.
*MACRO CALL - bootstr2 y  nc x1 to xp.
* NC = number of levels of the factor.
*NOTE: Standardize the independents.
*NOTE: YOU NEED THE DESIGN matrix -- use no intercept model put factors first.

boot_anc y1 2 p1 p2 p3 p4.

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

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

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

loop kk1 =  1 to ntreat.
compute samsize(kk1,1) = samsize1(1,kk1).
end loop.

print samsize.

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).
compute sumnn = -samsize(1,1).
compute j11 = 0.
loop j = 1 to ntreat.
compute nn = samsize(j,1).
compute sumnn = sumnn +nn.
loop kk1 = 1 to nn.
compute j11 = j11 +1.
compute id1 = trunc(nn*uniform(1,1)+1) + sumnn.
compute ystar1(j11,1) = y(id1,1).

loop kl = 1 to ncx.
compute xstar(j11,kl) = x(id1,kl).
end loop.

end loop.
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 j13 = 1 to ncx.
compute estimb(i,j13) = betastar(j13,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 = {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 .

!do !q !in (!3).

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



!enddefine.