PROGRAM ROBOLS ! ROBUST REGRESSION BY CAMPBELL METHOD
C     PROGRAM BY SK MISHRA, NEHU, SHILLONG (INDIA)
      PARAMETER(N=21,M=4,ITRN=50,NTYPE=2)! CHANGE THESE AS REQUIRED
C     N=NO. OF OBSERVATIONS, M=NO. OF VARIABLES INCLUDING Y
C     ITRN=NO. OF ITERATIONS (AT LEAST 50)
C     NTYPE =1 FOR CAMPBELL-I AND NTYPE = 2 FOR CAMPBELL-II REGRESSION
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION Z(N,M),WT(M),VX(M,M),VY(M),CF(M),AVZ(M)
      DIMENSION X(N,M),V(M,M),AV(M),W(N),XD(M),D(N),VV(M,M),DN(N)
      DATA B1,B2/2,1.25/ ! NOT TO BE CHANGED
C     READ DATA FROM FILE
C     -----------------------------------------------------------------
      OPEN(7,FILE='STACKLOSS.TXT')! INPUT FILE DATA [Y, X1, X2,..., XM]
C     -----------------------------------------------------------------
      D0=DSQRT(DFLOAT(M))+B1/DSQRT(2.D0)
      B22=B2**2
      DO I=1,N
      READ(7,*)(Z(I,J),J=1,M)
      ENDDO
      DO J=1,M
      AVZ(J)=0.D0
      DO I=1,N
      AVZ(J)=AVZ(J)+Z(I,J)
      ENDDO
      AVZ(J)=AVZ(J)/N
      DO I=1,N
      Z(I,J)=Z(I,J)-AVZ(J)
      ENDDO
      ENDDO
      DO I=1,N
      DO J=1,M
      X(I,J)=Z(I,J)
      ENDDO
      ENDDO
      CLOSE(7)
C     STANDARDIZE
      DO J=1,M
      AV(J)=0.D0
      XD(J)=0.D0
      DO I=1,N
      AV(J)=AV(J)+X(I,J)
      XD(J)=XD(J)+X(I,J)**2
      ENDDO
      AV(J)=AV(J)/N
      XD(J)=DSQRT(XD(J)/N-AV(J)**2)
      ENDDO
      DO J=1,M
      DO I=1,N
      X(I,J)=(X(I,J)-AV(J))/XD(J)
      ENDDO
      ENDDO
C     INITIALIZE WEIGHT VECTOR BY UNITY
      DO I=1,N
      W(I)=1.D0
      ENDDO
C     FIND SUM OF WEIGHTS
      DO ITER=1,ITRN

       SW=0.D0
       SSW=0.D0
       DO I=1,N
       SW=SW+W(I)
       SSW=SSW+W(I)**2
       ENDDO
       SSW=SSW-1.D0

C      COMPUTE MEAN VECTOR AND COVARIANCE MATRIX
       DO J=1,M
       AV(J)=0.D0
       DO I=1,N
       AV(J)=AV(J)+X(I,J)*W(I)
       ENDDO
       AV(J)=AV(J)/SW
       ENDDO
       DO J=1,M
       DO JJ=J,M
       V(J,JJ)=0.D0
       DO I=1,N
       V(J,JJ)=V(J,JJ)+(X(I,J)-AV(J))*(X(I,JJ)-AV(JJ))*W(I)**2
       ENDDO
       V(J,JJ)=V(J,JJ)/SSW
       IF(J.NE.JJ) V(JJ,J)=V(J,JJ)
       ENDDO
       ENDDO
       DO J=1,M
       DO JJ=1,M
       VV(J,JJ)=V(J,JJ)
       ENDDO
       ENDDO
C      INVERT V
       CALL MINV(V,M,DD) ! ON RETURN V IS INVERTED V
       DO I=1,N
       D(I)=0.D0
       DO J=1,M
       XD(J)=0.D0
       DO JJ=1,M
       XD(J)=XD(J)+(X(I,JJ)-AV(JJ))*V(JJ,J)
       ENDDO
       ENDDO
       DD=0.D0
       DO J=1,M
       DD=DD+XD(J)*(X(I,J)-AV(J))
       ENDDO
       DD=DSQRT(DD)
       D(I)=DD
       DN(I)=DD
       ENDDO
       IF(NTYPE.EQ.2) THEN
       CALL MEDIAN(DN,N,DNA,DNV)
       DO I=1,N
       DN(I)=DABS(DN(I)-DNA)
       ENDDO
       CALL MEDIAN(DN,N,DNAA,DNVV)
       ENDIF
       DNAA=DNAA/0.6745
       DO I=1,N
       IF(NTYPE.EQ.1) THEN
       IF(D(I).LE.D0)THEN
       WD= D(I)
       ELSE
       WD=D0*DEXP(-0.5D0*(D(I)-D0)**2/B22)
       ENDIF
       W(I)=1.D0
       IF(DABS(D(I)).GT.0.00001) W(I)=WD/D(I)
       ENDIF
       IF(NTYPE.EQ.2) THEN
       W(I)=0.D0
      DX=DABS(D(I)-DNA)
      IF(DX.LE.DNAA) W(I)=1.D0
      IF(DX.LE.2*DNAA.AND.DX.GT.DNAA) W(I)=.25D0
      IF(DX.LE.3*DNAA.AND.DX.GT.2*DNAA) W(I)=0.11D0
      IF(DX.LE.4*DNAA.AND.DX.GT.3*DNAA) W(I)=0.06D0
       ENDIF
       ENDDO
      ENDDO
      DO J=1,M
      DO JJ=1,M
      VX(J,JJ)=VV(J,JJ)
      V(J,JJ)=VV(J,JJ)/DSQRT(VV(J,J)*VV(JJ,JJ))
      ENDDO
      ENDDO
