C                                           *** STAT PACK ***
C     SUBROUTINE FOR KENDALL TAU CORRELATION.
C     CALLING SEQUENCE: CALL STKTAU(NV,NC,MV,MC,DATA,AV,NAMES)
C     WHERE NV - IS THE NUMBER OF COLUMNS ACTUALLY FILLED(VARIABLES)
C           NC - IS THE NUMBER OF ROWS ACTUALLY FILLED(OBSERVATIONS,CASES)
C           MV - IS MAXIMUM NUMBER OF COLUMNS, AS SPECIFIED IN MAIN.
C           MC - IS MAXIMUM NUMBER OF ROWS, AS SPECIFIED IN MAIN.
C           DATA - STORAGE FOR DATA DIMENSIONED FOR MAXIMUM MATRIX.
C           AV - IS A VECTOR DIMENSIONED AT LEAST NV.
C           NAMES - IS A VECTOR CONTAINING VARIABLE NAMES
C
C     PROGRAM OUTPUTS ENTIRE KENDALL TAU CORRELATION MATRIX FOR ALL DATA
C
      SUBROUTINE STKTAU(NV,NC,MV,MC,DATA,AV,NAMES)
      DIMENSION DATA(MC,MV),AV(1),NAMES(1)
      COMMON /DEV/ICC,IDATA,IOUT,IDLG,IDSK
      COMMON /PRNT/ LINPP,ICOPS,RUNPRG
      COMMON/EXTRA/HEDR(70),NSZ
      IF(IOUT.NE.21) WRITE(IOUT,5566)(HEDR(I),I=1,NSZ)
5566  FORMAT('1',70A1)
      IF(IOUT.EQ.21) CALL PRNTHD
      ISQ=7
      IF(IOUT.EQ.21) ISQ=15
      SWTCH=0
      WRITE(IOUT,10)
10    FORMAT('0',10X,'***** KENDALL TAU CORRELATION MATRIX *****')
      LINES=4
14    DO 1 I=1,NV
      DO 2 J=1,I
      IF(I.NE.J) GO TO 3
      AV(J)=1.0
      GO TO 2
3     AM=0.
      TX=0.
      TY=0.
      DO 4 K=1,NC-1
      SX=DATA(K,I)
      SY=DATA(K,J)
      DO 4 L=K+1,NC
      SM=(SX-DATA(L,I))*(SY-DATA(L,J))
      IF(SM.EQ.0) GO TO 5
      IF(SM.LT.0) AM=AM-1.
      IF(SM.GT.0) AM=AM+1.
      GO TO 4
5     IF((SX-DATA(L,I)).EQ.0)TX=TX+1.
      IF((SY-DATA(L,J)).EQ.0)TY=TY+1.
4     CONTINUE
      TC=NC
      XNGO=.5*(TC*(TC-1.))
      DEN=(XNGO-TX)*(XNGO-TY)
      IF(DEN.GT.0) GOTO 6
      SWTCH=1.
      AV(J)=100**4
      GO TO 2
6     AV(J)=AM/SQRT(DEN)
2     CONTINUE
      IF(IOUT.NE.21) GO TO 17
      KK=(I+ISQ-1)/ISQ+1
      LINES=LINES+KK
      IF(LINES.LE.(LINPP-KK)) GO TO 17
      WRITE(IOUT,15)
      DO 16 J=1,I-1,ISQ
      NEND=J-ISQ-1
      IF(NEND.GT.(I-1)) NEND=I-1
16    WRITE(IOUT,8)(NAMES(K),K=J,NEND)
      CALL PRNTHD
      LINES=2+KK
17    DO 1 J=1,I,ISQ
      NEND=J+ISQ-1
      IF(NEND.GT.I) NEND=I
      IF(J.EQ.1) WRITE(IOUT,7) NAMES(I),(AV(K),K=J,NEND)
      IF(J.NE.1) WRITE(IOUT,12)(AV(K),K=J,NEND)
7     FORMAT('0',A5,2X,15F8.4)
12    FORMAT(8X,15F8.4)
1     CONTINUE
      WRITE(IOUT,15)
15    FORMAT(1X)
      DO 13 J=1,NV,ISQ
      NEND=J+ISQ-1
      IF(NEND.GT.NV) NEND=NV
13    WRITE(IOUT,8)(NAMES(K),K=J,NEND)
8     FORMAT(8X,15(3X,A5))
      IF(SWTCH.NE.0) WRITE(IOUT,9)
9     FORMAT('0NOTE: ALL STARS COULD NOT BE CALCULATED BECAUSE '
     1,'ONE OF THE S.D.S WAS EQUAL TO ZERO')
      RETURN
      END
C                                           *** STAT PACK ***
C     SUBROUTINE FOR STEPWISE REGRESSION
C     CALLING SEQUENCE: CALL STSTRG(NV,NC,MV,MC,DATA,COR,VMN,STD,CO,NAMES)
C     WHERE NV - IS THE NUMBER OF COLUMNS ACTUALLY FILLED(VARIABLES)
C           NC - IS THE NUMBER OF ROWS ACTUALLY FILLED(OBSERVATIONS,CASES)
C           MV - IS MAXIMUM NUMBER OF COLUMNS, AS SPECIFIED IN MAIN.
C           MC - IS MAXIMUM NUMBER OF ROWS, AS SPECIFIED IN MAIN.
C           DATA - STORAGE FOR DATA DIMENSIONED FOR MAXIMUM MATRIX.
C           COR - IS THE CORRELATION MATRIX.
C           VMN - IS A VECTOR CONTAINING VARIABLE MEANS.
C           STD - IS A VECTOR CONTAINING VARIABLE STANDARD DEVIATIONS.
C           CO - IS A VECTOR DIMENSIONED AT LEAST MV
C           NAMES - IS A VECTOR CONTAINING VARIABLE NAMES
C
C     USER HAS MANY OPTIONS AVAILABLE IN STEPWISE REGRESSION,
C     AMONG WHICH ARE OUTPUT FOR ANALYSIS OF VARIANCE AND THE DURBIN-WATSON
C     TEST FOR AUTO CORRELATION.  LIKE THE MULTIPLE REGRESSION, OPTIONS
C     ARE AVAILABLE FOR STORING THE RESIDUALS. 
C     PROGRAM WAS ORIGINALLY OBTAINED FROM WAYNE STATE UNIVERSITY,
C     MODIFICATIONS MADE FOR W.M.U., AND FINAL MODIFICATION FOR STP.
C     ORIGINAL SOURCE DOCUMENT FOR STEPWISE REGRESSION WAS:
C     RALSTON AND WILF "MATHEMATICAL METHODS FOR DIGITAL COMPUTERS".
C
      SUBROUTINE STSTRG(NV,NC,MV,MC,DATA,COR,VMN,STD,CO,NAMES)
      DIMENSION DATA(MC,MV),COR(MV,MV),VMN(1),IV(20),DCOR(20,20)
      DIMENSION SP(20),SIGCO(20),STUDT(20),BETA(20),STD(1)
      DIMENSION OPT(6),AVG(20),INDEX(20),SIGMA(20),COEN(20)
      DIMENSION NAMES(1),IVA(20),CO(1)
      COMMON /DEV/ICC,IDATA,IOUT,IDLG,IDSK
      COMMON /PRNT/ LINPP,ICOPS,RUNPRG
      COMMON/EXTRA/HEDR(70),NSZ
