/*Example 10.2.1---*/
proc iml;
y1={21,21,20.5,23.5,21.5,20,21.5,23,20,16.5,24.5};
y2={20,21.5,24,24.5,23,21,22.5,23,21,19,25};
y3={21.5,24,24.5,25,22.5,21,23,23.5,22,19,28};
y4={23,25.5,26,26.5,23.5,22.5,25,24,21.5,19.5,28};
n=11; p=4;  f=n-1; m=2;
y=y1||y2||y3||y4;
yb=y`*J(n,1,1/n);
s= (y`-yb*J(1,n,1))*(y`-yb*J(1,n,1))`/(n-1);
B={ 1   1   1   1,
   -3  -1   1   3};
c={ 1  -1  -1   1,
  -1   3  -3   1};
F0=(f-p+m+1)*n/(f*(p-m))*yb`*c`*inv(c*s*c`)*c*yb;
Fa=finv(1-0.05, p-m,f-p+m+1);
print  F0 Fa;
if  F0>Fa then print "F0 > Fa: reject";
          else print "F0 < Fa: accept";
psi=inv(B*inv(s)*B`)*B*inv(s)*Yb;
Tpmsq=n*yb`*c`*inv(c*s*c`)*c*yb;
Tasq=f*m/(f-p+1)*finv(1-0.05,m,f-p+1);
a={1 2}`;
e1=sqrt(Tasq*(1+Tpmsq/f)*(a`*inv(B*inv(s)*B`)*a)/n);
low1=a`*psi-e1;
upp1=a`*psi+e1;
print low1 upp1;
U={0 1}; k=1; d=U*psi; uu=U`;
F00=(f-k-p+m+1)*n/(f*k)*d`*inv(U*inv(B*inv(s)*B`)*U`)*d/(1+Tpmsq/f);
Faa=Finv(1-0.05,k,f-k-p+m+1);
Taasq=f*k/(f-p+m-k+1)*finv(1-0.05,k,f-p+m-k+1);
e2=1/sqrt(n)*sqrt(Taasq*(1+Tpmsq/f)*uu`*inv(B*inv(s)*B`)*uu);
low2=uu`*psi-e2;
upp2=uu`*psi+e2;
print F00 Faa;
if  F00>Faa then print "F00 > Faa: reject";
          else print "F00 < Faa: accept";
print low2 upp2;
quit;
             F0            Fa
      0.0520224     4.2564947
             F0 < Fa: accept
           LOW1        Upp1
     (21.288957,  26.095369)
           F00           Faa
     44.304298     5.3176551
           F00 > Faa: reject
           LOW2         UPP2
     (0.3152223,   0.6494201)