C     DO J=1,M
C     WRITE(*,1)(V(J,JJ),JJ=1,M)
C     ENDDO
    1 FORMAT(8F9.3)
      WRITE(*,*)'-----------------'
C     WRITE(*,1)(AV(J),J=1,M)
C     FIND REGRESSION COEFFICIENTS
      DO J=1,M
      DO JJ=1,M
      VX(J,JJ)=0.D0
      DO I=1,N
        VX(J,JJ)=VX(J,JJ)+Z(I,J)*Z(I,JJ)*W(I)**2
      ENDDO
      VX(J,JJ)=VX(J,JJ)/N
      ENDDO
      ENDDO

      DO I=1,M
      VY(I)=VX(I,1)
      VX(1,I)=0.D0
      VX(I,1)=0.D0
      ENDDO
      VX(1,1)=1.D0
      VY(1)=0.D0
      CALL MINV(VX,M,DD)
      DO J=1,M
      CF(J)=0.D0
      DO JJ=1,M
      CF(J)=CF(J)+VX(J,JJ)*VY(JJ)
      ENDDO
      ENDDO
      SW=0.D0
      DO I=1,N
      SW=SW+W(I)
      ENDDO
      DO J=1,M
      AV(J)=0.D0
      DO I=1,N
      AV(J)=AV(J)+(Z(I,J)+AVZ(J))*W(I)
      ENDDO
      AV(J)=AV(J)/SW
      ENDDO
      CF(1)=AV(1)
      DO J=2,M
      CF(1)=CF(1)-CF(J)*AV(J)
      ENDDO
      WRITE(*,*)'REGRESSION COEFFICIENTS'
      WRITE(*,*)(CF(J),J=1,M)
      WRITE(*,*)'WEIGHTS. UNITY IS CLEAR INLIER. SMALLER IS THE WEIGHT'
      WRITE(*,*)'MORE STRONG IS THE OUTLIER. EXTREME OUTLIER IS ZERO'
      DO I=1,N
      WRITE(*,3) I,W(I)
      ENDDO
      OPEN(8,FILE='ROBOLSRESULTS.TXT')
      WRITE(8,*)'REGRESSION COEFFICIENTS'
      WRITE(8,*)(CF(J),J=1,M)
      WRITE(8,*)'WEIGHTS. UNITY IS CLEAR INLIER. SMALLER IS THE WEIGHT'
      WRITE(8,*)'MORE STRONG IS THE OUTLIER. EXTREME OUTLIER IS ZERO'
      DO I=1,N
      WRITE(8,3) I,W(I)
      ENDDO
      CLOSE(8)
    3 FORMAT(I5,F10.4)
      WRITE(*,*)'END OF THE PROGRAM'
      WRITE(*,*)'RESULTS ARE STORED IN FILE ROBOLSRESULTS.TXT'
      END
C     SUBROUTINE FOR MATRIX INVERSION
      SUBROUTINE MINV(A,N,D)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(N,N)
      U=1.D0
      D=U
      DO I=1,N
      D=D*A(I,I)
      A(I,I)=U/A(I,I)
         DO J=1,N
         IF(I.NE.J) A(J,I)=A(J,I)*A(I,I)
         ENDDO
         DO J=1,N
         DO K=1,N
         IF(I.NE.J.AND.K.NE.I) A(J,K)=A(J,K)-A(J,I)*A(I,K)
         ENDDO
         ENDDO
         DO J=1,N
         IF(J.NE.I) A(I,J)= -A(I,J)*A(I,I)
         ENDDO
      ENDDO
      RETURN
      END
C     -----------------------------------------------------------------
      SUBROUTINE MEDIAN(X,N,A,V) ! ------------------------------------
C     SUBROUTINE MEDIAN : FINDS MEDIAN (A) AND MEAN DEVIATION (V) OF A
C     GIVEN VARIATE, VARIATE X(N)
      PARAMETER (NMAX=1000)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(N),Z(NMAX)
C     STORE X IN Z
      DO I=1,N
      Z(I)=X(I)
      ENDDO
C     ARRANGE Z IN AN ASCENDING ORDER
      DO I=1,N-1
      DO J=I+1,N
      IF(Z(I).GT.Z(J)) THEN ! EXCHANGE
      TEMP=Z(I)
      Z(I)=Z(J)
      Z(J)=TEMP
      ENDIF
      ENDDO
      ENDDO
      K=(N+1)/2 ! K IS OBTAINED AS INT((N+1)/2.0D0)
      A=(Z(K)+Z(N+1-K))/2.D0 ! GIVES MEDIAN FOR ODD AS WELL AS EVEN N
C     FIND MEAN DEVIATION
      V=0.D0
      DO I=1,N
      V=V+DABS(Z(I)-A) ! A IS MEDIAN
      ENDDO
      V=V/N ! V IS MEAN DEVIATION FROM MEDIAN
C     WRITE(*,*)'MEDIAN =',A,'  MEAN DEVIATION =',V
      RETURN
      END