/**********************************************************************/


   /* Macro to Calculate ARL for Shewhart Chart */

options mprint nosymbolgen;

%macro shewhart(mu0=,sigma=,dl=,du=,d=,lcl=,ucl=);

data arls;
   mu0=&mu0; sigma=σ
   lcl=(&lcl-mu0)/sigma; ucl=(&ucl-mu0)/sigma;

   /* Compute ARL for various shifts in process mean       */
do shift = &dl to &du by &d;
   delta = shift/sigma;
   p=probnorm(ucl-delta) - probnorm(lcl-delta);
   arl = 1 / (1-p);
   output;
end;

proc print;
   var shift delta arl;
run;

%mend shewhart;
%shewhart(mu0=0,sigma=1,dl=0,du=3.5,d=.2,lcl=-3,ucl=3)


/**********************************************************************/


   /* Macro to Calculate ARL for Shewhart Chart with */
   /* Supplemental Run Rules                         */

options mprint ;

%macro shewsupp(side,mu0=,sigma=,dl=,du=,d=,c1=,r1=,r2=);

data arls;
   z11=0; z12=0; z13=0; z14=0;
   z21=0; z22=0; z23=0; z24=0;
   mu0=&mu0; sigma=σ

   /* Read Number of consecutive points */
   /* and region boundaries     */
%do i=1 %to 2;
   %do j=1 %to 4;
      %let z=%scan(&&r&i,&j,%str( ));
      %if %quote(&z)= %then %goto nexta;
      z&i&j=&z;
      z&i&j=(z&i&j - mu0)/sigma;
   %end;
%nexta: %end;

      /* Compute ARL for various shifts in process mean     */
   l=&c1;
   do shift = &dl to &du by &d;
      delta = shift/sigma;
      q=(probnorm(z12-delta) - probnorm(z11-delta))
         + (probnorm(z14-delta) - probnorm(z13-delta));
      p=(probnorm(z22-delta) - probnorm(z21-delta))
         + (probnorm(z24-delta) - probnorm(z23-delta));
      arlp = (1-q**l) / (1-q-(1-q**l)*p);
   
      /* Compute approximate ARL for Two-sided chart          */
      /* if necessary                                         */
   %if "&side" eq "ONESIDED" %then %do;
      q=(probnorm(z12+delta) - probnorm(z11+delta))
         + (probnorm(z14+delta) - probnorm(z13+delta));
      p=(probnorm(z22+delta) - probnorm(z21+delta))
         + (probnorm(z24+delta) - probnorm(z23+delta));
      arlm = (1-q**l) / (1-q-(1-q**l)*p);
      arl2 = (1/arlp) + (1/arlm);
      arlp = 1/arl2;
   %end;
arl=arlp;
output;
end;

proc print;
   var shift delta arl;
run;

%mend shewsupp;

%shewsupp(ONESIDED,mu0=0,sigma=1,dl=0,du=3.4,
           d=.2,c1=9,r1=0 3,r2=-1e8 0)

%shewsupp(TWOSIDED,mu0=0,sigma=1,dl=0,du=3.4,
           d=.2,c1=15,r1=-1 1,r2=-3 -1 1 3)

%shewsupp(TWOSIDED,mu0=0,sigma=1,dl=0,du=3.4,
           d=.2,c1=8,r1=-3 -1 1 3,r2=-1 1)

%shewsupp(ONESIDED,mu0=0,sigma=1,dl=0,du=3.4,
           d=.2,c1=3,r1=2 3,r2=-1e8 2)


/**********************************************************************/


   /* Example Code 1 */

   /* Example SAS code to Generate Cusum ARLs              */
data arls;
   mu0=0; sigma=1; h=8; k=.25;
   do shift=0 to 3.4 by .2;
      cusarl=cusumarl('TWOSIDED',shift,h,k);
      output;
   end;
run;

proc print;
run;


/**********************************************************************/


   /* Example Code 2 */

   /* SAS Code to Generate EWMA ARLs for Table 3           */
data ewma;
   lambda=.12;
   do shift=0,1;
      do h=2.35 to 2.95 by .10;
         ewmaarl = ewmaarl(shift,lambda,h); output;
      end;
   end;
run;

proc print;
   title 'EWMA Chart ARL Example';
run;


/**********************************************************************/


   /* Example Code 3 */

   /* Macro to compute approximate correlated data ARL     */
%macro approx(h=,k=,rho=,shift=,beta1=,beta2=);
data;
   h=&h; k=&k; shift=&shift; rho=ρ
   beta1=&beta1; beta2=&beta2;
   arl0 = cusumarl('TWOSIDED',shift,h,k);
   intercep = -probit(1/arl0);
   predpro = intercep + beta1*rho + beta2*rho**2;
   arlpred = 1 / probnorm(-predpro);
run;

proc print;
run;

%mend approx;
%approx(h=4,k=.5,rho=.3,shift=0,beta1=-1.92113,beta2=.60359)
%approx(h=4,k=.5,rho=.3,shift=1,beta1=-.52487,beta2=.087087)


/**********************************************************************/