10 REM     MULTREG1 (MULTREG Version 1.1)   7/9/85
12 REM with double precision when not compiled
15 REM COPYRIGHT (c) 1983, 1985 BY GREGORY C. DAHL
17 REM    * * *  ALL RIGHTS RESERVED  * * *
20 REM THIS PROGRAM RUNS A MULTIPLE REGRESSION
200 REM
210 REM DETAILED INSTRUCTIONS ARE GIVEN IN MULTREG.DOC
220 REM
230 ' Memory management under Microsoft Compiler has
235 ' been a problem.  Fixed data arrays for X and Y 
240 ' are inflexible and use up too much memory during
245 ' link loading.  The CLEAR command does not work
250 ' with OBSLIB.  An effort to store X and Y in
255 ' binary format in string variables produced an
260 ' unknown execution error.  Compilation has thus
263 ' not succeeded with OBSLIB.
265 '
266 ' The program can, however, be compiled with the
267 ' BRUN module, following instructions just below.
268 ' Data variables are poked into protected high
269 ' memory.  Current compiled dimension limits are
270 ' 17 indep. vbles and 300 observations.  Memory
280 ' limit also applies.
290 '
300 DEFINT A-C,F,I-N  : DEFDBL G-H,O-S,U-W 'SET MOST VARIABLES DOUBLE-PRECISION
310 ' CLEAR, &HB880 : BUFFSTART=&HB880 : TOPMEM%=CVI(CHR$(PEEK(6))+CHR$(PEEK(7)))
320 COMP = 0  ' FLAG FOR COMPILED VERSION
325 '
330 ' IF COMPILED: Change COMP flag       ! ! ! ! !
332 '              Make line 310 active
334 '              Delete lines 3005-3030, 3220-3347
335 '              Delete lines 1032-1035, 1037
337 '              Delete lines 1510, 1535, 1637-1700
338 '
340 IF COMP THEN NBLINK%=350 ELSE NBLINK%= 60 ' BLINK SPEED
350 OPTION BASE 1
360 WIDTH 255:FALSE%=0:TRUE%=NOT FALSE%
370 ' user supported software message
380 PRINT CHR$(&H1A);
390 PRINT TAB(16)"MULTREG (ver. 1.1)":PRINT
400 PRINT TAB(5)"COPYRIGHT (C) 1983,1985 BY GREGORY C. DAHL
410 IF COMP THEN PRINT TAB(6)CHR$(&H1B)CHR$(&H29)"PORTIONS COPYRIGHTED BY MICROSOFT, 1981"CHR$(&H1B)CHR$(&H28) ELSE PRINT
420 FX$="MULTREG.MSG":GOSUB 7200:IF NOT FLAG THEN PRINT:PRINT: PRINT CHR$(7)"Can't proceed without MULTREG.MSG.  Stopping.":STOP
430 OPEN "R",2,"MULTREG.MSG":FIELD#2,128 AS BUFF2$ 
440 FOR J=1 TO 7:GET#2,J:PRINT BUFF2$;:NEXT J
450 GOSUB 10000 'Press to return flashing message
460 '
470 SC% = 51 'SCREEN DISPLAY WIDTH
480 WID% = SC% : FOUT%=1
490 F1$="######.##" 'FORMAT FOR DISPLAYING DATA ON SCREEN
500 QUOTES$=CHR$(&H22) ' double "
510 '
520 ' 
530 '
1000 PRINT CHR$(&H1A):A$="MULTREG":GOSUB 12000
1010 PRINT:A$="Version 1.1":GOSUB 12000
1020 PRINT:PRINT:A$="COPYRIGHT (c) 1983,1985 BY GREGORY C. DAHL":GOSUB 12000
1030 PRINT:A$="* *   All rights reserved   * *":GOSUB 12000
1032 DOUBLE%=FALSE%:IF COMP THEN 1036
1033 PRINT:PRINT "  Single or double precision input data? (S/D)  S"CHR$(8);
1034 A$=INKEY$:IF A$="" THEN 1034 
1035 PRINT A$: IF A$="D" OR A$="d" THEN DOUBLE%=TRUE%
1036 IF COMP THEN I=(TOPMEM%-BUFFSTART-4)/4 ELSE I=INT(FRE(" ")/4-300)
1037 IF DOUBLE% THEN I=I/2 : DEFDBL X-Y
1040 PRINT:A$="[maximum memory size is about"+STR$(I)+" data points]":GOSUB 12000  
1050 FOUT%=1 '1=SCREEN,2=LPRINT,3=DISK FILE
1060 '	GET INPUT FILE NAME
1070 FOR J=8 TO 9:GET#2,J:PRINT BUFF2$;:NEXT J
1080 ' PRINT "FRE() = " FRE(" ") ' DEBUG * * * *
1090 PRINT:LINE INPUT"     ENTER DATA FILE NAME:  ";F$
1100 IF F$="" THEN FF1$="MULTREG.DTA":GOTO 1220 
1105 IF INSTR(F$,CHR$(3))>0 THEN STOP
1110 M=INSTR(F$,":"):IF M=0 THEN DRIVE$="B":GOTO 1150
1120 IF M=1 THEN 1200
1130 DRIVE$=MID$(F$,M-1,1):F$=RIGHT$(F$,LEN(F$)-M)
1140 IF DRIVE$<"A" OR DRIVE$>"D" THEN 1200
1150 M=INSTR(F$,"."):IF M=0 THEN EXT$=".DTA":GOTO 1170
1160 EXT$=MID$(F$,M,4):F$=LEFT$(F$,M-1)
1170 IF LEN(F$)=0 THEN F$="MULTREG"
1180 FF1$=DRIVE$+":"+F$+EXT$
1190 GOTO 1220
1200 PRINT:PRINT CHR$(7)"  Illegal file name " F$"  Try again.":GOSUB 10000:GOTO 1000
1210 ' 
1220 FX$=FF1$:GOSUB 7200:IF NOT FLAG THEN CLEAR:GOTO 300 'see if file exists
1230 OPEN "I",1,FF1$
1240 PRINT CHR$(&H1A);
1250 LINENUM=0:GOSUB 1730
1251 ' trim leading and trailing blanks and quotes
1252 FOR J=1 TO LEN(S$):IF LEFT$(S$,1)<>" " AND LEFT$(S$,1)<>QUOTES$ THEN 1256
1254   S$=MID$(S$,2,255):NEXT J 
1256 FOR J=1 TO LEN(S$):IF RIGHT$(S$,1)<>" " AND RIGHT$(S$,1)<>QUOTES$ THEN 1260
1258   S$=LEFT$(S$,LEN(S$)-1):NEXT J
1260 HEAD$=S$:GOSUB 1730:ALPHA=FALSE%:GOSUB 2030:V=XVAL:GOSUB 2030:N=XVAL 'get N and V from second line
1270 PRINT:PRINT:PRINT:PRINT
1280 PRINT "Header is:":PRINT:PRINT HEAD$:PRINT:PRINT
1290 PRINT "Number of independent variables is:"V
1300 PRINT:PRINT "Number of observations is:"N
1310 PRINT:PRINT:PRINT"Shall we proceed? (Y/N)  Y" CHR$(8);
1320 A$=INKEY$:IF A$="" THEN 1320
1330 PRINT A$
1340 IF A$="N" OR A$="n" THEN 6072
1350 M=V+1
1500 IF COMP THEN DIM XX#(18,18), XXINV#(18,18), A#(18,18),BETA#(18),C#(18),EST(300),ER(300),T(18),TPROB(18),VBLENAME$(19),VNAME$(19)
1510 IF NOT COMP THEN DIM XX#(M,M), XXINV#(M,M), A#(M,M), BETA#(M),C#(M),EST(N),ER(N),T(M),TPROB(M),VBLENAME$(M+1),VNAME$(M+1),X(N,M),Y(N)
1520 IF COMP THEN DEF FNY(L) = CVS(CHR$(PEEK(L*4+BUFFSTART))+CHR$(PEEK(L*4+BUFFSTART+1))+CHR$(PEEK(L*4+BUFFSTART+2))+CHR$(PEEK(L*4+BUFFSTART+3)))
1530 IF COMP THEN DEF FNX(L1,L2) = CVS(CHR$(PEEK((((L2)*N)+L1)*4+BUFFSTART))+CHR$(PEEK((((L2)*N)+L1)*4+BUFFSTART+1))+CHR$(PEEK((((L2)*N)+L1)*4+BUFFSTART+2))+CHR$(PEEK((((L2)*N)+L1)*4+BUFFSTART+3)))
1535 IF NOT COMP THEN DEF FNY(L)=Y(L) : DEF FNX(L1,L2)=X(L1,L2)
1540 ' X IS MATRIX OF DATA, ROWS ARE OBSERVATIONS
1550 '     COL 1 = 1'S, COLS 2 TO V+1 ARE INDEP. VBLES
1560 ' Y IS VECTOR OF DEPENDENT VBLE OBSERVATIONS
1570 ' (XX IS X'X, AND XXINV IS (X'X) INVERSE.  A IS A TEMPORARY WORKING MATRIX)
1580 PRINT CHR$(&H1A);:PRINT:PRINT "Reading input file:":PRINT
1590 ALPHA=TRUE%:GOSUB 1730:FOR I=1 TO V+1:GOSUB 2030:VBLENAME$(I)=A$
1600   IF A$="" AND I=1 THEN VBLENAME$(I)="DEP.VBLE"
1610   IF A$="" AND I>1 THEN VBLENAME$(I)="VBLE#"+MID$(STR$(I-1),2,LEN(STR$(I-1))-1)
1620 NEXT I
1630 ' INPUT Y,X(I,J)  
1635 ALPHA=FALSE%
1637 IF COMP THEN 1702
1640 FOR I=1 TO N
1650   GOSUB 1730 'Get input line 
1660   FOR J=0 TO V+1
1670   IF J=1 THEN X(I,J)=1:GOTO 1700' first col of X is vector of 1s
1680   GOSUB 2030
1690   IF J=0 THEN Y(I)=XVAL ELSE X(I,J)=XVAL
1700 NEXT J: NEXT I:CLOSE:GOTO 2320 ' Jump to input data display
1702 FOR I=1 TO N
1704   GOSUB 1730 'Get input line 
1706   FOR J=0 TO V+1
1708   IF J=1 THEN XV=1:LLOC=(((J)*N)+I)*4+BUFFSTART:GOSUB 20160:GOTO 1714
1710   GOSUB 2030: IF J=0 THEN LLOC=BUFFSTART+I*4:XV=XVAL:GOSUB 20160:GOTO 1714
1712   LLOC=(((J)*N)+I)*4+BUFFSTART:XV=XVAL:GOSUB 20160 'poke value in array buffer in high memory
1714 NEXT J: NEXT I:CLOSE:GOTO 2320
1720 ' READ A LINE FROM INPUT FILE:
1730 IF NOT EOF(1) THEN 2000
1740 ' End
1750 PRINT:PRINT CHR$(7)"END OF INPUT FILE REACHED":PRINT LINENUM "LINES READ":PRINT:PRINT: CLOSE:GOSUB 10000:RETURN
1760 ' 

	Input new line


2000 LINE INPUT#1, S$: PRINT S$
2005 IF INKEY$=CHR$(3) THEN STOP
2010 LINENUM=LINENUM+1:NP1=1:FLAGCOMMA=TRUE%:FLAGSTRING=FALSE%:RETURN
2020 ' 

	fetch a number XVAL from line S$ beginning at pointer NP1


2030 IF NP1> LEN(S$) THEN IF NOT ALPHA THEN GOTO 2290 ELSE A$="":RETURN
2040 LASTCHAR$=" "
2050 FOR NP=NP1 TO LEN(S$)
2060   L$=CHR$(ASC(MID$(S$,NP,1)) AND &H7F) 'strip high-order bit if any
2065 ' PRINT "L$="L$" LASTCHAR=";LASTCHAR$;" NP="NP" FLAGSTRING="FLAGSTRING" FLAGCOMMA="FLAGCOMMA ' DEBUG
2070   IF (LASTCHAR$="-" OR LASTCHAR$="+") AND L$=" " THEN 2210 ELSE LASTCHAR$=L$ ' ignore blanks after sign
2080   IF FLAGSTRING THEN 2190
2090   ' find beginning of string:
2100   IF L$=" " THEN 2210 'ignore leading blanks
2110   IF L$<>"," THEN 2140
2120     IF FLAGCOMMA THEN NP1=NP+1:NP2=NP1:GOTO 2230 ' output a zero if two commas in sequence
2130     FLAGCOMMA=TRUE% :GOTO 2210
2140   IF L$<>QUOTES$ THEN 2170
2150     NP2=INSTR(NP+1,S$,QUOTES$):IF NP2=0 THEN NP2=LEN(S$)+1'no close quote
2160     NP1=NP:GOTO 2230
2170   NP1=NP:FLAGSTRING=TRUE%:FLAGCOMMA=FALSE%:GOTO 2210 'beginning of string found
2180   ' find end of string:
2190   IF L$=" " OR L$="," THEN NP2=NP:GOTO 2230
2210 NEXT NP:NP2=LEN(S$)+1
2220 ' process string beginning at NP1 and ending at NP2-1 
2230 FLAGSTRING=FALSE%:FLAGCOMMA=FALSE%:IF L$<>QUOTES$ THEN 2250
2240   IF ALPHA THEN A$=MID$(S$,NP1+1,NP2-NP1-1):NP2=NP2+1:GOTO 2280 ELSE NP1=NP2+1:GOTO 2030 ' skip material in quotes if NOT ALPHA
2250 A$=MID$(S$,NP1,NP2-NP1):B$=LEFT$(A$,1)
2260 IF NOT ALPHA AND B$<>"-" AND B$<>"+" AND B$<>"." AND (B$<"0" OR B$>"9") AND LEN(A$)>0 THEN NP1=NP2+1:GOTO 2030 ' skip any string beginning with non-numeric
2270 XVAL = VAL(A$)
2280 NP1=NP2:RETURN 'begin at next location in S$
2290 PRINT CHR$(7)"AN ATTEMPT HAS BEEN MADE TO OBTAIN MORE NUMBERS
2300 PRINT "THAN ARE CONTAINED IN THIS LINE.  CHECK YOUR DATA 
2310 PRINT "FILE.":XVAL=0:A$="":RETURN
2320 GOSUB 15000 ' display input data
2330 '


	INVERT MATRIX X'X (ANSWER IS XXINV)


3000 PRINT:PRINT "

	INVERTING MATRIX.  PLEASE BE PATIENT....
3005 IF COMP THEN 3032
3010 FOR I=1 TO M:FOR J=1 TO M:XX#(I,J)=0
3020 FOR K=1 TO N:XX#(I,J)=XX#(I,J)+X(K,I)*X(K,J)
3030 NEXT K,J,I : GOTO 3040
3032 FOR I=1 TO M:FOR J=1 TO M:XX#(I,J)=0
3034 FOR K=1 TO N:XX#(I,J)=XX#(I,J)+FNX(K,I)*FNX(K,J)
3036 NEXT K,J,I : IF INKEY$=CHR$(3) THEN STOP
3040 FOR J=1 TO M: XXINV#(J,J)=1 : NEXT J
3050 FOR I=1 TO M:FOR J=1 TO M:A#(I,J)=XX#(I,J)::NEXT J,I 'MATRIX TO INVERT 
3060 FOR J=1 TO M:FOR I=J TO M
3065 ' PRINT "I,J,A#(I,J)=";I;J;A#(I,J) ' DEBUG * * * *
3070 IF A#(I,J)<>0 THEN 3090
3080 NEXT I : PRINT "SINGULAR MATRIX.  STOPPING":STOP
3090 FOR K=1 TO M
3100 S#=A#(J,K) : A#(J,K)=A#(I,K) : A#(I,K)=S#
3110 NEXT K
3120 T#=1/A#(J,J)
3130 FOR K=1 TO M : A#(J,K)=T#*A#(J,K) : XXINV#(J,K)=T#*XXINV#(J,K)
3140 NEXT K
3150 FOR L=1 TO M : IF L=J THEN 3210
3160 T#=-A#(L,J)
3170 FOR K=1 TO M
3180 A#(L,K)=A#(L,K)+T#*A#(J,K)
3190 XXINV#(L,K)=XXINV#(L,K)+T#*XXINV#(J,K)
3200 NEXT K
3210 NEXT L,J
3220 GOTO 3360 ' * * * * *
3230 '

	TEST INVERSE


3240 PRINT "TEST.  SHOULD BE IDENTITY MATRIX I:
3250 FOR I=1 TO M : FOR J=1 TO M : S#=0
3280 FOR K=1 TO M : S#=S#+XX#(I,K)*XXINV#(K,J) : NEXT K
3310 PRINT USING "#####.###############";S#;
3320 NEXT J:PRINT:NEXT I
3342 PRINT:PRINT"MATRIX XX#:"
3343 FOR I=1 TO M:FOR J=1 TO M:PRINT USING "#######.##### ";XX#(I,J);
3344 NEXT J:PRINT:NEXT I
3345 PRINT:PRINT"MATRIX XXINV#:"
3346 FOR I=1 TO M:FOR J=1 TO M:PRINT USING "#######.##### ";XXINV#(I,J);
3347 NEXT J:PRINT:NEXT I
3350 '

	COMPUTE ESTIMATED COEFFICIENTS BETA


3360 FOR I=1 TO M : C#(I)=0
3370 FOR J=1 TO N : C#(I)=C#(I)+FNX(J,I)*FNY(J)
3380 NEXT J,I
3390 FOR I=1 TO M : BETA#(I)=0
3400 FOR J=1 TO M : BETA#(I)=BETA#(I)+XXINV#(I,J)*C#(J)
3410 NEXT J,I
3420 '

	COMPUTE OTHER STATISTICS


3430 ' S# = Y'Y
3440 ' W# = (SUM Y)^2 / N
3450 W#=0:FOR I=1 TO N:W#=W#+FNY(I):NEXT I
3460 W#=W#*W#/N
3470 RSQR#=0:S#=0
3480 FOR I=1 TO M:RSQR#=RSQR#+BETA#(I)*C#(I):NEXT I
3490 FOR I=1 TO N:S#=S#+FNY(I)*FNY(I):NEXT I
3500 RSQR#=(RSQR#-W#)/(S#-W#)
3510 F! = RSQR# * (N-M) / ((1 - RSQR#)*(M-1))
3520 ' F DEGREES OF FREEDOM: M-1, N-M
3530 EMIN=1.7E+37 : EMAX = -EMIN 'INITIAL VALUES
3540 ERSQR# = 0
3550 FOR I=1 TO N
3560 EST(I)=0
3570 FOR J=1 TO M
3580 EST(I)=EST(I)+BETA#(J)*FNX(I,J)
3590 NEXT J
3600 ER(I)=FNY(I)-EST(I)
3610 ERSQR# = ERSQR#+ ER(I)*ER(I)
3620 IF ER(I)>EMAX THEN EMAX=ER(I)
3630 IF ER(I)<EMIN THEN EMIN=ER(I)
3640 NEXT I
3650 IF ABS(EMIN)>EMAX THEN EMAX=ABS(EMIN) ELSE EMIN=-EMAX ' MAKE ERROR LIMITS SYMMETRICAL FOR GRAPHING
3660 U# = SQR(ERSQR#/(N-M))
3670 FOR I=1 TO M: T(I) = BETA#(I)/(U# * SQR(XXINV#(I,I))) : GOSUB 21000 : NEXT I
3680 ' T HAS N-M DEGREES OF FREEDOM
3690 DW#=0: FOR I=2 TO N:DW#=DW#+(ER(I)-ER(I-1))*(ER(I)-ER(I-1)):NEXT I
3700 DW#=DW#/ERSQR# 'DURBAN-WATSON STATISTIC
3710 '

	DISPLAY RESULTS ON SCREEN


4000 PRINT CHR$(&H1A):PRINT HEAD$:PRINT
4010 PRINT "

ESTIMATED COEFFICIENTS:


4020 F3$ = "#######.###" 'FORMAT FOR t STATISTIC
4030 GOSUB 4040:GOTO 4170
4040 PRINT SPACE$(30) "t STATISTIC   t PROB"
4050 FOR I=1 TO M
4060 IF I=1 THEN PRINT "INTERCEPT:  "; ELSE PRINT LEFT$(VBLENAME$(I),10)SPACE$(10-LEN(VBLENAME$(I)))": ";
4070 IF ABS(BETA#(I))>100000! OR ABS(BETA#(I))<.001 THEN F$="#.######^^^^" ELSE F$="######.#####"
4080 PRINT USING F$;BETA#(I);
4082 A$="    ":IF TPROB(I)>.95 THEN A$="  * "
4084 IF TPROB(I)>.99 THEN A$=" ** "
4086 PRINT A$;:PRINT USING F3$;T(I);TPROB(I)
4090 NEXT I
4100 Z=RSQR#:PRINT "

R SQUARED:                "Z
4110 PRINT "

F STATISTIC:              "F!
4120 PRINT USING "(degrees of freedom: ##_,###)";M-1,N-M
4130 Z=DW#:PRINT "

DURBAN-WATSON STATISTIC:  "Z
4140 PRINT USING "

(t degrees of freedom: ###)";N-M
4150 RETURN
4160 '

	DISPLAY ERRORS


4170 MASK$="*          |          *"
4180 F$="###. #########.## #########.## #########.##"
4190 PRINT:GOSUB 10000
4200 PRINT CHR$(&H1A):A$=HEAD$:GOSUB 12000:A$="* *  ERRORS  * *":GOSUB 12000:PRINT:PRINT "OBS#     OBSERVED     ESTIMATED       ERROR":PRINT
4210 FOR I=1 TO N:PRINT USING F$;I,FNY(I),EST(I),ER(I)
4215   IF I MOD 19 = 0 THEN GOSUB 10000:PRINT CHR$(13)"                                            "CHR$(13); 'pause
4217   NEXT I
4220 PRINT:GOSUB 10000 'PRESS RETURN TO CONTINUE
4230 PRINT CHR$(&H1A):PRINT HEAD$:PRINT:PRINT "ERROR PATTERN:":PRINT
4240 M1$=STR$(CSNG(EMIN)):M2$=STR$(CSNG(EMAX)):L=LEN(M1$)+LEN(M2$):L=L/2:IF L<LEN(MASK$) THEN L=LEN(MASK$)-L ELSE L=2
4250 M1$=SPACE$(10-LEN(M1$)/2)+M1$+SPACE$(L)+M2$:PRINT M1$ 'MIN AND MAX VALUES
4260 PRINT TAB(11) STRING$(LEN(MASK$),"*")
4270 FOR I=1 TO N
4280 M$=MASK$:MID$(M$,CINT((ER(I)-EMIN)*(LEN(MASK$)-3)/(EMAX-EMIN))+2,1)="+"
4290 PRINT TAB(11) M$
4295 IF I MOD 19 = 0 THEN GOSUB 10000:PRINT CHR$(13)"                                            "CHR$(13); 'pause
4300 NEXT I
4310 PRINT TAB(11) STRING$(LEN(MASK$),"*")
4320 PRINT M1$
4330 GOSUB 10000
4340 ' 


	PRINT RESULTS TO PRINTER




5000 WID% = SC%
5010 INIT$ = "" 'INITIALIZATION STRING TO SEND TO PRINTER 
5020 F1$ = "########.#####" 'FORMAT FOR PRINTING DATA
5030 FF$=CHR$(&HC) 'FORM FEED
5040 FOUT%=1
5050 PRINT CHR$(&H1A):PRINT:PRINT:A$="YOU CAN NOW SEND DATA AND RESULTS":GOSUB 12000:A$="TO THE PRINTER":GOSUB 12000:PRINT:PRINT:A$="INPUT DATA":FOUT%=2
5060 WID%=80 'printer width in characters
5070 GOSUB 7000
5080 IF FLAG THEN 5130
5090 GOSUB 5500 'ASK FOR FORM FEED 
5100 '

	PRINT INPUT DATA


5110 GOSUB 15000
5120 '

	PRINT EST. COEFFICIENTS, ETC.


5130 A$="EST. COEFFS. AND STAT.":GOSUB 7000
5140 IF FLAG THEN 5200
5150 GOSUB 5500:IF FLAG THEN A$=HEAD$:GOSUB 12000
5160 LPRINT "



	ESTIMATED COEFFICIENTS:":LPRINT:LPRINT
5170 F3$ = "  #######.###" 'FORMAT FOR t STATISTIC
5180 GOSUB 20100:GOSUB 4040:GOSUB 20130 'print toggle, print statistics
5190 '

	PRINT ERRORS


5200 A$="ERRORS":GOSUB 7000
5210 IF FLAG THEN 5370
5220 GOSUB 5500
5230 MASK$="*          |          *"
5240 F$="###. #########.## #########.## #########.##"
5250 IF FLAG THEN A$=HEAD$:GOSUB 12000:A$=" * * *   ERRORS   * * *":GOSUB 12000
5260 M1$=STR$(CSNG(EMIN)):M2$=STR$(CSNG(EMAX)):L=LEN(M1$)+LEN(M2$):L=L/2:IF L<LEN(MASK$) THEN L=LEN(MASK$)-L ELSE L=2
5270 M1$=SPACE$(10-LEN(M1$)/2)+M1$+SPACE$(L)+M2$ 'MIN AND MAX VALUES
5280 LPRINT "


OBS#     OBSERVED     ESTIMATED       ERROR             ERROR PATTERN":LPRINT
5290 LPRINT:LPRINT SPACE$(LEN(F$)-3) M1$  
5300 LPRINT SPACE$(LEN(F$)+8) STRING$(LEN(MASK$),"*")
5310 FOR I=1 TO N
5320 M$=MASK$:MID$(M$,CINT((ER(I)-EMIN)*(LEN(MASK$)-3)/(EMAX-EMIN))+2,1)="+"
5330 LPRINT USING F$;I,FNY(I),EST(I),ER(I);:LPRINT SPACE$(8) M$
5340 NEXT I
5350 LPRINT SPACE$(LEN(F$)+8) STRING$(LEN(MASK$),"*")
5360 LPRINT SPACE$(LEN(F$)-3) M1$
5370 GOTO 6000
5380 '

	SEND FORM FEED


5500 PRINT "SEND A FORM FEED?  (Y/N)  N"CHR$(8);
5510 A$=INKEY$:IF A$="" THEN 5510
5520 FLAG=-1
5530 IF A$="Y" OR A$="y" THEN LPRINT FF$; ELSE A$="N":FLAG=0
5540 PRINT A$:RETURN
5550 ' 


	SEND RESULTS TO DISK FILE



6000 WID% = SC%
6010 INIT$ = "" 'INITIALIZATION STRING, IF ANY
6020 F1$ = "########.#####" 'FORMAT FOR DATA OUTPUT
6030 FF$="" 'FORM FEED
6040 FOUT%=1
6050 PRINT CHR$(&H1A):PRINT:PRINT:A$="YOU CAN NOW SEND DATA AND RESULTS":GOSUB 12000:A$="TO A DISK FILE":GOSUB 12000: WID%=80
6060 PRINT:PRINT:PRINT"ENTER A FILE NAME OR PRESS RETURN":PRINT:PRINT:LINE INPUT "    TO END OR REPEAT THE PROGRAM   "; F$
6070 IF F$<>"" THEN 6078 
6072 PRINT:INPUT "Run another regression? (Y/N)  ",A$
6074 IF A$<>"N" AND A$<>"n" THEN CLEAR:GOTO 300 'do again
6076 PRINT:IF COMP THEN PRINT "Free memory is: ";FRE(" "):PRINT
6077 END
6078 OPEN "O",#1,F$ : FOUT%=3 ' flag for disk out
6090 PRINT:PRINT:A$="THE INPUT DATA":GOSUB 7000
6100 IF FLAG THEN 6140
6110 '

	PRINT INPUT DATA


6120 GOSUB 15000
6130 '

	PRINT EST. COEFFICIENTS, ETC.


6140 A$="EST. COEFFS. AND STAT.":GOSUB 7000
6150 IF FLAG THEN 6310
6160 IF FLAG THEN A$=HEAD$:GOSUB 12000
6170 PRINT#1,:PRINT#1,:PRINT#1,:PRINT#1, "	ESTIMATED COEFFICIENTS:":PRINT#1,:PRINT#1,
6180 PRINT#1, SPACE$(30) "t STATISTIC   t PROB"
6190 F3$ = "#######.###" 'FORMAT FOR t STATISTIC
6200 FOR I=0 TO M-1
6210 IF I=0 THEN PRINT#1, "INTERCEPT:  "; ELSE PRINT#1, LEFT$(VBLENAME$(I+1),10)SPACE$(10-LEN(VBLENAME$(I+1)))": ";
6220 IF ABS(BETA#(I+1))>100000! OR ABS(BETA#(I+1))<.001 THEN F$="#.######^^^^" ELSE F$="######.#####"
6230 PRINT#1, USING F$;BETA#(I+1);
6232 A$="    ":IF TPROB(I+1)>.95 THEN A$="  * "
6234 IF TPROB(I+1)>.99 THEN A$=" ** "
6236 PRINT#1,A$;:PRINT#1, USING F3$;T(I+1);TPROB(I+1)
6240 NEXT I
6250 Z=RSQR#:PRINT#1, "

R SQUARED:                "Z
6260 PRINT#1, "

F STATISTIC:              "F!
6270 PRINT#1, USING "(degrees of freedom: ##_,###)";M-1,N-M
6280 Z=DW#:PRINT#1,"

DURBAN-WATSON STATISTIC:  "Z
6290 PRINT#1, USING "

(t degrees of freedom: ###)";N-M
6300 '

	PRINT ERRORS


6310 A$="ERRORS":GOSUB 7000
6320 IF FLAG THEN 6470
6330 MASK$="*          |          *"
6340 F$="###. #########.## #########.## #########.##"
6350 IF FLAG THEN A$=HEAD$:GOSUB 12000:A$=" * * *   ERRORS   * * *":GOSUB 12000
6360 M1$=STR$(CSNG(EMIN)):M2$=STR$(CSNG(EMAX)):L=LEN(M1$)+LEN(M2$):L=L/2:IF L<LEN(MASK$) THEN L=LEN(MASK$)-L ELSE L=2
6370 M1$=SPACE$(10-LEN(M1$)/2)+M1$+SPACE$(L)+M2$ 'MIN AND MAX VALUES
6380 PRINT#1,:PRINT#1,:PRINT#1, "OBS#     OBSERVED     ESTIMATED       ERROR             ERROR PATTERN":PRINT#1,
6390 PRINT#1,:PRINT#1, SPACE$(LEN(F$)-3) M1$  
6400 PRINT#1, SPACE$(LEN(F$)+8) STRING$(LEN(MASK$),"*")
6410 FOR I=1 TO N
6420 M$=MASK$:MID$(M$,CINT((ER(I)-EMIN)*(LEN(MASK$)-3)/(EMAX-EMIN))+2,1)="+"
6430 PRINT#1, USING F$;I,FNY(I),EST(I),ER(I);:PRINT#1, SPACE$(8) M$
6440 NEXT I
6450 PRINT#1, SPACE$(LEN(F$)+8) STRING$(LEN(MASK$),"*")
6460 PRINT#1, SPACE$(LEN(F$)-3) M1$
6470 PRINT:PRINT:A$="* * *  END OF PROGRAM  * * *":GOSUB 12000:PRINT:PRINT:CLOSE:STOP 'END OF PROGRAM
6480 '

	ASK FOR OUTPUT


7000 PRINT:IF FOUT%=2 THEN PRINT "PRINT "; ELSE PRINT "OUTPUT ";
7010 PRINT A$"?  (Y/N)  N"CHR$(8);
7020 A$=INKEY$:IF A$="" THEN 7020
7030 IF A$="Y" OR A$="y" THEN FLAG=FALSE% ELSE FLAG=TRUE%:A$="N"
7040 IF FLAG AND NOT FLAGINIT AND FOUT%=2 THEN LPRINT INIT$;:FLAGINIT=TRUE% 'PRINT INITIALIZATION STRING ONCE
7050 PRINT A$:RETURN 
7070 ' 


	 CHECK IF INPUT FILE EXISTS



7200 OPEN "R",1,FX$:GET#1:FLAG=TRUE%
7210 IF NOT EOF(1) THEN 7250
7220 CLOSE:PRINT:PRINT CHR$(7)"The input file "FX$" does not exist"
7230 PRINT"on the referenced drive.  Try again.":PRINT
7240 KILL FX$:GOSUB 10000:FLAG=FALSE% 
7250 CLOSE:RETURN
7260 ' 

	pause routine


10000 PRINT"           Press any key to continue...";
10010 FOR IB=1 TO NBLINK%:A$=INKEY$:IF A$=CHR$(3) THEN STOP
10015   IF A$<>""THEN RETURN ELSE NEXT IB
10020 PRINT CHR$(13)"                                       "CHR$(13);
10030 FOR IB=1 TO NBLINK%:A$=INKEY$:IF A$=CHR$(3) THEN STOP
10035   IF A$<>""THEN RETURN ELSE NEXT IB
10040 GOTO 10000
10990 ' 

	pause routine with exit


11000 PRINT"     Press any key to continue, X to exit...";
11010 FOR IB=1 TO NBLINK%:A$=INKEY$:IF A$<>""THEN GOTO 11050 ELSE NEXT IB
11020 PRINT CHR$(13)"                                            "CHR$(13);
11030 FOR IB=1 TO NBLINK%:A$=INKEY$:IF A$<>""THEN GOTO 11050 ELSE NEXT IB
11040 GOTO 11000
11050 IF A$="X" OR A$="x" THEN CLEAR:GOTO 300 ELSE RETURN
11150 '

	DISPLAY LINE, CENTERED


12000 IF FOUT%=1 THEN PRINT TAB((WID%-LEN(A$))/2) A$:RETURN
12010 IF FOUT%=2 THEN LPRINT TAB((WID%-LEN(A$))/2) A$:RETURN
12020 IF FOUT%=3 THEN PRINT#1, TAB((WID%-LEN(A$))/2) A$
12030 RETURN
12040 '

	DISPLAY INPUT DATA


15000 IF FOUT%=1 THEN PRINT CHR$(&H1A);
15010 A$=HEAD$:GOSUB 12000
15020 HEAD1$="* * * *  INPUT DATA  * * * *":A$=HEAD1$:GOSUB 12000
15030 FOR J=1 TO V+1 ' set up vble name headers
15040 VNAME$(J)=VBLENAME$(J):IF LEN(VNAME$(J))>LEN(F1$) THEN VNAME$(J)=LEFT$(VNAME$(J),LEN(F1$))
15050 I=(LEN(F1$)-LEN(VNAME$(J))-.4)/2:IF I<0 THEN I=0'spaces each side of name
15060 VNAME$(J)=SPACE$(I)+VNAME$(J)+SPACE$(I):IF LEN(VNAME$(J))<LEN(F1$) THEN VNAME$(J)=VNAME$(J)+" " ' take care of rounding error
15070 NEXT J
15080 I=(LEN(F1$)-8)/2:IF I<0 THEN I=0
15090 GOSUB 20050:GOSUB 20050:A$="OBS#  "+SPACE$(I)+"DEP. VBLE"+SPACE$(I+2)+"INDEPENDENT VARIABLES":GOSUB 20000
15100 GOSUB 20050:GOSUB 20050
15110 A$="       "+VNAME$(1)+"  ":GOSUB 20000
15120 F$="###.  "+F1$+" "
15130 K=INT((WID%-7)/LEN(F1$)): IF K>M THEN K=M 'NUMBER OF INDEP. VARIABLES THAT CAN BE DISPLAYED ON SCREEN
15140 FOR J=1 TO K-1:A$=VNAME$(J+1):GOSUB 20000:NEXT J
15150 IF FOUT%<>1 THEN 15310 
15160 PRINT
15170 FOR I=1 TO N
15180   PRINT USING F$;I,FNY(I);
15190   FOR J=2 TO K : PRINT USING F1$;FNX(I,J);
15200   NEXT J:PRINT
15205   IF I MOD 19 = 0 THEN GOSUB 11000:PRINT CHR$(&H1A) 'pause
15208   NEXT I:GOSUB 11000
15210 IF K=M THEN RETURN 'ALL INPUT DATA DISPLAYED
15220 K2=K
15230 K1=K2+1:K2=K1+(INT((WID%-6)/LEN(F1$)))-1:IF K2>M THEN K2=M
15240 PRINT CHR$(&H1A):A$=HEAD$:GOSUB 12000:A$=HEAD1$:GOSUB 12000:PRINT:PRINT"INDEP. VARIABLES, CONTINUED:":PRINT
15250 PRINT TAB(11);:FOR J=K1-1 TO K2-1:PRINT VNAME$(J+1);:NEXT J: PRINT
15260 FOR I=1 TO N:PRINT USING "###.    ";I;
15270   FOR J=K1 TO K2:PRINT USING F1$;FNX(I,J);
15280   NEXT J:PRINT
15285   IF I MOD 19 = 0 THEN GOSUB 11000:PRINT CHR$(&H1A) 'pause 
15287   NEXT I:GOSUB 11000
15290 IF K2<M THEN 15230
15300 RETURN
15310 IF FOUT%<>2 THEN 15470 
15320 LPRINT
15330 FOR I=1 TO N
15340   LPRINT USING F$;I,FNY(I);
15350   FOR J=2 TO K : LPRINT USING F1$;FNX(I,J);
15360   NEXT J:LPRINT:NEXT I:GOSUB 10000
15370 IF K=M THEN RETURN 'ALL INPUT DATA DISPLAYED
15380 K2=K
15390 K1=K2+1:K2=K1+(INT((WID%-6)/LEN(F1$)))-1:IF K2>M THEN K2=M
15400 LPRINT CHR$(&H1A):A$=HEAD$:GOSUB 12000:A$=HEAD1$:GOSUB 12000:LPRINT:LPRINT"INDEP. VARIABLES, CONTINUED:":LPRINT
15410 LPRINT TAB(11);:FOR J=K1-1 TO K2-1:LPRINT VNAME$(J+1);:NEXT J: LPRINT
15420 FOR I=1 TO N:LPRINT USING "###.      ";I;
15430 FOR J=K1 TO K2:LPRINT USING F1$;FNX(I,J);
15440 NEXT J:LPRINT:NEXT I:GOSUB 10000
15450 IF K2<M THEN 15230
15460 RETURN
15470 IF FOUT%<>3 THEN RETURN
15480 PRINT#1,
15490 FOR I=1 TO N
15500   PRINT#1, USING F$;I,FNY(I);
15510   FOR J=2 TO K : PRINT#1, USING F1$;FNX(I,J);
15520   NEXT J:PRINT#1,:NEXT I:GOSUB 10000
15530 IF K=M THEN RETURN 'ALL INPUT DATA DISPLAYED
15540 K2=K
15550 K1=K2+1:K2=K1+(INT((WID%-6)/LEN(F1$)))-1:IF K2>M THEN K2=M
15560 PRINT#1, CHR$(&H1A):A$=HEAD$:GOSUB 12000:A$=HEAD1$:GOSUB 12000:PRINT#1,:PRINT#1,"INDEP. VARIABLES, CONTINUED:":PRINT#1,
15570 PRINT#1, TAB(11);:FOR J=K1-1 TO K2-1:PRINT#1, VNAME$(J+1);:NEXT J: PRINT#1,
15580 FOR I=1 TO N:PRINT#1, USING "###.    ";I;
15590 FOR J=K1 TO K2:PRINT#1, USING F1$;FNX(I,J);
15600 NEXT J:PRINT#1,:NEXT I
15610 IF K2<M THEN 15230
15620 RETURN
15630 ' 

	send A$ to output device


20000 IF FOUT%=1 THEN PRINT A$;:RETURN
20010 IF FOUT%=2 THEN LPRINT A$;:RETURN
20020 IF FOUT%=3 THEN PRINT#1, A$;
20030 RETURN
20040 ' 

	send CR to output device


20050 IF FOUT%=1 THEN PRINT:RETURN
20060 IF FOUT%=2 THEN LPRINT:RETURN
20070 IF FOUT%=3 THEN PRINT#1,
20080 RETURN
20090 ' 

	print toggle on


20100 POKE(PEEK(2)*256+PEEK(1)+9),205
20110 RETURN
20120 ' 

	print toggle off


20130 POKE(PEEK(2)*256+PEEK(1)+9),195
20140 RETURN
20150 ' 

	poke single precision number XV into memory location LLOC


20160 POKE LLOC,ASC(MID$(MKS$(XV),1,1)):POKE LLOC+1,ASC(MID$(MKS$(XV),2,1)):POKE LLOC+2,ASC(MID$(MKS$(XV),3,1)):POKE LLOC+3,ASC(MID$(MKS$(XV),4,1))
20170 RETURN
20980 '


	COMPUTE TPROB      (4/23/85)



20990  ' COMPUTE T PROB FOR VALUE ABS(T(I)) AND N-M D.F.
21000 IF FLAGT THEN 21170
21005 INC! = .1 'step increment for integrals
21010 ' N-M = degrees of freedom
21030 '

	find coefficient:


21040 FLAGT=TRUE%:MM2! = (N-M-1)/2 
21050 ' scale factor to avoid overflow error:
21060 SCALE% = INT((MM2!+1)/10):SCALED!=MM2!
21070 GOSUB 21300
21080 ' PRINT "FACTORIAL OF ";MM2!;" = ";TEMP!;" SCALE = ";SCALE%
21090 COEFF! = (1/SQR((N-M)*3.14159))*TEMP!
21100 MM2!=(N-M-2)/2 : GOSUB 21300 : COEFF! = COEFF!/TEMP!
21110 ' PRINT "FACTORIAL OF ";MM2!;" = ";TEMP!
21120 ' PRINT "COEFF! = ";COEFF!
21130 ' RETURN ' DEBUG ****
21140 ' 

	find area under curve for XZ = 0 to ABS(T(I))


21150 '   (use area of trapezoids = average of two 
21160 '    sampling points * width/increment)
21170 MM=N-M : MM2!=-(N-M+1)/2
21180 TPROB(I) = (1^MM2!)/2 'opening value * 1/2, probability 0 to ABS(T(I))
21190 JTEST = ABS(T(I))/INC! ' integer value for limit T(I) in terms of increment
21200 FOR J = 1 TO 10000
21210   XZ!= J * INC!
21220   TPROB(I) = TPROB(I) + ((XZ!*XZ!/MM)+1)^MM2!
21230 '  IF J MOD 100 = 0 THEN PRINT XZ!, ((XZ!*XZ!/MM)+1)^MM2!, TPROB(I)
21240   IF J > JTEST THEN 21270
21250 NEXT J
21260 ' adjust for part-trapezoid at end:
21270 TPROB(I) = TPROB(I) + ((((ABS(T(I))-XZ!)/INC!)-1)/2)*((XZ!*XZ!/MM)+1)^MM2! + (((T(I)*T(I)/MM)+1)^MM2!)*(ABS(T(I))-XZ!)/(2*INC!)
21280 TPROB(I) = TPROB(I) * 2 * COEFF! * INC!: RETURN
21290 ' find (MM2!)!
21300 SCALE2%=SCALE%:IF MM2!<>INT(MM2!) THEN 21380 ' use gamma function
21310  TEMP!=1
21320   FOR J=2 TO MM2!:TEMP!=TEMP!*J
21330     IF TEMP!>SCALED!^10 THEN TEMP!=TEMP!/SCALED!^10:SCALE2%=SCALE2%-1
21340 '   PRINT "TEMP! =";TEMP! ' DEBUG *****
21350  NEXT J:TEMP!=TEMP!/SCALED!^(10*SCALE2%)
21360   RETURN
21370 'IF MM2!< 2 THEN 1070 'use integral to find gamma(MM2!)
21380 MM2!=(MM2!)+1
21390 TEMP!=1:FOR J=1 TO SCALE%:TEMP!=TEMP!*((MM2!^10)/SCALED!^10):NEXT J
21400 TEMP!=TEMP!*MM2!^(MM2!-10*SCALE%)
21410 TEMP! = TEMP!*EXP(-MM2!)*SQR(2*3.14159/MM2!)*(1+1/(12*MM2!)+1/(288*MM2!*MM2!)-139/(51840!*MM2!*MM2!*MM2!)-571/(2.48832E+06*MM2!*MM2!*MM2!*MM2!)) 'fromula from Stan. Math Tables p.316
21420 ' PRINT "SCALE 2 = " MM2!^-SCALE!
21430 RETURN
.48832E+06*MM2!*MM2!*MM2!*MM2!)) 'fromula from Stan. Math Tables p.316
21420 ' PRINT "SCALE 2 = " MM2!^-SCALE