*
* LQD(OPTION) DV  NBEG  NEND
* # <SUPPL. CARD> LIST OF EXPLANATORY VARIABLES--
* PLEASE DO NOT INCLUDE THE CONSTANT AS IT IS
* INCLUDED AUTOMATICALLY !!
*
* DV = THE DEPENDENT VARIABLE
* NBEG = THE FIRST OBSERVATION
* NEND = THE LAST OBSERVATION
*
* OPTION:
*   ITERATIONS = THE NUMBER OF SUBSAMPLES TO USE [2000]
*
* LQD.SRC COMPUTES A ROBUST LINEAR REGRESSION CALLED
* THE LEAST QUARTILE DIFFERENCE ESTIMATOR (LQD). IT WAS
* PROPOSED IN THE FOLLOWING PAPER:
*
*     CHRISTOPHE CROUX, PETER J. ROUSSEEUW, AND OLA
*     HOSSJER (1994), "GENERALIZED S-ESTIMATORS,"
*     JOURNAL OF THE AMERICAN STATISTICAL ASSOCIATION,
*     89, 1271-1281.
*
*  THE LQD IS CLOSELY RELATED TO OTHER HIGH-BREAKDOWN,
*  COMPUTER-INTENSIVE REGRESSION METHODS LIKE LEAST
*  MEDIAN OF SQUARES AND LEAST TRIMMED SQUARES. ALL THESE
*  METHODS ARE VERY EFFECTIVE AT DETECTING OUTLIERS AND
*  LEVERAGE POINTS IN DATA. WHEN THE RANDOM ERRORS ARE
*  IN FACT UNCONTAMINATED GAUSSIAN VARIABLES, LQD IS MORE
*  EFFICIENT THAN THESE OTHER HIGH-BREAKDOWN METHODS.
*  (ITS ASYMPTOTIC EFFICIENCY IS ABOUT 67 PERCENT
*  COMPARED TO LEAST SQUARES AT A GAUSSIAN DISTRIBUTION;
*  THE ASYMPTOTIC EFFICIENCY OF LMS IS ZERO AND THAT OF
*  LTS IS JUST 8 PERCENT.) MOREOVER, THE LQD REGRESSION 
*  COEFFICIENTS ARE ASYMPTOTICALLY GAUSSIAN WITH A COVARIANCE
*  MATRIX PROPORTIONAL TO THE INVERSE OF X'X. 
*
*  LIKE OTHER HIGH-BREAKDOWN REGRESSION METHODS, LQD CAN BE
*  COMPUTED BY A RESAMPLING SCHEME. THE USER CHOOSES THE
*  NUMBER OF SUBSAMPLES BY SETTING THE OPTION ITERATIONS. 
*  THE FOLLOWING GUIDELINES ARE SUGGESTED. 
* 
*           INDEPENDENT                  MINIMUM NUMBER
*            VARIABLES                   OF SUBSAMPLES
*                1                             500
*                2                           1,000
*                3                           1,500
*                4                           2,000
*                5                           2,500
*                6 OR MORE                   3,000
*
*  FOR MODERATE OR LARGE SAMPLES, LQD.SRC MAY RUN SLOWLY; USERS
*  MAY WANT TO USE LQD.SRC IN BATCH MODE OR DURING OFF-PEAK 
*  COMPUTER TIME.
*  OUTPUT OF LQD.SRC INCLUDES THE OLS REGRESSION, THE LQD
*  STANDARD ERROR OF ESTIMATE, THE LQD REGRESSION COEFFICIENTS
*  TOGETHER WITH THEIR ASYMPTOTIC STANDARD ERRORS AND T-RATIOS,
*  A LIST OF OBSERVATIONS WHICH ARE PROBABLE OUTLIERS, AND A
*  PLOT OF THE STANDARDIZED LQD AND OLS RESIDUALS.
PROCEDURE LQD  DV  NBEG  NEND
*
TYPE INTEGER NBEG NEND
TYPE SERIES DV
OPTION INTEGER ITERATIONS  2000
LOCAL VECTOR[REAL] CVEC WORK ACAND BROB
LOCAL REAL CR CRITMIN CRIT HARPO DN QN SD
LOCAL INTEGER J JJ N H KK KNEW JHELP NL NR SUMQ SUMP FOUND L ITR
LOCAL INTEGER WTOTAL WREST WLEFT WMID WRIGHT NN KCAND MGM NDF 
LOCAL VECTOR[INTEGER] IV LEFT RIGHT WEIGHT Q P IWCAND  
LOCAL SERIES[REAL] R ROLS RLQD
ENTER(VARYING) IV
INQUIRE(MATRIX=IV) K
COMPUTE K = K+1
COMPUTE N = NEND-NBEG+1
COMPUTE CRITMIN = 1000000.0
DIM WORK(N) ACAND(N) BROB(K)
DIM IV(K-1) LEFT(N) RIGHT(N) WEIGHT(N)
DIM  Q(N) P(N) IWCAND(N)
COMPUTE H=(N+K+1)/2
COMPUTE KK=H*(H-1)/2
COMPUTE ITR = ITERATIONS
*
* START THE ITERATIONS
*
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.

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

