PROGRAM ROB2SLS ! DEVELOPED BY SK MISHRA, NEHU, SHILLONG (INDIA)
      IMPLICIT DOUBLE PRECISION  (A-H,O-Z)
      CHARACTER *30 OFIL ! LENGTH OF INPUT DATA FILE NAME
      INTEGER *4 TIM1,TIM2,TIME
      PARAMETER (N=100, M=5, K=6, KMX=15, NITER=100, NOUT=60)
C     N=NO. OF OBSERVATIONS, M=NO. OF ENDOGENOUS VARIABLES
C     K=NO. OF EXOGENOUS VARIABLES, KMX=MAX DIMENSION FOR K+M
C     NITER=NO. OF ITERATIONS, NOUT=NO. OF OUTLIERS TO INTRODUCE
      PARAMETER (OMIN= 30.D0, OMAX=300.D0)! LIMITS ON OUTLIERS
      PARAMETER (OFIL='TSLDAT.TXT')! INPUT DATA FILE CONTAINING Y AND X
      COMMON /RNDM/IU,IV ! COMMON FOR RANDOM NUMBER GENERATOR
      DIMENSION Y(N,M),X(N,K),Z(N,KMX),YH(N,M),AT(M,M),BT(K,M),W(N)
      DIMENSION A(M,M),B(K,M),P(K,M),E(N,M),AA(M,M),BB(K,M),C(M+K)
      DIMENSION ZZ(M+K,M+K),XV(KMX,KMX),ARMS(M,M),BRMS(K,M),IOC(N)
      DIMENSION AMEAN(M,M),ASD(M,M),BMEAN(K,M),BSD(K,M)
C     ------------------ TRUE PARAMETERS --------------------------
      DATA ((AT(I,J),I=1,5),J=1,5)/-1,7,0,-6,0, 3,-1,5,0,0,  0,0,-1,3,0,
     & 6,0,0,-1,-3,  -11,0,9,0,-1/
      DATA ((BT(I,J),I=1,6),J=1,5)/0,5,0,-7,0,60,  3,0,-5,0,0,20,
     & 0,2,0,0,0,9,  0,4,0,0,-3,-8,  0,0,0,6,0, -11/
C     ---------MATRIX OF PRESENCE AND ABSENCE OF VARIABLES--------------
      DATA ((AA(I,J),I=1,5),J=1,5)/-1,1,0,1,0,  1,-1,1,0,0,  0,0,-1,1,0,
     & 1,0,0,-1,1,  1,0,1,0,-1/
      DATA ((BB(I,J),I=1,6),J=1,5)/0,1,0,1,0,1,  1,0,1,0,0,1,
     & 0,1,0,0,0,1,  0,1,0,0,1,1,  0,0,0,1,0, 1/
C     VALUE=1 FOR VARIABLES PRESENT, 0 FOR ABSENT, -1 FOR REGRESSAND Y
C     ------------------------------------------------------------------
      WRITE(*,*)'----- ROBUST WEIGHTED TWO-STAGE LEAST SQUARES --------'
      WRITE(*,*)'--- PROGRAM BY SK MISHRA, NEHU, SHILLONG (INDIA) -----'
      WRITE(*,*)'------------------------------------------------------'
      WRITE(*,*)'FEED THE CHOICE: 2SLS[0],W2SLS-OCP[1],W2SLS-MCP[2]'
      READ(*,*) NCH
      WRITE(*,*)'   '
      WRITE(*,*)'FEED RANDOM NUMBER SEED 4-DIGIT ODD INTEGER'
      READ(*,*) IU
C     READ Y AND THEN X  FROM INPUT FILE
      OPEN(7,FILE=OFIL)
      DO I=1,N
      READ(7,*)(Y(I,J),J=1,M)
      ENDDO
      DO I=1,N
      READ(7,*)(X(I,J),J=1,K)
      ENDDO
      CLOSE(7)
C     DATA READING OVER
      MK=M+K
C     INITIALIZE
      DO I=1,M
      DO J=1,M
      AMEAN(I,J)=0.D0
      ASD(I,J)=0.D0
      ARMS(I,J)=0.D0
      ENDDO
      ENDDO
      DO I=1,K
      DO J=1,M
      BMEAN(I,J)=0.D0
      BSD(I,J)=0.D0
      BRMS(I,J)=0.D0
      ENDDO
      ENDDO
      WRITE(*,*)'COMPUTING. PLEASE WAIT'
      TIM1=TIME()
      SIOC=0 ! TOTAL NO. OF POINTS DISTURBED
C     -----------------------------------------------------------------
      DO NIT=1,NITER ! BEGINNING OF THE OUTERMOST ITERATION LOOP
