## 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;

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;

/* 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;

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;
if ( arg < winsor ) then
phitab[ iadj ] = probnorm(arg - arg2);
end;
else do;
if ( iarg >= iwintg ) then
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;
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;

```