/* Example 17.5.1 test C*beta=d0  */
libname sriv'/u/mengd/sasuser/ssm';
/* n=24, k=3; B=100, alpha=95% */
/* C={1 2 -1, 3 7 4 }; d0={6,18} */
proc iml;
speed={38.8, 31.5, 10.6, 16.1,  7.7, 8.3,  8.5,  11.1,  8.6, 11.1,
        9.8,  7.8, 31.8, 31.6, 34.0, 28.9, 28.8, 10.5, 12.3, 13.2,
       11.4, 11.2, 10.3, 11.4};
density={20.4, 27.4, 106.2, 80.4, 141.3, 130.9, 121.7, 106.5, 
        130.5, 101.1, 123.9, 144.2, 20.9, 30.8, 26.5, 35.7, 30.0,
        106.2, 97.0, 90.1, 106.7, 99.3, 107.2, 109.1};
y=sqrt(speed);
x1=J(24,1,1);
x2=density;
x3=density#density;
x=x1||x2||x3;
e=(I(24)-x*inv(x`*x)*x`)*y;
create sriv.vehicle var{y x1 x2 x3 e };
append;
J=1/24*J(24,1,1);
eb=J`*e;
create sriv.tsq var{tsq};
b=inv(x`*x)*x`*y;
C={1 2 -1, 3 7 4};
d0={6, 18};
ssq=y`*(I(24)-x*inv(x`*x)*x`)*y/(24-3);
Tsq=(c*b-d0)`*inv(c*inv(x`*x)*c`)*(c*b-d0)/ssq;
append from Tsq;
hiimax=max(diag(x*inv(x`*x)*x`));
print c  b d0, eb hiimax;
quit;
data sriv.t;
input tistarsq;
cards;
;
run;
%macro ran(k= );
%do n=1 %to &k;
data rands&n;
retain seed &n&n;
do i=1 to 24;
call rantbl(seed,
 1/24, 1/24, 1/24, 1/24, 1/24, 1/24, 1/24, 1/24, 1/24, 1/24,
 1/24, 1/24, 1/24, 1/24, 1/24, 1/24, 1/24, 1/24, 1/24, 1/24,
 1/24, 1/24, 1/24, 1/24,iid);
output;
 end;
run;
/* Draw a bootstrap samples with replacement from e */
proc iml;
use rands&n;
read all var{iid} into ids;
use sriv.vehicle;
read all var{e} into e;
sume=0;
create sriv.et var{estar};
do i=1 to 24;
    estar=e[ids[i,],];
    append;
    sume=sume+estar;
end;
ebstar=sume/24;
/*  Calculate sistarsq */
use sriv.vehicle;
read all into z;
x=z[,2:4];
use sriv.et;
read all into et;
sistarsq=24/((24-1)*(24-3))*(et`*et-24*ebstar**2);
/* Calculate  Tistarsq */
C={1 2 -1, 3 7 4};
Tistarsq=et`*x*inv(x`*x)*c`*inv(c*inv(x`*x)*c`)*c*
inv(x`*x)*x`*et/sistarsq;
edit sriv.t ;
append from Tistarsq;
quit;
%end;
/* Order the values of Tistarsq  */
proc iml;
use sriv.t;
read all into tistarsq;
ord=rank(tistarsq);
do i=1 to 100;
   rr=ord[i,];
   if rr=96 then
   Tstarsq=tistarsq[i,];
end;
/* test */
use sriv.tsq;
read all into t00;
Tsq=t00[1,1];
print Tsq Tstarsq;
if Tsq >Tstarsq then print "reject";
             else print "accept";
quit;
%mend ran;
%ran(k=100)