|
|
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.
//* 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