C
C     ************************************************************
C
C      OPTIONS ARE ENTERED AND INTERUPTED
C
21    TOL=.0001
      COE=0
      EFIN=0
      EFOUT=0
      NOFRC=0
      OBSER=NC
      IF(ICC.NE.2) WRITE(IDLG,16)
16    FORMAT('0LIST THE OPTIONS YOU WISH SEPARATED BY COMMAS'/)
      READ (ICC,17)(SP(I),I=1,14)
17    FORMAT(14(A5,1X))
      IF(SP(1).EQ.'!') RETURN
      DO 22 I=1,6
22    OPT(I)=0
      DO 18 I=1,14
      EXIT=0
40    IF(SP(I).EQ.'    ') GO TO 41
      IF(SP(I).NE.'HELP') GO TO 19
      WRITE(IDLG,20)
20    FORMAT('0POSSIBLE OPTIONS ARE:'/
     1' "TOLER" - INDICATE A TOLERANCE OTHER THAN .0001'/
     2' "F-VAL" - INDICATE F-VALUES FOR ENTERING OR OMITTING',
     3' A VARIABLE'/
     4' "FORCE" - INDICATE VARIABLES TO BE FORCED INTO THE REGRESSION'
     5/' "DURWT" - OUTPUT DURBIN-WATSON TEST FOR AUTOCORRELATION'/
     6' "ANOVA" - ANALYSIS OF VARIANCE'/
     7' "RESID" - STORE RESIDUALS'/
     8'0IF NO OPTIONS ARE DESIRED TYPE A RETURN')
      GO TO 21
19    IF(SP(I).NE.'TOLER') GO TO 25
      IF(OPT(1).NE.0) WRITE(IDLG,23)
23    FORMAT('0OPTION "TOLER" LISTED MORE THAN ONCE')
      OPT(1)=1.
      GO TO 18
25    IF(SP(I).NE.'F-VAL') GO TO 27
      IF(OPT(2).NE.0) WRITE(IDLG,26)
26    FORMAT('0OPTION "F-VAL" LISTED MORE THAN ONCE')
      OPT(2)=1
      GO TO 18
27    IF(SP(I).NE.'FORCE') GO TO 30
      IF(OPT(3).NE.0) WRITE(IDLG,28)
28    FORMAT('0OPTION "FORCE" LISTED MORE THAN ONCE')
      OPT(3)=1
      GO TO 18
30    IF(SP(I).NE.'DURWT') GO TO 32
      IF(OPT(4).NE.0) WRITE(IDLG,31)
31    FORMAT('0OPTION "DURWT" LISTED MORE THAN ONCE')
      OPT(4)=1
      GO TO 18
32    IF(SP(I).NE.'ANOVA') GO TO 34
      IF(OPT(5).NE.0) WRITE(IDLG,33)
33    FORMAT('0OPTION "ANOVA" LISTED MORE THAN ONCE')
      OPT(5)=1
      GO TO 18
34    IF(SP(I).NE.'RESID') GO TO 36
      IF(OPT(6).NE.0) WRITE(IDLG,35)
35    FORMAT('0OPTION "RESID" LISTED MORE THAN ONCE')
      OPT(6)=1
      GO TO 18
36    WRITE(IDLG,37) SP(I)
37    FORMAT('0OPTION "',A5,'" DOES NOT EXIST')
      GO TO 21
18    CONTINUE
C
41    IF(OPT(1).NE.1) GO TO 44
      IF(ICC.NE.2) WRITE(IDLG,42)
42    FORMAT(' ENTER TOLERANCE ',$)
      READ(ICC,43)TOL
C
C     ************************************************************
C
43    FORMAT(F)
44    IF(OPT(2).NE.1) GO TO 60
      IF(ICC.NE.2) WRITE(IDLG,45)
45    FORMAT(' ENTER F-VALUE FOR ENTERING A VARIABLE? ',$)
      READ (ICC,43)EFIN
      IF(ICC.NE.2) WRITE(IDLG,46)
46    FORMAT(' ENTER F-VALUE FOR OMITTING A VARIABLE? ',$)
      READ(ICC,43) EFOUT
C
60    IF(OPT(6).NE.1) GOTO 3
57    IF(NV.EQ.MV) GO TO 56
62    IF(ICC.NE.2) WRITE(IDLG,61)
61    FORMAT('0WHICH VARIABLE DO YOU WANT THE RESIDUALS STORED UNDER? ',
     1$)
      GO TO 58
56    WRITE(IDLG,63)
63    FORMAT('0RESIDUALS CANNOT BE FIT IN WITHOUT OVERLAYING A VARIABLE'
     1/' PLEASE INDICATE WHICH ONE',$)
58    CALL ALPHA(INDXKK,1,J,IRET,IHELP,IERR,NAMES,MV)
      IF(IERR.EQ.1) GO TO 57
      IF(IRET.EQ.1) RETURN
      IF(IHELP.NE.1) GO TO 320
      WRITE(IDLG,59)
