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;