C     -----------------------------------------------------------------
      DO I=1,N
      IOC(I)=0 ! COUNTS THE NO. OF DATA POINTS UNDER PERTURATION
      DO J=1,M
      Z(I,J)=Y(I,J)
      ENDDO
      DO J=1,K
      Z(I,M+J)=X(I,J)
      ENDDO
      ENDDO
C     ADD NOUT NUMBER OF OUTLIERS AT RANDOM LOCATIONS
      IF(NOUT.EQ.0) THEN
      NOUTLIER=1
      MULT=0
      ELSE
      NOUTLIER=NOUT
      MULT=1
      ENDIF
      IOC(I)=0
      DO I=1,NOUTLIER
      CALL OUTLIER(N,M,OMIN,OMAX,IX,JX,OL)
      Z(IX,JX)=Z(IX,JX)+OL*MULT
      IOC(IX)=1
      ENDDO
C     ASSIGNMENT OF WEIGHTS -------------------------------------------
      DO I=1,N
      W(I)=1.D0
      ENDDO
      IF(NCH.EQ.1) CALL RCAMPBELL(Z,N,MK,W)
      IF(NCH.EQ.2) CALL MCAMPBELL(Z,N,MK,W)
C     -----------------------------------------------------------------
C     COMPUTE Z'Z
      DO J=1,MK
      DO JJ=1,MK
      ZZ(J,JJ)=0.D0
      DO I=1,N
      ZZ(J,JJ)=ZZ(J,JJ)+Z(I,J)*Z(I,JJ)*W(I)**2
      ENDDO
      ENDDO
      ENDDO
C     STORE X'X INTO XV
      DO J=1,K
      DO JJ=1,K
      XV(J,JJ)=ZZ(M+J,M+JJ)
      ENDDO
      ENDDO
C      WRITE(*,*)'XV MATRIX'
C      DO J=1,K
C      WRITE(*,1)(XV(J,JJ),JJ=1,K)
C      ENDDO
C     INVERT XX
      CALL MINV(XV,K,D)
C     PRE-MULTIPLY INVERTED XV BY X'Y
C      WRITE(*,*)'INVERTED MATRIX DET =',D
C      DO J=1,K
C      WRITE(*,2)(XV(J,JJ),JJ=1,K)
C      ENDDO
      DO J=1,K
      DO JJ=1,M
      P(J,JJ)=0.D0
      DO I=1,K
      P(J,JJ)=P(J,JJ)+XV(J,I)*ZZ(I+M,JJ)
      ENDDO
      ENDDO
      ENDDO
C     PRING COMPUTED P MATRIX
C      WRITE(*,*)'ESTIMATED P MATRIX'
C      DO I=1,K
C      WRITE(*,1)(P(I,J),J=1,M)
C      ENDDO
    1 FORMAT(6F12.4)
    2 FORMAT(6E12.5)
C     COMPUTE ESTIMATED Y; YH=XP
      DO I=1,N
      DO J=1,M
      YH(I,J)=0.D0
      DO JJ=1,K
      YH(I,J)=YH(I,J)+X(I,JJ)*P(JJ,J)
      ENDDO
      ENDDO
      ENDDO
C     SECOND STAGE OF 2-STAGE LEAST SQUARES
      DO J=1,M ! FOR M EQUATIONS
      DO I=1,N
      Z(I,1)=Y(I,J)
      ENDDO
      J1=0
      DO JJ=1,M
      IF(AA(JJ,J).GT.0) THEN
      J1=J1+1
      DO I=1,N
      Z(I,J1+1)=YH(I,JJ)
      ENDDO
      ENDIF
      ENDDO
      DO JJ=1,K
      IF(BB(JJ,J).GT.0) THEN
      J1=J1+1
      DO I=1,N
      Z(I,J1+1)=(-X(I,JJ))
C      ! ------ (-X DUE TO CHANGE IN SIGN OF Y)
      ENDDO
      ENDIF
      ENDDO
      J2=J1+1
      IF(NCH.EQ.1) CALL RCAMPBELL(Z,N,J2,W)
      IF(NCH.EQ.2) CALL MCAMPBELL(Z,N,J2,W)
      DO JA=1,J2
      DO JB=1,J2
      ZZ(JA,JB)=0.D0
      DO I=1,N
      ZZ(JA,JB)=ZZ(JA,JB)+Z(I,JA)*Z(I,JB)*W(I)**2
      ENDDO
      ENDDO
      ENDDO
      DO JA=1,J1
      DO JB=1,J1
      XV(JA,JB)=ZZ(JA+1,JB+1)
      ENDDO
      ENDDO
      CALL MINV(XV,J1,D)
      DO JA=1,J1
      C(JA)=0.D0
      DO JB=1,J1
      C(JA)=C(JA)+XV(JA,JB)*ZZ(JB+1,1)
      ENDDO
      ENDDO