59    FORMAT(' ENTER VARIABLE NAME OR NUMBER WHERE RESIDUALS CAN BE'/
     1' STORED.  NAME ASSIGNED TO RESIDUALS WILL BE "RESID"')
      GO TO 57
320    IF(INDXKK.GT.0) GO TO 3
      WRITE(IDLG,38)
      GO TO 57
3     IF(ICC.NE.2) WRITE(IDLG,1)
1     FORMAT('0LIST THE INDEPENDENT VARIABLES?'/)
      CALL ALPHA(IVA,19,NIND,IRET,IHELP,IERR,NAMES,NV)
      IF(IRET.EQ.1) RETURN
      IF(IERR.EQ.1) GO TO 3
      IF(IHELP.EQ.1) GO TO 3
4     K=1
      DO 5 I=1,NIND
      IV(I)=IVA(I)
      IF(IVA(I).GT.0) GO TO 5
      IV(I)=K
      K=K+1
5     CONTINUE
      DO 6 I=1,NIND-1
      IF(IVA(I).LT.0) GO TO 6
      DO 7 J=I+1,NIND
      IF(IVA(J).LT.0) GO TO 7
      IF(IVA(I).NE.IVA(J)) GO TO 7
      WRITE(IDLG,8)
8     FORMAT(' THE SAME INDEP. VARIABLE HAS BEEN LISTED TWICE')
      GO TO 3
7     CONTINUE
6     CONTINUE
      IF((NIND+1).LE.NV) GO TO 11
      WRITE(IDLG,9)
9     FORMAT(' YOU HAVE LISTED MORE VARIABLES THAN POSSIBLE'/
     1' WITH THE DATA SET AVAILABLE')
      GO TO 3
11    NIND=NIND+1
13    IF(ICC.NE.2) WRITE(IDLG,10)
10    FORMAT(' WHICH IS THE DEPENDENT VARIABLE? ',$)
      CALL ALPHA(J,1,I,IRET,IHELP,IERR,NAMES,NV)
      IF(IRET.EQ.1) RETURN
      IF(IERR.EQ.1) GOTO 13
      IF(IHELP.EQ.1) GOTO 13
      IVA(NIND)=J
      IV(NIND)=IVA(NIND)
      IF(IVA(NIND).LT.0) IV(NIND)=1
      IF(IVA(NIND).LT.0) GO TO 50
      DO 12 J=1,NIND-1
      IF(IVA(J).LT.0) GO TO 12
      IF(IVA(J).NE.IVA(NIND)) GO TO 12
      WRITE(IDLG,2)
2     FORMAT(' THE DEPENDENT VARIABLE  ALSO EXISTS AS AN INDEP. VAR.')
      GO TO 3
12    CONTINUE
50    IF(OPT(3).NE.1) GO TO 318
55    IF(ICC.NE.2) WRITE(IDLG,51)
51    FORMAT(' ENTER THE VARIABLES TO BE FORCED INTO',
     1' THE REGRESSION'/)
      CALL ALPHA(INDEX,20,NOFRC,IRET,IHELP,IERR,NAMES,NV)
      IF(IRET.EQ.1) RETURN
      IF(IERR.EQ.1) GO TO 55
      IF(IHELP.NE.1) GO TO 47
      WRITE(IDLG,48)
48    FORMAT(' ENTER THE VARIABLE NAME OR NUMBER TO BE FORCED INTO'/
     1' THE REGRESSION.  ONLY THOSE VARIABLES WHCIH WERE LISTED AS'/
     2' INDEPENDENT VARIABLES MAY BE USED.  *, ?, ALL ARE NOT'/
     3' ACCEPTABLE')
      GO TO 55
47    DO 49 K=1,NOFRC
      IF(INDEX(I).GT.0) GO TO 39
      WRITE(IDLG,38)
38    FORMAT(' ?, *, ALL MAY NOT BE USED HERE')
      GO TO 55
39    DO 53 J=1,NIND-1
      IF(INDEX(K).EQ.IVA(J)) GO TO 52
53    CONTINUE
      WRITE(IDLG,54) NAMES(INDEX(K))
54    FORMAT(' VARIABLE "',A5,'" WAS NOT LISTED AS A DEP. VAR.')
      GO TO 55
52    INDEX(K)=J
49    CONTINUE
      GO TO 318
C
C      *************************************************************
C
C     HERE THE DATA IS SET UP REPLACEING ALL *,?, AND  ALL WITH VALID
C     VARIABLES
C
300   J=NIND-1
314   IF(IVA(J).GT.0) GO TO 315
      IV(J)=IV(J)+1
      IF(IV(J).LE.NV) GO TO 316
315   J=J-1
      IF(J.GE.1) GO TO 314
      IF(IVA(NIND).GT.0) RETURN
      IV(NIND)=IV(NIND)+1
      IF(IV(NIND).GT.NV) RETURN
      K=1
      DO 312 I=1,NIND-1
      IF(IVA(I).GT.0) GO TO 312
      IV(I)=K
      K=K+1
312   CONTINUE
      J=NIND
      GO TO 318
316   K=IV(J)
      IF(J.EQ.NIND-1) GO TO 318
      DO 317 I=J+1,NIND-1
      IF(IVA(I).GT.0) GO TO 317
      K=K+1
      IF(K.GT.NV) GO TO 315
      IV(I)=K
317   CONTINUE
318   DO 313 I=1,NIND-1
      DO 313 K=I+1,NIND
      IF(IV(I).EQ.IV(K)) GO TO 300
313   CONTINUE
C
C     ********************************************************
C
14    DO 15 I=1,NIND
      K=IV(I)
      AVG(I)=VMN(K)
      SIGMA(I)=STD(K)*SQRT(NC-1.)
      DO 15 J=1,I
      L=IV(J)
      DCOR(I,J)=COR(K,L)
15    DCOR(J,I)=DCOR(I,J)
65    LDVAR=IV(NIND)
      ISTRE=O
      KDV1=0
      IDW=1
      SDEV2=0
      IDIF=0
      KTR=0
      SDIF2=0.0
      IF(IOUT.NE.21) WRITE(IOUT,5566)(HEDR(I),I=1,NSZ)
