*
* REDI(OPTION) DV NBEG NEND
* # <suppl. card> (list of explanatory variables--
* please do not list the CONSTANT as it is included
* automatically)
*
* DV = the dependent variable
* NBEG = the first observation
* NEND = the last observation
*
* Option:
*    ITERATIONS = number of subsamples to draw [3000]
*
* REDI computes a robust linear regression based on
* Rousseeuw's REsistant DIagnostic described in
*   P. J. Rousseeuw and A. M. Leroy, Robust Regression
*   and Outlier Detection, Wiley, 1987, pp. 238-245.
* Like other "high-breakdown" estimators (including least
* median of squares and least trimmed squares), REDI
* detects and downweights outliers and leverage points
* which may distort the regression estimate. REDI is
* computed by a resampling scheme. The user chooses
* the number of subsamples by setting the option
* ITERATIONS. Its default value, 3000 subsamples, should
* be adequate for most applications. The user should
* consult the Rousseeuw-Leroy book for additional details.
*
PROCEDURE REDI DV NBEG NEND
*
TYPE INTEGER NBEG NEND
TYPE SERIES DV
OPTION INTEGER ITERATIONS   3000
LOCAL VECTOR[INTEGER] IV M
LOCAL VECTOR[REAL] CVEC
LOCAL INTEGER I J K L N NH ITR
LOCAL REAL  CR CRITMIN CRIT RMED
LOCAL SERIES R RD ROLS RROB
*
ENTER(VARYING) IV
INQUIRE(MATRIX=IV) K
COMPUTE K = K+1
COMPUTE N = NEND-NBEG+1
COMPUTE NH = N/2 + K/2
COMPUTE CRITMIN = 1000000.0
COMPUTE ITR = ITERATIONS
SET RD 1 N = 0.0
DIM IV(K-1)
*
DO L = 1,ITR
DISPLAY(UNIT=SCREEN) 'ITERATION ' L

* AT EACH ITERATION L, DRAW A RANDOM SUBSAMPLE OF K OBSERVATIONS.

:10 BOOT M 1 K 1 N
SET R 1 N = 0.0
DO I = 1,K
COMPUTE R(M(I)) = 1.0
DO J = 1,K
IF I .EQ. J
BEGIN
BRANCH 20
END
IF M(I) .EQ. M(J)
BEGIN
BRANCH 10
END
:20 END DO J
END DO I

* SKIP ANY SUBSAMPLE THAT IS SINGULAR OR NEARLY SO.

CMOM(SMPL=R,CORR)
# IV
EIGEN %CMOM CVEC
COMPUTE CR = CVEC(1)/CVEC(K-1)
COMPUTE CR = SQRT(CR)
IF CR .GT. 25.0
BEGIN
BRANCH 10
END

* USE THE SUBSAMPLE TO COMPUTE A SET OF REGRESSION COEFFICIENTS.
* THEN COMPUTE THE RESIDUALS OVER THE ENTIRE DATA SET.

LINREG(NOPRINT,SMPL=R) DV
# CONSTANT IV
PRJ ROLS 1 N
SET R 1 N = DV-ROLS

* COMPUTE THE RESISTANT DIAGNOSTIC FOR EACH OBSERVATION

SET ROLS 1 N = ABS(R(T))
ORDER ROLS 1 N
COMPUTE RMED = ROLS(NH)
SET RD 1 N = %MAX(RD(T),ABS(R(T))/RMED)

END DO L

* HAVING TRIED ITR SUBSAMPLES, COMPUTE
* A SET OF RESISTANT DIAGNOSTICS

STATISTICS(NOPRINT,FRACTILES) RD 1 N
SET RD 1 N = RD / %MEDIAN

* USE THE RD'S AS WEIGHTS IN A REGRESSION

SET ROLS 1 N = RD**2
LINREG(NOPRINT,SPREAD=ROLS) DV 1 N R
# CONSTANT IV

* STANDARDIZE THE RESIDUALS

SET R 1 N = R/SQRT(%SEESQ)

* NEXT, RUN AN OLS REGRESSION, OMITTING THE LARGE
* ROBUST RESIDUALS.

DISPLAY 'ROBUST REDI REGRESSION OMITTING OUTLIERS AND'
DISPLAY 'LEVERAGE POINTS'

SET RROB 1 N = ABS(R(T)) .LE. 2.5
LINREG(SMPL=RROB) DV
# CONSTANT IV
PRJ ROLS 1 N
SET R 1 N = DV-ROLS
SET R 1 N = R/SQRT(%SEESQ)


* FINALLY, COMPUTE AN OLS REGRESSION WITH NO
* ADJUSTMENT FOR ROBUSTNESS

DISPLAY 'ORDINARY LEAST SQUARES REGRESSION'

LINREG DV 1 N ROLS
# CONSTANT IV
SET ROLS 1 N = ROLS/SQRT(%SEESQ)


* REPORT ANY LARGE ROBUST RESIDUALS.

DISPLAY 'LARGE REDI RESIDUALS AND CORRESPONDING'
DISPLAY 'OLS RESIDUALS'
DISPLAY 'OBS.      REDI         OLS'

DO I = 1,N
IF ABS(R(I)) .GT. 2.5
BEGIN
DISPLAY I R(I) ROLS(I)
END
END DO I

SET RROB 1 N = R
SET R 1 N = T
SCATTER(STYLE=SYMBOL,PATTERNS,KEY=BELOW) 2
# R RROB 1 N
# R ROLS 1 N
END REDI