C      WRITE(*,*)'EQUATION NUMBER =', J
C      WRITE(*,1)(C(JA),JA=1,J1)
C     PLACE COEFFICIENTS IN PROPER MATRIX CELL
      J1=0
      DO JA=1,M
      IF(AA(JA,J).GT.0) THEN
      J1=J1+1
      A(JA,J)=C(J1)
      ELSE
      IF(JA.EQ.J) A(JA,J)=-1
      ENDIF
      ENDDO
      DO JA=1,K
      IF(BB(JA,J).GT.0) THEN
      J1=J1+1
      B(JA,J)=C(J1)
      ENDIF
      ENDDO
C      WRITE(*,*)'EQUATION NUMBER =', J
C      WRITE(*,1)(A(JA,J),JA=1,M)
C      WRITE(*,1)(B(JA,J),JA=1,K)
      ENDDO
C     ------------------------------------------------------------------
      DO I=1,M
      DO J=1,M
      AMEAN(I,J)=AMEAN(I,J)+A(I,J)
      ASD(I,J)=ASD(I,J)+A(I,J)**2
      ARMS(I,J)=ARMS(I,J)+(A(I,J)-AT(I,J))**2
      ENDDO
      ENDDO
      DO I=1,K
      DO J=1,M
      BMEAN(I,J)=BMEAN(I,J)+B(I,J)
      BSD(I,J)=BSD(I,J)+B(I,J)**2
      BRMS(I,J)=BRMS(I,J)+(B(I,J)-BT(I,J))**2
      ENDDO
      ENDDO
      DO I=1,N
      SIOC=SIOC+IOC(I)
      ENDDO
      ENDDO ! END OF THE OUTERMOST ITERATION LOOP
      WRITE(*,*)' '
      WRITE(*,*)'----------- RESULTS OF MONTE CARLO EXPERIMENT --------'

      WRITE(*,*)' '
      WRITE(*,*)'MEAN OF TRANSPOSE(A) COEFFICIENTS'
      DO J=1,M
      WRITE(*,1)(AMEAN(I,J)/NITER,I=1,M)
      ENDDO
      WRITE(*,*)'SD OF TRANSPOSE(A) COEFFICIENTS'
      DO J=1,M
      WRITE(*,1)(DSQRT(ASD(I,J)/NITER-(AMEAN(I,J)/NITER)**2),I=1,M)
      ENDDO
      WRITE(*,*)'RMS OF TRANSPOSE(A) COEFFICIENTS'
      DO J=1,M
      WRITE(*,1)(DSQRT(ARMS(I,J)/NITER),I=1,M)
      ENDDO
      WRITE(*,*)' '
      WRITE(*,*)'MEAN OF TRANSPOSE(B) COEFFICIENTS'
      DO J=1,M
      WRITE(*,1)(BMEAN(I,J)/NITER,I=1,K)
      ENDDO
      WRITE(*,*)'SD OF TRANSPOSE(B) COEFFICIENTS'
      DO J=1,M
      WRITE(*,1)(DSQRT(BSD(I,J)/NITER-(BMEAN(I,J)/NITER)**2),I=1,K)
      ENDDO
      WRITE(*,*)'RMS OF TRANSPOSE(B) COEFFICIENTS'
      DO J=1,M
      WRITE(*,1)(DSQRT(BRMS(I,J)/NITER),I=1,K)
      ENDDO
      WRITE(*,*)'-----------------------------------------------------'
      WRITE(*,*)' '
      IF(NOUT.NE.0) THEN
      WRITE(*,*)'POINTS (%) PURTURBED IN THE EXPERIMENT=',SIOC/NITER
      ENDIF
      TIM2=TIME()
      WRITE(*,*)'END OF THE EXPERIMENT. TIME TAKEN(SECONDS)=',TIM2-TIM1
      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=100)
      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
C     ------------------------------------------------------------------
C     CAMPBELL COVARIANCE MATRIX  (TYPE-I)
      SUBROUTINE RCAMPBELL(Z,N,MM,W)
      PARAMETER(NV=100,MV=15,ITRN=200)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(NV,MV),V(MV,MV),AV(MV),W(NV),XD(MV)
      DIMENSION D(NV),VV(MV,MV),DN(NV),Z(NV,MV)
      DATA B1,B2/2, 1.5/