5566  FORMAT('1',70A1)
      IF(IOUT.EQ.21) CALL PRNTHD
      WRITE(IOUT,67)
67    FORMAT('0',10X,'***** STEPWISE REGRESSION *****')
      WRITE(IOUT,66)NIND,NAMES(IV(NIND))
66    FORMAT('0',I2,' VARIABLES; VARIABLE: ',A5,' IS DEPENDENT')
      LINES=6
      DEFR=NC-1.
      IF(OPT(3).NE.1) GO TO 76
      LINES=LINES+(NOFRC+4)/5
      WRITE(IOUT,70)NOFRC,(NAMES(IV(INDEX(I))),I=1,NOFRC)
70    FORMAT('0',I2,' VARIABLES FORCED INTO REGRESSION: ',
     15(A5,1X)/(38X,5(A5,1X)))
C
C     SOLVE SYSTEM OF EQUATIONS FOR FORCED VARIABLES
C
      DO 71 L=1,NOFRC
      K=INDEX(L)
      DIAG=DCOR(K,K)
      DO 72 I=1,NIND
      IF(I.EQ.K) GO TO 72
      FACTR=DCOR(I,K)/DIAG
      DO 73 J=1,NIND
      IF(J.EQ.K) GO TO 73
      DCOR(I,J)=DCOR(I,J)-(FACTR*DCOR(K,J))
73    CONTINUE
72    CONTINUE
      DO 74 I=1,NIND
74    DCOR(I,K)=-DCOR(I,K)/DIAG
      DO 75 J=1,NIND
75    DCOR(K,J)=DCOR(K,J)/DIAG
71    DCOR(K,K)=1.0/DIAG
      FACTR=NOFRC
      DEFR=DEFR-FACTR
      FLEVL=(DCOR(K,NIND)**2/DIAG)*DEFR/DCOR(NIND,NIND)
      NOENT=K
      NOSTP=NOFRC-1
      GO TO 77
76    NOSTP=-1
C
C     SELECT VARIABLES BY ANAL OF VAR AND SOLVE EQUATIONS
C
77    NOSTP=NOSTP+1
      IF(DCOR(NIND,NIND).LT.0.0) WRITE(IOUT,78)
78    FORMAT('0LAST DIAGONAL ELEMENT IS NEGATIVE')
      IF(DCOR(NIND,NIND).EQ.0.0) WRITE(IOUT,79)
79    FORMAT('0LAST DIAGONAL ELEMENT IS ZERO')
      IF(DCOR(NIND,NIND).LE.0.0) GO TO  300
      SIGY=SIGMA(NIND)*SQRT(DCOR(NIND,NIND)/DEFR)
      DEFR=DEFR-1.0
      IF(DEFR.GT.0.0) GO TO 80
      WRITE(IDLG,81)DEFR
81    FORMAT('0DEGREES OF FREEDOM IS',F8.2)
      GO TO 300
80    VMIN=-10000.0
      VMAX=0.0
      NOIN=0
      DO 85 I=1,NIND-1
      DIAG=DCOR(I,I)
      IF(DIAG.GT.0.0) GO TO 83
      IF(DIAG+TOL.GE.0) GO TO 85
      WRITE(IDLG,82)NAMES(IV(I))
82    FORMAT('0NEGATIVE DIAGONAL ELEMENT BELOW TOLERANCE VARIABLE: ',A5)
      GO TO 300
83    IF((DIAG-TOL).LT.0.0) GO TO 85
      VAR=DCOR(I,NIND)*DCOR(NIND,I)/DIAG
      SIG=SIGMA(I)
      IF(VAR.EQ.0.0) GO TO85
      IF(VAR.GT.0.0) GO TO 89
      IF(SIG.LT.0.0)GO TO 86
      IF(SIG.GT.0.0) GO TO 87
90    WRITE(IDLG,91)SIG,NAMES(IV(I))
91    FORMAT('0STANDARD DEVIATION IS',G,' FOR VAR ',A5)
      GO TO 300
86    SIG=-SIG
      GO TO 88
87    IF(VAR.LE.VMIN) GO TO 88
      VMIN=VAR
      NOMIN=I
88    NOIN=NOIN+1
      INDEX(NOIN)=I
      COEN(NOIN)=DCOR(I,NIND)*SIGMA(NIND)/SIG
      SIGCO(NOIN)=(SIGY/SIG)*SQRT(DIAG)
      GO TO 85
89    IF(SIG.LE.0.0) GO TO 90
      IF(VAR.LE.VMAX) GO TO 85
      VMAX=VAR
      NOMAX=I
85    CONTINUE
      IF(NOIN.LT.0) GO TO 110
      IF(NOIN.GT.0) GO TO 94
      KTR=KTR+1
      IF(IOUT.NE.21) GO TO 400
      LINES=LINES+2
      IF(LINES.LE.LINPP) GO TO 400
      CALL PRNTHD
      LINES=4
400   WRITE(IOUT,92) SIGY
92    FORMAT('0STANDARD ERROR OF Y =',G)
      IF(KTR.LE.1) GO TO 96
113   IF(IOUT.NE.21) GO TO 401
      LINES=LINES+2
      IF(LINES.LE.LINPP) GO TO 401
      CALL PRNTHD
      LINES=4
401   WRITE(IOUT,93)
93    FORMAT('0ERROR, NO VARIABLE CAN ENTER THE EQUATION')
      RETURN
94    CNST=AVG(NIND)
      DO 95 J=1,NOIN
      L=INDEX(J)
95    CNST=CNST-(COEN(J)*AVG(L))
      FLEV=FLEVL
      NOEN=NOENT
      K1=K
      IF(VMIN.EQ.-10000.0) GO TO 96
      FLEVL=-VMIN*(DEFR+1.0)/DCOR(NIND,NIND)
      IF(EFOUT.LE.FLEVL) GO TO 96
      K=NOMIN
      DEFR=DEFR+2.0
      NOENT=0
      GO TO 97
96    FLEVL=VMAX*DEFR/(DCOR(NIND,NIND)-VMAX)
      IF(EFIN.GE.FLEVL) GO TO 112
      K=NOMAX
      NOENT=K
