FLS-TVLR: Flexible Least Squares for
Time-Varying Linear Regression

Last Updated: 6 January 2016

Site maintained by:
Leigh Tesfatsion
Professor of Econ, Math, and ECpE
Heady Hall 375
Iowa State University
Ames, Iowa 50011-1070
FAX: 515-294-0221
http://www2.econ.iastate.edu/tesfatsi/
tesfatsi AT iastate.edu

FLS Homepage

About FLS-TVLR:

The Fortran program FLS-TVLR implements the flexible least squares (FLS) approach to time-varying linear regression developed by Robert E. Kalaba and Leigh Tesfatsion in "Time-Varying Linear Regression Via Flexible Least Squares" (pdf,2.4M), Computers and Mathematics with Applications Vol. 17 (1989), pp. 1215-1245. The published article is available from Science Direct.

In analogy to the determination of Pareto efficiency frontiers for multicriteria decision problems, the basic FLS-TVLR objective is to determine the Residual Efficiency Frontier (REF) for time-varying linear regression problems. The REF is the collection of all coefficient trajectory estimates that yield vector-minimal sums of squared residual measurement and dynamic errors conditional on a given set of observations.

The FLS-TVLR program has been incorporated into the statistical packages SHAZAM and GAUSS - TSM (Time Series Methods). See the FLS Homepage for more detailed information about the general FLS approach to model estimation.

The FLS-TVLR Program:

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-TVLR:                                                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 © 2013 Leigh Tesfatsion. All Rights Reserved.