*SET SEED 1234321.
* Implement Sherman's bootstrap for GLM.

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

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

*Get X.
get x /varaibles = !2/file = */missing = omit.
*print x.
compute nr = nrow(x).
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 = 200.
compute b1star = make(bign,1,1).
compute jbig = make(bign,1,1).
compute fstar = make(bign,1,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 bb(i,1) = betastar(2,1).
*end if.
end loop.


*compute msebb.
*compute meanbb = t(jbig)*bb/bign.
*compute msebb = t(bb-meanbb*jbig)*(bb-meanbb*jbig)/(bign-1).
*compute fcal_bb = mst_str/msebb.
*save bb/outfile = 'c:\bb.sav'.
* sort fstar.
*print fstar.
loop iloop = 1 to bign.
loop jloop = 2 to bign - iloop+1.

do if (fstar(jloop-1,1) > fstar(jloop,1)).
compute atemp = fstar(jloop-1,1) .
compute btemp = fstar(jloop,1) .
compute fstar(jloop-1,1) = btemp.
compute fstar(jloop,1) = atemp.
end if.

end loop.
end loop.
*print fstar.
*print fstar(5,1).
*print fstar(95,1).
* compute cnt = cnt/bign.



*print cnt.
save fstar/outfile = 'c:\fstar.sav'.
* *print fcal.
compute idfg = bign*.975.
compute idfg1 = bign*.95.
*print idfg.
*print idfg1.
compute resul = {fcal,df_h,dfdenom,fstar(idfg1,1),fstar(idfg,1)}.
save resul /outfile = 'c:\resul.sav'.

end if.


end matrix.

!enddefine.



define  bigrun (!positional !cmdend).


get
file ='c:\resul0.sav'.
SAVE OUTFILE='C:\resul1.sav'
 /COMPRESSED.
SAVE OUTFILE='C:\concord3.sav'
 /COMPRESSED.



!do !i = 1 !to !1.
GET
  FILE='c:\downloads\concord1.sav'.
compute id11 = rv.uniform(0,1) -0.1.
execute.
compute x1 = rv.normal(0,1).
compute x2 = rv.normal(0,1).
compute x3 = rv.normal(0,1).
do if   (id11  GT  0)  .
*compute water80 = rv.normal(0,1).
 COMPUTE water81 = 6+ x1 + 3*x2 + x3+rv.normal(0,1).
 else .
compute   x1 =3*(rv.chi(1) -1).
compute   x2 =3*(rv.chi(1) -1).
compute   x3 =3*(rv.chi(1) -1).
COMPUTE water81 = 6+ x1 + 3*x2 + x3+ 3*(rv.chi(1) -1).

end if. 

*compute water80 = rv.normal(0,1).
*compute x2 = rv.normal(0,1).
*COMPUTE water81 = 6+  3*(rv.chi(1)-1).





compute c = 1.
EXECUTE .
save outfile= 'c:\temp3.sav'
 /COMPRESSED.


*ADD FILES /FILE=*
 /FILE='C:\concord3.sav'.
*EXECUTE.


*SAVE OUTFILE='C:\concord3.sav'
 /COMPRESSED.

*GET
  FILE='C:\temp1.sav'.


 bootstr2 water81  c x1 x2 x3 .

GET
  FILE='C:\resul1.sav'.
ADD FILES /FILE=*
 /FILE='C:\resul.sav'.
EXECUTE.
save outfile = 'c:\temp1.sav'
 /COMPRESSED.
SAVE OUTFILE='C:\resul1.sav'
 /COMPRESSED.

!doend

!enddefine.

bigrun     100.

COMPUTE id = $casenum .
EXECUTE .


FILTER OFF.
USE ALL.
SELECT IF(id>1).
EXECUTE .




compute idf1 = idf.f(.95,col2,col3).
EXECUTE .



compute ct1 = 0.
if (Col1 > col5) ct1 =1.
execute.

compute ct2= 0.
IF (col1 > idf1) ct2 = 1 .
EXECUTE .


DESCRIPTIVES
  VARIABLES=ct1 ct2
  /STATISTICS=MEAN   .