Recursive Residuals, CUSUM, and CUSS Plots

 /****************************************************************/
 /*          S A S   S A M P L E   L I B R A R Y                 */
 /*                                                              */
 /*    NAME: CUSUM                                               */
 /*   TITLE: Recursive Residuals, CUSUM, and CUSS Plots          */
 /* PRODUCT: IML                                                 */
 /*    KEYS: MATRIX REGR                                         */
 /*   PROCS: IML PLOT                                            */
 /*    DATA:                                                     */
 /*                                                              */
 /* SUPPORT: LWB                         UPDATE:                 */
 /*     REF:                                                     */
 /*    MISC: CONVERTED FROM MATRIX TO IML USING MATIML           */
 /*                                                              */
 /****************************************************************/
*---------------------------CUSUM------------------------------*
| THIS IML ROUTINE COMPUTES RECURSIVE RESIDUALS AND OUTPUTS    |
| THEIR CUMULATIVE SUM AND SS (CALLED CUSUM AND CUSUMSS) TO A  |
| SAS DATA SET TO BE PLOTTED.  THE PLOT SHOULD SHOW IF THE     |
| REGRESSION PARAMETERS CHANGE OVER TIME.                      |
| FROM: BROWN,DURBIN, AND EVANS, "TECHNIQUES FOR TESTING THE   |
|       CONSTANCY OF REGRESSION RELATIONSHIPS OVER TIME",      |
|       JRSS B V37, P149 (1975).                               |
*------------------------------------------------JPS--25OCT79--*;

TITLE 'CUSUM for Data Simulated with Change At 25th Observation';
DATA A;
      DO N=1 TO 40;
         IF N>25 THEN B=5;
         ELSE B=10;
         X0=1;    X1=RANUNI(198721);    X2=RANUNI(198721)*5;
         YDAT= X0+B*X1+X2 + RANNOR(198721);
         OUTPUT;
      END;

PROC IML;
      USE A;
      READ ALL VAR {X0 X1 X2} INTO X;
      READ ALL VAR {YDAT} INTO Y;
      CLOSE A;
*--------CUSUM----------;
      N=NROW(X);
      P=NCOL(X);
      P1=P+1;

      *---FIRST P+1 OBSERVATIONS---;
      Z=X[1:P1,];
      YY=Y[1:P1,];
      XY=Z`*YY;
      IXPX=INV(Z`*Z);
      B=IXPX*XY;
      R=Y[P1,]-Z[P1,]*B;

      *---RECURSION UNTIL N---;
      START RECURSE;
         DO I=P+2 TO N;
            XI=X[I,]; YI=Y[I,];
            XY=XY+XI`*YI;
            IXPX=INVUPDT(IXPX,XI);
            B=IXPX*XY;
            R=R//(YI-XI*B);
         END;
      FINISH;
      RUN RECURSE;

      *---PRINT ESTIMATES---;
      PRINT B[L="FINAL PARAMETER ESTIMATES"];
      PRINT IXPX[L="INVERSE XPX MATRIX"];

      *---FORM LOWER TRIANGLE OF ONES---;
      M=N-P;
      TRI= (1:M)`*REPEAT(1,1,M)  >= REPEAT(1,M,1)*(1:M) ;

      *---FORM CSUM AND CSUMSS---;
      CUSUM=TRI*R;
      CUSUMSS=TRI*(R#R);

      *---OUTPUT TO DATASET---;
      RNAMES={ N   RESIDUAL   CUSUM   CUSUMSS  };
      X= (P1:N)` || R || CUSUM || CUSUMSS;
      PRINT  X [COLNAME=RNAMES];
      CREATE B FROM X [COLNAME=RNAMES];
      APPEND FROM X;
QUIT;


PROC SGSCATTER DATA=B;
LABEL RESIDUAL='Recursive Residual';
compare x=N y=(RESIDUAL CUSUM CUSUMSS);
run;
title;