C	WESTERN MICHIGAN UNIVERSITY
C	STEPR.F4 (FILE NAME ON LIBRARY DECTAPE)
C	STEPR, 1.3.2 (CALLING NAME, SUBLST #)
C	STEPWISE REGRESSION
C	ORIGINAL PROGRAM CAME FROM WAYNE STATE UNIVERSITY.  LATER
C	 IT WAS MODIFIED BY SAM ANEMA AND R.R. BARR.
C	LIBRARY DECTAPE PROG. USED:  USAGE.MAC
C	FORWMU PROGS. USED:  TTYPTY, EXISTS, DEVCHG, PRINTS, XPRODH
C	APLIB PROGS. USED:  IO, GETFOR, FISHER
C	INTERNAL SUBR. USED:  LANG2, TRANS, LANG1, REGGR, TNUM
C	ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS PUT IN BY WG
C
C
C
      DIMENSION X(70),SX(70),NN(70),SXY(70,70),XMEAN(70),IDEN(16)
      DIMENSION FMD(70),IFT(96),IT(70),KTRAN(50)
      DIMENSION KVAR1(50),KVAR2(50),KVAR3(50),CONST(50),IFOPT(15)
      DIMENSION VAR(70),IRROR(2)
      INP=2
      IOUT=3
      IND=-4
      NOUTD = -1
      WRITE(NOUTD,2021)
2021  FORMAT('1',20X,'WMU STEPWISE REGRESSION'//)
C	CALL USAGE('STEPR ')
C---------------TTYPTY RETURNS ZERO - TTY JOB, MINUS ONE - BATCH JOB;
C--------------- ICODE USED AS ARG. IN CALL IO.
      CALL TTYPTY(ICODE)
C---------------IDV1, IDV2 ARE RETURNED.  OTHER ARGS. ARE INPUT.
C--------------- 1 MEANS OUTPUT? PRINTS, 0 MEANS INPUT? PRINTS.
      CALL IO(IOUT,IDV1,NOUTD,IND,1,ICODE)
1111  NTR=0
C---------------ITYPE REFERS TO MISSING DATA REPLACEMENT METHOD
C--------------- SEE ST. 110.
      ITYPE=3
      CALL IO(INP,IDV2,NOUTD,IND,0,ICODE)
C---------------IFT, ISTD ARE RETURNED.  96=NO. OF OBJ. TIME FORMAT
C---------------WORDS. (6 LINES)  2 MEANS F-TYPE FORMAT ONLY.
      CALL GETFOR(NOUTD,IND,IFT,ISTD,96,2)
      IF(ISTD.EQ.1)IFT(1)='(10F)'
1     WRITE(NOUTD,2)
2     FORMAT(' ENTER NUMBER OF VARIABLES'/)
      READ(IND,3)NVAR
3     FORMAT(I)
      IF(NVAR.LT.71.AND.NVAR.GT.1)GO TO 5
      CALL TTYPTY(ICC)
      IF(ICC.EQ.-1)CALL EXIT
      GO TO 1
5     WRITE(NOUTD,4)
4     FORMAT(' ENTER IDENTIFICATION IF DESIRED OTHERWISE RETURN'/)
      READ(IND,11)(IDEN(I),I=1,16)
11    FORMAT(16A5)
      REWIND 1
121      WRITE(NOUTD,6)
6      FORMAT(' ENTER OPTIONS OR TYPE "HELP"'/)
      READ(IND,7)(IT(I),I=1,70)
7     FORMAT(70A1)
	IF(IT(1).NE.'H'.OR.IT(2).NE.'E'.OR.IT(3).NE.'L')GO TO 91
	WRITE(NOUTD,9991)
9991	FORMAT(/,' OPTION TABLE',//' CODE		DISCRIPTION',/,
     1	' ----		-----------',/,' MIS*		MISSING DATA',/,
     2	' TRA*		TRANSFORMATION OF VARIABLES',/,
     3	' MVS		MEANS STANDARD DEV., AND VARIABLES',/,
     4	' XPR		RAW SUMS OF SQUARES AND CROSS PRODUCTS',/,
     5	' COV		VARIANCE - COVARIANCE MATRIX',/,
     6	' COR		CORRELATION, FISHER Z, AND T MATRICES',/,
     7	' FVL		F - VALUES',/,
     8	' FOR		FORCE VARIABLES',/,
     9	' ELM		ELIMINATE VARIABLES',/,
     A	' RES		RESIDUALS',/,
     B	' ANA		ANALYSIS OF VARIANCE AT FINAL STAGE',/,
     C	' DUR		DURBIN - WATSON STATISTICS',/,
     D	' ZER		FORCE REGRESSION THROUGH ZERO',//,
     E	' * READ PERTINATE REMARKS IN SECTION 3.0',/,
     F	' OF LIBRARY PROGRAM WRITEUP #1.3.2 BEFORE USING.',//)
	GO TO 121
91      DO 8 I=1,15
8     IFOPT(I)=0
      CALL LANG1(IT,IFOPT)
C---------------IFOPT(I)=1 MEANS USER CHOOSES THIS OPTION.
C--------------- IFOPT(1), ... , IFOPT(15) MEANS MVS, COR, COV,
C--------------- RES, ANA, DUR, ZER, TRA, MIS, OBJ, FVL, FOR, ELM,
C--------------- XPR, TOL RESPECTIVELY - OBJ AND TOL DO NOT APPEAR
C--------------- IN WRITE UP IN HELP LIST.
      IF(IFOPT(9).NE.1)GO TO 120
      WRITE(NOUTD,102)
102   FORMAT(' IS THERE MORE THAN ONE MISSING DATA SYMBOL ?'/)
      READ(IND,103)MDS
103   FORMAT(A3)
      IF(MDS.EQ.'YES') GO TO 104
      WRITE(NOUTD,106)
106   FORMAT(' ENTER MISSING DATA SYMBOL'/)
      READ(IND,107)FMD(1)
107   FORMAT(10F)
      DO 108 I=2,NVAR
108   FMD(I)=FMD(1)
      GO TO 111
104   WRITE(NOUTD,109)
109   FORMAT(' ENTER MISSING DATA SYMBOLS'/)
      READ(IND,107)(FMD(I),I=1,NVAR)
111   WRITE(NOUTD,110)
110   FORMAT(' HOW WOULD YOU LIKE TO COMPENSATE FOR MISSING DATA? TYPE:'
     1/'     1 - TO REPLACE MISSING DATA BY MEAN VALUE'/
     2'     2 - TO DELETE THE OBSERVATION'/)
      READ(IND,3)ITYPE
      IF(ITYPE .EQ.1.OR.ITYPE.EQ.2)GO TO 120
      CALL TTYPTY(ICC)
      IF(ICC.EQ.-1)CALL EXIT
      GO TO 111
120   MAX=0
C---------------IFOPT(8)=1 MEANS CHOOSE TRANS
      IF(IFOPT(8).NE.1)GO TO 30
C---------------INDICES OF VARS. I, J, K AS PRESCIBED BY USER
C--------------- (SEE WRITE UP) FOR THE TRANSFORMATIONS ARE RETURNED
C--------------- IN KVAR1, KVAR2, KVAR3 RESPECTIVELY FOR LATER
C--------------- USE BY SUBR. TRANS.
C---------------KTRAN(I)=ONE OF 15 POSSIBLE TRANSF. LISTED IN SECTION
C--------------- 3.0 OF WRITE UP UNDER OPTION 2;  CONST=A AS
C--------------- DESCRIBED IN THE SAME SECTION.
      CALL LANG2(KTRAN,CONST,KVAR1,KVAR2,KVAR3,IT,NTR)
C---------------NTR=NO. OF TRANSF. AND IS RETURNED BY SUBR. LANG2.
      DO 9 I=1,NTR
      IF(KVAR1(I).GT.MAX)MAX=KVAR1(I)
9     CONTINUE
C---------------NVAR=NO. OF VARS, MAX=MAXIMUM INDEX OF VARS.
C--------------- GENERATED BY TRANSF.
30    NVTR=MAX0(NVAR,MAX)
40    REWIND 20
      REWIND 1
      DO 299 J=1,NVTR
      SX(J)=0.0
      NN(J)=0
      DO 300 I=1,NVTR
      SXY(I,J)=0.0
300   CONTINUE
299   CONTINUE
      K=0
C---------------IDV2 IS RETURNED BY IO.
      IF(IDV2.EQ.'TTY')GO TO 13
      WRITE(NOUTD,14)
14    FORMAT(' YOUR DATA IS BEING READ'/)
      GO TO 17
13    IF(ISTD.EQ.1)WRITE(NOUTD,15)
15    FORMAT(' ENTER DATA (AT MOST 10 PER LINE)'/)
      IF(ISTD.EQ.0)WRITE(NOUTD,16)
16    FORMAT(' ENTER DATA'/)
17    IRROR(1)=0
      IRROR(2)=0
      GO TO (301,351,391),ITYPE
C---------------HERE IF WE REPLACE MISSING DATA BY MEAN VALUE,
C--------------- SEE ST. 110+1.
301   READ(INP,IFT,END=340)(X(I),I=1,NVAR)
303   K=K+1
      DO 304 I=1,NVAR
C---------------FIND(I)=MISSING DATA SYMBOLS SEE ST. 109+1.
      IF(X(I).EQ.FMD(I)) GO TO 304
      SX(I)=SX(I)+X(I)
      NN(I)=NN(I)+1
304   CONTINUE
      WRITE(20)(X(I),I=1,NVAR)
      GO TO 301
340   REWIND 20
      DO 341 I=1,NVAR
      XMEAN(I)=SX(I)/FLOAT(NN(I))
341   SX(I)=0.0
C---------------K=NO. OF ABS., SEE ST. 303.
      DO 343 J=1,K
      READ(20)(X(I),I=1,NVAR)
      DO 342 I=1,NVAR
      IF(X(I).EQ.FMD(I))X(I)=XMEAN(I)
342   CONTINUE
      CALL TRANS(KTRAN,CONST,KVAR1,KVAR2,KVAR3,X,NTR,IRROR,KAR)
C---------------KAR RETURNED BY SUBR. TRANS AND KAR INDICATES
C--------------- ERRORS DETECTED BY SUBR. TRANS.
C---------------X(I) CONTAINS NEW VALUES GENERATED BY SUBR. TRANS.
      IF(KAR.EQ.1)GO TO 343
      WRITE(1)(X(I),I=1,NVTR)
      CALL XPRODH(X,SX,SXY,NVTR,70)
343   CONTINUE
      K=K-IRROR(1)-IRROR(2)
      GO TO 2020
C---------------HERE WE DELETE ENTIRE OBS. IF ONE VAR. HAS 
C--------------- MISSING DATA.
351   READ(INP,IFT,END=2020)(X(I),I=1,NVAR)
360   DO 354 I=1,NVAR
      IF(X(I).EQ.FMD(I))GO TO 351
354   CONTINUE
      CALL TRANS(KTRAN,CONST,KVAR1,KVAR2,KVAR3,X,NTR,IRROR,KAR)
      IF(KAR.EQ.1)GO TO 351
      K=K+1
      WRITE(1)(X(I),I=1,NVTR)
      CALL XPRODH(X,SX,SXY,NVTR,70)
      GO TO 351
C---------------HERE IF NO MISSING DATA
391   READ(INP,IFT,END=2020)(X(I),I=1,NVAR)
395   CALL TRANS(KTRAN,CONST,KVAR1,KVAR2,KVAR3,X,NTR,IRROR,KAR)
      IF(KAR.EQ.1)GO TO 391
      K=K+1
      WRITE(1)(X(I),I=1,NVTR)
      CALL XPRODH(X,SX,SXY,NVTR,70)
      GO TO 391
2020	IF(K.GT.1)GO TO 20210
      IF(IRROR(1).GT.0)WRITE(NOUTD,7654)IRROR(1)
      IF(IRROR(2).GT.0)WRITE(NOUTD,7655)IRROR(2)
	WRITE(NOUTD,7656)K
7656	FORMAT(' ONLY ',I2,' OBSERVATIONS REMAIN',/)
	IF(ICODE.EQ.-1)CALL EXIT
	GO TO 1111
C---------------VAR(I)=VARIANCE OF ITH VAR., SX(I)=MEAN OF ITH VAR.
C---------------SEE ST. 2023+3.
20210  DO 2022 I=1,NVTR
      SX(I)=SX(I)/FLOAT(K)
C---------------THIS IS COMP. FORMULA FOR VARIANCE TAKING NOTE
C--------------- OF THE FACT THAT ON THE LINE ABOVE
C--------------- SX(I)=SX(I)/FLOAT(K).  SEE STEPR FOLDER.
      VAR(I)=(SXY(I,I)-FLOAT(K)*SX(I)**2)/FLOAT(K-1)
      IF(ABS(VAR(I)).LT..0000001)GO TO 19
2022  CONTINUE
      NVAR=NVTR
      WRITE(IOUT,12)(IDEN(I),I=1,16)
12    FORMAT('1',16A5/)
      WRITE(IOUT,2026)K
2026  FORMAT(/' THE NUMBER OF OBSERVATIONS IS',I5//)
      IF(IRROR(1).GT.0)WRITE(IOUT,7654)IRROR(1)
7654  FORMAT(' NUMBER OF OBS. REJECTED BECAUSE OF ILLEGAL TRANS. LOG IS'
     1,I5/)
      IF(IRROR(2).GT.0)WRITE(IOUT,7655)IRROR(2)
7655  FORMAT(' NUMBER OF OBS. REJECTED BECAUSE OF ILLEGAL TRANS.',
     1' DIVISION IS',I5/)
      IF(IFOPT(1).NE.1)GO TO 50
      WRITE(IOUT,2023)
2023  FORMAT(4X,'VARIABLE',4X,'MEAN',8X,'VARIANCE',6X,'STD. DEV.'/)
      DO 2024 I=1,NVAR
      SD=SQRT(VAR(I))
      WRITE(IOUT,2025)I,SX(I),VAR(I),SD
2025  FORMAT(6X,I2,2X,3G)
2024  CONTINUE
50    IF(IFOPT(14).NE.1)GO TO 55
      WRITE(IOUT,2030)
2030  FORMAT('1RAW SUMS OF SQUARES AND CROSS PRODUCTS'/)
      KSIZE=5
C---------------RETURNED BY TTYPTY.  ICODE=-1 MEANS BATCH JOB
      IF(ICODE.EQ.-1)KSIZE=10
      KF=1
2032  KL=MIN0(KF+KSIZE-1,NVAR)
      WRITE(IOUT,2043)(I,I=KF,KL)
      DO 2034 I=KF,NVAR
      KLL=MIN0(KL,I)
C---------------RETURNED BY XPRODH
      WRITE(IOUT,2035)I,(SXY(I,J),J=KF,KLL)
2035  FORMAT(/I3,1X,12G12.5)
2034  CONTINUE
      IF(NVAR.EQ.KL)GO TO 55
      KF=KF+KSIZE
      WRITE(IOUT,2066)
      GO TO 2032
C---------------IFOPT(7)=1 MEANS FORCE EQ. THRU ZERO
55    IF(IFOPT(7).EQ.1)GO TO 3220
C---------------CALCULATE AND WRITE OUT COVARIANCE MATRIX.
      DO 2040 I=1,NVAR
      DO 2040 J=1,I
C---------------SEE COMMENTS WITH ST. 20210+2 ABOVE
2040  SXY(I,J)=(SXY(I,J)-FLOAT(K)*SX(I)*SX(J))/FLOAT(K-1)
C---------------IFOPT(3)=1 MEANS PRINT COV. MATRIX
      IF(IFOPT(3).NE.1)GO TO 60
      WRITE(IOUT,2041)
2041  FORMAT('1COVARIANCES'/)
      KSIZE=5
      IF(ICODE.EQ.-1)KSIZE=10
      KF=1
2042  KL=MIN0(KF+KSIZE-1,NVAR)
      WRITE(IOUT,2043)(I,I=KF,KL)
2043  FORMAT(4X,12(5X,I2,5X))
      DO 2044 I=KF,NVAR
      KLL=MIN0(KL,I)
      WRITE(IOUT,2045)I,(SXY(I,J),J=KF,KLL)
2045  FORMAT(/I3,1X,12F12.5)
2044  CONTINUE
      IF(NVAR.EQ.KL)GO TO 60
      KF=KF+KSIZE
      WRITE(IOUT,2066)
      GO TO 2042
C---------------HERE FROM ST. 55 I.E. USER WANTS FORCE THRU ZERO;
C--------------- MATRIX SUBMITTED TO SUBROUTINE REGGR IS SIMILAR TO
C--------------- CORRELATION MATRIX BUT MEANS ARE NOT SUBTRACTED
C--------------- I.E., WE USE RAW SUMS OF SQUARES AND CROSS PRODUCTS
C--------------- SEE STEPR FOLDER.
3220  REWIND 20
      DO 3221 I=1,NVAR
3221  X(I)=SQRT(SXY(I,I))
      DO 3222 I=1,NVAR
      DO 3222 J=1,I
      SXY(I,J)=SXY(I,J)/(X(I)*X(J))
3222  SXY(J,I)=SXY(I,J)
      WRITE(20)(SX(I),I=1,NVAR)
      WRITE(20)(X(I),I=1,NVAR)
      DO 3223 I=1,NVAR
3223  WRITE(20)(SXY(I,J),J=I,NVAR)
C---------------CALL REGGR
      GO TO 70
60    REWIND 20
      DO 2050 I=1,NVAR
      X(I)=SQRT(VAR(I)*FLOAT(K-1))
C---------------FORM CORRELATION MATRIX
      DO 2050 J=1,I
      SXY(I,J)=SXY(I,J)/SQRT(VAR(I)*VAR(J))
      IF(I.EQ.J)SXY(I,J)=1.0
2050  SXY(J,I)=SXY(I,J)
      WRITE(20)(SX(I),I=1,NVAR)
      WRITE(20)(X(I),I=1,NVAR)
      DO 2049  I=1,NVAR
 2049  WRITE(20)(SXY(I,J),J=I,NVAR)
C---------------IFOPT(2)=1 MEANS PRINT CORRELATION MATRIX
      IF(IFOPT(2).NE.1)GO TO 70
      WRITE(IOUT,2051)
2051  FORMAT('1CORRELATIONS'/)
      KSIZE=7
      IF(ICODE.EQ.-1)KSIZE=14
      KF=1
2052  KL=MIN0(KF+KSIZE-1,NVAR)
      WRITE(IOUT,2053)(I,I=KF,KL)
2053  FORMAT(4X,14(4X,I2,3X))
      DO 2054 I=KF,NVAR
      KLL=MIN0(KL,I)
      WRITE(IOUT,2055)I,(SXY(I,J),J=KF,KLL)
2055  FORMAT(/I3,1X,14F9.5)
2054  CONTINUE
      IF(NVAR.EQ.KL)GO TO 2058
      KF=KF+KSIZE
      WRITE(IOUT,2066)
      GO TO 2052
2058  WRITE(IOUT,2059)
2059  FORMAT('1T-VALUES'/)
      KF=1
      N=NVAR
2060  KL=MIN0(KF+KSIZE-1,N)
      WRITE(IOUT,2061)(I,I=KF,KL)
2061  FORMAT(4X,14(3X,I2,4X))
      DO 2062 I=KF,N
      KLL=MIN0(KL,I)
      DO 2063 J=KF,KLL
      IF(I.EQ.J)GO TO 2064
      SXY2=1.0-SXY(I,J)**2
      IF(SXY2.LE.0.0)GO TO 2064
      X(J)=SXY(I,J)*SQRT(FLOAT(K-2)/SXY2)
      GO TO 2063
2064  X(J)=0.0
2063  CONTINUE
      WRITE(IOUT,2065)I,(X(J),J=KF,KLL)
2065  FORMAT(/I3,1X,14F9.4)
2062  CONTINUE
      IF(NVAR.EQ.KL)GO TO 2068
      KF=KF+KSIZE
      WRITE(IOUT,2066)
2066  FORMAT(///)
      GO TO 2060
2068  WRITE(IOUT,2069)
2069  FORMAT('1Z-VALUES'/)
      KF=1
2070  KL=MIN0(KF+KSIZE-1,N)
      WRITE(IOUT,2061)(I,I=KF,KL)
      DO 2072 I=KF,N
      KLL=MIN0(KL,I)
      DO 2073 J=KF,KLL
      IF(I.EQ.J)GO TO 2074
      IF(1.0-SXY(I,J))2071,2074,2071
2071  SXY2=(1.0+SXY(I,J))/(1.0-SXY(I,J))
      IF(SXY2.LE.0.0)GO TO 2074
      X(J)=SQRT(FLOAT(K-3))*.5*ALOG(SXY2)
      GO TO 2073
2074  X(J)=0.0
2073  CONTINUE
      WRITE(IOUT,2065)I,(X(J),J=KF,KLL)
2072  CONTINUE
      IF(NVAR.EQ.KL)GO TO 70
      KF=KF+KSIZE
      WRITE(IOUT,2066)
      GO TO 2070
C---------------K=NO. OF OBS. (SEE ST. 303)
70    CALL REGGR(SXY,K,VAR,SX,NVAR,IFOPT)
      WRITE(NOUTD,71)
71    FORMAT(' DO YOU WISH TO REANALYZE THE SAME DATA?'/)
      READ(IND,18)KWHER
18    FORMAT(A3)
      IF(KWHER.EQ.'YES')GO TO 70
C---------------ST. 1111 IS ON FIRST PAGE.  NEXT PROMPTING
C--------------- IS INPUT?
      GO TO 1111
19    WRITE(NOUTD,21)I
21    FORMAT(' VARIABLE ',I3,' HAS ZERO VARIANCE'/)
      IF(ICODE.EQ.-1)CALL EXIT
      GO TO 1111
      END
C---------------KTRAN, CONST, KVAR1, KVAR2, KVAR3, NTR ARE
C--------------- INPUT.  X, KAR ARE OUTPUT.
      SUBROUTINE TRANS(KTRAN,CONST,KVAR1,KVAR2,KVAR3,X,NTR,IRROR,KAR)
      DIMENSION X(70),IRROR(2)
      DIMENSION KTRAN(50),KVAR1(50),KVAR2(50),KVAR3(50),CONST(50)
      KAR=0
      IF(NTR.EQ.0)RETURN
      DO 10 I=1,NTR
      GO TO (20,21,22,23,24,25,26,27,28,29,30,31,32,33,34),KTRAN(I)
C---------------SEE OPTION 2 PAGE 5 OF WRITE UP
20    X(KVAR1(I))=-X(KVAR2(I))
      GO TO 10
21    X(KVAR1(I))=X(KVAR2(I))
      GO TO 10
22    X(KVAR1(I))=X(KVAR2(I))**CONST(I)
      GO TO 10
23    X(KVAR1(I))=X(KVAR2(I))+X(KVAR3(I))
      GO TO 10
24    X(KVAR1(I))=X(KVAR2(I))*X(KVAR3(I))
      GO TO 10
25    IF(X(KVAR2(I)).GT.0.0)GO TO 40
      IRROR(1)=IRROR(1)+1
      KAR=1
      RETURN
40    X(KVAR1(I))=ALOG(X(KVAR2(I)))
      GO TO 10
26    IF(X(KVAR2(I)).GT.0.0)GO TO 41
      IRROR(1)=IRROR(1)+1
      KAR=1
      RETURN
41    X(KVAR1(I))=ALOG10(X(KVAR2(I)))
      GO TO 10
27    X(KVAR1(I))=EXP(X(KVAR2(I)))
      GO TO 10
28    X(KVAR1(I))=10.0**X(KVAR2(I))
      GO TO 10
29    X(KVAR1(I))=X(KVAR2(I))+CONST(I)
      GO TO 10
30    X(KVAR1(I))=X(KVAR2(I))*CONST(I)
      GO TO 10
31    X(KVAR1(I))=SIN(X(KVAR2(I)))
      GO TO 10
32    X(KVAR1(I))=COS(X(KVAR2(I)))
      GO TO 10
33    IF(X(KVAR3(I)).NE.0.0)GO TO 42
      IRROR(2)=IRROR(2)+1
      KAR=1
      RETURN
42    X(KVAR1(I))=X(KVAR2(I))/X(KVAR3(I))
      GO TO 10
34    X(KVAR1(I))=X(KVAR2(I))-X(KVAR3(I))
10    CONTINUE
      RETURN
      END
C---------------KTRAN, CONST, KVAR1, KVAR2, KVAR3, NTR ARE OUTPUT;
C--------------- IT IS INPUT.  NTR=NO. OF TRANSF.,
      SUBROUTINE LANG2(KTRAN,CONST,KVAR1,KVAR2,KVAR3,IT,NTR)
      DIMENSION KTRAN(50),CONST(50),KVAR1(50),KVAR2(50)
      DIMENSION KVAR3(50),IT(70)
      NOUTD=-1
      IND=-4
      NTR=0
      WRITE(NOUTD,2001)
2001  FORMAT(' ENTER TRANSFORMATIONS'/)
2000  READ(IND,100)(IT(I),I=1,70)
100   FORMAT(70A1)
      IER=0
      K=1
      DO 1 I=1,70
      IF(IT(I).EQ.' ')GO TO 1
      IT(K)=IT(I)
      K=K+1
1     CONTINUE
      NK=K-1
      DO 2 I=K,70
2     IT(I)=' '
      IF(IT(1).EQ.'E'.AND.IT(2).EQ.'N'.AND.IT(3).EQ.'D')GO TO 3000
      NTR=NTR+1
      IF(IT(1).EQ.'(')GO TO 500
3     WRITE(NOUTD,1100)
1100  FORMAT(' ERROR IN TRANSFORMATION'/)
      CALL TTYPTY(ICC)
      IF(ICC.EQ.-1)CALL EXIT
      NTR=NTR-1
      GO TO 2000
500   I=2
      CALL TNUM(I,X,IT,IER)
      IF(IER.NE.0) GO TO 3
      KVAR1(NTR)=X+.01
      I=I+1
      IF(IT(I).EQ.'=')GO TO 501
      GO TO 3
501   I=I+1
      IF(IT(I).EQ.'(')GO TO 300
      IF(IT(I).EQ.'E')GO TO 400
      IF(IT(I).EQ.'L')GO TO 600
      IF(IT(I).EQ.'1')GO TO 700
      IF(IT(I).EQ.'S')GO TO 800
      IF(IT(I).EQ.'C')GO TO 900
      IF(IT(I).EQ.'-')GO TO 1000
      GO TO 3
300   I=I+1
      CALL TNUM(I,X,IT,IER)
      IF(IER.NE.0)GO TO 3
      KVAR2(NTR)=X+.01
      I=I+1
      IF(IT(I).EQ.'*')GO TO 310
      IF(IT(I).EQ.'+')GO TO 320
      IF(IT(I).EQ.'-')GO TO 330
      IF(IT(I).EQ.'/')GO TO 340
      IF(IT(I).EQ.' ')GO TO 350
      GO TO 3
310   I=I+1
      IF(IT(I).EQ.'*')GO TO 304
      IF(IT(I).EQ.'(')GO TO 306
      IF(IT(I).EQ.'+'.OR.IT(I).EQ.'-'.OR.IT(I).EQ.'.')GO TO 308
      IF(IT(I).LE.'9'.AND.IT(I).GE.'0')GO TO 308
      GO TO 3
304   I=I+1
      CALL TNUM(I,X,IT,IER)
      IF(IER.NE.0)GO TO 3
      CONST(NTR)=X
      KTRAN(NTR)=3
      GO TO 2000
306   I=I+1
      CALL TNUM(I,X,IT,IER)
      IF(IER.NE.0)GO TO 3
      KVAR3(NTR)=X+.01
      KTRAN(NTR)=5
      GO TO 2000
308   CALL TNUM(I,X,IT,IER)
      IF(IER.NE.0)GO TO 3
      CONST(NTR)=X
      KTRAN(NTR)=11
      GO TO 2000
320   I=I+1
      IF(IT(I).EQ.'(')GO TO 322
      IF(IT(I).EQ.'.')GO TO 326
      IF(IT(I).GE.'0'.AND.IT(I).LE.'9')GO TO 326
      GO TO 3
322   I=I+1
      CALL TNUM(I,X,IT,IER)
      IF(IER.NE.0)GO TO 3
      KVAR3(NTR)=X+.01
      KTRAN(NTR)=4
      GO TO 2000
326   CALL TNUM(I,X,IT,IER)
      IF(IER.NE.0)GO TO 3
      CONST(NTR)=X
      KTRAN(NTR)=10
      GO TO 2000
330   I=I+1
      IF(IT(I).EQ.'(')GO TO 332
      IF(IT(I).EQ.'.')GO TO 336
      IF(IT(I).GE.'0'.AND.IT(I).LE.'9')GO TO 336
      GO TO 3
332   I=I+1
      CALL TNUM(I,X,IT,IER)
      IF(IER.NE.0)GO TO 3
      KVAR3(NTR)=X+.01
      KTRAN(NTR)=15
      GO TO 2000
336   I=I-1
      CALL TNUM(I,X,IT,IER)
      IF(IER.NE.0)GO TO 3
      CONST(NTR)=X
      KTRAN(NTR)=10
      GO TO 2000
340   I=I+1
      IF(IT(I).EQ.'(')GO TO 341
      GO TO 3
341   I=I+1
      CALL TNUM(I,X,IT,IER)
      IF(IER.NE.0)GO TO 3
      KVAR3(NTR)=X+.01
      KTRAN(NTR)=14
      GO TO 2000
350   KTRAN(NTR)=2
      GO TO 2000
400   I=I+4
      IF(IT(I-3).EQ.'*'.AND.IT(I-2).EQ.'*'.AND.IT(I-1).EQ.'(')GO TO 450
      GO TO 3
450   CALL TNUM(I,X,IT,IER)
      IF(IER.NE.0)GO TO 3
      KVAR2(NTR)=X+.01
      KTRAN(NTR)=8
      GO TO 2000
600   I=I+1
      IF(IT(I).EQ.'N')GO TO 620
      IF(IT(I).EQ.'O')GO TO 640
      GO TO 3
620   I=I+1
      IF(IT(I).EQ. '(')GO TO 625
      GO TO 3
625   I=I+1
      CALL TNUM(I,X,IT,IER)
      KVAR2(NTR)=X+.01
      KTRAN(NTR)=6
      GO TO 2000
640   I=I+2
      IF(IT(I-1).EQ.'G'.AND.IT(I).EQ.'(')GO TO 645
      GO TO 3
645   I=I+1
      CALL TNUM(I,X,IT,IER)
      IF(IER.NE.0)GO TO 3
      KVAR2(NTR)=X+.01
      KTRAN(NTR)=7
      GO TO 2000
700   I=I+4
      IF(IT(I-3).EQ.'0'.AND.IT(I-2).EQ.'*'.AND.IT(I-1).EQ.'*'.AND.IT(I).
     1   EQ.'(')GO TO 750
      GO TO 3
750   I=I+1
      CALL TNUM(I,X,IT,IER)
      IF(IER.NE.0)GO TO 3
      KVAR2(NTR)=X+.01
      KTRAN(NTR)=9
      GO TO 2000
800   I=I+3
      IF(IT(I-2).EQ.'I'.AND.IT(I-1).EQ.'N'.AND.IT(I).EQ.'(')GO TO 850
      GO TO 3
850   I=I+1
      CALL TNUM(I,X,IT,IER)
      IF(IER.NE.0)GO TO 3
      KVAR2(NTR)=X+.01
      KTRAN(NTR)=12
      GO TO 2000
900   I=I+3
      IF(IT(I-2).EQ.'O'.AND.IT(I-1).EQ.'S'.AND.IT(I).EQ.'(')GO TO 950
      GO TO 3
950   I=I+1
      CALL TNUM(I,X,IT,IER)
      IF(IER.NE.0)GO TO 3
      KVAR2(NTR)=X+.01
      KTRAN(NTR)=13
      GO TO 2000
1000  I=I+1
      IF(IT(I).EQ.'(')GO TO 1050
      GO TO 3
1050  I=I+1
      CALL TNUM(I,X,IT,IER)
      IF(IER.NE.0)GO TO 3
      KVAR2(NTR)=X+.01
      KTRAN(NTR)=1
      GO TO 2000
3000  RETURN
      END
C---------------IT INPUT, IFOPT OUTPUT.
      SUBROUTINE LANG1(IT,IFOPT)
      DIMENSION ICOM(15,3),IT(70),IFOPT(15)
      DATA (ICOM(I,1),I=1,10)/'M','C','C','R','A','D','Z','T','M','O'/
      DATA (ICOM(I,2),I=1,10)/'V','O','O','E','N','U','E','R','I','B'/
      DATA (ICOM(I,3),I=1,10)/'S','R','V','S','A','R','R','A','S','J'/
      DATA (ICOM(I,1),I=11,15)/'F','F','E','X','T'/
      DATA (ICOM(I,2),I=11,15)/'V','O','L','P','O'/
      DATA (ICOM(I,3),I=11,15)/'L','R','M','R','L'/
      K=1
C**** WMU AM: 1.3.2-2, MTO, 1-FEB-78
      DO 1I=1,69
C**** END = LANG1, STMT #1-4
      IF(IT(I).EQ.' ')GO TO 1
      IT(K)=IT(I)
      K=K+1
1     CONTINUE
      NK=K
      DO 2 I=K,70
2     IT(I)=' '
C**** WMU AM: 1.3.2-2, MTO, 1-FEB-78
      DO 3 I=4,NK
      IF(IT(I).NE.','.AND.IT(I).NE.' ')GO TO 3
      DO 3 J=1,15
      IF(IT(I-3).EQ.ICOM(J,1).AND.IT(I-2).EQ.ICOM(J,2).AND.IT(I-1).EQ.
     1ICOM(J,3))IFOPT(J)=1
C**** END = LANG1, STMT #2+5
3     CONTINUE
      RETURN
      END
C---------------I, X, IER ARE OUTPUT;  IT IS INPUT
      SUBROUTINE TNUM(I,X,IT,IER)
      DIMENSION IT(70),IN(15),ITP(3)
      DO 10 J=I,70
      IF(IT(J).EQ.'.')GO TO 10
      IF(IT(J).GE.'0'.AND.IT(J).LE.'9')GO TO 10
      IF(IT(J).EQ.' '.OR.IT(J).EQ.')')GO TO 20
      IF(J.NE.I)IER=1
      IF(IT(J).EQ.'+')GO TO 10
      IF(IT(J).EQ.'-')GO TO 10
      IER=1
10    CONTINUE
20    NPT=0
      KF=J-1
      DO 30 J=I,KF
      IF(IT(J).EQ.'.')NPT=NPT+1
30    CONTINUE
      IF(NPT.GT.1)IER=1
      IF(IER.NE.0)RETURN
      NCH=15-(KF-I)
      DO 40 J=1,15
40    IN(J)=' '
      K=I
      DO 50 J=NCH,15
      IN(J)=IT(K)
50    K=K+1
      ENCODE(15,60,ITP)IN
60    FORMAT(15A1)
      DECODE(15,70,ITP)X
70    FORMAT(F15.0)
      I=KF+1
      RETURN
      END
C  PDP-10  VERSION OF ESSO STEPWISE REGRESSION ANALYSIS - PART 2
C     BLOCK E - READ MATRIX, ELIMINATE VARIABLES
C---------------ALL ARGS. ARE INPUT.  ARRAY=EITHER CORRELATION
C--------------- MATRIX OR A MATRIX DERIVED FROM CORRELATION COEFFS.
C--------------- WHERE MEANS ARE NOT SUBTRACTED.  SEE ST. 3220 AND
C--------------- ST. 2049+1 IN MAIN PROG.;  NDATA=NO. OF OBS.;
C--------------- SIGMA=VARIANCES;  AVG=MEANS, NTERS=NO. OF VAR.,
C--------------- IFOPT=USER SELECTED OPTIONS.
      SUBROUTINE REGGR(ARRAY,NDATA,SIGMA,AVG,NTERS,IFOPT)
C**L.M.#1.3.2-1 RRB 26SEP77
C**ADDED ARRAY "IFORCE" TO REPLACE "INDEX" BEFORE DO LOOP 104
C**SO LIST OF FORCED VARIABLES IS NOT OVER-WRITTEN.
      DIMENSION COEN(70),SIGCO(70),INDEX(70),SIGMA(70),STUDT(70),
     1DATA(70),AVG(70),BETA(70),ARRAY(70,70),INDX2(70),NAME(70),
     2LAME(70),KAME(70),KFR(70),IFOPT(15),IFORCE(70)
      IND=-4
      NOUTD=-1
      IOUT=3
      TOL=.0001
      REWIND 20
      REWIND 1
      INP=5
      COE=0.0
      EFIN=0.0
      EFOUT=0.0
      NTERM=NTERS
      NOELM=0
      NOFRC=0
      OBSER=FLOAT(NDATA)
      WRITE(NOUTD,2000)
2000  FORMAT(///' WHICH IS THE DEPENDENT VARIABLE?'/)
      READ(IND,2001)NDVAR
2001  FORMAT(I)
      IF(IFOPT(15).NE.1)GO TO 2019
      WRITE(NOUTD,2022)
2022  FORMAT(' ENTER TOLERANCE'/)
      READ(IND,2003)TOL
2019  IF(IFOPT(11).NE.1)GO TO 3862
      WRITE(NOUTD,2002)
2002  FORMAT(' ','ENTER F-VALUE FOR ENTERING A VARIABLE.'/)
      READ(IND,2003)EFIN
2003  FORMAT(F)
      WRITE(NOUTD,2004)
2004  FORMAT(' ','ENTER F-VALUE FOR OMITTING A VARIABLE.'/)
      READ(IND,2003)EFOUT
3862  IF(IFOPT(13).NE.1)GO TO 2006
      WRITE(NOUTD,2005)
2005  FORMAT(' ','HOW MANY VARIABLES WOULD YOU LIKE TO ELIMINATE?'/)
      READ(IND,2001)NOELM
      IF(NOELM)2006,2006,2007
2007  WRITE(NOUTD,2008)
2008  FORMAT(' ','WHICH ARE THEY?(MAX: 20 PER LINE)'/)
      READ(IND,2009)(INDX2(I),I=1,NOELM)
2009  FORMAT(20I)
2006  NOVAR=NDVAR
      IF(IFOPT(12).NE.1)GO TO 715
      WRITE(NOUTD,2010)
2010  FORMAT(' ','ENTER NUMBER OF VARIABLES TO BE FORCED INTO THE ',
     1'REGRESSION.'/)
      READ(IND,2001)NOFRC
      IF(NOFRC)2011,715,2011
2011  WRITE(NOUTD,2008)
C**L.M.#1.3.2-1 RRB 26SEP77 - CHANGED "INDEX" TO "IFORCE"
      READ(IND,2009)(IFORCE(I),I=1,NOFRC)
715   LDVAR = NDVAR
      ISTRE = 0
      KDV1=0
      IDW=1
      SDEV2=0.0
      IDIF=0
      KTR = 0
      SDIF2 = 0.0
      IF(NOELM)28,28,302
28    NOELM=0
      GO TO 70
302    IF(NOELM)994,994,304
994   NPCRD = 2
      NCOL=6
995   WRITE(IOUT,10)NPCRD,NCOL
10    FORMAT(24H ERROR ON PARAMETERCARD,I3,7X,4HCOL.,I2)
      GO TO 1004
304   IF(NOELM - 1)994,303,72
72    L15=NOELM - 1
      DO 451 K8 = 1,L15
      K9=K8+1
C---------------SORT VARS. TO BE ELEMINATED.  INDX2(I) ARE
C--------------- INDICES OF VARS. TO BE ELIM.
      DO 451 K10=K9,NOELM
      IF(INDX2(K8) - INDX2(K10))451,451,450
450   LN=INDX2(K8)
      INDX2(K8)=INDX2(K10)
      INDX2(K10)=LN
451   CONTINUE
303   L=NOELM+1
      INDX2(L)=0
C---------------OBSER=NO. OF OBS. SEE ST. 2000-2.
70    IF(OBSER)71,73,73
71    IDW=0
      IFDEV=0
      OBSER= -OBSER
C---------------NTERM=NO. OF VARS. SEE ST. 2000-5.  COEN(I)=
C--------------- SQRT. OF ROW SUMS OF SQUARES.  (SEE ST. 3222+2
C--------------- AND 3221 IN MAIN PROG.)  IF USER CHOOSES FORCE
C--------------- THRU ZERO OPTION;  OTHERWISE COEN(I)=SQRT. OF
C--------------- QUANTITY VARIANCE TIMES NO. OF OBS. MINUS 1 (SEE
C--------------- ST. 60+2 IN MAIN PROG.) IF USER DOES NOT CHOOSE FORCE
C--------------- THRU ZERO OPTION.  (ALSO SEE ST. 2050+2 IN MAIN PROG.)
73    READ(20)(AVG(I),I=1,NTERM)
      READ(20)(COEN(I),I=1,NTERM)
        LTERM=NTERM
      DO 1050 K=1,NTERM
1050  NAME(K)=K
518   WRITE(IOUT,502)NDVAR
502   FORMAT(13H/VARIABLE NO.,I4,14H  IS DEPENDENT///)
      IF(NOELM) 700,703,700
700   WRITE(IOUT,701)(INDX2(I),I=1,NOELM)
701   FORMAT(1H0,18X,'ELIMINATE ', 30I3/(30X,30I3))
      GO TO 1703
703   WRITE(IOUT,1701)NTERM
1701  FORMAT(1H ,I4,29H VARIABLES    NONE ELIMINATED)
1703  J=1
      L = 1
      N=INDX2(1)
C---------------COMPRESS VARS. NOT TO BE ELIMINATED.  KFR(I) WILL
C---------------GIVE INDICES FOR NON-ELIMINATED VARS. IN SUCH A WAY THAT
C--------------- THE KTH POSITION GIVES ROW AND COL. NO. OF MATRIX ARRAY
C--------------- FOR THE KTH VARIABLE.
      DO 87 I=1,NTERM
      IF(NOELM)994,85,86
86    IF(I-N)88,89,88
88    AVG(J)=AVG(I)
85    SIGMA(J)=COEN(I)
      KAME(J) = NAME(I)
      K6 = NAME(I)
      KFR(K6) = J
      IF(I - NDVAR)851,850,851
850   KDV = J
 851  J=J+1
      GO TO 87
89    L=L+1
      N=INDX2(L)
87    CONTINUE
852   J1 = J - 1
869   IF(KDV - J1)870,871,3380
C---------------PROG. SHOULD NEVER GET TO 3380 SINCE BY ST. 850
C--------------- KDV IS ONE OF THE NON-ELIMINATED VARS.
3380  WRITE(IOUT,3381)
3381  FORMAT(66H NO. OF DEPENDENT VARIABLE IS GREATER THAN THE NUMBER OF
     1 VARIABLES)
C---------------GO TO RETURN
      GO TO 1004
870   R1=AVG(KDV)
      R2=SIGMA(KDV)
      N1=KAME(KDV)
C---------------PUT VALUES FOR DEP. VAR. IN LAST POSITION OF
C---------------NON-ELIMINATED VARS.
      DO 872 I=KDV,J1-1
      AVG(I)=AVG(I+1)
      SIGMA(I)=SIGMA(I+1)
      KAME(I)=KAME(I+1)
872   CONTINUE
      AVG(J1)=R1
      SIGMA(J1)=R2
      KAME(J1)=N1
      DO 875 I=1,NTERM
      IF(KFR(I).EQ.KDV)GO TO 873
875   CONTINUE
873   DO 874 J=I+1,NTERM
874   KFR(J)=KFR(J)-1
871   M = 1
      LSTRT = 1
      DO 331 I=1,NTERM
      READ(20)(COEN(J),J=I,NTERM)
      IF(NOELM)994,333,332
332   L = LSTRT
      N=INDX2(L)
      IF(I -  N)339,330,339
330   LSTRT = LSTRT+1
      GO TO 331
339   K=M
      L=LSTRT
      N=INDX2(L)
      DO 336 J=I,NTERM
      IF(J-N) 335,337,335                                  
 335   ARRAY(M,K)=COEN(J)
3352  K=K+1
      GO TO 336
337   L=L+1
      N=INDX2(L)
336   CONTINUE
      M=M+1
      GO TO 331
333   DO 338  K=I,NTERM
338   ARRAY(I,K)=COEN(K)
331   CONTINUE
      KDV1=KDV
      DEFR=OBSER-1.0
C---------------NOVAR=NO. OF TERMS NOT TO BE ELIMINATED
      NOVAR = NTERM - NOELM
      NOVMI=NOVAR-1
      DO 300 I=2,NOVAR
      IM1=I-1
      DO 300 J=1,IM1
300   ARRAY(I,J)=ARRAY(J,I)
      IF(NOELM) 994,3312,3310
C---------------KDV1=INDEX OF DEP. VAR. SEE ST. 331+1.
3310  IF(KDV1 - NOVAR)3311,3049,3380
3311  NDVAR = KDV1
      NTERM = NOVAR
      GO TO 3382
3312  IF(NDVAR - NTERM)3382,3049,3380
3382  DO 3000 I=1,NTERM
      T2 = ARRAY(I,NDVAR)
      DO 3001 J=NDVAR,NTERM-1
3001  ARRAY(I,J)=ARRAY(I,J+1)
3000  ARRAY(I,NTERM)=T2
      DO 3002 I=1,NTERM
      T2=ARRAY(NDVAR,I)
      DO 3003 J=NDVAR,NTERM-1
3003  ARRAY(J,I)=ARRAY(J+1,I)
3002  ARRAY(NTERM,I)=T2
C---------------NOFRC=NO. OF VARS. TO BE FORCED.
3049  IF(NOFRC)305,350,305
C**L.M.#1.3.2-1 RRB 26SEP77 - CHANGED "INDEX" TO "IFORCE"
305   WRITE(IOUT,702)(IFORCE(I),I=1,NOFRC)
702   FORMAT(23X, 'FORCE ',  30I3/(30X,30I3))
C
C     BLOCK F - SOLVE SYSTEM OF EQUATIONS FOR COEFFICIENTS
313   DO 306 L=1,NOFRC
C**L.M.#1.3.2-1 RRB 26SEP77 - CHANGED "INDEX" TO "IFORCE"
      K=IFORCE(L)
      K= KFR(K)
      DIAG=ARRAY(K,K)
      DO 307 I=1,NOVAR
      IF(I-K)308,307,308
308   FACTR=ARRAY(I,K)/DIAG
      DO 309 J=1,NOVAR
      IF(J-K)310,309,310
310   ARRAY(I,J)=ARRAY(I,J)-(FACTR*ARRAY(K,J))
309   CONTINUE
307   CONTINUE
      DO 311 I=1,NOVAR
311   ARRAY(I,K)=-ARRAY(I,K)/DIAG
      DO 312 J=1,NOVAR
312   ARRAY(K,J)=ARRAY(K,J)/DIAG
306   ARRAY(K,K)=1.0/DIAG
      FACTR = NOFRC
      DEFR = DEFR-FACTR
      FLEVL=(ARRAY(K,NOVAR)**2/DIAG)*DEFR/ARRAY(NOVAR,NOVAR)
      NOENT = K
      NOSTP=NOFRC-1
6004  GO TO 1006
350   NOSTP=-1
6005  GO TO 1006
C
C     BLOCK G-SELECT VARIABLES BY ANAL OF VARIANCE AND SOLVE EQUATIONS
 1006 NOSTP=NOSTP+1
      IF(ARRAY(NOVAR,NOVAR))147,14,101
147   WRITE(IOUT,1014)
1014  FORMAT(34H LAST DIAGONAL ELEMENT IS NEGATIVE)
      GO TO 1004
14    WRITE(IOUT,15)
15    FORMAT(30H LAST DIAGONAL ELEMENT IS ZERO)
      GO TO 1004
101   IF(ARRAY(NOVAR,NOVAR).LT.0.0)WRITE(NOUTD,2222)
2222  FORMAT(' 1R')
      SIGY=SIGMA(NOVAR)*SQRT( ARRAY(NOVAR,NOVAR)/DEFR )
      DEFR=DEFR-1.0
      IF(DEFR) 16,103,103
16    WRITE(IOUT,17)DEFR
17    FORMAT(22H DEGREES OF FREEDOM IS,F8.2)
      GO TO 1004
 103  VMIN=-10000.0
      VMAX=0.0
      NOIN=0
      DO 104 I=1,NOVMI
	DIAG=ARRAY(I,I)
C**L.M.#1.3.2-1 RRB 26SEP77
C**EXCLUDE FORCED VARIABLES FROM ENTER/REMOVE TESTING
	IF(DIAG.EQ.0)GO TO 104
	IF(NOFRC.EQ.0)GO TO 105
	DO 1991 JJ=1,NOFRC
1991	IF(KAME(I).EQ.IFORCE(JJ))GO TO 1992
	GO TO 105
1992	SIG=ABS(SIGMA(I))
	IF(SIG.EQ.0)GO TO 18
	GO TO 134
C**END OF L.M.#1.3.2-1
105   IF(DIAG)107,104,108
107   IF(DIAG+TOL)109,104,104
109   WRITE(IOUT,12)KAME(I)
12    FORMAT( 43H NEGATIVE DIAGONAL ELEMENT BELOW TOLERANCE,,8HVARIABLE,
     1I4)
      GO TO 1004
108   IF(DIAG-TOL)104,110,110
110   VAR=ARRAY(I,NOVAR)*ARRAY(NOVAR,I)/DIAG
      SIG=SIGMA(I)
      IF(VAR)106,104,112
106   IF(SIG)129,18,122
18    WRITE(IOUT,19)SIG,KAME(I)
19    FORMAT(22H STANDARD DEVIATION IS,F14.7,13HFOR VARIABLE ,I3)
      GO TO 1004
129   SIG=-SIG
      GO TO 134
122   IF(VAR-VMIN)134,134,113
113   VMIN=VAR
C---------------INDEX OF VAR. CORRESPONDING TO MIN. VAR.
      NOMIN=I
C---------------NOIN IS NO. OF EQS. IN REG. EQ. (SEE ST. 133+2)
134   NOIN=NOIN+1
      INDEX(NOIN)=I
      COEN(NOIN)=ARRAY(I,NOVAR)*SIGMA(NOVAR)/SIG
      IF(DIAG.LT.0.0)WRITE(NOUTD,2223)
2223  FORMAT(' 2R')
      SIGCO(NOIN)=(SIGY/SIG)*SQRT(DIAG)
1341  LAME(NOIN)=KAME(I)
      IF(KAME(I)-LDVAR)104,1342,104
1342  LAME(NOIN)=J1
      GO TO 104
112   IF(SIG) 18,18,146
146   IF(VAR-VMAX)104,104,114
114   VMAX=VAR
C---------------NOMAX=INDEX OF VAR. CORRESPONDING TO MAX. VAR.
       NOMAX=I
104   CONTINUE
      IF(NOIN) 99,115,116
115   KTR=KTR+1
      WRITE(IOUT,400)SIGY
400   FORMAT(25H0  STANDARD ERROR OF Y = G)
      IF(KTR - 1)123,123,455
455   WRITE(IOUT,456)
456   FORMAT(38H ERROR, NO VARIABLE CAN ENTER EQUATION)
      RETURN
116   IF(IFOPT(7).NE.1)GO TO 117
118   CNST=0.0
      SALPHA=0.0
      TALPHA=0.0
C**** WMU AM: 1.3.2-2, MTO, 30-JAN-78
      IDFALP=NDATA-NOIN
      DO 6031 I=1,NOIN
6031  SIGCO(I)=SIGCO(I)*SQRT(FLOAT(IDFALP-1)/IDFALP)
      GO TO 120
117   CNST=AVG(NOVAR)
      SALPHA=1.0/FLOAT(NDATA)
      DO 119 J=1,NOIN
      L=INDEX(J)
      DO 8022 I=1,NOIN
      LL=INDEX(I)
      SALPHA=SALPHA+AVG(L)*AVG(LL)*ARRAY(L,LL)/(SIGMA(L)*SIGMA(LL))
8022  CONTINUE
119   CNST=CNST-(COEN(J)*AVG(L))
      SALPHA=SQRT(SALPHA)*SIGY
      TALPHA=CNST/SALPHA
      IDFALP=NDATA-NOIN-1
120   FLEV=FLEVL
      IF (IFOPT(7).EQ.1) FLEV=IDFALP*FLEV/(IDFALP-1)
      NOEN=NOENT
      K1=K
      IF(VMIN+10000.0)121,123,121
121   FLEVL=-VMIN*(DEFR+1.0)/ARRAY(NOVAR,NOVAR)
      IF(EFOUT-FLEVL)123,123,124
124   K=NOMIN
      DEFR=DEFR+2.0
      NOENT=0
      GO TO 126
123   FLEVL=VMAX*DEFR/(ARRAY(NOVAR,NOVAR)-VMAX)
      IF(EFIN - FLEVL)125,1007,1007
125   K=NOMAX
      NOENT=K
126   CONTINUE
      COEFDT=1.0-ARRAY(NOVAR,NOVAR)
      COD=COEFDT-COE
      COE=COEFDT
      IF(COEFDT.LT.0.0)WRITE(NOUTD,2224)
2224  FORMAT(' 3R')
      RRR=SQRT(COEFDT)
127   IF(NOIN)99,128,130
130   IF(NOEN)131,131,132
132   WRITE(IOUT,401)NOSTP,KAME(K1)
401   FORMAT(///9H STEP NO. I3/20H   VARIABLE ENTERING I6)
      ISTRE = K1
      GO TO 133
131   IF(K1 - ISTRE)8037,8038,8037
8038  IF(K1)8033,8037,8033
8037  WRITE(IOUT,402)NOSTP,KAME(K1)
402   FORMAT(///9H STEP NO. I3/19H   VARIABLE REMOVED I6)
133	FPRB=FISHER(1,IDFALP,FLEV)
C****** END = REGGR, STMT #133
C**L.M.#1.3.2-1 RRB 26SEP77 - ADD TPROB TO STEPS
	TPRB=FISHER(1,IDFALP,TALPHA*TALPHA)
      WRITE(IOUT,403)FLEV,FPRB,SIGY,COEFDT,RRR,COD,IDFALP,CNST,SALPHA,
     1TALPHA,TPRB,(LAME(J),COEN(J),SIGCO(J),J=1,NOIN)
403   FORMAT(3X,'F LEVEL',2X,F12.4,5X,'F-PROB =',2X,F8.5,
     1/3X,'STANDARD ERROR OF ESTIMATE = ',
     1F12.4/'   COEFFICIENT OF DETERMINATION = ',F8.5/
     1'   COEFFICIENT OF MULTIPLE REGRESSION = ',F8.5/
     1 '   INCREASE IN COEFFICIENT OF DETERMINATION = ',F8.5/
     13X,'DEGREES OF FREEDOM =',I6,/,
     15X,'CONSTANT',5X,'STD. ERR.',13X,'T',9X,'T-PROB'/1X,3G,F8.5/13X,
     1'VARIABLE   COEFFICIENT   STD ERRO',
     2'R OF COEF'  /(16H              X=I3,F15.5,F15.5))
C**END OF L.M.#1.3.2-1
128   DIAG=ARRAY(K,K)
      DO 140 I=1,NOVAR
      IF(I-K)141,140,141
141   FACTR=  ARRAY(I,K)/DIAG
      DO 143 J=1,NOVAR
      IF(J-K)142,143,142
142   ARRAY(I,J)=ARRAY(I,J)-(FACTR*ARRAY(K,J))
143   CONTINUE
140   CONTINUE
      DO 144 I=1,NOVAR
144   ARRAY(I,K)=-ARRAY(I,K)/DIAG
      DO 145 J=1,NOVAR
145   ARRAY(K,J)=ARRAY(K,J)/DIAG
      ARRAY(K,K)=1.0/DIAG
      IF(NOSTP-250)1006,1450,1450
1450   WRITE(IOUT,1451)NOSTP
1451  FORMAT(25H ERROR, NUMBER OF STEPS =,I4,15H CHECK F-VALUES)
      GO TO 1004
99    WRITE(IOUT,13)NOIN
13    FORMAT(15H ERROR, NOIN IS,I7)
      GO TO 1004
C
C     BLOCK H-WRITE FINAL OUTPUT
C     COMMON ARRAY,STUDT,BETA,INDEX,NOVAR,DEFR,NOSTP,
C    1 NOIN,NOEN,NUMBR,YPRED,DEV,I,J,K,
C    2 NDATA,NTERM,IFINV,IFDEV,NOELM,IFEOF,INDX2,
C    3 COEN,SIGCO,K1,CNST,SIGY,FLIV
1007  IF( NOIN.EQ.0)GO TO 455
      DO 201 I=1,NOIN
      J=INDEX(I)
      STUDT(I)=COEN(I)/SIGCO(I)
201   BETA(I)=ARRAY(J,NOVAR)
      IF(NOEN)997,202,203
997   WRITE(IOUT,22)NOEN
22    FORMAT(14H ERROR, NOEN =,I6)
      GO TO 1004
203   WRITE(IOUT,601)NOSTP,KAME(K1)
601   FORMAT(///9H STEP NO. I3/20H   VARIABLE ENTERING I6)
      ISTRE = K1
      GO TO 204
202   IF(K1-ISTRE)8036,8035,8036
8035  IF(K1)8033,8036,8033
8033  WRITE(IOUT,8034)
8034  FORMAT (  42H VARIABLE REMOVED SAME AS VARIABLE ENTERED )
      GO TO 1004
8036  WRITE(IOUT,602)NOSTP,KAME(K1)
602   FORMAT(///9H STEP NO. I3/19H   VARIABLE REMOVEDI6)
204   COEFDT=1.0-ARRAY(NOVAR,NOVAR)
      COD=COEFDT-COE
      IF(COEFDT.LT.0.0)TYPE 2225
2225  FORMAT(' 4R')
      RRR=SQRT(COEFDT)
C****** WMU AM: 1.3.2-2, MTO, 30-JAN-78
	FPRB=FISHER(1,IDFALP,FLEV)
	TPRB=FISHER(1,IDFALP,TALPHA*TALPHA)
      WRITE(IOUT,600)FLEV,FPRB,SIGY,COEFDT,RRR,COD,IDFALP,CNST,SALPHA,
     1TALPHA,TPRB
600   FORMAT(12H   F LEVEL =,F12.4,5X,'F-PROB =',2X,F8.5,
     1/'   STANDARD ERROR OF ESTIMATE = ',
     1F12.4/'   COEFFICIENT OF DETERMINATION = ',F8.5/
     1'   COEFFICIENT OF MULTIPLE REGRESSION = ',F8.5/
     1'   INCREASE IN COEFFICIENT OF DETERMINATION = ',F8.5/
     13X,'DEGREES OF FREEDOM =',I6,/,
     15X,'CONSTANT',5X,'STD. ERR.',13X,'T',9X,'T-PROB'/1X,3G,F8.5//
     1' VARIABLE',4X,'COEFF.  STD ERR OF COEFF.',
     #17X,'T-VALUE',3X,'T-PROB',/37X,'STANDARDIZED COEFF.')
	DO 6002 J=1,NOIN
	TPRB=FISHER(1,IDFALP,STUDT(J)*STUDT(J))
C****** END = REGGR, STMT #6002-1
6002	WRITE(IOUT,6001)LAME(J),COEN(J),SIGCO(J),BETA(J),STUDT(J),TPRB
6001	FORMAT(2X,'X(',I2,')',F13.5,F13.6,F20.6,F9.3,F9.5)
	IF(DEFR.EQ.0)WRITE(IOUT,8002)
8002	FORMAT(/,' NO MORE VARIABLES ARE INCLUDED BECAUSE',/,
     1	' DEGREES OF FREEDOM ARE ZERO',/)
8003  IF(IFOPT(4).NE.1.AND.IFOPT(5).NE.1.AND.IFOPT(6).NE.1)GO TO 1000
      IF(IFOPT(4).NE.1)GO TO 2027
205   WRITE(IOUT,604)
604   FORMAT(///20X,'PREDICTED VS ACTUAL RESULTS' /
     111X,'OBS. NO.     ACTUAL      PREDICTED     DEVIATION' //)
2027  SSS=0.
      SST=0.
      DO 2070  J=1,NDATA
      NUMBR=J
      READ(1)(DATA(L),L=1,LTERM)
8     YPRED=CNST
      SSS=DATA(LDVAR)+SSS
      SST=SST+DATA(LDVAR)**2
11    DO 208I=1,NOIN
      K=LAME(I)
208   YPRED=YPRED+COEN(I)*DATA(K)
      DEV=DATA(LDVAR)-YPRED
51    SDEV2 = SDEV2 + DEV*DEV
33    IF(IDIF)35,36,35
35    DIF=DEV-DEVP
      SDIF2 = SDIF2 + DIF*DIF
      GO TO 37
36    IDIF = 1
37    DEVP = DEV
      IF(IFOPT(4).NE.1)GO TO 2070
207   WRITE(IOUT,605)NUMBR,DATA(LDVAR),YPRED,DEV
2070  CONTINUE
605   FORMAT(8H        I8,3F14.5)
      IF(IFOPT(5).NE.1)GO TO 2472
C**** WMU AM: 1.3.2-2, MTO, 30-JAN-78
      IF (IFOPT(7).NE.1) SST=SST-SSS*SSS/FLOAT(NDATA)
      SSR=SST-SDEV2
      NERR=IDFALP
      NTOTL=IDFALP+NOIN
      ERMS=SDEV2/FLOAT(NERR)
      RMS=SSR/FLOAT(NOIN)
      F=RMS/ERMS
	FPRB=FISHER(NOIN,IDFALP,F)
C**** END = REGGR, STMT #2432-2
      WRITE(IOUT,2432)
2432  FORMAT(//20X,'ANALYSIS OF VARIANCE')
      WRITE(IOUT,2020)SSR,NOIN,RMS,F,FPRB,SDEV2,NERR,ERMS,SST,NTOTL
2020  FORMAT(/4X,'SOURCE',6X,'SUM OF SQ.',6X,'DF',5X,'MEAN SQ.',
     16X,'F',6X,'F-PROB'/' REGRESSION',F17.6,I6,F15.5,F9.3,F9.5/,
     16X,'ERROR',F17.6,
     2I6,F15.5/6X,'TOTAL',F17.6,I6)
2472  IF(IFOPT(6).NE.1)GO TO 1000
      DWS=SDIF2/SDEV2
      DW1=4.-DWS
      WRITE(IOUT,40)SDIF2,DWS,DW1
40    FORMAT(//14X,38HDURBIN-WATSON TEST FOR AUTOCORRELATION//
     110X,28HSUM OF SQUARED DIFFERENCES =,F14.4/10X,5HDWS =,
     2F12.8/10X,9H4 - DWS =,F12.8)
1000  RETURN
1004	CALL TTYPTY(IERCOD)
	IF(IERCOD.NE.0)CALL EXIT
	RETURN
      END

