Cusum ARLs Using Markov Chain

 /****************************************************************/
 /*       S A S   S A M P L E   L I B R A R Y                    */
 /*                                                              */
 /*    NAME: CUSARLM                                             */
 /*   TITLE: Cusum ARLs Using Markov Chain                       */
 /* PRODUCT: QC                                                  */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: SQC                                                 */
 /*   PROCS: IML                                                 */
 /*    DATA:                                                     */
 /*                                                              */
 /*     REF: Hawkins, D. M. (1992).  "Evaluation of Average Run  */
 /*          Lengths of Cumulative Sum Charts for an Arbitrary   */
 /*          Data Distribution."  Communications in Statistics-- */
 /*          Simulation.  21(4) 1001-1020.                       */
 /*                                                              */
 /*   NOTES: This example provides a SAS/IML program to compute  */
 /*          the ARL of cusum schemes using a Markov chain       */
 /*          method.  The program is a translation of a program  */
 /*          due to Hawkins (1992).                              */
 /*                                                              */
 /*          To support distributions other than the normal,     */
 /*          replace the two calls to probnorm and make          */
 /*          modifications to denrat as described by Hawkins.    */
 /*                                                              */
 /****************************************************************/

  proc iml;

  start arlcomp;

  /* input parameters */
  di      =  5.0;     /* decision interval h    (H=)     */
  ref     =  0.5;     /* reference value k      (K=)     */
  denrat  =  0.0;     /* common denominator              */
  arg2    =  0.0;     /* distribution parameter (DELTA=) */
  winsor  =  10.0;    /* winsorization value    (W=)     */
  eps     =  0.001;

  /* initialize output parameters */
  regarl  = 0.0;
  firarl  = 0.0;

  /* constants */
  esterr  = 0.0;
  maxhlf  = 7;
  maxpow  = 2##maxhlf; /* 2**maxhlf                */
  toler   = 1.0e-12;   /* subtractive cancellation */

  /* working matrices */
  tran    = j( maxpow + 2, maxpow + 2, 0 );
  answer  = j( maxpow + 2, 1, 0 );
  copy    = j( maxpow + 2, maxpow + 2, 0 );
  phitab  = j( 4#maxpow + 5, 1, 0 );
  richar1 = j( maxhlf, maxhlf, . );
  richar2 = j( maxhlf, maxhlf, . );

  /* check parameters for feasibility */
  if ( di <= 0.0 )
  then do;
     print , "ERROR: Nonpositive decision interval";
     stop;
     end;
  if ( winsor <= ref )
  then do;
     print , "ERROR: Winsorization value less than reference value";
     stop;
     end;
  if ( denrat < 0 )
  then do;
     print , "ERROR: Common denominator is less than zero.";
     stop;
     end;

  /* are the data continuous? */
  if ( denrat = 0 ) then realx  = 1; else realx = 0;

  /* other initialization */
  mbig   = 1;
  idenrn = 0;

  /* subscript adjustments */
  phadj   = 2 # maxpow + 3;

  if ( realx = 0 )
  then do;
     if ( ref <= 0.0 )
        then knumer = int( ref*denrat - 0.5 );
     else
        knumer = int( ref#denrat + 0.5 );
     idenrn = int( denrat );
     mbig   = int( di#denrat - 0.5 );
     iwintg = int( winsor#denrat + 0.5 );
     mb1    = mbig + 1;
     if ( mbig >= maxpow )
     then do;
        print , "ERROR: di#denrat exceeds maxpow";
        stop;
        end;
     deviat  = abs( knumer - ref#denrat    ) +
               abs( mbig - di#denrat + 1   ) +
               abs( iwintg - winsor#denrat ) ;
     if ( deviat > 0.01 )
     then do;
        print , "ERROR: di, ref or winsor not a multiple of 1/denrat";
        stop;
        end;
     end;

  answer[ 1 ] = 0;

  /* set up coefficient matrix of linear equations */
  do mloop = 1 to maxhlf;

     if ( realx = 1 )
     then do;
        mbig    = 2 # mbig;
        mb1     = mbig + 1;
        delta   = di / mbig;
        mtwo    = 2 # mbig;
        lastnx  = mtwo + 2;
        fact    = 4.0;
        end;
     else lastnx = mbig;

     do i= -lastnx - idenrn to lastnx;

        iadj = i + phadj;

        if ( realx = 1 )
        then do;
           fact = 5.0 - fact;
           arg  = 0.5 # delta # i + ref;
           if ( arg > winsor ) then
              phitab[ iadj ] = fact / 6.0;
           else
              phitab[ iadj ] = probnorm(arg - arg2) # fact / 6.0;
           end;

        else do;
           iarg = i + knumer;
           if ( mod( iarg, idenrn ) = 0 )
           then do;
              arg = iarg / idenrn;
              phitab[ iadj ] = 1.0;
              if ( arg < winsor ) then
                 phitab[ iadj ] = probnorm(arg - arg2);
              end;
           else do;
              phitab[ iadj ] = phitab[ iadj - 1 ];
              if ( iarg >= iwintg ) then
                 phitab[ iadj ] = 1.0;
              end;
           end;
        end;

     do j = 1 to mbig;

        if ( realx = 1 )
        then do;

           inx = 2 # j;

           tran[ j+1, 1 ] = phitab[ -inx +     phadj ] +
                            phitab[ -inx + 1 + phadj ] +
                            phitab[ -inx + 2 + phadj ] ;

           tran[ 1, j+1 ] = 6.0 # ( phitab[ inx +     phadj ] -
                                    phitab[ inx - 2 + phadj ]   );

           tran[ mb1 + 1, j + 1 ] = 0.0;

           do i = 1 to mbig;
              inx = 2 # ( i - j );
              tran[ j+1, i+1 ] = ( phitab[ inx + 2 + phadj ] -
                                   phitab[ inx     + phadj ]   ) +
                                 ( phitab[ inx + 1 + phadj ] -
                                   phitab[ inx - 1 + phadj ]   ) +
                                 ( phitab[ inx     + phadj ] -
                                   phitab[ inx - 2 + phadj ]   ) ;
              end;
           end;

        else do;

           tran[ j+1, 1 ] = phitab[ -j + phadj ];
           tran[ 1, j+1 ] = phitab[ j + phadj ] - phitab[ j-1 + phadj ];

           do i = 1 to mbig;
              tran[ j+1, i+1 ] = phitab[ i-j     + phadj ] -
                                 phitab[ i-j - 1 + phadj ] ;
              end;
           end;
        end;

     tran[ 1, 1 ] = phitab[ phadj ];
     if ( realx = 1 ) then tran[ 1, 1 ] = 6.0 # phitab[ phadj ];
     tran[ mb1 + 1, mb1 + 1 ] = 1.0;
     tran[ mb1 + 1, 1       ] = 0.0;

     do j = 0 to mb1;

        copy[ j + 1, mb1 + 1 ] = 1.0;
        do i = 0 to mbig;

           copy[ j + 1, i + 1 ] = -tran[ j + 1, i + 1 ];
           if ( i = j ) then
              copy[ j + 1, i + 1 ] = copy[ j + 1, i + 1 ] + 1.0;

           end;
        end;

     /* solve the equations  */
     result = sweep( copy, 1 : mbig + 1 );
     answer = result[ , mb1 + 1 ];

     /* exit loop if data are integer */
     if ( realx = 0 )
     then do;
        regarl = answer[ 1 ];
        firarl = answer[ ( mbig + 1 ) / 2 + 1 ];
        print , "Algorithm terminated for integer case";
        print , regarl;
        print , firarl;
        stop;
        end;

     /* for a head start a*di other than di/2, set the subscript */
     /* above to the integer part of a*(mbig+1)                  */

     richar1[ mloop, 1 ] = answer[ 1 ];
     richar2[ mloop, 1 ] = ( answer[ mbig / 2 + 1 ] +
                             answer[ mbig / 2 + 2 ]   ) / 2.0 ;

     /* for a head start a*di other than di/2, use linear inter- */
     /* polation of the table entries bracketing a*mbig          */

     /* perform Richardson extrapolation */
     pow = 1.0;
     do jl = 2 to mloop;
        pow = pow # 4;
        richar1[ mloop, jl ] = ( pow # richar1[ mloop, jl-1 ] -
                               richar1[ mloop-1, jl-1 ] ) / ( pow - 1 );
        richar2[ mloop, jl ] = ( pow # richar2[ mloop, jl-1 ] -
                               richar2[ mloop-1, jl-1 ] ) / ( pow - 1 );
        end;

     if ( mloop ^= 1 )
     then do;

        oacce1 = richar1[ mloop - 1, mloop - 1 ];
        oacce2 = richar2[ mloop - 1, mloop - 1 ];
        acce1  = richar1[ mloop, mloop ];
        acce2  = richar2[ mloop, mloop ];
        esterr = max( abs( acce1 - oacce1 ),
                      abs( acce2 - oacce2 )  ) / abs( acce1 );

        /* test for convergence */
        if ( esterr < eps )
        then do;

           /* diagnostic */
           print , richar1, richar2;

           regarl = acce1;
           firarl = acce2;
           print , "Algorithm converged";
           print , regarl;
           print , firarl;
           stop;
           end;
        end;

     end;

  /* diagnostic */
  print , richar1, richart2;
  print , "ERROR: Convergence to requested precision not obtained";

  finish arlcomp;
  run arlcomp;

quit;
run;