Trailing-Edge
-
PDP-10 Archives
-
decuslib10-09
-
43,50466/tscd.for
There are 2 other files named tscd.for in the archive. Click here to see a list.
C WESTERN MICHIGAN UNIVERSITY
C TSCD.FOR (FILE NAME ON LIBRARY DECTAPE)
C TSCD, 1.13.1 (CALLING NAME, SUBLST NO.)
C TIME SERIES CHANGE DETECTION
C THE TEST PROCEDURES UTILIZED IN THIS PROGRAM ARE
C MODIFICATIONS OF THE METHODS AND PROGRAMMING SHOWN IN
C "CHANGE DETECTION MODEL FOR SERIALLY CORRELATED DATA"
C BY JONES, CROWELL, AND KAPUNIANI IN
C PSYCHOLOGICAL BULLETIN, 1969, VOL. 71, NO. 5, 352-358.
C PROGRAMMED BY B. HOUCHARD, ACADEMIC COMPUTER CENTER, WMU.
C LIBRARY DECTAPE PROGS. USED: USAGE.MAC
C FORWMU PROGS. USED: TTYPTY, EXISTS, PRINTS, DEVCHG
C BNKLIB PROGS. USED: GETID, GETFR1, FISHER
C INTERNAL SUBR. USED: TCHAN, PLOT, CHIPRB, CUNO, IOTS
C ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS PUT IN BY WG
C
C LIMITATIONS:
C
C (1) # OF POST DATA IS BETWEEN 1 AND 2000
C (2) # OF POST DATA + # OF PRE DATA SHOULD BE BETWEEN 1 AND 4000
C (3) ONLY 1 LINE ALLOTED FOR OBJECT TIME FORMAT
C (4) ONLY F-TYPE FORMAT ALLOWED
C
C***********************************************************************
C
DIMENSION X(4000),POST(2000),IDATE(2),AVE(2),S(2),B(2),SP(2),
1 SS(2),PP(2),T(2),NOTF(16)
COMMON/IOBLK/IDLG,ICC,INP,IOUT,IO2,IO3,ICODE,IPLT,NAMI(2),ITYCH
COMMON/IOBLKA/NAMO(2),IPJ,IPG,NCOPYS
COMMON/SGETFR/ISTD,ITYPE
COMMON/SID/ID(16),ISTOP
COMMON/BLOCK1/N,MA,X,POST
DOUBLE PRECISION NAMM
EQUIVALENCE (NAMM,NAMO)
DATA PP/'PRE ','POST'/
C
C*********************************************************************
C DEVICES USED:
C IDLG--DEVICE USED TO COMMUNICATE WITH USER
C IT IS ALWAYS SET TO -1
C ICC---DEVICE USED TO ACCEPT USER'S RESPONSES
C IT IS ALWAYS SET TO -4
C INP---DEVICE USED TO READ IN THE DATA
C ITS LOGICAL NUMBER IS ASSIGNED BY SUBRUTINE IOTS
C IOUT--DEVICE USED TO WRITE OUT THE REPORT
C ITS LOGICAL NUMBER IS ASSIGNED BY SUBROUTINE IOTS
C NOUT--DEVICE USED TO PLOT THE GRAPH
C NOUT= 1 IF IO3 IS NOT 'TTY'
C = IOUT OTHERWISE
C***********************************************************************
C
IDLG=-1
ICC=-4
INP=2
IOUT=3
C
C***********************************************************************
C ADD ONE TO LIBRARY PROGRAM USAGE COUNTER
C***********************************************************************
C
C CALL USAGE('TSCD')
C
C***********************************************************************
C CHECK IF JOB IS FROM TELETYPE OR PSEUDO-TELETYPE
C IF ICODE = 0 JOB IS FROM TELETYPE
C = -1 JOB IS FROM PSUEDO-TELETYPE
C***********************************************************************
C
CALL TTYPTY(ICODE)
C
C***********************************************************************
C GATHER ALL INPUT/OUTPUT INFORMATION
C OUTPUT OPTION IS AVAILABLE ONLY ONCE IN THE PROGRAM
C***********************************************************************
C
CALL IOTS(1)
IPLT=0
IDUM=0
NOUT=IOUT
10 CALL IOTS(0)
C
C***********************************************************************
C FORMAT SUBROUTINE, ITYPE = 0 MEANS F-TYPE ONLY
C***********************************************************************
C
ITYPE=0
C---------------NOTF IS RETURNED.
CALL GETFR2(0,16,NOTF)
C
C***********************************************************************
C GATHER ALL OTHER PROBLEM INFORMATON
C***********************************************************************
C
C****** WMU AM: 1.13.1-1, MTO, 2-FEB-78
WRITE (IDLG,36)
36 FORMAT (' ENTER HEADER')
READ (ICC,37) ID
37 FORMAT (16A5)
DO 38 ISTOP = 16, 1, -1
38 IF (ID(ISTOP).NE.' ') GOTO 39
39 CONTINUE
C****** END = MAIN, STMT #30-1
C
C
C
30 WRITE(IDLG,31)
31 FORMAT(' # OF PRE DATA AND # OF POST DATA, SEPARATED BY COMMA--'
1,$)
READ(ICC,32) N,MA
32 FORMAT(2I)
IF ((MA.GE.1).AND.(MA.LE.2000)) GO TO 330
WRITE(IDLG,331) MA
331 FORMAT(' PROGRAM RESTRICTION: # OF POST DATA SHOULD BE
1 BETWEEN 1 AND 2000'/23X,'INCLUSIVE. USER ENTERED ',I7,
2 ' TRY AGAIN.'/)
GO TO 332
330 NMA=N+MA
IF ((NMA.GT.0).AND.(NMA.LE.4000)) GO TO 40
WRITE(IDLG,33) NMA
33 FORMAT(' PROGRAM RESTRICTION: # OF PRE DATA + # OF POST DATA
1 SHOULD BE'/23X,'BETWEEN 1 AND 4000 INCLUSIVE. USER ENTERED'/
2 23X,I7,', TRY AGAIN.'/)
332 IF (ICODE.GE.0) GO TO 30
34 WRITE(IDLG,35)
35 FORMAT('+JOB TERMINATES'/)
CALL EXIT
40 WRITE(IDLG,41)
41 FORMAT(' LINEAR TRENDS IN THE PRE AND POST DATA? (YES OR NO)--',
1 $)
READ(ICC,42) ITREND
42 FORMAT(A3)
WRITE(IDLG,44)
44 FORMAT(' MODEL ADEQUACY CHECK? (YES OR NO)--',$)
READ(ICC,42) IMAC
WRITE(IDLG,43)
43 FORMAT(' PLOT OF THE PRE, POST AND PREDICTED POST DATA? (YES
1 OR NO)--',$)
READ(ICC,42) IPLOT
IF (IPLOT.NE.'YES') GO TO 500
IF (IO3.NE.'TTY') GO TO 500
IDUM=IDUM+1
IF (IDUM.GT.1) GO TO 500
WRITE(IDLG,432)
432 FORMAT('-THE GRAPH WILL BE AT THE OUTPUT WINDOW AFTER THIS
1 RUN'/)
NOUT=1
IPLT=1
NAMO(1)='TSGRP'
NAMO(2)='.DAT'
CALL DEFINE FILE(NOUT,0,NV,NAMM,0,0)
C
C************************************************************************
C ADJUST FORMAT IF NECESSARY
C START READING DATA
C***********************************************************************
C
500 IF (ISTD.EQ.0) GO TO 51
NOTF(1)='(10F)'
DO 50 I=2,16
50 NOTF(I)=' '
51 IF (IO2.EQ.'TTY') GO TO 520
WRITE(IDLG,521)
521 FORMAT(' PLEASE WAIT, YOUR DATA IS BEING PROCESSED'/)
GO TO 54
520 WRITE(IDLG,52)
52 FORMAT(' ENTER PRE DATA'/)
IF (ISTD.NE.0) WRITE(IDLG,53)
53 FORMAT('+10 NUMBERS PER LINE, SEPARATED BY COMMAS'/)
54 READ(INP,NOTF)(X(I),I=1,N)
IF (IO2.NE.'TTY') GO TO 56
WRITE(IDLG,55)
55 FORMAT(' ENTER POST DATA'/)
IF (ISTD.NE.0) WRITE(IDLG,53)
56 READ(INP,NOTF)(X(I),I=N+1,NMA)
C
C***********************************************************************
C START WRITING OUT THE REPORT
C***********************************************************************
C
IF (IO3.EQ.'TTY') GO TO 65
CALL DATE (IDATE)
WRITE(IOUT,60) IDATE,ID,NOTF,N,MA,ITREND,IMAC,IPLOT
60 FORMAT('1'/'3WESTERN MICHIGAN UNIVERSITY'/' TIME SERIES
1 TEST FOR CHANGE DETECTION PROGRAM'/'-CALLING NAME: TSCD'/
2 ' DATE RUN: ',2A5/'-TITLE: ',16A5/' FORMAT: ',16A5//
3 '-NUMBER OF PRE DATA',25('.'),I5/' NUMBER OF POST DATA',
4 24('.'),I5/' LINEAR TRENDS?',29('.'),2X,A3/
5 ' MODEL ADEQUACY CHECKS?',21('.'),2X,A3/
6 ' PLOT OF THE PRE, POST AND PREDICTED DATA?.. ',A3)
WRITE(IOUT,610)
610 FORMAT('1')
GO TO 70
65 WRITE(IOUT,66) (ID(I),I=1,ISTOP)
66 FORMAT('1',16A5)
WRITE(IOUT,660) N,MA
660 FORMAT(' NUMBER OF PRE DATA =',I6/' NUMBER OF POST DATA =',I6)
70 CALL TCHAN(IMAC,AVE(1),S(1))
C
C***********************************************************************
C LINEAR TRENDS TESTS OPTION
C***********************************************************************
C
IF (ITREND.NE.'YES') GO TO 90
S(2)=0
AVE(2)=0
DO 81 I=N+1,NMA
81 AVE(2)=AVE(2)+X(I)
AVE(2)=AVE(2)/MA
K=0
IEND=N
SS(1)=FLOAT((N-1)*N*(N+1))/12
SS(2)=FLOAT((MA-1)*MA*(MA+1))/12
DO 83 I=1,2
SP(I)=0
DO 84 J=1,IEND
84 SP(I)=SP(I)+(X(J+K)-AVE(I))*J
B(I)=SP(I)/SS(I)
WRITE(IOUT,85) PP(I),AVE(I),PP(I),B(I)
85 FORMAT('-AVERAGE OF ',A4,'-DATA =',F12.4/' ESTIMATED LINEAR
1 SLOPE OF ',A4,'-DATA =',G12.4)
K=N
IEND=MA
83 CONTINUE
WRITE(IOUT,88)
88 FORMAT('-NOTE: IF EITHER ONE OF THE TRENDS IS DIFFERENT FROM
1 ZERO, THEN THE'/8X,'DATA IS PERHAPS NOT STATIONARY, AND THE MODEL
2 ADEQUACY OPTION'/8X,'SHOULD BE USED TO CHECK THE VALIDITY OF THE
3 FIRST ORDER MARKOV'/8X,'MODEL.')
C
C*************************************************************************
C PLOT OPTION
C CREATE A DISK FILE FOR TELETYPE OUTPUT JOBS
C***********************************************************************
C
90 IF (IPLOT.EQ.'YES') CALL PLOT(NOUT,ID)
WRITE(IDLG,91)
91 FORMAT('-')
GO TO 10
END
C
C SUBROUTINE TCHAN TEST FOR CHANGE USING FIRST-ORDER
C AUTOREGRESSIVE MODEL
C
C
C THIS SUBROUTINE WAS ORIGINALLY TAKEN FROM THE ARTICLE
C "CHANGE DETECTION MODEL FOR SERIALLY CORRELATED DATA"
C BY RICHARD H. JONES, DAVID H. CROWELL AND LINDA KAPUNIAI,
C PSYCHOLOGICAL BULLETINE 1969, VOL. 71, NO.5, PAGES 352-358.
C
C
C THE ABOVE SUBROUTINE IS COMPLETELY REWRITTEN. THE CHANGES ARE:
C
C (1) REPLACE AND ADD MORE VECTORS,
C (2) PRESERVE PRE AND POST DATA VECTOR FOR FUTURE USE,
C (3) MORE CALCULATIONS ADDED,
C (4) REPORT ON RESULTS ARE WRITTEN IN THE SUBROUTINE.
C
C SOME OF THE ABOVE CHANGES ARE RECOMMENDED BY DR. MICHAEL
C STOLINE, STATISTICAL CONSULTANT,WMU COMPUTER CENTER.
C
C WRITTEN BY BERENICE G. HOUCHARD
C WMU COMPUTER CENTER
C MAY, 1973
C
C
C---------------IMAC IS INPUT. OTHER ARGS. ARE RETURNED. IOUT
C--------------- IS INUT THRU COMMON /IOBLK/. W, MA, X ARE INPUT THRU
C--------------- /BLOCK1/. POST IS RETURNED THRU COMMON /BLOCK1/.
C SUBROUTINE TCHAN(IMAC,AVE,R1)
C ARGUMENTS RETURNED TO THE CALLING PROGRAM:
C AVE-------MEAN OF PRE-DATA
C R(1)-----
C
C
C
SUBROUTINE TCHAN(IMAC,AVE,R1)
DIMENSION X(4000),POST(2000),RR(2000),T(2000),TR(1898),
1 FRR(2000),R(2000),E(2000),EX(1998),M(3),MD(3),CH(3),P(3)
C
C
C TR(N-2)--T TEST ON PARTIAL SERIAL CORRELATION COEFFICIENTS LAGS
C 2 TO N-1
C T(MA)----T TEST ON POST DIFFERENCES
C EX(N-2)--PARTIAL SERIAL CORRELATION COEFFICIENTS
C RR(N)----SERIAL CORRELATIONS
C R(N)-----WORKING STORAGE VECTOR
C
COMMON/IOBLK/IDLG,ICC,INP,IOUT,IO2,IO3,ICODE,IPLT,NAMI(2),ITYCH
COMMON/BLOCK1/N,MA,X,POST
EQUIVALENCE (RR,T,TR),(FRR,R)
DATA STAR/'*'/
C
C CALCULATION OF PRE-MEAN
C
AVE=0
DO 10 I=1,N
10 AVE=AVE+X(I)
AVE=AVE/N
C
C ESTIMATE COVARIANCE FUNCTION OF PRE
C
DO 20 I=1,N
R(I)=0
L=I-1
DO 19 J=1,N-L
19 R(I)=R(I)+(X(J+L)-AVE)*(X(J)-AVE)
20 RR(I)=R(I)/R(1)
C
C ESTIMATE PARAMETERS
C
R1=R(1)
A=R(2)/R(1)
V=R(1)-A*R(2)
N2=N-2
SD=SQRT(V/N2)
NEND=82
IF (N.LT.5) NEND=N*11+27
WRITE(IOUT,50) AVE,A,SD,(STAR,I=1,NEND)
50 FORMAT('-MEAN OF PRE DATA =',F12.4/
1 ' ESTIMATED AUTOREGRESSION COEFFICIENTS =',F12.4/
2 ' ESTIMATED ONE STEP PREDICTION STANDARD DEVIATION =',F12.4//
3 '-',13A1/' * L A G *'/
4 ' * FROM- TO *',9X,'S E R I A L',5X,'CORRELATIONS'/
5 ' ',69A1)
IEND=1
51 IST=IEND+1
IEND=IST+4
IF (IEND.GT.N) IEND=N
IST1=IST-1
IEND1=IEND-1
WRITE(IOUT,52) IST1,IEND1,(RR(I),I=IST,IEND)
52 FORMAT(' *',I4,' -',I4,' *',5F11.4)
IF (IEND.NE.N) GO TO 51
WRITE(IOUT,53)
53 FORMAT(' *',11X,'*')
IF (IMAC.NE.'YES') GO TO 40
C
C CALCULATE PARTIAL CORRELATION FUNCTION T-TEST
C
M(1)=N/3
M(2)=N/2
M(3)=2*N/3
DO 100 I=3,M(3)
FN=N+1-I
EX(I-2)=(R(I)-A*R(I-1))/V
EE=EX(I-2)*EX(I-2)
100 TR(I-2)=EX(I-2)*SQRT(FN/(1-EE))
C
C CALCULATE RESIDUES OF PRE-DATA FIRST ORDER AUTOREGRESSIVE MODEL
C AND CALCULATE AUTOCORRELATION OF FITTED RESIDUALS AND BOX-PIERCE
C TEST
C
E(1)=X(1)-AVE
AVV=E(1)
DO 101 I=2,N
E(I)=(X(I)-AVE)-A*(X(I-1)-AVE)
101 AVV=AVV+E(I)
AVV=AVV/N
FR1=0
DO 102 I=1,N
102 FR1=FR1+(E(I)-AVV)**2
DO 104 I=2,N
FR=0
L=I-1
DO 103 J=1,N-L
103 FR=FR+(E(J+L)-AVV)*(E(J)-AVV)
104 FRR(I)=FR/FR1
DO 200 I=1,3
MD(I)=M(I)-2
200 CH(I)=0
DO 201 I=2,M(1)+1
201 CH(1)=CH(1)+FRR(I)**2
CH(1)=N*CH(1)
DO 202 I=M(1)+2,M(2)+1
202 CH(2)=CH(2)+FRR(I)**2
CH(2)=N*CH(2)+CH(1)
DO 203 I=M(2)+2,M(3)+1
203 CH(3)=CH(3)+FRR(I)**2
CH(3)=N*CH(3)+CH(2)
DO 210 I=1,3
CALL CHIPRB(MD(I),CH(I),P(I),IERR)
IF (IERR.EQ.1) P(I)=9999E-20
210 CONTINUE
C
C STANDARD ERROR AND CONFIDENCE LIMITS FOR AUTOREGRESSIVE
C COEFFICIENT
C
SEE=SQRT((1-A*A)/N)
CLL=A-2*SEE
CLU=A+2*SEE
C
C
C
WRITE(IOUT,300) (STAR,I=1,NEND)
300 FORMAT('-',13A1/' * L A G *'/' * FROM- TO *',9X,'SERIAL
1 CORRELATIONS OF RESIDUALS'/1X,69A1)
IEND=1
301 IST=IEND+1
IEND=IST+4
IF (IEND.GT.N) IEND=N
IST1=IST-1
IEND1=IEND-1
WRITE(IOUT,52) IST1,IEND1, (FRR(I),I=IST,IEND)
IF (IEND.NE.N) GO TO 301
WRITE(IOUT,302)
302 FORMAT(' *',11X,'*'/'-BOX-PIERCE CHI-SQUARE TESTS FOR RESIDUAL
1 CORRELATION IN THE FITTED MODEL'//2X,'CHI-SQUARE',3X,'DEGREES
2 OF'/4X,'VALUES',6X,'FREEDOM',5X,'PROBABILITY'/2X,10('-'),
3 3X,10('-'),3X,11('-')/)
DO 303 I=1,3
303 WRITE(IOUT,304) CH(I),MD(I),P(I)
304 FORMAT(1X,F11.3,I8,9X,F6.4)
WRITE(IOUT,305) A,SEE,CLL,CLU,N2
305 FORMAT('-ESTIMATED AUTOREGRESSION COEFFICIENT =',F12.3/
1 ' STANDARD ERROR OF AUTOREGRESSION COEFFICIENT =',F12.3/
2 '-95% CONFIDENCE LIMITS FOR THE AUTOREGRESSION COEFFICIENT'/
3 ' LOWER LIMIT =',F12.3,8X,'UPPER LIMIT =',F12.3/
4 '-T-TESTS FOR PARTIAL SERIAL CORRELATION COEFFICIENT'/
5' (EACH T-TEST IS BASED ON ',I7,' DEGREES OF FREEDOM)'/
6 '-',5X,'* PARTIAL SERIAL'/2X,'LAG *',3X,'CORR COEFF',8X,'T VALUE'
7 /1X,35('*'))
DO 306 I=2,M(3)
306 WRITE(IOUT,307) I,EX(I-1),TR(I-1)
307 FORMAT(1X,I4,' *',2F13.4)
WRITE(IOUT,308)
308 FORMAT(6X,'*'/'-NOTE: IF THE BOX-PIERCE CHI-SQUARE TEST(S) ARE
1 SIGNIFICANT (LARGE)'/8X,'THEN THE FIRST ORDER AUTOREGRESSIVE
2 MODEL DOES NOT FIT THE'/8X,'PRE-TEST DATA AND SOME OTHER MODEL
3 SHOULD BE USED.'//8X,'IF THE 95% CONFIDENCE INTERVAL FOR THE
4 AUTOREGRESSIVE'/8X,'COEFFICIENT CONTAINS ONE, THEN THIS IS
5 EVIDENCE THAT THE'/8X,'PRE-TEST DATA IS NON-STATIONARY. THE
6 NEW SPECIFIED MODEL'/8X,'SHOULD BE BASED ON AT LEAST FIRST
7 DIFFERENCES OF THE PRE-TEST'/8X,'DATA.'//8X,'IF SOME OF T-TESTS
8 FOR PARTIAL SERIAL CORRELATIONS ARE'/8X,'SIGNIFICANTLY
9 DIFFERENT FROM ZERO, THEN THIS IS EVIDENCE THAT'/8X,'THE NEW
1 MODEL BE BASED ON EITHER MORE AUTOREGRESSIVE PARAMETERS'/8X,
2 'OR A MOVING AVERAGE STRUCTURE.'/)
C
40 F=0
DO 41 J=1,MA
K=MA+1-J
POST(K)=AVE+A*(X(N+K-1)-AVE)
T(K)=(X(N+K)-POST(K))/SD
41 F=F+T(K)**2
F=F/MA
MEND=82
IF (MA.LT.5) MEND=MA*11+27
WRITE(IOUT,60) (STAR,I=1,MEND)
60 FORMAT(//'-',13A1/
1 ' *OBSERVATION*'/
2 ' * FROM- TO *',9X,'PREDICTED POST DATA'/
3 ' ',69A1)
IEND=0
61 IST=IEND+1
IEND=IST+4
IF (IEND.GT.MA) IEND=MA
WRITE(IOUT,52) IST,IEND,(POST(I),I=IST,IEND)
IF (IEND.NE.MA) GO TO 61
WRITE(IOUT,53)
N2=N-2
PROB=FISHER(MA,N2,F)
WRITE(IOUT,70) F,MA,N2,PROB,N2,(STAR,I=1,MEND)
70 FORMAT(//'-F TEST ON ENTIRE POST INTERVAL HAS F-VALUE =',F12.4/
1 ' WITH ',I6,' AND',I6,' DEGREES OF FREEDOM'/
2 ' PROBABILITY VALUE =',F10.5/
3 '-NOTE: A SIGNIFICANT F-VALUE IS EVIDENCE THAT A CHANGE
4 IN MEAN'/8X,'RESPONSE HAS OCCURRED BETWEEN PRE AND POST
5 CONDITIONS.'/8X,'THE T-VALUES BELOW CAN BE USED TO IDENTIFY
6 SPECIFIC SIGNIFICANT'/8X,'DEVIANT POST OBSERVATIONS'///
7 '-T TESTS ON POSTDIFFERENCES'/
8 ' (EACH T-VALUE HAS',I6,' DEGREES OF FREEDOM)'/
9 '-',13A1/' *POST OBSER.*'/
1 ' * FROM- TO *',9X,'T',5X,'V A L U E S'/1X,69A1)
IEND=0
71 IST=IEND+1
IEND=IST+4
IF (IEND.GT.MA) IEND=MA
WRITE(IOUT,52) IST,IEND, (T(I),I=IST,IEND)
IF (IEND.NE.MA) GO TO 71
WRITE(IOUT,53)
RETURN
END
C
C---------------BOTH ARGS. ARE INPUT. N, MA, X, POST ARE INPUT
C--------------- THRU COMMON /BLOCK1/.
C SUBROUTINE PLOT(NOUT,ID)
C
C THIS SUBROUTINE PLOTS THE PRE, POST AND PREDICTED POST DATA.
C
C ALL ARGUMENTS ARE PASSED FROM THE MAIN PROGRAM:
C
C NOUT--OUTPUT DEVICE
C ID----PROBLEM IDENTICATION VECTOR OR HEADER VECTOR
C
SUBROUTINE PLOT(NOUT,ID)
DIMENSION X(4000),POST(1),ID(1),CUTPT(11),ROW(0/100)
COMMON/BLOCK1/N,MA,X,POST
DATA DOT,CROSS,CIRCLE,STAR/'.','X','O','*'/
C
C***********************************************************************
C DETERMINE THE MAXIMUM AND MINIMUM OF THE DATA POINTS
C*************************************************************************
C
NMA=N+MA
XMIN=X(1)
XMAX=X(1)
DO 10 I=2,NMA
IF (X(I).LT.XMIN) XMIN=X(I)
IF (X(I).GT.XMAX) XMAX=X(I)
10 CONTINUE
DO 11 I=1,MA
IF (POST(I).LT.XMIN) XMIN=POST(I)
IF (POST(I).GT.XMAX) XMAX=POST(I)
11 CONTINUE
DELTA=(XMAX-XMIN)/100.
C
C***********************************************************************
C WRITE OUT HEADING FOR THE GRAPH
C***********************************************************************
C
DO 20 I=1,11
20 CUTPT(I)=XMIN+DELTA*(I-1)*10
WRITE(NOUT,21) (ID(I),I=1,16),DELTA,(CUTPT(I),I=1,11,2),
1 (CUTPT(I),I=2,11,2)
21 FORMAT('1TIME SERIES GRAPH ON PRE, POST AND PREDICTED POST
1 DATA'//1X,16A5/'-GRAPH INTERVAL =',G12.4/' NOTE: . = PRE
2 DATA'/8X,'X = POST DATA'/8X,'O = PREDICTED POST DATA'/
3 8X,'* = INTERSECTION OF POST AND PREDICTED POST DATA'/
4 '-',6(G10.4,10X)/7X,6('^',19X)/7X,5('^',3X,G10.4,6X),'^ ',
5 'PRE & POST PREDICTED'/7X,10('^',9X),'^ D A T A',5X,
6 'POST DATA'/7X,10('+---------'),'+')
C
C***********************************************************************
C CALCULATE THE INDICES AND PRINT OUT ONE ROW AT A TIME
C***********************************************************************
C
ROUND=DELTA/2
DO 30 I=1,N
Y=X(I)
DO 31 J=0,100
31 ROW(J)=' '
I1=(Y-XMIN)/DELTA+ROUND
ROW(I1)=DOT
30 WRITE(NOUT,32) I,ROW,Y
32 FORMAT(1X,I4,2X,101A1,2X,G10.4,2X,G10.4)
DO 40 I=N+1,NMA
Y=X(I)
Z=POST(I-N)
DO 41 J=0,100
41 ROW(J)=' '
I1=(Y-XMIN)/DELTA+ROUND
I2=(Z-XMIN)/DELTA+ROUND
ROW(I1)=CROSS
ROW(I2)=CIRCLE
IF (I1.EQ.I2) ROW(I1)=STAR
I1=I-N
40 WRITE(NOUT,32) I1,ROW,Y,Z
RETURN
END
C THIS SUBROUTINE AND SUBROUTINE CUNO ARE TAKEN FROM STP.
C
C
C *** STAT PACK ***
C SUBROUTINE USED FOR DETERMINING CHI SQUARE PROBABILITIES
C CALLING SEQUENCE: CALL CHIPRB(K,X,Y,IERR)
C WHERE K - NUMBER OF DEGREES OF FREEDOM
C X - CHI SQUARE VALUE
C Y - PROBABILITY ASSOCIATED WITH CHI SQUARE
C IERR - ERROR WAS ENCOUNTERED WHEN ATTEMPTING TO CALCULATE
C PROBABILITY
C
C ROUTINE WAS WRITTEN BY CHARLES NAGY OF WESTERN, AND ADAPTED
C TO STP BECAUSE THE FISHER ROUTINE WOULD NOT GIVE ENOUGH
C ACCURACY FOR THE PSYCHOLOGISTS. CALLS SUBROUTINE CUNO.
C
C---------------K, X ARE INPUT. Y, IERR ARE OUTPUT.
SUBROUTINE CHIPRB(K,X,Y,IERR)
DIMENSION F(25)
F(1)=.5
F(2)=.598706326
F(3)=.691462461
F(4)=.773372648
F(5)=.841344746
F(6)=.894350226
F(7)=.933192799
F(8)=.959940843
F(9)=.977249868
F(10)=.987775527
F(11)=.993790335
F(12)=.997020237
F(13)=.998650102
F(14)=.999422975
F(15)=.999767371
F(16)=.999911583
F(17)=.999968329
F(18)=.999989311
F(19)=.999996602
F(20)=.999998983
F(21)=.999999713
F(22)=.999999924
F(23)=.999999981
F(24)=.999999996
F(25)=.999999999
IERR=0
Y=0
IF((K.LE.0).OR.(K.GT.100)) IERR=1
IF(X.GT.141) IERR=1
IF(IERR.EQ.1) RETURN
IF(X.LE.0) GO TO 13
IF(K.GE.4) GO TO 4
GO TO (1,2,3),K
1 P=SQRT(X)
CALL CUNO(P,S,IERR,F)
IF(IERR.EQ.1) RETURN
Y=2.*S-1
GO TO 13
2 Y=1.-(1./EXP(X/2.))
GO TO 13
3 P=SQRT(X)
CALL CUNO(P,S,IERR,F)
IF(IERR.EQ.1) RETURN
Y=(2.*S-1)-P/(1.25331414*EXP(X/2.))
GO TO 13
4 M=K/2
IF(K.EQ.2*M) GO TO 6
P=SQRT(X)
CALL CUNO(P,S,IERR,F)
IF(IERR.EQ.1) RETURN
Y=2.*S-1.
S=X/2.
C=1./(.62665707*P*EXP(S))
P=S
T=.5
GO TO 7
6 C=1./EXP(X/2.)
Y=1.-C
S=0
P=1
T=0
7 DO 8 I=1,M-1
T=T+1
P=P*(X/(T*2.))
8 S=S+P
Y=Y-C*S
13 Y=1.-Y
END
C *** STAT PACK ***
C SUBROUTINE USED IN FINDING PROB FOR CHI SQUARE
C CALLING SEQUENCE: CALL CUNO(X,Y,IERR,F)
C ORIGINALLY WRITTEN BY CHARLES NAGY OF WMU.
C
C---------------X, F ARE INPUT. Y, IERR ARE RETURNED.
SUBROUTINE CUNO(X,Y,IERR,F)
DIMENSION F(1)
W=X
IF(W.LT.0) W=-W
IF(W.LE.6.125) GO TO 2
1 Z=1
GO TO 7
2 K=INT(4.*W)+1
A=.25*(K-1.)
IF(W-A) 10,3,4
3 Z=F(K)
GO TO 7
4 IF(W-(A+.125))6,6,5
5 K=K+1
A=A+.25
6 H=W-A
ASQ=A*A
C1=((-ASQ+10.)*ASQ-15.)*A
C2=(6.*ASQ-36.)*ASQ+18.
C3=(-30.*ASQ+90.)*A
C4=120.*(ASQ-1.)
C5=-360.*A
C6=(((((C1*H+C2)*H+C3)*H+C4)*H+C5)*H+720.)*H
Z=F(K)+C6/(720.*SQRT(6.28318531*EXP(ASQ)))
7 Y=Z
IF(X.LT.0.) Y=1.-Z
RETURN
10 IERR=1
RETURN
END
C****** WMU AM: 1.13.1, #1, MTO, 25-JAN-78
C------ FISHER ROUTINE REMOVED ------
C****** END
C---------------IORO IS INPUT. IDLG, NDEVI, NDEVO, IDVI, IDVO, IPLT,
C--------------- ITYCH ARE INPUT THRU COMMON /IOBLK/. NDEVI, IDVI,
C--------------- IDVO ARE MODIFIED. NAMI IS RETURNED THRU COMMON /IOBLK/.
C---------------NAMO, IPJ, IPG, NCOPYS ARE RETURNED THRU COMMON
C--------------- /IOBLKA/.
SUBROUTINE IOTS(IORO)
DIMENSION IN(50),INAME(2),B(10),NAM(2)
COMMON/IOBLK/IDLG,IRSP,NDEVI,NDEVO,IDVI,IDVO,ICODE,IPLT,NAMI(2),
1 ITYCH
COMMON/IOBLKA/NAMO(2),IPJ,IPG,NCOPYS
DOUBLE PRECISION JNAME
EQUIVALENCE (JNAME,INAME)
L1="555004020100
L2="565004020100
IF(JONCE.EQ.0)ITMP=NDEVI
NDEVI=ITMP
IF ((IORO.AND.1).EQ.0)IDV=IDVI
IF ((IORO.AND.1).EQ.1)IDV=IDVO
1 GO TO(401,403,402,404),IORO+1
401 WRITE(IDLG,310)
310 FORMAT(' INPUT? (TYPE HELP IF NEEDED)--',$)
402 IDEV=NDEVI
GO TO 405
403 WRITE(IDLG,311)
311 FORMAT(' OUTPUT? (TYPE HELP IF NEEDED)--',$)
404 IDEV=NDEVO
405 READ(IRSP,10,END=201)IN
10 FORMAT(50A1)
IF ((IN(1).EQ.'F').AND.(IN(2).EQ.'I').AND.(IN(3).EQ.'N').AND.
1(IN(4).EQ.'I')) GO TO 201
IF ((IN(1).EQ.'S').AND.(IN(2).EQ.'A').AND.(IN(3).EQ.'M').AND.
1(IN(4).EQ.'E').AND.(IN(5).EQ.' ')) GO TO 212
IF ((IN(1).EQ.'H').AND.(IN(2).EQ.'E').AND.(IN(3).EQ.'L').AND.
1 (IN(4).EQ.'P')) GO TO (500,600),IORO+1
ITYFLG=0
CALL RELEAS(IDEV)
NEVER=0
ICOLN=0
ILBR=0
ISL=0
IPROJ=0
IPROG=0
INAME(1)=' '
INAME(2)=' '
IDV=' '
K=0
IDP=0
12 K=K+1
IF(K.GT.50)GO TO 15
IF (IN(K).EQ.'.') IDP=1
IF(IN(K).EQ.':')GO TO 13
IF(IN(K).EQ."555004020100)GO TO 14
IF(IN(K).EQ.'/')GO TO 23
GO TO 12
13 ICOLN=K+4
DO 20 I=50,K+4,-1
20 IN(I)=IN(I-4)
DO 27 I=0,3
27 IN(K+I)=' '
K=K+4
GO TO 12
14 ILBR=K+9
DO 21 I=50,K+9,-1
21 IN(I)=IN(I-9)
DO 22 I=K,K+8
22 IN(I)=' '
K=K+9
GO TO 12
23 ISL=K
GO TO 12
15 IF(ILBR.EQ.0)GO TO 31
30 ENCODE(12,40,B)(IN(I),I=ILBR+1,ILBR+12)
40 FORMAT(12A1)
DECODE(12,41,B)IPROJ,IPROG
41 FORMAT(2O)
31 IF (IDP.NE.0) GO TO 32
DO 33 I=ICOLN+9,ICOLN+1,-2
IF (IN(I).NE.' ') GO TO 34
33 CONTINUE
I=6
34 IN(I+1)='.'
32 ENCODE(10,42,INAME)(IN(I),I=ICOLN+1,ICOLN+10)
42 FORMAT(10A1)
IF ((INAME(1).EQ.'FINIS').OR.(INAME(1).EQ.'FINI')) GO TO 201
IF(ICOLN.EQ.0)GO TO 101
100 ENCODE(5,44,IDV)(IN(I),I=1,5)
44 FORMAT(5A1)
101 IF(ISL.EQ.0)GO TO 24
ENCODE(5,44,B)(IN(I),I=ISL+1,ISL+5)
DECODE(5,46,B)NCOPYS
46 FORMAT(I)
24 IF(IDV.NE.' ')GO TO 124
IF(INAME(1).EQ.' ')GO TO 28
IDV='DSK'
GO TO 124
28 IF(ICODE.EQ.-1)GO TO 125
IDV='TTY'
GO TO 124
125 IF((IORO.AND.1).EQ.0)IDV='CDR'
IF((IORO.AND.1).EQ.1)IDV='LPT'
124 CALL DEVCHG(IDV,IDEV)
IF(IDV.EQ.'DSK')GO TO 102
IF(IDV.EQ.'LPT')GO TO 104
IF(IDV.LE."422510134500.AND.IDV.GE."422510130100)GO TO 102
213 IF(IDV.EQ.'TTY'.AND.(IORO.AND.1).EQ.0)GO TO 214
GO TO 410
104 INAME(1)='OUTAA'
INAME(2)='A.AAA'
IPR=1
LPT=IDEV
CALL DEVCHG('DSK',IDEV)
105 CALL EXISTS(IDEV,INAME,MRK)
IF(MRK.EQ.1)GO TO 211
INAME(2)=INAME(2)+2
GO TO 105
211 NAM(1)=INAME(1)
NAM(2)=INAME(2)
102 IF(INAME(1).NE.' ')GO TO 302
IF((IORO.AND.1).EQ.0)INAME(1)='INPUT'
IF((IORO.AND.1).EQ.1)INAME(1)='OUTPT'
INAME(2)='.DAT'
302 IF((IORO.AND.1).EQ.1)GO TO 303
CALL EXISTS(IDEV,INAME,MRK,IPROJ,IPROG)
IF(MRK.EQ.0)GO TO 303
WRITE(IDLG,305)
305 FORMAT(' FILE DOES NOT EXIST'/)
IF(ICODE.EQ.-1)CALL EXIT
GO TO 1
303 CALL DEFINE FILE(IDEV,0,NEVER,JNAME,IPROJ,IPROG)
GO TO 213
201 IF (IPR.EQ.0) GO TO 2010
CALL RELEAS(LPT)
CALL PRINTS(NAM,2,1,NCOPYS)
CALL EXIT
2010 IF (IPLT.NE.1) CALL EXIT
ENDFILE 1
CALL PRINTS(NAMO,2,1,1)
CALL EXIT
212 IF(ITYFLG.EQ.1)GO TO 215
IF ((IORO.AND.1).EQ.O)REWIND IDEV
GO TO 410
C NO TTY: SAME OPTION IF NO CHANNEL PROVIDED IN ITYCH
214 IF(ITYCH.LT.1)GO TO 410
IF(IONCE.NE.1)CALL DEVCHG('DSK',ITYCH)
IONCE=1
IF(ITYFLG.EQ.1)GO TO 215
ITYFLG=1
CALL RELEAS(ITYCH)
CALL DEFINE FILE(ITYCH,0,NV,'TTYDAT.TMP',0,0)
410 IOROA=IORO.AND.1
IF(IOROA.EQ.1)GO TO 411
IPJ=IPROJ
IPG=IPROG
IDVI=IDV
NDEVI=IDEV
NAMI(1)=INAME(1)
NAMI(2)=INAME(2)
GO TO 412
411 NAMO(1)=INAME(1)
NAMO(2)=INAME(2)
IDVO=IDV
412 CONTINUE
JONCE=1
RETURN
215 REWIND ITYCH
IDEV=ITYCH
GO TO 410
500 WRITE(IDLG,501)
501 FORMAT('-THIS ANSWER DEFINES WHERE THE PROGRAM IS TO FIND THE
1 INPUT DATA. IT'/' USUALLY CONSISTS OF A DEVICE, POSSIBLY A
2 FILENAME WITH OR WITHOUT AN'/' EXTENSION, AND A PROJECT-
3PROGRAMMER NUMBER.'//' POSSIBLE DEVICES ARE:'//6X,'DEVICES',3X,
4 'DESCRIPTION'/6X,7('-'),3X,11('-')/6X,'TTY:',6X,'TERMINAL'/
5 6X,'DSK:',6X,'DISK (FILENAME AND EXTENSION, PROJECT-PROGRAMMER
6 NUMBER'/22X,'MAY BE USED)'/6X,'CDR:',6X,'CARD READER (THIS
7 DEVICE IS NOT APPLICABLE ON TERMINAL'/30X,'JOBS)'/6X,'DTA#:',5X,
8 'DECTAPE UNIT (USER''S DECTAPE SHOULD ALREADY BE MOUNTED)'/6X,
9 'MTA#:',5X,'MAGTAPE UNIT (USER''S MAGTAPE SHOULD ALREADY BE
1 MOUNTED'/30X,'AND POSITIONED)'///' DEFAULTS:'//' (1) IF NO INPUT
2 DEVICE IS SPECIFIED BUT A FILENAME IS GIVEN, THE'/6X,'DEFAULT
3 DEVICE WILL BE DSK:'//' (2) IF A DEVICE WHICH REQUIRES A
4 FILENAME AND EXTENSION IS SPECIFIED,'/6X,'BUT NO FILENAME IS
5 GIVEN, THE DEFAULT NAME WILL BE INPUT.DAT'//' (3) IF NO RESPONSE
6 IS GIVEN, I.E. A CARRIAGE RETURN <CR> IS ENTERED,'/6X,'THE
7 DEFAULT DEVICE IS TTY: ON JOBS RUN FROM TERMINALS; AND'/28X,'CDR:
8 ON BATCH JOBS'//' (4) IF DSK: IS SPECIFIED AS THE INPUT DEVICE
9 AND NO PROJECT-PROGRAMMER'/6X,'NUMBER IS GIVEN, THE USER''S
1 PROJECT-PROGRAMMER NUMBER WILL BE'/6X,'ASSUMED.'///)
WRITE(IDLG,502) L1,L2
502 FORMAT(' EXAMPLES: DATA.DAT'/14X,'TEST.DAT',A1,'420,420',A1/
2 14X,'MTA0:'/14X,'DTA2:FILE1'//' NOTE: THE FOLLOWING RESPONSES
3 ARE VALID AFTER THE FIRST "INPUT?"'//' (1) SAME COMMAND. IF THE
4 DATA FILE TO BE USED IS THE SAME AS THE'/6X,'PRECEEDING ONE, THE
5 USER MAY SIMPLY ENTER "SAME"'//' (2) FINISH COMMAND. THE USER
6 MUST ENTER "FINISH" TO EXIT FROM THE'/6X, 'PROGRAM. THIS ENSURES
7 THAT OUTPUT ASSIGNED TO LPT: WILL BE'/6X,'PRINTED. FAILURE TO
8 USE THE "FINISH" COMMAND MAY RESULT IN'/6X,'LOSING THE ENTIRE
9 OUTPUT FILE.'//' (3) A ^Z (CONTROL Z) WILL RESULT IN THE SAME
1 ACTION AS THE "FINISH"'/6X,'COMMAND.'///)
503 CALL RELEAS(IDLG)
GO TO (401,403), IORO+1
600 WRITE(IDLG,601)
601 FORMAT('-THE ANSWER DEFINES WHERE THE OUTPUT FROM THE PROGRAM
1 IS TO BE PLACED.'/' IT USUALLY CONSISTS OF A DEVICE AND POSSIBLY
2 A FILENAME WITH OR WITH-'/' OUT AN EXTENSION.'//' POSSIBLE
3 DEVICES ARE:'//6X,'DEVICE',3X,'DESCRIPTION'/6X,6('-'),3X,
4 11('-')/6X,'TTY:',5X,'TERMINAL'/6X,'DSK:',5X,'DISK (FILENAME
5 AND EXTENSION MAY BE USED)'/6X,'LPT:',5X,'LINEPRINTER (MULTIPLE
6 COPIES MAY BE REQUESTED BY USE OF'/29X,'THE "/COPIES" COMMAND)'/
7 6X,'DTA#:',4X,'DECTAPE UNIT (USER''S DECTAPE SHOULD ALREADY
8 BE MOUNTED;'/29X,'FILENAME AND EXTENSION MAY BE USED.)'/
9 6X,'MTA#:',4X,'MAGTAPE UNIT (USER''S MAGTAPE SHOULD ALREADY
1 BE MOUNTED'/29X,'AND POSITIONED)'///' DEFAULTS:'//' (1) IF NO
2 OUTPUT DEVICE IS SPECIFIED BUT A FILENAME IS GIVEN, THE'/6X,
3 'DEFAULT DEVICE WILL BE DSK:'//' (2) IF A DEVICE WHICH REQUIRES
4 A FILENAME AND EXTENSION IS SPECIFIED,'/6X,'BUT NO FILENAME IS
5 GIVEN, THE DEFAULT NAME WILL BE OUTPT.DAT'//' (3) IF NO RESPONSE
6 IS GIVEN, I.E. A CARRIAGE RETURN <CR> IS ENTERED,'/6X,'THE
7 DEFAULT DEVICE IS TTY: ON JOBS RUN FROM TERMINALS; AND'/28X,'LPT:
8 ON BATCH JOBS'//' (4) IF LPT: IS LISTED AS THE OUTPUT DEVICE,
9 THE NUMBER OF COPIES WILL'/6X,'DEFAULT TO 1.'///
1' EXAMPLES: LPT:/2'/14X,'RPT.DAT'/14X,'DTA0:NAME.DAT'///)
GO TO 503
END
SUBROUTINE GETFR2(ISEL,NCHAR,IFMT)
DIMENSION IFMT(1),ITYPE(0/2),NAM(2),IN(80)
C****** WMU AM: 1.13.1-1, MTO, 2-FEB-78
COMMON /IOBLK/ IDLG
C****** END = GETFR2, STMT #99-5
COMMON /SGETFR/ISTD,LTYPE
DOUBLE PRECISION NAME
EQUIVALENCE (NAME,NAM)
DATA ITYPE/'F','A','I'/
C
C
99 IBEG=1
ISTD=0
NPT=0
IPAREN=0
IF(NCHAR.EQ.0) PAUSE 'SYSTEM ERROR CONTACT COMPUTER CENTER'
C
C WRITE HEADER
C
IF((LTYPE.EQ.3).OR.(ISEL.EQ.1)) WRITE(IDLG,100)
100 FORMAT(' ENTER FORMAT ENCLOSED IN PARENTHESES'/)
IF((LTYPE.LT.3).AND.(ISEL.EQ.0)) WRITE(IDLG,101) ITYPE(LTYPE)
101 FORMAT(' ENTER FORMAT ENCLOSED IN PARENTHESES: (',A1,
. '-TYPE ONLY)'/)
98 IF(NPT.LT.2) CALL GES(IN,80,IERR)
IF(NPT.EQ.2) READ(29,102) IN
C
C FIND # OF CHARACTERS
C
DO 10 I=80,1,-1
IF(IN(I).NE.' ') GOTO 20
10 CONTINUE
ISTD=1
RETURN
C
C STANDARD FORMAT REQUESTED; RETURN
C
20 LAST=I
C
C GET RID OF EXTRA SPACES
C
I=0
25 I=I+1
IF(I.GT.LAST) GOTO 30
IF(IN(I).NE.' ') GOTO 25
DO 15 J=I,LAST
15 IN(J)=IN(J+1)
LAST=LAST-1
GOTO 25
C
C
30 IF(NPT.NE.0) GOTO 300
IF(IN(1).EQ.'H'.AND.IN(2).EQ.'E'.AND.IN(3).EQ.'L'.AND.IN(4)
. .EQ.'P') GOTO 999
IF(IN(1).EQ.'S'.AND.IN(2).EQ.'A'.AND.IN(3).EQ.'M'.AND.IN(4).EQ.
. 'E') RETURN
IF((IN(1).NE.'(').AND.(IN(1).NE."401004020100)) GOTO 9999
NPT=1
DO 35 I=1,NCHAR
35 IFMT(I)=' '
IF(IN(1).EQ.'(') GOTO 300
C
C COMMAND FILE
C
NPT=2
DO 40 I=2,LAST
IF(IN(I).EQ.'.') GOTO 45
40 CONTINUE
LAST=LAST+1
IN(LAST)='.'
45 J=LAST-1
ENCODE(J,102,NAM(1)) (IN(I),I=2,LAST)
102 FORMAT(80A1)
CLOSE(UNIT=29)
OPEN(UNIT=29,DEVICE='DSK',FILE=NAME,ACCESS='SEQIN')
GOTO 98
C
C READ FORMAT
C
300 DO 50 I=1,LAST
IF(IN(I).EQ.'(') IPAREN=IPAREN+1
IF(IN(I).EQ.')') IPAREN=IPAREN-1
50 CONTINUE
IF(IBEG+((LAST+4)/5).GT.NCHAR) GOTO 9999
ENCODE (LAST,102,IFMT(IBEG)) (IN(I),I=1,LAST)
IBEG=IBEG+(LAST+4)/5
IF(IPAREN.LT.1) GOTO (200,201),NPT
GOTO 98
C
C RETURN
C
201 CLOSE (UNIT=29)
200 RETURN
C
C ERROR AND HELP
C
9999 WRITE(IDLG,103)
103 FORMAT('-ERROR: Format incorrectly specified'/)
GOTO 99
999 I=NCHAR*5
WRITE(IDLG,104) I
104 FORMAT(' Any FORMAT specification must comply with the FORTRAN-1
.0 Format'/' requirement. The FORMAT must also be enclosed
.in parentheses'/' and be no more than ',I3,' characters in length'
. //' Example: ENTER FORMAT ENCLOSED IN PARENTHESES'
. /11X,'(I2,F3.0,1X,F2.0,I1)'/)
GOTO 99
END