1C -------------------------------------------------------------------- 2C KS IMPLEMENTS THE KALMAN SMMOOTHER RECURSIONS in 3C Koopman (1997), JASA, 92, 440, 1630-38 4C Developed by A.Rossi, C.Planas and G.Fiorentini 5C 6C XS = E[x(t)|y(1),...,y(nobs)] 7C PS = V[x(t)|y(1),...,y(nobs)], t = 1,2,...,nobs 8C 9C State-space format: y(t) = c(t)z(t) + H(t)x(t) + G(t)u(t) 10C x(t) = a(t) + F(t)x(t-1) + R(t)u(t) 11C 12C y(t) (ny x 1) ny = # of endogenous series 13C z(t) (nz x 1) nz = # of exogenous series 14C x(t) (nx x 1) nx = # of continous states 15C u(t) (nu x 1) nu = # of shocks 16C c(t) (ny x nz x ns1) ns1 = # of states for c(t) 17C H(t) (ny x nx x ns2) ns2 = # of states for S2(t) 18C G(t) (ny x nu x ns3) ns3 = # of states for S3(t) 19C a(t) (nx x ns4) ns4 = # of states for S4(t) 20C F(t) (nx x nx x ns5) ns5 = # of states for S5(t) 21C R(t) (nx x nu x ns6) ns6 = # of states for S6(t) 22C 23C d(1): order of integration of the system 24C d(2): number of non-stationary elements 25C 26C Copyright (C) 2010-2014 European Commission 27C 28C This file is part of Program DMM 29C 30C DMM is free software developed at the Joint Research Centre of the 31C European Commission: you can redistribute it and/or modify it under 32C the terms of the GNU General Public License as published by 33C the Free Software Foundation, either version 3 of the License, or 34C (at your option) any later version. 35C 36C DMM is distributed in the hope that it will be useful, 37C but WITHOUT ANY WARRANTY; without even the implied warranty of 38C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 39C GNU General Public License for more details. 40C 41C You should have received a copy of the GNU General Public License 42C along with DMM. If not, see <http://www.gnu.org/licenses/>. 43C -------------------------------------------------------------------- 44 SUBROUTINE KS(nobs,d,ny,nz,nx,nu,ns,S,yk,IYK,c,H,G,a,F,R,XS,PS) 45C INPUT 46 INTEGER nobs,d(2),ny,nz,nx,nu,ns(6),S(nobs,6),IYK(nobs,ny+1) 47 DOUBLE PRECISION yk(nobs,ny+nz),c(ny,max(nz,1),ns(1)), 48 1 H(ny,nx,ns(2)),G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)), 49 2 R(nx,nu,ns(6)) 50 51C OUTPUT 52 DOUBLE PRECISION XS(nobs,nx),PS(nobs,nx,nx) 53 54C LOCALS 55 INTEGER imain,iny,I,J,IFAIL,FiRANK,ITIME 56 INTEGER IPIV(nx) 57 DOUBLE PRECISION TOL,SUMW1 58 DOUBLE PRECISION,ALLOCATABLE:: Pi(:,:,:),HPs(:,:),HPi(:,:), 59 1 Fi(:,:),Fs(:,:),Fim(:,:),Fsm(:,:),PHFs(:,:),PHFi(:,:),FFF(:,:), 60 2 Mi(:,:),Ms(:,:),Ci(:,:),Kst(:,:,:),Kit(:,:,:),W1(:),WORK(:), 61 3 WORK1(:),PFP(:,:),APPO(:,:),APPO1(:,:),COM(:,:),RG(:,:), 62 4 FP(:,:),HP1(:,:),V(:,:),CC(:,:),HPV(:,:),X0(:),P0(:,:), 63 5 RECR(:),RECRI(:),RECN(:,:),XT(:,:),PT(:,:,:),INN(:,:), 64 1 Vinv(:,:,:),Vis(:,:,:),Vii(:,:,:) 65 66C EXTERNAL SUBROUTINES 67 EXTERNAL INVF,INVFBIS,LYAP,DSYEV,DPOTRF,DPOTRI,DGETRF,DGETRI 68 69 ALLOCATE(Pi(d(1),nx,nx),HPs(ny,nx),HPi(ny,nx), 70 1 Fi(ny,ny),Fs(ny,ny),Fim(ny,ny),Fsm(ny,ny), 71 2 PHFs(nx,ny),PHFi(nx,ny),FFF(ny,ny),Mi(nx,ny),Ms(nx,ny),Ci(nx,nx), 72 3 Kst(d(1),nx,ny),Kit(d(1),nx,ny)) 73 74 ALLOCATE(W1(ny),WORK(64*nx),WORK1(64*ny), 75 1 PFP(nx,nx),APPO(nx,nx),APPO1(nx,ny),COM(ny+1,ny),RG(nx,ny)) 76 77 ALLOCATE(FP(nx,nx),HP1(ny,nx),V(ny,ny),CC(nx,nx), 78 1 HPV(nx,ny),X0(nx),P0(nx,nx),RECR(nx),RECRI(nx),RECN(nx,nx)) 79 80 ALLOCATE(XT(nobs,nx),PT(nobs,nx,nx),INN(nobs,ny), 81 1 Vinv(nobs,ny,ny),Vis(d(1),ny,ny),Vii(d(1),ny,ny)) 82 83 TOL = 1.D-3 84C Unconditional mean and variance 85 IF (d(1).EQ.0) THEN ! stationary models 86 IF(SUM(ABS(a(:,S(1,4)))).EQ.0.D0) THEN 87 XT(1,:) = 0.D0 ! X(1|0) 88 ELSE 89 APPO = -F(:,:,S(1,5)) 90 DO 1 I = 1,nx 911 APPO(I,I) = 1.D0+APPO(I,I) 92C CALL F07ADF(nx,nx,APPO,nx,IPIV,IFAIL) 93C CALL F07AJF(nx,APPO,nx,IPIV,WORK,64*nx,IFAIL) 94 CALL DGETRF(nx,nx,APPO,nx,IPIV,IFAIL) 95 CALL DGETRI(nx,APPO,nx,IPIV,WORK,64*nx,IFAIL) 96 97 DO 2 I =1,nx 982 XT(1,I) = SUM(APPO(I,:)*a(:,S(1,4))) ! inv(I-F)*a 99 ENDIF 100 101C P(1|0) - F*P(1|0)*F' = R*R' 102 CALL LYAP(nx,nu,TOL,F(:,:,S(1,5)),R(:,:,S(1,6)),PT(1,:,:)) 103 ELSE 104C ----------------------------------------------------------- 105C Non-stationary models 106C Define X(1) = aa + A*eta + B*delta (A*B' = 0) 107C eta~N(0,I), delta~N(0,k*I) k -> +inf 108C X(1)~N(aa,P), P=Ps+k*Pi, Ps=AA', Pi=BB'. 109C CARE!! aa (uncond. mean),Ps, and Pi to be filled by users 110C ----------------------------------------------------------- 111 XT(1,1:nx) = 0.D0 ! X(1|0) 112 PT(1,1:nx,1:nx) = 0.D0 ! P(1|0) 113 IF (d(2).LT.nx) THEN 114 IF(SUM(ABS(a(d(2)+1:nx,S(1,4)))).NE.0.D0) THEN 115 APPO(d(2)+1:nx,d(2)+1:nx) = -F(d(2)+1:nx,d(2)+1:nx,S(1,5)) 116 DO 3 I = d(2)+1,nx 1173 APPO(I,I) = 1.D0+APPO(I,I) 118C CALL F07ADF(nx-d(2),nx-d(2),APPO(d(2)+1:nx,d(2)+1:nx),nx-d(2), 119C 1 IPIV(d(2)+1:nx),IFAIL) 120C CALL F07AJF(nx-d(2),APPO(d(2)+1:nx,d(2)+1:nx),nx-d(2), 121C 1 IPIV(d(2)+1:nx),WORK,64*nx,IFAIL) 122 CALL DGETRF(nx-d(2),nx-d(2),APPO(d(2)+1:nx,d(2)+1:nx),nx-d(2), 123 1 IPIV(d(2)+1:nx),IFAIL) 124 CALL DGETRI(nx-d(2),APPO(d(2)+1:nx,d(2)+1:nx),nx-d(2), 125 1 IPIV(d(2)+1:nx),WORK,64*nx,IFAIL) 126 127 DO 4 I = d(2)+1,nx 1284 XT(1,I) = SUM(APPO(I,d(2)+1:nx)*a(d(2)+1:nx,S(1,4))) ! inv(I-F)*a 129 ENDIF 130 131C Lyapunov eqn 132 CALL LYAP(nx-d(2),nu,TOL,F(d(2)+1:nx,d(2)+1:nx,S(1,5)), 133 1 R(d(2)+1:nx,1:nu,S(1,6)),PT(1,d(2)+1:nx,d(2)+1:nx)) 134 ENDIF 135 136 Pi(:,:,:) = 0.D0 137 DO 5 I = 1,d(2) 1385 Pi(1,I,I) = 1.D0 139 140 DO 200 imain = 1,d(1) 141 iny = IYK(imain,ny+1) 142 DO 30 I=1,iny 143 DO 30 J=1,nx 14430 HPs(I,J) = SUM(H(IYK(imain,I),:,S(imain,2))*PT(imain,:,J)) 145 146 DO 40 I=1,iny 147 Fs(I,I) = SUM(HPs(I,:)*H(IYK(imain,I),:,S(imain,2))) 148 + +SUM(G(IYK(imain,I),:,S(imain,3))*G(IYK(imain,I),:,S(imain,3))) 149 DO 40 J=1,I-1 150 Fs(I,J) = SUM(HPs(I,:)*H(IYK(imain,J),:,S(imain,2))) 151 + +SUM(G(IYK(imain,I),:,S(imain,3))*G(IYK(imain,J),:,S(imain,3))) 15240 Fs(J,I) = Fs(I,J) 153 154 DO 50 I=1,iny 155 DO 50 J=1,nx 15650 HPi(I,J) = SUM(H(IYK(imain,I),:,S(imain,2))*Pi(imain,:,J)) 157 158 DO 60 I=1,iny 159 Fi(I,I) = SUM(HPi(I,:)*H(IYK(imain,I),:,S(imain,2))) 160 DO 60 J=1,I-1 161 Fi(I,J) = SUM(HPi(I,:)*H(IYK(imain,J),:,S(imain,2))) 16260 Fi(J,I) = Fi(I,J) 163 164C -------------------------------------------------------------------------- 165C Computes inverse of the innovation variance matrix 166C Cases: ny = 1, Fi is scalar >0 (or 0 not considered) 167C ny > 1, Fi is full rank or singular (or 0 matrix not considered) 168C -------------------------------------------------------------------------- 169 IF (iny.EQ.1) THEN 170 Fsm = 0.D0 171 Fim = 1.D0/Fi 172 FFF = Fim*Fs*Fim 173 ELSE 174 175 IFAIL = -1 176 COM(1:iny,1:iny) = Fi(1:iny,1:iny) 177C CALL F02FAF('N','U',iny,COM(1:iny,1:iny),iny,W1(1:iny), 178C 1 WORK1,64*iny,IFAIL) 179 CALL DSYEV('N','U',iny,COM(1:iny,1:iny),iny,W1(1:iny), 180 1 WORK1,64*iny,IFAIL) 181 182 FiRANK = 0 183 SUMW1 = SUM(ABS(W1(1:iny))) 184 DO 70 I=1,iny 185 W1(I) = W1(I)/SUMW1 18670 IF (W1(I).GT.1.D-10) FiRANK=FiRANK+1 187 FiRANK = min(FiRANK,d(2)) 188 189 IF(FiRANK.EQ.iny) THEN 190 Fsm = 0.D0 191 COM(1:iny,1:iny) = Fi(1:iny,1:iny) 192 IFAIL = -1 193C CALL F01ADF(iny,COM(1:iny+1,1:iny),iny+1,IFAIL) 194 CALL DPOTRF('L',iny,COM(1:iny,1:iny),iny,IFAIL) ! COM = L*L' 195 CALL DPOTRI('L',iny,COM(1:iny,1:iny),iny,IFAIL) ! COM = VV^-1 196 197 DO 80 I=1,iny 198 Fim(I,I) = COM(I,I) 199 DO 80 J=1,I-1 200 Fim(I,J) = COM(I,J) 20180 Fim(J,I) = Fim(I,J) 202 203 DO 81 I=1,iny 204 DO 81 J=1,iny 20581 COM(I,J) = SUM(Fim(I,1:iny)*Fs(1:iny,J)) ! Fim x Fs 206 207 DO 82 I=1,iny 208 FFF(I,I) = SUM(COM(I,1:iny)*Fim(1:iny,I)) 209 DO 82 J=1,I-1 210 FFF(I,J) = SUM(COM(I,1:iny)*Fim(1:iny,J)) ! Fim x Fs x Fim 21182 FFF(J,I) = FFF(I,J) 212 ELSE 213 SUMW1=0.D0 214 DO I=Firank+1,iny 215 SUMW1 = SUMW1 + Fi(I,I) 216 ENDDO 217 IF (SUMW1.GT.0.D0) THEN 218 CALL INVFBIS(Fs(1:iny,1:iny),Fi(1:iny,1:iny),iny,FiRANK, 219 1 Fsm(1:iny,1:iny),Fim(1:iny,1:iny),FFF(1:iny,1:iny)) 220 ELSE 221 CALL INVF(Fs(1:iny,1:iny),Fi(1:iny,1:iny),iny,FiRANK, 222 1 Fsm(1:iny,1:iny),Fim(1:iny,1:iny),FFF(1:iny,1:iny)) 223 ENDIF 224 ENDIF 225 ENDIF 226 Vis(imain,1:iny,1:iny) = Fsm(1:iny,1:iny) 227 Vii(imain,1:iny,1:iny) = Fim(1:iny,1:iny) 228 229C ------------------------------------------------------------------ 230C X(d|d) = X(d|d-1)+((Ps*H'+R*G')*Fsm+Pi*H'*Fim)*(Y(d)-H*X(d|d-1)-c) 231C ------------------------------------------------------------------ 232 DO 85 I = 1,nx 233 DO 85 J = 1,iny 234 RG(I,J) = 235 # SUM(R(I,1:nu,S(imain,6))*G(IYK(imain,J),1:nu,S(imain,3))) 23685 HPs(J,I) = HPs(J,I) + RG(I,J) ! HPs = (Ps*H'+R*G')' 237 238 DO 90 I = 1,nx 239 DO 90 J = 1,iny 240 PHFs(I,J) = SUM(HPs(1:iny,I)*Fsm(1:iny,J)) 24190 PHFi(I,J) = SUM(HPi(1:iny,I)*Fim(1:iny,J)) 242 243C Innovations 244 DO 100 I=1,iny 245100 INN(imain,I) = yk(imain,IYK(imain,I)) 246 + - SUM(H(IYK(imain,I),1:nx,S(imain,2))*XT(imain,1:nx)) 247 + - SUM(c(IYK(imain,I),1:nz,S(imain,1))*yk(imain,ny+1:ny+nz)) 248 249 DO 110 I=1,nx 250110 X0(I) = XT(imain,I) 251 + + SUM((PHFs(I,1:iny)+PHFi(I,1:iny))*INN(imain,1:iny)) 252 253C P(d|d) = P(d|d-1) + Pi*H'*Fim*Fs*Fim*H*Pi - Ps*H'*Fsm*H*Ps - Ps*H'*Fim*H*Pi - (Ps*H'*Fim*H*Pi)' 254C - Ps*H'*Fsm*H*Ps 255 DO 120 I = 1,nx 256 APPO(I,I) = -SUM(PHFs(I,1:iny)*HPs(1:iny,I)) 257 DO 120 J = 1,I-1 258 APPO(I,J) = -SUM(PHFs(I,1:iny)*HPs(1:iny,J)) 259120 APPO(J,I) = APPO(I,J) 260 261C - Ps*H'*Fim*H*Pi - (Ps*H'*Fim*H*Pi)' 262 DO 130 I = 1,nx 263 APPO(I,I) = APPO(I,I) - SUM(HPs(1:iny,I)*PHFi(I,1:iny)) 264 + - SUM(PHFi(I,1:iny)*HPs(1:iny,I)) 265 DO 130 J = 1,I-1 266 APPO(I,J) = APPO(I,J) - SUM(HPs(1:iny,I)*PHFi(J,1:iny)) 267 + - SUM(PHFi(I,1:iny)*HPs(1:iny,J)) 268130 APPO(J,I) = APPO(I,J) 269 270C Pi*H'*Fim*Fs*Fim*H*Pi 271 DO 140 I = 1,nx 272 DO 140 J = 1,iny 273140 APPO1(I,J) = SUM(HPi(1:iny,I)*FFF(1:iny,J)) 274 275 DO 150 I = 1,nx 276 PFP(I,I) = SUM(APPO1(I,1:iny)*HPi(1:iny,I)) 277 DO 150 J = 1,I-1 278 PFP(I,J) = SUM(APPO1(I,1:iny)*HPi(1:iny,J)) 279150 PFP(J,I) = PFP(I,J) 280 281 P0(:,:) = PT(imain,:,:) + PFP(:,:) + APPO(:,:) 282 283C Mi = F*Pi*H' 284 DO 151 I = 1,nx 285 DO 151 J = 1,nx 286151 PFP(I,J) = SUM(F(I,1:nx,S(imain,5))*Pi(imain,1:nx,J)) ! F*Pi 287 288 DO 152 I = 1,nx 289 DO 152 J = 1,iny 290152 Mi(I,J) = SUM(PFP(I,1:nx)*H(IYK(imain,J),1:nx,S(imain,2))) 291 292C Ms = F*Ps*H' + R*G' 293 DO 153 I = 1,nx 294 DO 153 J = 1,nx 295153 FP(I,J) = SUM(F(I,:,S(imain,5))*PT(imain,:,J)) ! F*Ps 296 297 DO 154 I = 1,nx 298 DO 154 J = 1,iny 299154 Ms(I,J)=RG(I,J)+SUM(FP(I,1:nx)*H(IYK(imain,J),1:nx,S(imain,2))) 300 301C Kit = Ms*Fim - Mi*FFF 302C Kst = Ms*Fsm + Mi*Fim 303 DO 155 I = 1,nx 304 DO 155 J = 1,iny 305 Kit(imain,I,J) = SUM(Ms(I,1:iny)*Fim(1:iny,J)) 306 + - SUM(Mi(I,1:iny)*FFF(1:iny,J)) 307155 Kst(imain,I,J) = SUM(Ms(I,1:iny)*Fsm(1:iny,J)) 308 + + SUM(Mi(I,1:iny)*Fim(1:iny,J)) 309 310C ---------------------------------- 311C Predictions X(t+1|t) and P(t+1|t) 312C ---------------------------------- 313 IF (imain.LT.d(1)) THEN 314C Pi+1 = F*Pi*F'-Ci 315 316C Ci = Mi*Fim*Mi' 317 DO 164 I = 1,nx 318 DO 164 J = 1,iny 319164 RG(I,J) = SUM(Mi(I,1:iny)*Fim(1:iny,J)) ! Mi*Fim 320 321 DO 166 I = 1,nx 322 Ci(I,I) = SUM(RG(I,1:iny)*Mi(I,1:iny)) 323 DO 166 J = 1,I-1 324166 Ci(I,J) = SUM(RG(I,1:iny)*Mi(J,1:iny)) 325 326 DO 168 I = 1,nx 327 Pi(imain+1,I,I)=SUM(PFP(I,1:nx)*F(I,1:nx,S(imain+1,5)))-Ci(I,I) 328 DO 168 J = 1,I-1 329 Pi(imain+1,I,J)=SUM(PFP(I,1:nx)*F(J,1:nx,S(imain+1,5)))-Ci(I,J) 330168 Pi(imain+1,J,I) = Pi(imain+1,I,J) 331 332 ENDIF 333 334C X(t+1|t) = a + F*X(t|t) 335 DO 170 I=1,nx 336170 XT(imain+1,I) = a(I,S(imain+1,4)) 337 + + SUM(F(I,1:nx,S(imain+1,5))*X0(1:nx)) 338 339C P(t+1|t) = F*PddF' + R*R' 340 DO 172 I = 1,nx 341 DO 172 J = 1,nx 342172 APPO(I,J) = SUM(F(I,:,S(imain+1,5))*P0(:,J)) ! F*Pdd 343 344 DO 180 I = 1,nx 345 PT(imain+1,I,I) = SUM(APPO(I,1:nx)*F(I,1:nx,S(imain+1,5))) 346 + + SUM(R(I,1:nu,S(imain+1,6))*R(I,1:nu,S(imain+1,6))) 347 DO 180 J = 1,I-1 348 PT(imain+1,I,J) = SUM(APPO(I,1:nx)*F(J,1:nx,S(imain+1,5))) 349 + + SUM(R(I,1:nu,S(imain+1,6))*R(J,1:nu,S(imain+1,6))) 350180 PT(imain+1,J,I) = PT(imain+1,I,J) 351 352200 CONTINUE 353 ENDIF 354 355 DO 400 imain = d(1)+1,nobs 356 iny = IYK(imain,ny+1) 357C ------------------------------- 358C Innovations: INN = yk-H*X1-c*z 359C ------------------------------- 360 DO 210 I=1,iny 361210 INN(imain,I) = yk(imain,IYK(imain,I)) 362 + - SUM(H(IYK(imain,I),1:nx,S(imain,2))*XT(imain,:)) 363 + - SUM(c(IYK(imain,I),1:nz,S(imain,1))*yk(imain,ny+1:ny+nz)) 364 365C ---------------------------------------------------------- 366C Innovation variance V = H*P1*H' + G*G' + H*R*G' + G*R'*H' 367C ---------------------------------------------------------- 368 DO 220 I=1,iny 369 DO 220 J=1,nx 370220 HP1(I,J) = SUM(H(IYK(imain,I),1:nx,S(imain,2))*PT(imain,1:nx,J)) 371 372 DO 221 I=1,nx 373 DO 221 J=1,iny 374221 RG(I,J)=SUM(R(I,1:nu,S(imain,6))*G(IYK(imain,J),1:nu,S(imain,3))) 375 376 DO 222 I=1,iny 377 DO 222 J=1,iny 378222 COM(I,J)=SUM(H(IYK(imain,I),1:nx,S(imain,2))*RG(1:nx,J)) ! H*R*G' 379 380 DO 230 I=1,iny 381 V(I,I) = SUM(HP1(I,1:nx)*H(IYK(imain,I),1:nx,S(imain,2))) 382 # + SUM(G(IYK(imain,I),1:nu,S(imain,3))* 383 # G(IYK(imain,I),1:nu,S(imain,3))) + 2.*COM(I,I) 384 DO 230 J=1,I-1 385 V(I,J) = SUM(HP1(I,1:nx)*H(IYK(imain,J),1:nx,S(imain,2)))+ 386 # SUM(G(IYK(imain,I),1:nu,S(imain,3))* 387 # G(IYK(imain,J),1:nu,S(imain,3)))+COM(I,J)+COM(J,I) 388230 V(J,I) = V(I,J) 389 390C ------------------------------------------------------------------- 391C Updating equations: 392C x0 = x1 + (P1*H'+R*G')*Vinv*INN 393C p0 = p1 - (P1*H'+R*G')*Vinv*(P1*H'+R*G')' 394C ------------------------------------------------------------------- 395 IF (iny.GT.0) THEN 396 COM(1:iny,1:iny) = V(1:iny,1:iny) 397 IFAIL = -1 398C CALL F01ADF(iny,COM(1:iny+1,1:iny),iny+1,IFAIL) 399 CALL DPOTRF('L',iny,COM(1:iny,1:iny),iny,IFAIL) ! COM = L*L' 400 CALL DPOTRI('L',iny,COM(1:iny,1:iny),iny,IFAIL) ! COM = VV^-1 401 402 DO 240 I=1,iny 403 Vinv(imain,I,I) = COM(I,I) 404 DO 240 J=1,I-1 405 Vinv(imain,I,J) = COM(I,J) 406240 Vinv(imain,J,I) = Vinv(imain,I,J) 407 408 DO 260 I=1,nx 409 DO 260 J=1,iny 410260 HPV(I,J) = SUM((HP1(1:iny,I)+RG(I,1:iny))*Vinv(imain,1:iny,J)) 411 412 DO 270 I=1,nx 413270 X0(I) = XT(imain,I)+SUM(HPV(I,1:iny)*INN(imain,1:iny)) 414 415 DO 280 I=1,nx 416 P0(I,I) = PT(imain,I,I) 417 + - SUM(HPV(I,1:iny)*(HP1(1:iny,I)+RG(I,1:iny))) 418 DO 280 J=1,I-1 419 P0(I,J) = PT(imain,I,J) 420 + - SUM(HPV(I,1:iny)*(HP1(1:iny,J)+RG(J,1:iny))) 421280 P0(J,I) = P0(I,J) 422 423 ELSE 424 425 X0(1:nx) = XT(imain,1:nx) 426 P0(1:nx,1:nx) = PT(imain,1:nx,1:nx) 427 428 ENDIF 429 430C ------------------------------------ 431C Prediction x1 = c+F*x0 432C Prediction var. P1 = F*p0*F'+ R*R' 433C ------------------------------------ 434 IF (imain.LT.nobs) THEN 435 DO 290 I=1,nx 436290 XT(imain+1,I) = a(I,S(imain+1,4)) 437 + + SUM(F(I,1:nx,S(imain+1,5))*X0(1:nx)) 438 439 DO 300 I=1,nx 440 DO 300 J=1,nx 441300 FP(I,J) = SUM(F(I,1:nx,S(imain+1,5))*P0(1:nx,J)) 442 443 DO 310 I=1,nx 444 PT(imain+1,I,I) = SUM(FP(I,:)*F(I,:,S(imain+1,5))) 445 + + SUM(R(I,1:nu,S(imain+1,6))*R(I,1:nu,S(imain+1,6))) 446 DO 310 J=1,I-1 447 PT(imain+1,I,J) = SUM(FP(I,:)*F(J,:,S(imain+1,5))) 448 + + SUM(R(I,1:nu,S(imain+1,6))*R(J,1:nu,S(imain+1,6))) 449310 PT(imain+1,J,I) = PT(imain+1,I,J) 450 ENDIF 451400 CONTINUE 452 453C **** SMOOTHING BAKWARD RECURSIONS **** 454 RECR(1:nx) = 0.D0 455 RECN(1:nx,1:nx) = 0.D0 456 DO 600 ITIME = nobs,d(1)+1,-1 457 iny = IYK(ITIME,ny+1) 458 459C R*G' and H'*Vinv 460 DO 420 J=1,iny 461 DO 420 I=1,nx 462 RG(I,J)=SUM(R(I,1:nu,S(ITIME,6))*G(IYK(ITIME,J),1:nu,S(ITIME,3))) 463420 PHFs(I,J) = 464 # SUM(H(IYK(ITIME,1:iny),I,S(ITIME,2))*Vinv(ITIME,1:iny,J)) 465 466C F(t+1)*P(t|t-1) 467 DO 430 I=1,nx 468 DO 430 J=1,nx 469430 FP(I,J) = SUM(F(I,1:nx,S(min(nobs,ITIME+1),5))*PT(ITIME,1:nx,J)) 470 471C H'*Vinv*H 472 DO 440 I=1,nx 473 APPO(I,I) = SUM(PHFs(I,1:iny)*H(IYK(ITIME,1:iny),I,S(ITIME,2))) 474 DO 440 J=1,I-1 475 APPO(I,J) = SUM(PHFs(I,1:iny)*H(IYK(ITIME,1:iny),J,S(ITIME,2))) 476440 APPO(J,I) = APPO(I,J) 477 478C L(t) = F(t+1)-(F(t+1)*P(t|t-1)*H(t)'+R(t)*G(t)')*Vinv(t)*H(t) 479 DO 450 I=1,nx 480 DO 450 J=1,nx 481450 PFP(I,J) = F(I,J,S(min(nobs,ITIME+1),5)) 482 + - SUM(FP(I,1:nx)*APPO(1:nx,J)) 483 + - SUM(RG(I,1:iny)*PHFs(J,1:iny)) 484 485C r(t-1) = H(t)'*Vinv(t)*INN(t)+L(t)'*r(t) 486 DO 460 I=1,nx 487460 WORK(I) = SUM(PFP(1:nx,I)*RECR(1:nx)) 488 489 DO 470 I=1,nx 490470 RECR(I) = WORK(I) + SUM(PHFs(I,1:iny)*INN(ITIME,1:iny)) 491 492C N(t-1) = H(t)'*Vinv(t)*H(t)+L(t)'*N(t)*L(t) 493 DO 480 I=1,nx 494 DO 480 J=1,nx 495480 CC(I,J) = SUM(PFP(1:nx,I)*RECN(1:nx,J)) ! L(t)'*N(t) 496 497 DO 490 I=1,nx 498 RECN(I,I) = APPO(I,I) + SUM(CC(I,1:nx)*PFP(1:nx,I)) 499 DO 490 J=1,I-1 500 RECN(I,J) = APPO(I,J) + SUM(CC(I,1:nx)*PFP(1:nx,J)) 501490 RECN(J,I) = RECN(I,J) 502 503C X(t|T) = X(t|t-1) + P(t|t-1)*r(t-1) 504 DO 500 I = 1,nx 505500 XS(ITIME,I) = XT(ITIME,I) + SUM(PT(ITIME,I,1:nx)*RECR(1:nx)) 506 507C P(t|T) = P(t|t-1) - P(t|t-1)*N(t-1)*P(t|t-1) 508 DO 510 I=1,nx 509 DO 510 J=1,nx 510510 CC(I,J) = SUM(PT(ITIME,I,1:nx)*RECN(1:nx,J)) ! P(t|t-1)*N(t-1) 511 512 DO 520 I=1,nx 513 PS(ITIME,I,I) = PT(ITIME,I,I) - SUM(CC(I,1:nx)*PT(ITIME,1:nx,I)) 514 DO 520 J=1,I-1 515 PS(ITIME,I,J) = PT(ITIME,I,J) - SUM(CC(I,1:nx)*PT(ITIME,1:nx,J)) 516520 PS(ITIME,J,I) = PS(ITIME,I,J) 517 518600 CONTINUE 519 520C INITIAL KALMAN SAMOOTING 521 RECRI(1:nx) = 0.D0 522 DO 800 ITIME = d(1),1,-1 523 iny = IYK(ITIME,ny+1) 524C L(t) = F(t)-Kst(t)*H(t) 525 DO 610 I=1,nx 526 DO 610 J=1,nx 527610 PFP(I,J) = F(I,J,S(ITIME,5)) 528 + - SUM(Kst(ITIME,I,1:iny)*H(IYK(ITIME,1:iny),J,S(ITIME,2))) 529 530C r(t-1) = H(t)'*Fsm*INN(t) + L(t)'*r(t) 531C ri(t-1) = H(t)'*Fim*INN(t) + L(t)'*ri(t) + Li(t)'*r(t) 532 DO 620 I=1,nx 533 WORK(I) = SUM(PFP(1:nx,I)*RECR(1:nx)) ! L(t)'*r(t) 534620 WORK(nx+I) = SUM(PFP(1:nx,I)*RECRI(1:nx)) ! L(t)'*ri(t) 535 536C Li(t) = -Kit(t)*H(t) 537 DO 621 I=1,nx 538 DO 621 J=1,nx 539621 PFP(I,J) = 540 # - SUM(Kit(ITIME,I,1:iny)*H(IYK(ITIME,1:iny),J,S(ITIME,2))) 541 542C L(t)'*ri(t) + Li(t)'*r(t) 543 DO 622 I=1,nx 544622 WORK(nx+I) = WORK(nx+I)+SUM(PFP(1:nx,I)*RECR(1:nx)) 545 546C H'*Fsm 547 DO 625 I=1,nx 548 DO 625 J=1,iny 549625 PHFs(I,J) = 550 # SUM(H(IYK(ITIME,1:iny),I,S(ITIME,2))*Vis(ITIME,1:iny,J)) 551 552 DO 630 I=1,nx 553630 RECR(I) = WORK(I) + SUM(PHFs(I,1:iny)*INN(ITIME,1:iny)) 554 555C H'*Fim 556 DO 631 I=1,nx 557 DO 631 J=1,iny 558631 PHFs(I,J) = 559 # SUM(H(IYK(ITIME,1:iny),I,S(ITIME,2))*Vii(ITIME,1:iny,J)) 560 561 DO 632 I=1,nx 562632 RECRI(I) = WORK(nx+I) + SUM(PHFs(I,1:iny)*INN(ITIME,1:iny)) 563 564C X(d|T) = X(d|d-1) + Psd*r(t-1) + Pid*ri(t-1) 565 DO 660 I = 1,nx 566660 XS(ITIME,I) = XT(ITIME,I) + SUM(PT(ITIME,I,1:nx)*RECR(1:nx)) 567 + + SUM(Pi(ITIME,I,1:nx)*RECRI(1:nx)) 568 569800 CONTINUE 570 571 DEALLOCATE(Pi,HPs,HPi,Fi,Fs,Fim,Fsm,PHFs,PHFi,FFF,Mi,Ms,Ci, 572 1 Kst,Kit,W1,WORK,WORK1,PFP,APPO,APPO1,COM,RG,FP,HP1,V,CC,HPV, 573 2 X0,P0,RECR,RECRI,RECN,XT,PT,INN,Vinv,Vis,Vii) 574 575 RETURN 576 END 577 578C This is the variance and must be completed!! 579C N(t-1) = H(t)'*Vinv(t)*H(t)+L(t)'*N(t)*L(t) 580c DO 640 I=1,nx 581c DO 640 J=1,nx 582c640 CC(I,J) = SUM(PFP(1:nx,I)*RECN(1:nx,J)) ! L(t)'*N(t) 583c DO 650 I=1,nx 584c RECN(I,I) = APPO(I,I) + SUM(CC(I,1:nx)*PFP(1:nx,I)) 585c DO 650 J=1,I-1 586c RECN(I,J) = APPO(I,J) + SUM(CC(I,1:nx)*PFP(1:nx,J)) 587c650 RECN(J,I) = RECN(I,J) 588C P(t|T) = P(t|t-1) - P(t|t-1)*N(t-1)*P(t|t-1) 589c DO 670 I=1,nx 590c DO 670 J=1,nx 591c670 CC(I,J) = SUM(PT(ITIME,I,1:nx)*RECN(1:nx,J)) ! P(t|t-1)*N(t-1) 592c 593c DO 680 I=1,nx 594c PS(ITIME,I,I) = PT(ITIME,I,I) - SUM(CC(I,1:nx)*PT(ITIME,1:nx,I)) 595c DO 680 J=1,I-1 596c PS(ITIME,I,J) = PT(ITIME,I,J) - SUM(CC(I,1:nx)*PT(ITIME,1:nx,J)) 597c680 PS(ITIME,J,I) = PS(ITIME,I,J) 598 599