* COMPUTE THE ROBUST SCALE QN

DO I = 1,N
COMPUTE LEFT(I)=N-I+2
COMPUTE RIGHT(I)=N
END DO I
COMPUTE JHELP=N*(N+1)/2
COMPUTE KNEW=KK+JHELP
COMPUTE NL=JHELP
COMPUTE NR=N*N
COMPUTE FOUND = 0
DO MGM = 1,1
:200 IF((NR-NL.LE.N) .OR. (FOUND.EQ.1))
BRANCH 1010
COMPUTE J=1
DO I =2,N
IF (LEFT(I).LE.RIGHT(I))
   {
   COMPUTE WEIGHT(J)=RIGHT(I)-LEFT(I)+1
   COMPUTE JHELP=LEFT(I)+WEIGHT(J)/2
   COMPUTE WORK(J)=ROLS(I)-ROLS(N+1-JHELP)
   COMPUTE J=J+1
   }
END DO I

* START WHIMED

COMPUTE NN=J-1
COMPUTE WTOTAL=0
DO I=1,NN
COMPUTE WTOTAL=WTOTAL+WEIGHT(I)
END DO I
COMPUTE WREST=0
:100 SET WORK1 1 NN = WORK(T)
ORDER WORK1 1 NN
COMPUTE TRY = WORK1(NN/2+1)
COMPUTE WLEFT=0
COMPUTE WMID=0
COMPUTE WRIGHT=0
DO I = 1,NN
IF (WORK(I) .LT. TRY)
COMPUTE WLEFT=WLEFT+WEIGHT(I)
IF (WORK(I) .GT. TRY)
COMPUTE WRIGHT=WRIGHT+WEIGHT(I)
IF (WORK(I) .EQ. TRY)
COMPUTE WMID=WMID+WEIGHT(I)
END DO I
IF ((2*WREST+2*WLEFT).LE.WTOTAL)
BRANCH 5050
COMPUTE KCAND=0
DO I = 1,NN
IF (WORK(I).LT.TRY)
   {
   COMPUTE KCAND=KCAND+1
   COMPUTE ACAND(KCAND)=WORK(I)
   COMPUTE IWCAND(KCAND)=WEIGHT(I)
   }
END DO I
COMPUTE NN=KCAND
BRANCH 6060
:5050 IF ((2*WREST+2*WLEFT+2*WMID).GT.WTOTAL)
   {
   COMPUTE WHIMED=TRY
   BRANCH 900
   }
COMPUTE KCAND=0
DO I=1,NN
IF (WORK(I).GT.TRY)
   {
   COMPUTE KCAND=KCAND+1
   COMPUTE ACAND(KCAND)=WORK(I)
   COMPUTE IWCAND(KCAND)=WEIGHT(I)
   }
END DO I
COMPUTE NN=KCAND
COMPUTE WREST=WREST+WLEFT+WMID
:6060 DO I = 1,NN
COMPUTE WORK(I)=ACAND(I)
COMPUTE WEIGHT(I) = IWCAND(I)
END DO I
BRANCH 100

* END WHIMED

:900 COMPUTE TRIAL= WHIMED
COMPUTE J=0
DO I=N,1,-1
:45 IF ((J.LT.N) .AND. ((ROLS(I)-ROLS(N-J)) .LT. TRIAL))
   {
   COMPUTE J=J+1
   BRANCH 45
   }
COMPUTE P(I)=J
END DO I
COMPUTE J=N+1
DO I =1,N
:55 IF((ROLS(I)-ROLS(N-J+2)).GT.TRIAL)
  {
  COMPUTE J=J-1
  BRANCH 55
  }