97    COEFDT=1.0-DCOR(NIND,NIND)
      COD=COEFDT-COE
      COE=COEFDT
      RRR=SQRT(COEFDT)
      IF(NOIN.LT.0) GO TO 110
      IF(NOIN.EQ.0) GO TO 104
      IF(NOEN.LE.0) GO TO 99
      IF(IOUT.NE.21) GO TO 402
      LINES=LINES+6
      IF(LINES.LE.LINPP) GO TO 402
      CALL PRNTHD
      LINES=8
402   WRITE(IOUT,98)NOSTP,NAMES(IV(K1))
98    FORMAT(///'0STEP NO.',I3/' ENTERING VARIABLE: ',A5)
      ISTRE=K1
      GO TO 102
99    IF(K1.NE.ISTRE) GO TO 100
      IF(K1.NE.0) GO TO 119
100   IF(IOUT.NE.21) GO TO 403
      LINES=LINES+6
      IF(LINES.LE.LINPP) GO TO 403
      CALL PRNTHD
      LINES=8
403   WRITE(IOUT,101)NOSTP,NAMES(IV(K1))
101   FORMAT(///'0STEP NO.',I3/' VARIABLE: ',A5,' REMOVED')
102   NDGN=1
      NDGD=NC-NOIN+1
      PROB=FISHER(NDGN,NDGD,FLEV)
      IF(IOUT.NE.21) GO TO 404
      LINES=LINES+6
      IF(LINES.LE.LINPP) GO TO 404
      CALL PRNTHD
      LINES=8
404   WRITE(IOUT,103)FLEV,PROB,SIGY,COEFDT,RRR,COD,CNST
103   FORMAT(3X,'F-LEVEL',2X,G,'  WITH PROB. ',F7.4/
     13X,'STANDARD ERROR OF',' ESTIMATE',G/
     23X,'COEFFICIENT OF DETERMINATION =',G/
     23X,'COEFFICIENT OF MULTIPLE REGRESSION =',G/
     33X,'INCREASE IN COEFFICIENT OF DETERMINATION =',G/
     43X,'CONSTANT',G13.5)
      IF(IOUT.NE.21) GO TO 407
      LINES=LINES+2
      IF(LINES.LE.(LINPP-5)) GO TO 407
      CALL PRNTHD
      LINES=4
407   WRITE(IOUT,501)
501   FORMAT('0VARIABLE',3X,'COEFFICIENT',5X,'STD ERROR OF COEFF.')
      DO 502 J=1,NOIN
      IF(IOUT.NE.21) GO TO 408
      LINES=LINES+1
      IF(LINES.LE.LINPP) GO TO 408
      CALL PRNTHD
      WRITE(IOUT,501)
      LINES=5
408   WRITE(IOUT,503) NAMES(IV(INDEX(J))),COEN(J),SIGCO(J)
503   FORMAT(2X,A5,5X,G15.7,1X,G15.7)
502   CONTINUE
104   DIAG=DCOR(K,K)
      DO 105 I=1,NIND
      IF(I.EQ.K) GO TO 105
      FACTR=DCOR(I,K)/DIAG
      DO 106 J=1,NIND
      IF(J.EQ.K) GO TO 106
      DCOR(I,J)=DCOR(I,J)-(FACTR*DCOR(K,J))
106   CONTINUE
105   CONTINUE
      DO 107 I=1,NIND
107   DCOR(I,K)=-DCOR(I,K)/DIAG
      DO 108 J=1,NIND
108   DCOR(K,J)=DCOR(K,J)/DIAG
      DCOR(K,K)=1.0/DIAG
      IF(NOSTP.LT.250) GO TO 77
      WRITE(IDLG,109) NOSTP
109   FORMAT('0ERROR, NUMBER OF STEPS =',I4,' CHECK F-VALUES')
      GO TO 300
110   WRITE(IDLG,111)NOIN
111   FORMAT('0ERROR, NOIN IS',I7)
      GO TO 300
C
C     FINAL OUTPUT
C
112   IF(NOIN.EQ.0) GO TO 113
      DO 114 I=1,NOIN
      J=INDEX(I)
      STUDT(I)=COEN(I)/SIGCO(I)
114   BETA(I)=DCOR(J,NIND)
      IF(NOEN.GE.0) GO TO 116
      WRITE(IDLG,115)NOEN
115   FORMAT('0ERROR, NOEN=',I6)
      GO TO 300
116   IF(NOEN.EQ.0) GO TO 118
      IF(IOUT.NE.21) GO TO 405
      LINES=LINES+6
      IF(LINES.LE.LINPP) GO TO 405
      CALL PRNTHD
      LINES=8
405   WRITE(IOUT,98)NOSTP,NAMES(IV(K1))
      ISTRE=K1
      GO TO 123
118   IF(K1.NE.ISTRE) GO TO 121
      IF(L1.EQ.0) GO TO 121
119   WRITE(IDLG,120)
120   FORMAT(' VARIABLE REMOVED SAME AS VARIABLE ENTERED')
      GO TO 300
121   IF(IOUT.NE.21) GO TO 406
      LINES=LINES+6
      IF(LINES.LE.LINPP) GO TO 406
      CALL PRNTHD
      LINES=8
406   WRITE(IOUT,101) NOSTP,NAMES(IV(K1))
123   COEFDT=1.0-DCOR(NIND,NIND)
      COD=COEFDT-COE
      RRR=SQRT(COEFDT)
      NDGN=1
      NDGD=NC-NOIN+1
      PROB=FISHER(NDGN,NDGD,FLEV)
      IF(IOUT.NE.21) GO TO 409
      LINES=LINES+6
      IF(LINES.LE.LINPP) GO TO 409
      CALL PRNTHD
      LINES=8
409   WRITE(IOUT,103)FLEV,PROB,SIGY,COEFDT,RRR,COD,CNST
      IF(IOUT.NE.21) GO TO 410
      LINES=LINES+3
      IF(LINES.LE.(LINPP-5)) GO TO 410
      CALL PRNTHD
      LINES=5
410   WRITE(IOUT,504)
504   FORMAT('0',24X,'STD.ERR.',5X,'STANDARDIZED'/
     1' VARIABLE',4X,'COEFF.',6X,'OF COEFF.',4X,'COEFFICIENT',
     25X,'T-VALUE',5X,'PROB.')
      DO 1100 I=1,NOIN
      NDGN=1
      NDGD=NC-NOIN-1
      PROB=STUDT(I)**2
      PROB=FISHER(NDGN,NDGD,PROB)
      IF(IOUT.NE.21) GO TO 1100
      LINES=LINES+1
      IF(LINES.LE.LINPP) GO TO 1100
      CALL PRNTHD
      WRITE(IOUT,504)
      LINES=6
1100  WRITE(IOUT,1101)NAMES(IV(INDEX(I))),COEN(I),SIGCO(I),BETA(I),
     1STUDT(I),PROB
1101  FORMAT(2X,A5,4X,G12.4,1X,G12.4,1X,G,1X,G11.4,1X,F6.3)
      IF((OPT(4).NE.1).AND.(OPT(5).NE.1).AND.(OPT(6).NE.1)) GO TO 300
      IF(OPT(6).NE.1) GO TO 135
      DO 200 I=1,NV
200   CO(I)=0
      IF(INDXKK.LE.NV) GO TO 135
      INDXKK=NV+1
      NV=INDXKK
135   NAMES(INDXKK)='RESID'
      SSS=0
      SST=0
      DO 125 J=1,NC
      YPRED=CNST
      SSS=DATA(J,IV(NIND))+SSS
      SST=SST+DATA(J,IV(NIND))**2
      DO 126 I=1,NOIN
      K=IV(INDEX(I))
126   YPRED=YPRED+COEN(I)*DATA(J,K)
      DEV=DATA(J,IV(NIND))-YPRED
      SDEV2=SDEV2+DEV**2
      IF(IDIF.EQ.0) GO TO 127
      DIF=DEV-DEVP
      SDIF2=SDIF2+DIF*DIF
      GO TO 128
127   IDIF=1
128   DEVP=DEV
      IF(OPT(6).NE.1)GO TO 125
      DATA(J,INDXKK)=DEV
      DO 201 I=1,NV
201   CO(I)=DATA(J,I)*DEV+CO(I)
125   CONTINUE
      IF(OPT(6).NE.1) GO TO 202
      VMN(INDXKK)=0
      STD(INDXKK)=SQRT(CO(INDXKK)/(NC-1.))
      DO 203 I=1,NV
      IF(I.EQ.INDXKK) GO TO 205
      IF((STD(I)*STD(INDXKK)).EQ.0) GO TO 204
      COR(I,INDXKK)=CO(I)/((NC-1.)*STD(I)*STD(INDXKK))
      GO TO 203
204   COR(I,INDXKK)=0
      GO TO 203
205   COR(I,INDXKK)=1.
203   COR(INDXKK,I)=COR(I,INDXKK)
202   IF(OPT(5).NE.1) GO TO 132
      SST=SST-SSS*SSS/NC
      SSR=SST-SDEV2
      NERR=NC-(NOIN+1)
      NTOTL=NERR+NOIN
      ERMS=SDEV2/NERR
      RMS=SSR/NOIN
      F=RMS/ERMS
      PROB=FISHER(NOIN,NERR,F)
      IF(IOUT.NE.21) GO TO 411
      LINES=LINES+6
      IF(LINES.LE.LINPP) GO TO 411
      CALL PRNTHD
      LINES=8
411   WRITE(IOUT,130)
130   FORMAT('0ANALYSIS OF VARIANCE')
      WRITE(IOUT,131)SSR,NOIN,RMS,F,PROB,SDEV2,NERR,ERMS,SST,NTOTL
131   FORMAT(4X,'SOURCE',6X,'SUM OF SQ.',6X,'DF',5X,
     1'MEAN SQ.',6X,'F',7X,'PROB.'/' REGRESSION',F17.6,I6,F15.5,F9.3
     2,2X,F7.4/6X,'ERROR',F17.6,I6,F15.5/6X,'TOTAL',F17.6,I6)
132   IF(OPT(4).NE.1) GO TO 300
      DWS=SDIF2/SDEV2
      DW1=4.-DWS
      IF(IOUT.NE.21) GO TO 412
      LINES=LINES+7
      IF(LINES.LE.LINPP) GO TO 412
      CALL PRNTHD
      LINES=9
412   WRITE(IOUT,133)SDIF2,DWS,DW1
133   FORMAT(//14X,'DURBIN-WATSON TEST FOR AUTOCORRRELATION'//
     110X,'SUM OF SQUARED DIFFERENCES =',G/10X,'DWS=',G/
     210X,'4 - DWS =',G)
      GO TO 300
      END
C                                           *** STAT PACK ***
C     SUBROUTINE FOR PARTIAL CORRELATIONS
C     CALLING SEQUENCE: CALL PCORR(NV,NC,MV,MC,COR,SP,NAMES)
C     WHERE NV - IS THE NUMBER OF COLUMNS ACTUALLY FILLED(VARIABLES)
C           NC - IS THE NUMBER OF ROWS ACTUALLY FILLED(OBSERVATIONS,CASES)
C           MV - IS MAXIMUM NUMBER OF COLUMNS, AS SPECIFIED IN MAIN.
C           MC - IS MAXIMUM NUMBER OF ROWS, AS SPECIFIED IN MAIN.
C           COR - IS THE CORRELATION MATRIX.
C           SP - IS A VECTOR DIMENSIONED AT LEAST NV.
C           NAMES - IS A VECTOR CONTAINING VARIABLE NAMES
C
C     GIVES PARTIAL CORRELATIONS MATRIX FOR SPECIFIED VARIABLES.
C
      SUBROUTINE PCORR(NV,NC,MV,MC,COR,SP,NAMES)
      COMMON /DEV/ICC,IDATA,IOUT,IDLG,IDSK
      COMMON /PRNT/ LINPP,ICOPS,RUNPRG
      COMMON/EXTRA/HEDR(70),NSZ
      DIMENSION COR(MV,MV),SP(1),DCOR(20,20),XVER(20,20)
      DIMENSION IV(20),NAMES(1),IVA(20)
      ISQ=7
      IF(IOUT.EQ.21) ISQ=15
1     IF(ICC.NE.2) WRITE(IDLG,2)
2     FORMAT('0WHICH VARIABLES? ',$)
      CALL ALPHA(IVA,20,NN,IRET,IHELP,IERR,NAMES,NV)
      IF(IRET.EQ.1) RETURN
      IF(IERR.EQ.1) GO TO 1
      IF(IHELP.EQ.1) GO TO 1
      K=1
      DO 3 I=1,NN
      IV(I)=IVA(I)
      IF(IVA(I).GT.0) GO TO 3
      IV(I)=K
      K=K+1
3     CONTINUE
      DO 4 I=1,NN-1
      IF(IVA(I).LT.0) GOTO 4
      DO 5 J=I+1,NN
      IF(IVA(I).LT.0) GO TO 5
      IF(IVA(I).NE.IVA(J)) GO TO 5
      WRITE(IDLG,6)
6     FORMAT(' THE SAME VARIABLE HAS BEEN LISTED TWICE')
      GO TO 1
5     CONTINUE
4     CONTINUE
      IF(NN.LE.NV) GO TO 58
      WRITE(IDLG,7)
7     FORMAT(' YOU HAVE LISTED MORE VARIABLES THAN IS POSSIBLE'/
     1' WITH THE DATA SET AVAILABLE')
      GO TO 1
50    J=NN
54    IF(IVA(J).GT.0) GOTO 55
      IV(J)=IV(J)+1
      IF(IV(J).LE.NV) GO TO 56
55    J=J-1
      IF(J.GE.1) GO TO 54
      RETURN
56    K=IV(J)
      IF(J.EQ.NN) GOTO 58
      DO 57 I=J+1,NN
      IF(IVA(I).GT.0) GO TO 57
      K=K+1
      IF(K.GT.NV) GO TO 55
      IV(I)=K
57    CONTINUE
58    DO 53 I=1,NN-1
      DO 53 K=I+1,NN
      IF(IV(I).EQ.IV(K)) GO TO 50
53    CONTINUE
      DO 8 I=1,NN
      DO 8 J=1,NN
8     DCOR(I,J)=COR(IV(I),IV(J))
C
C     *****************************************************************
C     INVERSE BY LINEAR ROW TRANSFORMATION
C
      DO 9 I=1,NN
      DO 10 J=1,NN
10    XVER(I,J)=0
9     XVER(I,I)=1.0
      DO 11 I=1,NN
      IF(((DCOR(I,I)+100.)-100.).NE.0.0) GO TO 15
      IF(I.EQ.NN) GO TO 30
      DO 12 J=I+1,NN
      IF(((DCOR(J,I)+100.)-100.).NE.0.0) GO TO 13
12    CONTINUE
30     WRITE(IDLG,31)
31     FORMAT('0THE INVERSE DOES NOT EXIST')
      GO TO 50
13    DO 14 K=1,NN
      DCOR(I,K)=DCOR(I,K)+DCOR(J,K)
14    XVER(I,K)=XVER(I,K)+XVER(J,K)
15    G=DCOR(I,I)
      DO 16 J=1,NN
      DCOR(I,J)=DCOR(I,J)/G
16    XVER(I,J)=XVER(I,J)/G
      DO 17 L=1,NN
      IF(L.EQ.I) GO TO 17
      G=DCOR(L,I)
      DO 18 J=1,NN
      DCOR(L,J)=DCOR(L,J)-G*DCOR(I,J)
18    XVER(L,J)=XVER(L,J)-G*XVER(I,J)
17    CONTINUE
11    CONTINUE
C
C     ***************************************************************
C
      IF(IOUT.NE.21) WRITE(IOUT,5566)(HEDR(I),I=1,NSZ)
5566  FORMAT('1',70A1)
      IF(IOUT.EQ.21) CALL PRNTHD
      WRITE(IOUT,19)
19    FORMAT('0',10X,'***** PARTIAL CORRELATION MATRIX *****')
      LINES=4
      DO 20 I=2,NN
      DO 21 J=1,I-1
21    SP(J)=-XVER(I,J)/SQRT(XVER(I,I)*XVER(J,J))
      IF(IOUT.NE.21) GO TO 27
      L=(I+ISQ-2)/ISQ+1
      LINES=LINES+L
      IF(LINES.LE.(LINPP-L-1)) GO TO 27
      WRITE(IOUT,25)
      DO 28 J=1,I-2,ISQ
      NEND=J+ISQ-1
      IF(NEND.GT.(I-2)) NEND=I-2
28    WRITE(IOUT,23)(NAMES(IV(J)),K=J,NEND)
      CALL PRNTHD
      LINES=2+L
27    DO 20 K=1,I-1,ISQ
      NEND=K+ISQ-1
      IF(NEND.GT.I-1) NEND=I-1
      IF(K.EQ.1)WRITE(IOUT,22) NAMES(IV(I)),(SP(J),J=K,NEND)
22    FORMAT('0',A5,1X,15F8.4)
      IF(K.NE.1) WRITE(IOUT,26) (SP(J),J=K,NEND)
26    FORMAT(7X,15F8.4)
20    CONTINUE
      WRITE(IOUT,25)
25    FORMAT(1X)
      DO 24 I=1,NN,ISQ
      NEND=I+ISQ-1
      IF(NEND.GT.NN-1) NEND=NN-1
24    WRITE(IOUT,23)(NAMES(IV(J)),J=I,NEND)
23    FORMAT(7X,15(3X,A5))
      GO TO 50
      END
C     SUBROUTINE FOR WILCOXON MATCHED PAIRS SIGNED-RANK TEST.
C     CALLING SEQUENCE: CALL WILCX(NV,NC,MV,MC,DATA,IV,H,NAMES)
C     WHERE NV - IS THE NUMBER OF COLUMNS ACTUALLY FILLED(VARIABLES)
C           NC - IS THE NUMBER OF ROWS ACTUALLY FILLED(OBSERVATIONS,CASES)
C           MV - IS MAXIMUM NUMBER OF COLUMNS, AS SPECIFIED IN MAIN.
C           MC - IS MAXIMUM NUMBER OF ROWS, AS SPECIFIED IN MAIN.
C           DATA - STORAGE FOR DATA DIMENSIONED FOR MAXIMUM MATRIX.
C           IV - IS A VECTOR DIMENSIONED AT LEAST NC.
C           H - IS A VECTOR DIMENSIONED FOR AT LEAST NC.
C           NAMES - IS A VECTOR CONTAINING VARIABLE NAMES
C
C     PROGRAM AUTOMATICALLY OUTPUTS ENTIRE MATRIX, FORMAING IT ONE LINE
C     AT A TIME.  MAKES USE OF A PARTITIONING SORT.
C
      SUBROUTINE WILCX(NV,NC,MV,MC,DATA,IV,H,NAMES)
      DIMENSION DATA(MC,MV),IV(1),H(1),NAMES(1)
C     FOLLOWING DIMENSION IS FOR USE IN THE SORT
      DIMENSION IU(16),IL(16)
      COMMON/DEV/ICC,IDATA,IOUT,IDLG,IDSK
      COMMON /PRNT/ LINPP,ICOPS,RUNPRG
      COMMON/EXTRA/HEDR(70),NSZ
      INTEGER BASE
      IF(IOUT.NE.21) WRITE(IOUT,5566) (HEDR(I),I=1,NSZ)
5566  FORMAT('1',70A1)
      IF(IOUT.EQ.21) CALL PRNTHD
      WRITE(IOUT,1)
1     FORMAT('0',10X,'***** WILCOXON MATCHED-PAIRS SIGNED-RANK TEST',
     1' *****'/'0VAR. VS VAR.',7X,'MEAN',10X,'S.D.',4X,'N',6X
     2,'Z',10X,'T')
      LINES=6
      DO 2 I=1,NV-1
      DO 3 JJ=I+1,NV
      L=1
      DO 4 K=1,NC
      IF(DATA(K,I)-DATA(K,JJ))5,4,6
5     IV(L)=1
      H(L)=DATA(K,JJ)-DATA(K,I)
      GO TO 7
6     IV(L)=0
      H(L)=DATA(K,I)-DATA(K,JJ)
7     L=L+1
4     CONTINUE
C 
C     SORT BY PARTITION - ACM
C
      NN=L-1
      IF(NN.LT.2) GO TO 27
      M=1
      II=1
      J=NN
71    IF(II.GE.J) GO TO 78
72    K=II
      IJ=(J+II)/2
      T=H(IJ)
      IF(H(II).LE.T) GO TO 73
      H(IJ)=H(II)
      H(II)=T
      T=H(IJ)
      ISAV=IV(IJ)
      IV(IJ)=IV(II)
      IV(II)=ISAV
73    LL=J
      IF(H(J).GE.T) GO TO 75
      H(IJ)=H(J)
      H(J)=T
      T=H(IJ)
      ISAV=IV(IJ)
      IV(IJ)=IV(J)
      IV(J)=ISAV
      IF(H(II).LE.T) GO TO 75
      H(IJ)=H(II)
      H(II)=T
      T=H(IJ)
      ISAV=IV(IJ)
      IV(IJ)=IV(II)
      IV(II)=ISAV
      GO TO 75
74    H(LL)=H(K)
      H(K)=TT
      ISAV=IV(LL)
      IV(LL)=IV(K)
      IV(K)=ISAV
75    LL=LL-1
      IF(H(LL).GT.T) GO TO 75
      TT=H(LL)
76    K=K+1
      IF(H(K).LT.T) GO TO 76
      IF(K.LE.LL) GO TO 74
      IF((LL-II).LE.(J-K)) GO TO 77
      IL(M)=II
      IU(M)=LL
      II=K
      M=M+1
      GO TO 79
77    IL(M)=K
      IU(M)=J
      J=LL
      M=M+1
      GO TO 79
78    M=M-1
      IF(M.EQ.0) GO TO 20
      II=IL(M)
      J=IU(M)
79    IF((J-II).GE.11) GO TO 72
      IF(II.EQ.1) GO TO 71
      II=II-1
80    II=II+1
      IF(II.EQ.J) GO TO 78
      T=H(II+1)
      ISAV=IV(II+1)
      IF(H(II).LE.T) GO TO 80
      K=II
81    H(K+1)=H(K)
      IV(K+1)=IV(K)
      K=K-1
      IF(T.LT.H(K)) GO TO 81
      H(K+1)=T
      IV(K+1)=ISAV
      GO TO 80
C
C     END SORT BEGIN RANKING
C
20    SAV=H(1)
      SUM=1
      IVV=1
      DO 21 K=2,NN
      IF(H(K).EQ.SAV) GO TO 23
      H(K-1)=SUM/IVV
      DO 22 L=K-IVV,K-1
22    H(L)=H(K-1)
      SAV=H(K)
      SUM=K
      IVV=1
      GO TO 21
23    IVV=IVV+1
      SUM=SUM+K
21    CONTINUE
      H(NN)=SUM/IVV
      DO 24 L=NN-IVV+1,NN
24    H(L)=H(NN)
C
C     RANKING COMPLETE
C
C     FIND NUMBER OF ORIGINAL NEGATIVES
C
27    SD=0
      Z=0
      SAV=0
      IVV=0
      SUM=0
      IF(NN.LT.2) GO TO 28
      DO 25 L=1,NN
      IF(IV(L).EQ.0) GO TO 25
      SUM=SUM+H(L)
      IVV=IVV+1
25    CONTINUE
      SAV=(NN*(NN+1))/2
      IF((IVV*2).LT.NN) GO TO 26
C
C     IF IVE COUNTED THE WROND SIGN ADJUST THE SUM BY KNOWING 
C     THE TOTAL OF THE RANKS IS N(N+1)/2
C
      SUM=SAV-SUM
      IVV=NN-IVV
26    SAV=SAV/2
      SD=SQRT(((2.*NN+1.)*SAV)/6.)
      Z=(SUM-SAV)/SD
28    IF(IOUT.NE.21) GO TO 31
      LINES=LINES+1
      IF(LINES.LE.LINPP) GO TO 31
      CALL PRNTHD
      WRITE(IOUT,1)
      LINES=7
31    WRITE(IOUT,30) NAMES(I),NAMES(JJ),SAV,SD,NN,Z,SUM
30    FORMAT(1X,A5,2X,A5,F11.2,F14.2,I5,F8.2,F11.2)
3     CONTINUE
2     CONTINUE
      RETURN
      END

