A Flexible Least Squares Program for
Time-Varying Linear Regression

Last Updated: 26 January 2024

Site Maintained By:
  Leigh Tesfatsion
  Prof. Emerita of Economics
  Courtesy Research Professor
   of Electrical & Comp. Eng.
  Heady Hall 260
  Iowa State University
  Ames, Iowa 50011-1054
  tesfatsi AT iastate.edu

FLS Homepage
FLS Intro Article
FLS-TVLR Article
GFLS Intro Article
FLS Residual Efficiency Frontier Figure

FLS-TVLR Program: Overview

The FLS-TVLR Program is an open-source Fortran program that implements the Flexible Least Squares (FLS) approach to Time-Varying Linear Regression (TVLR) developed in the following article:

Robert E. Kalaba and Leigh Tesfatsion, "Time-Varying Linear Regression Via Flexible Least Squares" (WPPreprint,pdf,2.4MB), Computers and Mathematics with Applications Vol. 17 (1989), pp. 1215-1245.

Note: The published version of this article is available from Science Direct.

In analogy to the determination of Pareto efficiency frontiers for multicriteria decision problems, the basic FLS-TVLR objective for a given TVLR model, conditional on a given time-series data set T, is to determine REF(T), the Residual Efficiency Frontier for T. As depicted at the top of this website, REF(T) is the collection of all estimated sequences for the TVLR model coefficients that yield vector-minimal sums of squared residual measurement and dynamic errors, conditional on T.

The FLS-TVLR Program has been incorporated into the statistical packages SHAZAM and GAUSS - TSM (Time Series Methods).

See the FLS Homepage for more extensive conceptual and code-availability information for FLS-TVLR as well as for FLS approaches developed for more general model formulations.

FLS-TVLR Program: Source Code (Fortran):

Software Release Disclaimer:

The software below is provided as-is, without warranty of any kind. It is released here as freeware by Robert E. Kalaba and Leigh Tesfatsion (the copyright holders) under the terms of the Artistic License Agreement.
 