COMPUTE Q(I)=J
END DO I
COMPUTE SUMP=0
COMPUTE SUMQ=0
DO I=1,N
COMPUTE SUMP=SUMP+P(I)
COMPUTE SUMQ=SUMQ+Q(I)-1
END DO I
IF (KNEW.GT.SUMP)
BRANCH 2020
DO I=1,N
COMPUTE RIGHT(I)=P(I)
END DO I
COMPUTE NR=SUMP
BRANCH 3030
:2020 IF (KNEW.LE.SUMQ)
BRANCH 4040
DO I=1,N
COMPUTE LEFT(I)=Q(I)
END DO I
COMPUTE NL=SUMQ
BRANCH 3030
:4040 COMPUTE QN=TRIAL
COMPUTE FOUND=1
:3030 BRANCH 200
:1010 IF FOUND .EQ. 0
   {
   COMPUTE J=1
   DO I =2,N
   IF (LEFT(I) .LE. RIGHT(I))
           {
           DO JJ=LEFT(I),RIGHT(I)
           COMPUTE WORK(J)=ROLS(I)-ROLS(N-JJ+1)
           COMPUTE J=J+1
           END DO JJ
           }
     END DO I
     SET WORK2 1 J-1 = WORK(T)
     ORDER WORK2 1 J-1
     COMPUTE QN = WORK2(KNEW-NL)
     }
END MGM

* IS QN IN THIS SUBSAMPLE THE SMALLEST SO FAR ? IF SO, SAVE
* THE RESIDUALS AND THE REGRESSION COEFFICIENTS. THEN DRAW
* ANOTHER SUBSAMPLE.

IF QN .LT. CRITMIN
BEGIN
COMPUTE CRITMIN = QN
SET RLQD 1 N = R
COMPUTE BROB = %BETA
END
END DO L

* HAVING TRIED ITR SUBSAMPLES, REPORT THE ROBUST SCALE ESTIMATE
* QN AND THE ROBUST REGRESSION COEFFICIENTS

DO MGM = 1,1
COMPUTE NDF = N-K
IF (NDF .GT. 9)
BRANCH 7070
IF (NDF.EQ.2)
COMPUTE DN=0.399
IF (NDF.EQ.3)
COMPUTE DN=0.994
IF (NDF.EQ.4)
COMPUTE DN=0.512
IF (NDF.EQ.5)
COMPUTE DN=0.844
IF (NDF.EQ.6)
COMPUTE DN=0.611
IF (NDF.EQ.7)
COMPUTE DN=0.857
IF (NDF.EQ.8)
COMPUTE DN=0.669
IF (NDF.EQ.9)
COMPUTE DN=0.872
BRANCH 8080
:7070 COMPUTE HARPO = %FRAC(NDF/2.0)
IF HARPO .NE. 0.000
COMPUTE DN=NDF/(NDF+1.4)
IF HARPO .EQ. 0.000
COMPUTE DN=NDF/(NDF+3.8)
:8080 COMPUTE QN = DN*2.2219*CRITMIN
END DO MGM

* COMPUTE THE LEAST-SQUARES REGRESSION FOR THE FULL SAMPLE

DISPLAY '     '
DISPLAY 'ORDINARY LEAST SQUARES REGRESSION'
DISPLAY '     '
LINREG DV 1 N ROLS
# IV CONSTANT
SET ROLS 1 N = ROLS/SQRT(%SEESQ)
DISPLAY  '     '
DISPLAY 'LQD STANDARD ERROR OF ESTIMATE =  '  QN
DISPLAY '    '
DISPLAY 'LQD REGRESSION ESTIMATES (THE CONSTANT IS LISTED LAST)'
DISPLAY  '     '
DISPLAY 'VAR     LQD COEF  ASYMP. STD DEV    T-RATIO'
DO J = 1,K
COMPUTE SD = QN*SQRT(%XX(J,J))
DISPLAY J BROB(J) SD BROB(J)/SD
END DO J
DISPLAY '   '
DISPLAY 'LARGE ROBUST STANDARDIZED RESIDUALS'
DISPLAY 'OBS         LQD           OLS'   
STATISTICS(NOPRINT,FRACTILES) RLQD 1 N
SET RLQD 1 N = (RLQD-%MEDIAN)/QN
DO T = 1,N
IF ABS(RLQD(T)) .GT. 2.5
DISPLAY T RLQD(T) ROLS(T)
END DO T

* GRAPH STANDARDIZED RESIDUALS

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