C     SOME DEFINITIONS
      D0=DSQRT(DFLOAT(M))+B1/DSQRT(2.D0)
      B22=B2**2
      M=MM-1
      DO I=1,N
      DO J=1,M
      X(I,J)=Z(I,J)
      ENDDO
      ENDDO
      NSKIP=1
      IF(NSKIP.NE.1) THEN  ! DO NOT STANDARDIZE THE VARIABLES
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
      ENDIF
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
       D(I)=DD
       DN(I)=DD
       ENDDO
       DO I=1,N
       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.1.0D-05) W(I)=WD/D(I)
       ENDDO
      ENDDO
      DO J=1,M
      DO JJ=1,M
      V(J,JJ)=VV(J,JJ)
      ENDDO
      ENDDO
C      DO J=1,M
C      WRITE(*,1)(V(J,JJ),JJ=1,M)
C      ENDDO
    1 FORMAT(10F7.3)
C      WRITE(*,*)'------------------------------------------------------'
C     WEIGTING Z MATRIX
      DO J=1,MM
      DO I=1,N
      Z(I,J)=Z(I,J)*W(I)
      ENDDO
      ENDDO
      RETURN
      END
C     ------------------------------------------------------------------
C     CAMPBELL COVARIANCE MATRIX  (TYPE-II)
      SUBROUTINE MCAMPBELL(Z,N,MM,W)
      PARAMETER(NV=100,MV=15,ITRN=200)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(NV,MV),V(MV,MV),AV(MV),W(NV),XD(MV)
      DIMENSION D(NV),VV(MV,MV),DN(NV),Z(NV,MV)
      M=MM-1
      DO I=1,N
      DO J=1,M
      X(I,J)=Z(I,J)
      ENDDO
      ENDDO
      NSKIP=1
      IF(NSKIP.NE.1) THEN  ! DO NOT STANDARDIZE THE VARIABLES
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
      ENDIF
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
       D(I)=DD
       DN(I)=DD
       ENDDO
       CALL MEDIAN(DN,N,DNA,DNV)
       DO I=1,N
       DN(I)=DABS(DN(I)-DNA)
       ENDDO
       CALL MEDIAN(DN,N,DNAA,DNVV)
      DNAA=DNAA/0.6745D0
      DO I=1,N
      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)=0.25D0
      IF(DX.LE.3*DNAA.AND.DX.GT.2*DNAA) W(I)=0.1111111111D0
      IF(DX.LE.4*DNAA.AND.DX.GT.3*DNAA) W(I)=0.0625D0
      ENDDO
      ENDDO
      DO J=1,M
      DO JJ=1,M
      V(J,JJ)=VV(J,JJ)
      ENDDO
      ENDDO
C      DO J=1,M
C      WRITE(*,1)(V(J,JJ),JJ=1,M)
C      ENDDO
   1  FORMAT(10F7.3)
C      WRITE(*,*)'------------------------------------------------------'
C     WEIGHTING Z MATRIX
      DO J=1,MM
      DO I=1,N
      Z(I,J)=Z(I,J)*W(I)
      ENDDO
      ENDDO
      RETURN
      END
C     -----------------------------------------------------------------
C     SUBROUTINE FOR MATRIX INVERSION
      SUBROUTINE MINV(A,N,D)
      PARAMETER(KMX=15)! MUST BE CONSISTEN WITH KMX IN THE MAIN PROGRAM
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(KMX,KMX)
      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     -----------------------------------------------------------------
C     RANDOM NUMBER GENERATOR (UNIFORM BETWEEN 0 AND 1 - BOTH EXCLUSIVE)
       SUBROUTINE RANDOM(RAND)
       IMPLICIT DOUBLE PRECISION  (A-H,O-Z)
       COMMON /RNDM/IU,IV
       IV=IU*65539
       IF(IV.LT.0) THEN
       IV=IV+2147483647+1
       ENDIF
       RAND=IV
       IU=IV
       RAND=RAND*0.4656613E-09
       RETURN
       END
C     -----------------------------------------------------------------
      SUBROUTINE OUTLIER(N,M,OMIN,OMAX,IX,JX,OL)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /RNDM/IU,IV
    1 CALL RANDOM(RAND)
      KK=INT(RAND*N*M)+1
      JX=MOD(KK,M)
      IF(JX.EQ.0) GOTO 1
      IX=KK/M+1
      CALL RANDOM(RAND)
      OL=OMIN+(OMAX-OMIN)*(RAND-0.5D0)
C     OL IS THE SIZE (QUANTUM) OF PERTURBATION AND IX AND JX POINT
C     TO THE LOCATION WHERE OL WOULD BE ADDED TO DATA, Z.
      RETURN
      END