//*                                                                     00000010
// EXEC FORT1CLG,REGION=512K                                            00000020
//*                                                                     00000030
/*JOBPARM COPIES=6                                                      00000040
C                                                                       00000050
C      PROGRAM FLS:                                                     00000060
C         EQUATION NUMBERS IN COMMENT STATEMENTS REFER TO               00000070
C         ''TIME VARYING LINEAR REGRESSION VIA FLEXIBLE LEAST SQUARES'' 00000080
C         BY R. KALABA AND L. TESFATSION, COMPUTERS AND MATHEMATICS     00000090
C         WITH APPLICATIONS 17, (1990), PAGES 1215-1245.                00000100
C      LAST UPDATED: 14 JUNE 1994                                       00000110
C                                                                       00000120
C      MAIN PROGRAM                                                     00000130
C                                                                       00000140
       IMPLICIT REAL*8(A-H,O-Z)                                         00000150
       REAL*8M                                                          00000160
C      THIS PROGRAM IS CURRENTLY DIMENSIONED FOR A MAXIMUM OF           00000170
C      110 OBSERVATIONS, WITH REGRESSOR VECTORS OF A MAXIMUM            00000180
C      DIMENSION 10.                                                    00000190
       DIMENSION X(10,110),Y(110),R(10,10)                              00000200
       DIMENSION XN(10),Z(10,10),M(10,10,110),V(10)                     00000210
       DIMENSION U(10),E(10,110),W(10),Q(10,10),A(10,10)                00000220
       DIMENSION C(10,10),D(10),B(10,110),DIFVEC(10)                    00000230
       DIMENSION VVV(10),ZZ(10,10),ZZINV(10,10),BOLSE(10)               00000240
       DIMENSION BOLS(10),XY(10),XXT(10,10),XXTINV(10,10)               00000250
       DIMENSION XTBOLS(110),OLSR(110)                                  00000260
       DIMENSION TRUEB(10,110)                                          00000270
C       THE FOLLOWING SUBROUTINE INITIALIZES THE PENALTY WEIGHT         00000280
C       AMU, THE DIMENSION K OF THE REGRESSOR VECTORS, AND              00000290
C       THE NUMBER OF OBSERVATIONS NCAP.  IT ALSO GENERATES             00000300
C       THE K BY NCAP MATRIX X OF REGRESSOR VALUES AND THE              00000310
C       NCAP BY 1 VECTOR Y OF SCALAR OBSERVATIONS.                      00000320
        CALL INPUT(AMU,K,NCAP,X,Y,TRUEB)                                00000330
C       INITIALIZATION FOR THE AUXILIARY MATRIX RN = QN-1 + AMU*I       00000340
C       AT TIME N = 1                                                   00000350
        DO 30 I=1,K                                                     00000360
        DO 40 J=1,K                                                     00000370
        R(I,J) = 0.0D+00                                                00000380
        IF(I.EQ.J) R(I,J) = AMU                                         00000390
   40   CONTINUE                                                        00000400
   30   CONTINUE                                                        00000410
        DO 50 N=1,NCAP                                                  00000420
C       FORM THE REGRESSOR VECTOR XN AND THE SCALAR                     00000430
C       OBSERVATION YN                                                  00000440
        DO 60 I=1,K                                                     00000450
        XN(I)=X(I,N)                                                    00000460
   60   CONTINUE                                                        00000470
        YN=Y(N)                                                         00000480
C       CALCULATE THE INVERSE Z OF THE MATRIX RN+XNXNT VIA              00000490
C       THE WOODBURY MATRIX INVERSION LEMMA                             00000500
        CALL WOOD(K,R,XN,Z)                                             00000510
C       CALCULATE & STORE THE K BY K MATRICES MN AND THE K BY 1         00000520
C       VECTORS EN APPEARING IN EQUATIONS (5.7B) AND (5.7C)             00000530
        DO 70 I=1,K                                                     00000540
        DO 80 J=1,K                                                     00000550
        M(I,J,N)=AMU*Z(I,J)                                             00000560
   80   CONTINUE                                                        00000570
   70   CONTINUE                                                        00000580
        DO 90 I=1,K                                                     00000590
        V(I)=XN(I)*YN                                                   00000600
   90   CONTINUE                                                        00000610
        DO 100 I=1,K                                                    00000620
        IF(N.EQ.1) U(I)=0.0D+00                                         00000630
        IF(N.GT.1) U(I)=AMU*E(I,N-1)                                    00000640
  100   CONTINUE                                                        00000650
        DO 110 I=1,K                                                    00000660
        W(I)=U(I)+V(I)                                                  00000670
  110   CONTINUE                                                        00000680
        DO 120 I=1,K                                                    00000690
        E(I,N)=0.0D+00                                                  00000700
        DO 130 J=1,K                                                    00000710
        E(I,N)=E(I,N)+Z(I,J)*W(J)                                       00000720
  130   CONTINUE                                                        00000730
  120   CONTINUE                                                        00000740
        DO 140 I=1,K                                                    00000750
        DO 150 J=1,K                                                    00000760
        R(I,J)=-AMU*AMU*Z(I,J)                                          00000770
        IF(I.EQ.J)R(I,J)=(2.0D+00)*AMU+R(I,J)                           00000780
  150   CONTINUE                                                        00000790
  140   CONTINUE                                                        00000800
   50   CONTINUE                                                        00000810
C       CALCULATE THE (5.15) FLS ESTIMATE BN FOR THE                    00000820
C       FINAL TIME N = NCAP                                             00000830
        DO 160 I=1,K                                                    00000840
        DO 170 J=1,K                                                    00000850
        Q(I,J)=-AMU*M(I,J,NCAP-1)                                       00000860
        IF(I.EQ.J) Q(I,J)=Q(I,J)+AMU                                    00000870
  170   CONTINUE                                                        00000880
  160   CONTINUE                                                        00000890
C       OBTAIN THE INVERSE C OF THE MATRIX A=(QN-1 + XNXNT)             00000900
C       IN EQUATION (5.15)                                              00000910
        DO 180 I=1,K                                                    00000920
        DO 190 J=1,K                                                    00000930
        A(I,J)=Q(I,J)+X(I,NCAP)*X(J,NCAP)                               00000940
  190   CONTINUE                                                        00000950
  180   CONTINUE                                                        00000960
        CALL INV(K,A,C)                                                 00000970
C       FORM THE VECTOR D=(PN-1 + XNYN) IN EQUATION (5.15)              00000980
        DO 200 I=1,K                                                    00000990
        D(I)=AMU*E(I,NCAP-1)+X(I,NCAP)*Y(NCAP)                          00001000
  200   CONTINUE                                                        00001010
C       POSTMULTIPLY C BY D TO OBTAIN BNCAP                             00001020
        DO 210 I=1,K                                                    00001030
        B(I,NCAP)=0.0D+00                                               00001040
        DO 220 J=1,K                                                    00001050
        B(I,NCAP)=B(I,NCAP)+C(I,J)*D(J)                                 00001060
  220   CONTINUE                                                        00001070
  210   CONTINUE                                                        00001080
C       USE EQUATIONS (5.16) TO OBTAIN SMOOTHED FLS ESTIMATES           00001090
C       FOR B1,...,BNCAP-1                                              00001100
        NCAP1=NCAP-1                                                    00001110
        DO 230 N=1,NCAP1                                                00001120
        L=NCAP-N                                                        00001130
        DO 240 I=1,K                                                    00001140
        B(I,L)=E(I,L)                                                   00001150
        DO 250 J=1,K                                                    00001160
        B(I,L)=B(I,L)+M(I,J,L)*B(J,L+1)                                 00001170
  250   CONTINUE                                                        00001180
  240   CONTINUE                                                        00001190
  230   CONTINUE                                                        00001200
        WRITE(6,2020)                                                   00001210
 2020   FORMAT(1X,'HERE ARE THE FLS ESTIMATES FOR B1 AND THE            00001220
     &  TRUE B1')                                                       00001230
        WRITE(6,37)(B(1,N),TRUEB(1,N),N=1,NCAP)                         00001240
        WRITE(6,2030)                                                   00001250
 2030   FORMAT(1X,'HERE ARE THE FLS ESTIMATES FOR B2 AND THE            00001260
     &  TRUE B2')                                                       00001270
        WRITE(6,37)(B(2,N),TRUEB(2,N),N=1,NCAP)                         00001280
   37   FORMAT(1X,2D25.10)                                              00001290
C       CALCULATING RSUBM FROM EQUATION (3.1)                           00001300
        SUM =0.0D+00                                                    00001310
        DO 500 N=1,NCAP                                                 00001320
        SUM1=0.0D+00                                                    00001330
        DO 510 J=1,K                                                    00001340
        SUM1=SUM1+X(J,N)*B(J,N)                                         00001350
  510   CONTINUE                                                        00001360
        DIF=Y(N)-SUM1                                                   00001370
        DIFSQ=DIF*DIF                                                   00001380
        SUM=SUM+DIFSQ                                                   00001390
  500   CONTINUE                                                        00001400
        RSUBM=SUM                                                       00001410
C       CALCULATING RSUBD FROM EQUATION (3.2)                           00001420
        SUM=0.0D+00                                                     00001430
        DO 520 N=1,NCAP1                                                00001440
        DO 530 J=1,K                                                    00001450
        DIFVEC(J)=B(J,N+1)-B(J,N)                                       00001460
  530   CONTINUE                                                        00001470
        SUM1=0.0D+00                                                    00001480
        DO 540 J=1,K                                                    00001490
        SUM1=SUM1+DIFVEC(J)*DIFVEC(J)                                   00001500
  540   CONTINUE                                                        00001510
        SUM=SUM+SUM1                                                    00001520
  520   CONTINUE                                                        00001530
        RSUBD=SUM                                                       00001540
C       CALCULATING THE INCOMPATIBILITY COST AMU*RSUBD + RSUBM          00001550
        COST=AMU*RSUBD+RSUBM                                            00001560
        WRITE(6,580) RSUBM,RSUBD,COST                                   00001570
  580   FORMAT(1H0,'HERE ARE RSUBM,RSUBD,COST'/1X,3D20.10)              00001580
C       FIRST VALIDATION CHECK:                                         00001590
C       CALCULATION OF THE ESTIMATE BOLSE FOR OLS FROM THE MATRIX       00001600
C       AVERAGE OF THE FLS ESTIMATES GIVEN BY EQUATION (6.2)            00001610
        DO 810 I=1,K                                                    00001620
        VVV(I)=0.0D+00                                                  00001630
        DO 820 N=1,NCAP                                                 00001640
        SUM1=0.0D+00                                                    00001650
        DO 830 J=1,K                                                    00001660
  830   SUM1=SUM1+X(J,N)*B(J,N)                                         00001670
  820   VVV(I)=VVV(I)+X(I,N)*SUM1                                       00001680
  810   CONTINUE                                                        00001690
        DO 840  J=1,K                                                   00001700
        DO 850 I=1,K                                                    00001710
  850   ZZ(I,J)=0.0D+00                                                 00001720
  840   CONTINUE                                                        00001730
        DO 860 N=1,NCAP                                                 00001740
        DO 870 I=1,K                                                    00001750
        DO 880 J=1,K                                                    00001760
  880   ZZ(I,J)=X(I,N)*X(J,N)+ZZ(I,J)                                   00001770
  870   CONTINUE                                                        00001780
  860   CONTINUE                                                        00001790
        CALL INV(K,ZZ,ZZINV)                                            00001800
        DO 890 I=1,K                                                    00001810
        BOLSE(I)=0.0D+00                                                00001820
        DO 900 J=1,K                                                    00001830
  900   BOLSE(I)=BOLSE(I)+ZZINV(I,J)*VVV(J)                             00001840
  890   CONTINUE                                                        00001850
        WRITE(6,910)                                                    00001860
  910   FORMAT(1H0,'COMPONENTS OF BOLSE')                               00001870
        WRITE(6,920) (BOLSE(I),I=1,K)                                   00001880
  920   FORMAT(1X,D35.10)                                               00001890
C       CALCULATING THE OLS ESTIMATE BOLS FROM THE USUAL                00001900
C       FORMULA BOLS=(XXT)-1*XY                                         00001910
        DO 1600 I=1,K                                                   00001920
        SUM=0.0D+00                                                     00001930
        DO 1610 N=1,NCAP                                                00001940
 1610   SUM=SUM+X(I,N)*Y(N)                                             00001950
 1600   XY(I)=SUM                                                       00001960
        DO 1620 I=1,K                                                   00001970
        DO 1630 J=1,K                                                   00001980
        XXT(I,J)=0.0D+00                                                00001990
        DO 1640 N=1,NCAP                                                00002000
 1640   XXT(I,J)=XXT(I,J)+X(I,N)*X(J,N)                                 00002010
 1630   CONTINUE                                                        00002020
 1620   CONTINUE                                                        00002030
        CALL INV(K,XXT,XXTINV)                                          00002040
        DO 1650 I=1,K                                                   00002050
        BOLS(I)=0.0D+00                                                 00002060
        DO 1660 J=1,K                                                   00002070
 1660   BOLS(I)=BOLS(I)+XXTINV(I,J)*XY(J)                               00002080
 1650   CONTINUE                                                        00002090
        WRITE(6,1670)                                                   00002100
 1670   FORMAT(1H0,'COMPONENTS OF BOLS')                                00002110
        WRITE(6,1671) (BOLS(I),I=1,K)                                   00002120
 1671   FORMAT(1X,D35.10)                                               00002130
C       CALCULATING THE SUM OF SQUARED RESIDUAL MEASUREMENT             00002140
C       ERRORS OLSRM FOR OLS                                            00002150
        DO 1800 N=1,NCAP                                                00002160
        XTBOLS(N)=0.0D+00                                               00002170
        DO 1810 I=1,K                                                   00002180
 1810   XTBOLS(N)=XTBOLS(N)+X(I,N)*BOLS(I)                              00002190
 1800   CONTINUE                                                        00002200
        DO 1820 N=1,NCAP                                                00002210
 1820   OLSR(N)=Y(N)-XTBOLS(N)                                          00002220
        OLSRM=0.0D+00                                                   00002230
        DO 1830 N=1,NCAP                                                00002240
 1830   OLSRM=OLSRM+OLSR(N)*OLSR(N)                                     00002250
        WRITE(6,1840)                                                   00002260
 1840   FORMAT(1H0,'SUM OF SQUARED RESIDUAL MEASUREMENT ERRORS          00002270
     &  FOR OLS')                                                       00002280
        WRITE(6,1841) OLSRM                                             00002290
 1841   FORMAT(1X,D27.10)                                               00002300
C       SECOND VALIDATION CHECK: FIRST-ORDER CONDITION TEST             00002310
C       HOW WELL DO THE FLS ESTS SATISFY THE FOC CONDITIONS (A.14)      00002320
        CALL FOCTST(AMU,K,NCAP,X,Y,B)                                   00002330
        STOP                                                            00002340
        END                                                             00002350
C                                                                       00002360
        SUBROUTINE WOOD(K,R,X,Z)                                        00002370
C       CALCULATES THE INVERSE Z OF A MATRIX OF THE FORM R+XXT          00002380
C       VIA THE WOODBURY MATRIX INVERSION LEMMA                         00002390
        IMPLICIT REAL*8(A-H,O-Z)                                        00002400
        DIMENSION R(10,10),X(10),Z(10,10),S(10,10),V(10)                00002410
        DIMENSION XNUM(10,10),U(10,10)                                  00002420
C       CALCULATE THE INVERSE S OF THE K BY K MATRIX R                  00002430
        CALL INV(K,R,S)                                                 00002440
C       CALCULATE THE K BY 1 VECTOR V=SX=R-1*X                          00002450
        DO 10 I=1,K                                                     00002460
        V(I)=0.0D+00                                                    00002470
        DO 20 J=1,K                                                     00002480
        V(I)=V(I)+S(I,J)*X(J)                                           00002490
   20   CONTINUE                                                        00002500
   10   CONTINUE                                                        00002510
C       CALCULATE XNUM=VVT=R-1*XXT*R-1                                  00002520
        DO 30 I=1,K                                                     00002530
        DO 40 J=1,K                                                     00002540
        XNUM(I,J)=V(I)*V(J)                                             00002550
   40   CONTINUE                                                        00002560
   30   CONTINUE                                                        00002570
C       CALCULATE Y=(1+VTV)=(1+XT*R-1*X)                                00002580
        Y=1.0D+00                                                       00002590
        DO 50 I=1,K                                                     00002600
        Y=Y+X(I)*V(I)                                                   00002610
   50   CONTINUE                                                        00002620
C       CALCULATE U=XNUM/Y                                              00002630
        DO 60 I=1,K                                                     00002640
        DO 70 J=1,K                                                     00002650
        U(I,J)=XNUM(I,J)/Y                                              00002660
   70   CONTINUE                                                        00002670
   60   CONTINUE                                                        00002680
C       CALCULATE Z=S-U                                                 00002690
        DO 80 I=1,K                                                     00002700
        DO 90 J=1,K                                                     00002710
        Z(I,J)=S(I,J)-U(I,J)                                            00002720
   90   CONTINUE                                                        00002730
   80   CONTINUE                                                        00002740
        RETURN                                                          00002750
        END                                                             00002760
C                                                                       00002770
        SUBROUTINE INV(K,A,C)                                           00002780
C       CALCULATES THE INVERSE C OF A K BY K MATRIX A                   00002790
        IMPLICIT REAL*8(A-H,O-Z)                                        00002800
        DIMENSION A(10,10),C(10,10),B(10,20)                            00002810
        DO 5 J=1,K                                                      00002820
        DO 6 I=1,K                                                      00002830
    6   B(I,J)=A(I,J)                                                   00002840
    5   CONTINUE                                                        00002850
        K2=K*2                                                          00002860
        DO 7 J=1,K                                                      00002870
        DO 8 I=1,K                                                      00002880
        B(I,K+J)=0.0D+00                                                00002890
        IF(I.EQ.J) B(I,K+J)=1.0D+00                                     00002900
    8   CONTINUE                                                        00002910
    7   CONTINUE                                                        00002920
C       THE PIVOT OPERATION STARTS HERE                                 00002930
        DO 9 L=1,K                                                      00002940
        PIVOT=B(L,L)                                                    00002950
        DO 13 J=L,K2                                                    00002960
   13   B(L,J)=B(L,J)/PIVOT                                             00002970
C       TO IMPROVE THE ROWS                                             00002980
        DO 14 I=1,K                                                     00002990
        IF(I.EQ.L) GO TO 14                                             00003000
        AIL=B(I,L)                                                      00003010
        DO 15 J=L,K2                                                    00003020
   15   B(I,J)=B(I,J)-AIL*B(L,J)                                        00003030
   14   CONTINUE                                                        00003040
    9   CONTINUE                                                        00003050
        DO 45 I=1,K                                                     00003060
        DO 46 J=1,K                                                     00003070
   46   C(I,J)=B(I,K+J)                                                 00003080
   45   CONTINUE                                                        00003090
        RETURN                                                          00003100
        END                                                             00003110
C                                                                       00003120
        SUBROUTINE INPUT(AMU,K,NCAP,X,Y,TRUEB)                          00003130
        IMPLICIT REAL*8(A-H,O-Z)                                        00003140
        DIMENSION X(10,110),Y(110),TRUEB(10,110)                        00003150
C                                                                       00003160
C       RUN FOR ELLIPTICAL TRUE B WITH NOISY OBSERVATIONS               00003170
C                                                                       00003180
        K=2                                                             00003190
        AMU=1.0D+00                                                     00003200
        NCAP=30                                                         00003210
        SIGMA=0.00D+00                                                  00003220
        DO 3030 I=1,NCAP                                                00003230
        AI=DFLOAT(I)                                                    00003240
        PI=(DATAN(1.0D+00))*4.0D+00                                     00003250
        TRUEB(1,I)=.5D+00*DSIN((2.0D+00*PI/30.0D+00)*AI)                00003260
        TRUEB(2,I)=DCOS((2.0D+00*PI/30.0D+00)*AI)                       00003270
 3030   CONTINUE                                                        00003280
        X(1,1)=1.0D+00                                                  00003290
        X(2,1)=1.0D+00                                                  00003300
        DO 3010 I=2,NCAP                                                00003310
        AI=DFLOAT(I)                                                    00003320
        X(1,I)=DSIN(10.0D+00+(AI))+.01D+00                              00003330
        X(2,I)=DCOS(10.0D+00+(AI))                                      00003340
 3010   CONTINUE                                                        00003350
 4020   CONTINUE                                                        00003360
        DO 3020 I=1,NCAP                                                00003370
C       THE DISCREPANCY TERM UU IN THE NEXT LINE IS SET TO ZERO; IT
C       CAN INSTEAD BE SET EQUAL TO A PSEUDO-RANDOM NUMBER TO TEST FLS
C       IN THE PRESENCE OF NOISE IN THE OBSERVATIONS.
        UU = 0.0D+00                                                    00003380
        Y(I)=X(1,I)*TRUEB(1,I)+X(2,I)*TRUEB(2,I)+UU                     00003390
 3020   CONTINUE                                                        00003400
        RETURN                                                          00003410
        END                                                             00003420
C                                                                       00003430
        SUBROUTINE FOCTST(AMU,K,NCAP,X,Y,B)                             00003440
        IMPLICIT REAL*8(A-H,O-Z)                                        00003450
        DIMENSION X(10,110),Y(110),B(10,110),DIF(10)                    00003460
        WRITE(6,100)                                                    00003470
  100   FORMAT(1H0,'HERE ARE THE FOC TEST RESULTS FOR EQUATIONS (A.14)')00003480
        DO 1 N=1,NCAP                                                   00003490
        IF (N.NE.1) GO TO 9000                                          00003500
        SUM=0.0D+00                                                     00003510
        DO 2 J1=1,K                                                     00003520
        SUM=SUM+X(J1,N)*B(J1,N)                                         00003530
    2   CONTINUE                                                        00003540
        SUM=SUM-Y(N)                                                    00003550
        DO 3 J1=1,K                                                     00003560
        DIF(J1)=SUM*X(J1,N)-AMU*(B(J1,N+1)-B(J1,N))                     00003570
    3   CONTINUE                                                        00003580
        WRITE (6,200)  N,(DIF(J1),J1=1,K)                               00003590
  200   FORMAT(1X,'FOR N EQUAL TO',I5/1X,6D12.3)                        00003600
        GO TO 1                                                         00003610
 9000   IF (N.EQ.NCAP) GO TO 9001                                       00003620
        SUM=0.0D+00                                                     00003630
        DO 4 J1=1,K                                                     00003640
        SUM=SUM+X(J1,N)*B(J1,N)                                         00003650
    4   CONTINUE                                                        00003660
        SUM=SUM-Y(N)                                                    00003670
        DO 5 J1=1,K                                                     00003680
        DIF(J1)=SUM*X(J1,N)-AMU*(B(J1,N+1)-B(J1,N))                     00003690
     &  +AMU*(B(J1,N)-B(J1,N-1))                                        00003700
    5   CONTINUE                                                        00003710
        WRITE(6,200) N,(DIF(J1),J1=1,K)                                 00003720
        GO TO 1                                                         00003730
 9001   SUM=0.0D+00                                                     00003740
        DO 6 J1=1,K                                                     00003750
        SUM=SUM+X(J1,N)*B(J1,N)                                         00003760
    6   CONTINUE                                                        00003770
        SUM=SUM-Y(N)                                                    00003780
        DO 7 J1=1,K                                                     00003790
        DIF(J1)=SUM*X(J1,N)+AMU*(B(J1,N)-B(J1,N-1))                     00003800
    7   CONTINUE                                                        00003810
        WRITE(6,200) N,(DIF(J1),J1=1,K)                                 00003820
    1   CONTINUE                                                        00003830
        RETURN                                                          00003840
        END                                                             00003850

Copyright © Leigh Tesfatsion. All Rights Reserved.