|
 |
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:
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.
Software Release Disclaimer:
//* 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