---------------------------------------------------------------------
****INVOKE THE IML PROGRAM AND DEFINE THE MODULE WJGLM****;
PROC IML;
*RESET NONAME;

****DEFINE MODULES TO PERFORM ALL CALCULATIONS****;
**DEFINE MODULE FOR INITIAL SPECIFICATIONS****;
START INITIAL(C,U,Y,R,X) GLOBAL(NX,NTOT,BOBS,WOBS,WOBS1);
IF NROW(U)=0 THEN U=I(NCOL(Y));
IF NROW(C)>NCOL(C) THEN PRINT
  'ERROR:NUMBER OF ROWS OF C EXCEEDS NUMBER OF COLUMNS';
IF NCOL(U)>NROW(U) THEN PRINT
  'ERROR: NUMBER OF COLUMNS OF U EXCEEDS NUMBER OF ROWS';
DO I=1 TO NCOL(NX);
  X1=J(NX[I],1,I);
  IF I=1 THEN X=X1;
  ELSE X=X//X1;
END;
X=DESIGN(X);
NTOT=NROW(Y);
WOBS=NCOL(Y);
BOBS=NCOL(X);
WOBS1=WOBS-1;
R=C@U`;*print r;
FINISH;

****DEFINE MODULE TO COMPUTE LEAST SQUARES OR TRIMMED MEANS****;
START MNMOD (Y,OPT1,BHAT,BHATW,MUHAT,YT,DF) GLOBAL(BOBS,WOBS,WOBS1,NTOT,NX,PER,X);
IF OPT1=0 THEN DO;
  BHAT=INV(X`*X)*X`*Y;*print bhat;
  BHATW=BHAT;
  YT=Y;
  DF=NX-1;
END;
IF OPT1=1 THEN DO;
  BHAT=J(BOBS,WOBS,0);
  BHATW=BHAT;
  YT=J(NTOT,WOBS,0);
  DF=J(1,BOBS,0);
  F=1;
  M=0;
  DO J=1 TO NCOL(NX);
    SAMP=NX[J];
    L=M+SAMP;
    G=INT(PER#SAMP);
    DF[J]=SAMP-2#G-1;
    DO K=1 TO NCOL(Y);
    TEMP=Y[F:L,K];
    NV=TEMP;
    TEMP[RANK(NV),]=NV;
    TRIMY=TEMP[G+1:SAMP-G,];
    TRIMMN=SUM(TRIMY)/(DF[J]+1);
    BHAT[J,K]=TRIMMN;
    MINT=MIN(TRIMY);
    MAXT=MAX(TRIMY);
    DO P=1 TO NROW(NV);
          IF NV[P]<=MINT THEN NV[P]=MINT;
          IF NV[P]>=MAXT THEN NV[P]=MAXT;
        END;
    YT[F:L,K]=NV;
    WINMN=SUM(NV)/SAMP;
    BHATW[J,K]=WINMN;
   END;
   M=L;
   F=F+NX[J];
  END;
END;
MUHAT=SHAPE(BHAT,BOBS#WOBS);
FINISH;

***DEFINE MODULE TO COMPUTE SIGMA MATRIX****;
START SIGMOD (YT,BHATW,DF,SIGMA) GLOBAL(BOBS,WOBS,WOBS1,NTOT,NX,PER,X);
SIGMA=J(WOBS#BOBS,WOBS#BOBS,0);
DO I=1 TO BOBS;
  SIGB=(YT#X[,I]-X[,I]*BHATW[I,])`*(YT#X[,I]-X[,I]*BHATW[I,])/((DF[I]+1)#DF[I]);
  F=I#WOBS-WOBS1;
  L=I#WOBS;
  SIGMA[F:L,F:L]=SIGB;
END;
FINISH;

****DEFINE MODULE TO COMPUTE TEST STATISTIC****;
START TESTMOD(SIGMA,MUHAT,R,DF,FSTAT,DF1,DF2,C,U) GLOBAL(BOBS,WOBS,WOBS1,NTOT,NX,PER,X);
T=(R*MUHAT)`*INV(R*SIGMA*R`)*(R*MUHAT);*print r, muhat, t;
A=0;
IMAT=I(WOBS);
DO I=1 TO BOBS;
  QMAT=J(BOBS#WOBS,BOBS#WOBS,0);
  F=I#WOBS-WOBS1;
  L=I#WOBS;
  QMAT[F:L,F:L]=IMAT;
  PROD=(SIGMA*R`)*INV(R*SIGMA*R`)*R*QMAT;
  A=A+(TRACE(PROD*PROD)+TRACE(PROD)**2)/DF[I];
END;
A=A/2;
DF1=NROW(R);
DF2=DF1#(DF1+2)/(3#A);
CVAL=DF1+2#A-6#A/(DF1+2);
FSTAT=T/CVAL;
if NCOL(NX)=1 | (NROW(C)=1 & C[1,+]=1) then do;
  pp=NCOL(u);
  df=c*df`;
  FSTAT=(df-pp+1)*T/(df*pp);*print fstat;
  df2=(df-pp+1);
end;
FINISH;
****DEFINE MODULES TO PERFORM BOOTSTRAP****;
***DEFINE MODULE TO GENERATE BOOTSTRAP DATA AND CENTRE DATA****;
START BOOTDAT (Y,BHAT,YB)   GLOBAL(BOBS,WOBS,WOBS1,NTOT,NX,PER,X,SEED);
  F=1;
  M=0;
  DO J=1 TO BOBS;
   L=M+NX[J];
   TEMP=Y[F:L,];
   BVAL=TEMP;
   DO P=1 TO NROW(TEMP);
     RVAL=UNIFORM(SEED);
         BVAL[P,]=TEMP[CEIL(NROW(TEMP)#RVAL),];
   END;
   IF J=1 THEN YB=BVAL;
   ELSE YB=YB//BVAL;
   M=L;
   F=F+NX[J];
 END;
****CENTRE THE BOOTSTRAP DATA****;
M=0;
F=1;
DO I=1 TO BOBS;
  L=M+NX[I];
  MVAL=BHAT[I,];
  DO Q=F TO L BY 1;
    YB[Q,]=YB[Q,]-MVAL;
  END;
  M=L;
  F=F+NX[I];
END;
FINISH;

****DEFINE MODULE TO COMPUTE BOOTSTRAP STATISTIC****;
START BOOTSTAT(YB,OPT1,R,FSTATB,C,U) GLOBAL(BOBS,WOBS,WOBS1,NTOT,NX,PER,X,SEED);
  CALL MNMOD(YB,OPT1,BHATB,BHATBW,MUHATB,YTB,DFB);
  CALL SIGMOD(YTB,BHATBW,DFB,SIGMAB);
  CALL TESTMOD(SIGMAB,MUHATB,R,DFB,FSTATB,DF1B,DF2B,C,U);
FINISH;

****COMPUTE WELCH-JAMES STATISTIC****;
START WJGLM;
CALL INITIAL(C,U,Y,R,X);
CALL MNMOD(Y,OPT1,BHAT,BHATW,MUHAT,YT,DF);
CALL SIGMOD(YT,BHATW,DF,SIGMA);
CALL TESTMOD(SIGMA,MUHAT,R,DF,FSTAT,DF1,DF2,C,U);
IF OPT2=1 THEN DO;
   DO SIMLOOP=1 TO NUMSIM;
      CALL BOOTDAT(Y,BHAT,YB);
      CALL BOOTSTAT(YB,OPT1,R,FSTATB,C,U);
          IF SIMLOOP=1 THEN FMAT=FSTATB;
        ELSE FMAT=FMAT//FSTATB;
   END;
TEMPB=FMAT;
FMAT[RANK(FMAT)]=TEMPB;
END;
***CALCULATE SIGNIFICANCE LEVEL****;
RESULTS=J(4,1,0);
RESULTS[1]=FSTAT;
RESULTS[2]=DF1;
RESULTS[3]=DF2;
IF OPT2=0 THEN RESULTS[4]=1-PROBF(RESULTS[1],DF1,DF2);
IF OPT2=1 THEN DO;
  AVEC=(FSTAT<=FMAT);
  PVAL=SUM(AVEC)/NUMSIM;
  RESULTS[4]=PVAL;
END;
****PRINT WELCH-JAMES RESULTS****;
PRINT'WELCH-JAMES APPROXIMATE DF SOLUTION';
IF OPT1=0 THEN PRINT'LEAST SQUARES MEANS & VARIANCES';
IF OPT1=1 THEN PRINT'TRIMMED MEANS & WINSORIZED VARIANCES';
IF OPT1=1 THEN DO;
  PRINT 'PERCENTAGE OF TRIMMING:';
  PRINT PER[FORMAT=4.2];
END;
IF OPT2=0 THEN PRINT'F DISTRIBUTION CRITICAL VALUE';
IF OPT2=1 THEN PRINT'BOOTSTRAP CRITICAL VALUE FOR SINGLE TEST STATISTIC';
IF OPT2=1 THEN DO;
  PRINT 'NUMBER OF BOOTSTRAP SAMPLES:';
  PRINT NUMSIM[FORMAT=4.0],;
  PRINT 'STARTING SEED:';
  PRINT SEED[FORMAT=15.0],;
END;
PRINT'CONTRAST MATRIX:';
PRINT R[FORMAT=4.1],;
MUHAT=MUHAT`;
PRINT 'MEAN VECTOR:';
PRINT MUHAT[FORMAT=10.4],;
PRINT 'SIGMA MATRIX:';
PRINT SIGMA[FORMAT=10.4],;
RESLAB={"TEST STATISTIC" "NUMERATOR DF" "DENOMINATOR DF" "P-VALUE"};
PRINT 'SIGNIFICANCE TEST RESULTS:';
PRINT RESULTS[ROWNAME=RESLAB FORMAT=10.4]/;
FINISH;

****DEFINE MODULE TO COMPUTE BOOTSTRAP RESULTS WITH FWR CONTROL****;
START BOOTCOM;
CALL INITIAL(C,U,Y,R,X);
CALL MNMOD(Y,OPT1,BHAT,BHATW,MUHAT,YT,DF);
CALL SIGMOD(YT,BHATW,DF,SIGMA);
DO I=1 TO NROW(R);
CM=R[I,];
CALL TESTMOD(SIGMA,MUHAT,CM,DF,FSTAT,DF1,DF2,C,U);
IF I=1 THEN DO;
   CMMAT=FSTAT;
   DF1MAT=DF1;
   DF2MAT=DF2;
END;
IF I>1 THEN DO;
    CMMAT=CMMAT||FSTAT;
    DF1MAT=DF1MAT||DF1;
    DF2MAT=DF2MAT||DF2;
END;
END;
DO SIMLOOP=1 TO NUMSIM;
CALL BOOTDAT(Y,BHAT,YB);
DO Q=1 TO NROW(R);
  CM=R[Q,];
  CALL BOOTSTAT(YB,OPT1,CM,FSTATB,C);
  IF Q=1 THEN FROW=FSTATB;
    ELSE FROW=FROW||FSTATB;
END;
IF SIMLOOP=1 THEN FMAT=FROW;
  ELSE FMAT=FMAT//FROW;
END;
FMAX=J(NUMSIM,1,0);
DO K=1 TO NUMSIM;
 FMAX[K]=MAX(FMAT[K,]);
END;
TEMPB=FMAX;
FMAX[RANK(FMAX)]=TEMPB;
RESULTS=J(3,NROW(R),0);
RESULTS[1,]=CMMAT;
RESULTS[2,]=DF1MAT;
RESULTS[3,]=DF2MAT;
QCRIT=ROUND((1-ALPHA)#NUMSIM);
CRITV=FMAX[QCRIT];
****PRINT RESULTS****;
PRINT'WELCH-JAMES APPROXIMATE DF SOLUTION';
IF OPT1=0 THEN PRINT'LEAST SQUARES MEANS & VARIANCES';
IF OPT1=1 THEN PRINT'TRIMMED MEANS & WINSORIZED VARIANCES';
IF OPT1=1 THEN DO;
  PRINT 'PERCENTAGE OF TRIMMING:';
  PRINT PER[FORMAT=4.2];
END;
PRINT'BOOTSTRAP CRITICAL VALUE FOR FWR CONTROL';
PRINT 'NUMBER OF BOOTSTRAP SAMPLES:';
PRINT NUMSIM[FORMAT=4.0];
PRINT 'STARTING SEED:';
PRINT SEED[FORMAT=15.0],;
PRINT 'SIGNIFICANCE LEVEL:';
PRINT ALPHA[FORMAT=3.2],;
PRINT'CONTRAST MATRIX:';
PRINT R[FORMAT=4.1],;
MUHAT=MUHAT`;
PRINT 'MEAN VECTOR:';
PRINT MUHAT[FORMAT=10.4],;
PRINT 'SIGMA MATRIX:';
PRINT SIGMA[FORMAT=10.4],;
RESLAB={"TEST STATISTIC" "NUMERATOR DF" "DENOMINATOR DF" "SIGNIFICANCE"};
PRINT 'SIGNIFICANCE TEST RESULTS:';
PRINT RESULTS[ROWNAME=RESLAB FORMAT=10.4],;
PRINT 'CRITICAL VALUE:';
PRINT CRITV[FORMAT=5.2]/;
FINISH;
Y={41, 38.4, 24.4, 25.9, 21.9, 18.3, 13.1, 27.3, 28.5, -16.9,
   26.0, 17.4, 21.8, 15.4, 27.4, 19.2, 22.4, 17.7, 26.0, 29.4, 21.4, 26.6,
   22.7, 10.1, 6.1, 20.4, 7.3, 14.3, 15.5, -9.9, 6.8, 28.2, 17.9, -9.0,
  -12.9, 14.0, 6.6, 12.1, 15.7, 39.9, -15.9, 54.6, -14.7, 44.1, -9.0};
nx={23 22};
c={1 -1};
opt1=1;
per=.20;
opt2=1;
seed=38190;
numsim=599;
run wjglm;
quit;