1C ---------------------------------------------------------------------- 2C GCK Implements the SINGLE-MOVE Sampler of 3C Gerlach, Carter and Kohn (2000): Efficient Bayesian Inference 4C for Dynamic Mixture Models, JASA, 95,451, pp.819-28 5C Developed by A.Rossi, C.Planas and G.Fiorentini 6C 7C Pr[Z(t)|Z(\t),Y] pto Pr[Y^(t+1,T)|Y^t,Z] x Pr[Y(t)|Y^(t-1),Z^t] 8C x Pr[Z(t)|Z(\t)] 9C 10C State-space format: y(t) = c(t)z(t) + H(t)x(t) + G(t)u(t) 11C x(t) = a(t) + F(t)x(t-1) + R(t)u(t) 12C 13C y(t) (ny x 1) ny = # of endogenous series 14C z(t) (nz x 1) nz = # of exogenous series 15C x(t) (nx x 1) nx = # of continous states 16C u(t) (nu x 1) nu = # of shocks 17C c(t) (ny x nz x ns1) ns1 = # of states for c(t) 18C H(t) (ny x nx x ns2) ns2 = # of states for H(t) 19C G(t) (ny x nu x ns3) ns3 = # of states for G(t) 20C a(t) (nx x ns4) ns4 = # of states for a(t) 21C F(t) (nx x nx x ns5) ns5 = # of states for F(t) 22C R(t) (nx x nu x ns6) ns6 = # of states for R(t) 23C 24C FURTHER INPUT: 25C 26C nobs: # of observatios 27C d(1): order of integration of the system 28C d(2): number of non-stationary elements 29C nv: # of discrete latent variables (S1,S2,...) 30C ns: ns1,ns2,... 31C nstot: total # of states (states of S1 x S2 x ...x Snv) 32C nt: dimension of theta 33C np: dimension of psi 34C PMAT: (nstot x nstot) one-step transition probabilities 35C PE: ergodic distribution of S1 x S2 x ...x Snv 36C INFOS: (9 x 6) set latent variables 37C nstot: total # of states i.e. ns1 x ns2 x ...x nsv 38C 39C OUTPUT: 40C 41C Z:(nobs x 1) takes values {1,2,...,nstot} 42C 43C Copyright (C) 2010-2014 European Commission 44C 45C This file is part of Program DMM 46C 47C DMM is free software developed at the Joint Research Centre of the 48C European Commission: you can redistribute it and/or modify it under 49C the terms of the GNU General Public License as published by 50C the Free Software Foundation, either version 3 of the License, or 51C (at your option) any later version. 52C 53C DMM is distributed in the hope that it will be useful, 54C but WITHOUT ANY WARRANTY; without even the implied warranty of 55C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 56C GNU General Public License for more details. 57C 58C You should have received a copy of the GNU General Public License 59C along with DMM. If not, see <http://www.gnu.org/licenses/>. 60C ---------------------------------------------------------------------- 61 SUBROUTINE GCK(nobs,d,ny,nz,nx,nu,nv,ns,nstot,nt,np,yk,IYK, 62 1 theta,psi,INFOS,pdll,Z,S) 63 64 USE dfwin 65 INTERFACE 66 SUBROUTINE DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R) 67 INTEGER ny,nz,nx,nu,ns(6),nt 68 DOUBLE PRECISION theta(nt) 69 DOUBLE PRECISION c(ny,max(1,nz),ns(1)),H(ny,nx,ns(2)), 70 1 G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)),R(nx,nu,ns(6)) 71 END SUBROUTINE 72 END INTERFACE 73 POINTER (pdll,fittizia) ! ASSOCIATE pointer P alla DLL ad una varibile fittizia 74 POINTER (pdesign,DESIGN) ! IMPORTANT associo il puntatore pdesign alla Interface definita 75 76C INPUT 77 INTEGER nobs,d(2),ny,nz,nx,nu,nv,ns(6),nstot,nt,np,IYK(nobs,ny+1), 78 1 INFOS(9,6) 79 DOUBLE PRECISION yk(nobs,ny+nz),theta(nt),psi(np) 80 81C INPUT/OUTPUT 82 INTEGER Z(nobs),S(nobs,6) 83 84C LOCALS 85 INTEGER IT,I,J,K,IFAIL,ISEQ,IMAX(1),iny,NIFS,KKK 86 INTEGER IS(6),SH(3),IND(max(1,d(1)),6),SEQ(nv),IFS(nstot),dc(2) 87 DOUBLE PRECISION c(ny,max(nz,1),ns(1)),H(ny,nx,ns(2)), 88 1 G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)),R(nx,nu,ns(6)) 89 DOUBLE PRECISION PMAT(nstot,nstot),PE(nstot), 90 1 P1(MAX(1,INFOS(8,1)*MIN(INFOS(9,1),1)), 91 1 MAX(1,INFOS(8,1)*MIN(INFOS(9,1),1))), 92 1 P2(MAX(1,INFOS(8,2)*MIN(INFOS(9,2),1)), 93 1 MAX(1,INFOS(8,2)*MIN(INFOS(9,2),1))), 94 1 P3(MAX(1,INFOS(8,3)*MIN(INFOS(9,3),1)), 95 1 MAX(1,INFOS(8,3)*MIN(INFOS(9,3),1))), 96 1 P4(MAX(1,INFOS(8,4)*MIN(INFOS(9,4),1)), 97 1 MAX(1,INFOS(8,4)*MIN(INFOS(9,4),1))), 98 1 P5(MAX(1,INFOS(8,5)*MIN(INFOS(9,5),1)), 99 1 MAX(1,INFOS(8,5)*MIN(INFOS(9,5),1))), 100 1 P6(MAX(1,INFOS(8,6)*MIN(INFOS(9,6),1)), 101 1 MAX(1,INFOS(8,6)*MIN(INFOS(9,6),1))) 102 DOUBLE PRECISION OM(nobs,nx,nx),MU(nobs,nx) 103 DOUBLE PRECISION HRG(ny,nu),VV(ny,ny),HF(ny,nx),COM(ny+1,ny), 104 1 HRGV(nu,ny),B(nx,ny),RR(nx,nx),BH(nx,nx),BHRG(nx,nu),DD(nx,nx), 105 2 CS(nx,nx),CC(nx,nx),AA(nx,nx),DI(nx,nx),COM1(nx+1,nx),Ha(ny), 106 3 OMC(nx,nx),OMCDIC(nx,nx),AOMCDIC(nx,nx),AOMCDICOM(nx,nx), 107 4 VVHF(ny,nx),Xdd(0:max(1,d(1)),nx),Pdd(0:max(1,d(1)),nx,nx), 108 5 WORK(3*nx),LAM(nx),AUX,U 109 INTEGER, ALLOCATABLE:: IRANK(:) 110 DOUBLE PRECISION, ALLOCATABLE:: DLL(:),PROB(:),PROBL(:), 111 1 XT(:,:),PT(:,:,:),PMUL(:,:,:) 112 DOUBLE PRECISION EPS,ONE,ZERO 113 DATA EPS/1.D-14/,ONE/1.0D0/,ZERO/0.0D0/ 114 DOUBLE PRECISION genunf,LEMMA4,MARKOVP 115 116 pdesign = getprocaddress(pdll, "design_"C) 117 CALL DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R) 118 CALL DESIGNZ(nv,np,psi,INFOS,P1,P2,P3,P4,P5,P6) 119C PALL(i,j) = Pr[Z(t+1)=i|Z(t)=j], Z = S1 x S2 x ... x Snv 120 CALL PPROD(nv,nstot,INFOS,P1,P2,P3,P4,P5,P6,PMAT) 121C ERGODIC solves PE: PE*(I-P') = 0 122 CALL ERGODIC(nstot,PMAT,PE) 123 124C OMEGA and MU RECURSIONS 125 OM(:,:,:)= ZERO 126 MU(:,:) = ZERO 127 DO 250 IT = nobs-1,1,-1 128C INT2SEQ map from Z(IT+1) to IS = (k1 k2 k3 k4 k5 k6) 129 CALL INT2SEQ(Z(IT+1),nv,INFOS,SEQ,IS) 130 iny = IYK(IT+1,ny+1) 131 132 DO 10 I=1,iny 133 Ha(I) = SUM(H(IYK(IT+1,I),:,IS(2))*a(:,IS(4))) ! H*a (ny x 1) 134 DO 10 J=1,nu 13510 HRG(I,J) = SUM(H(IYK(IT+1,I),:,IS(2))*R(:,J,IS(6))) 136 + + G(IYK(IT+1,I),J,IS(3)) ! HR+G (ny x nu) 137 138 DO 20 I=1,iny 139 DO 20 J=1,iny 14020 VV(I,J) = SUM(HRG(I,1:nu)*HRG(J,1:nu)) ! (HR+G)*(HR+G)' (ny x ny) 141 142 DO 30 I=1,iny 143 DO 30 J=1,nx 14430 HF(I,J)=SUM(H(IYK(IT+1,I),:,IS(2))*F(:,J,IS(5))) ! HF(ny x nx) 145 146 COM(1:iny,1:iny) = VV(1:iny,1:iny) 147 IFAIL = -1 148C CALL F01ADF(iny,COM(1:iny+1,1:iny), iny+1, IFAIL) 149 CALL DPOTRF('L',iny,COM(1:iny,1:iny),iny,IFAIL) ! COM = L*L' 150 CALL DPOTRI('L',iny,COM(1:iny,1:iny),iny,IFAIL) ! COM = VV^-1 151 DO 40 I=1,iny 152 DO 40 J=1,I 153 VV(I,J) = COM(I,J) 15440 VV(J,I) = VV(I,J) ! inv[(HR+G)*(HR+G)'] (ny x ny) 155 156C B = R*(H*R+G)'*VV (nx x ny) 157 DO 50 I=1,nu 158 DO 50 J=1,iny 15950 HRGV(I,J) = SUM(HRG(1:iny,I)*VV(1:iny,J)) ! (H*R+G)'*VV (nu x ny) 160 161 DO 60 I=1,nx 162 DO 60 J=1,iny 16360 B(I,J) = SUM(R(I,1:nu,IS(6))*HRGV(1:nu,J)) ! B (nx x ny) 164 165 DO 70 I=1,nx 166 DO 70 J=1,nx 167 RR(I,J) = SUM(R(I,:,IS(6))*R(J,:,IS(6))) ! RR' (nx x nx) 16870 BH(I,J) = SUM(B(I,1:iny)*H(IYK(IT+1,1:iny),J,IS(2))) ! BH (nx x nx) 169 170C FIND CS such that CS*CS' = RR'-B*HRG*R' (nx x nx) 171 DO 80 I=1,nx 172 DO 80 J=1,nu 17380 BHRG(I,J) = SUM(B(I,1:iny)*HRG(1:iny,J)) 174 175 DO 90 I=1,nx 176 DO 90 J=1,I 17790 CC(I,J) = RR(I,J) - SUM(BHRG(I,1:nu)*R(J,1:nu,IS(6))) 178 179 IFAIL=-1 180C CALL F02FAF('V','L',nx,CC,nx,LAM,WORK,3*nx,IFAIL) 181 CALL DSYEV( 'V','L',nx,CC,nx,LAM,WORK,3*nx,IFAIL) 182 DO 100 I=1,nx 183 IF (LAM(I).LE.EPS) LAM(I)= ZERO 184100 CS(:,I) = CC(1:nx,I)*DSQRT(LAM(I)) 185 186C AA = F - B*HF (nx x nx) 187 DO 110 I=1,nx 188 DO 110 J=1,nx 189110 AA(I,J) = F(I,J,IS(5)) - SUM(B(I,1:iny)*HF(1:iny,J)) 190 191C OMC = OM(+1)*CS (nx x nx) 192 DO 120 I=1,nx 193 DO 120 J=1,nx 194120 OMC(I,J) = SUM(OM(IT+1,I,:)*CS(:,J)) 195 196C DD = I + CS'*OM(+1)*CS (nx x nx) 197 DD(:,:) = ZERO 198 DO 130 I=1,nx 199 DD(I,I) = ONE 200 DO 130 J=1,I 201 DD(I,J) = DD(I,J) + SUM(CS(:,I)*OMC(:,J)) 202130 DD(J,I) = DD(I,J) 203 204C DI = inv(DD) (nx x nx) 205 COM1(1:nx,1:nx) = DD(1:nx,1:nx) 206 IFAIL = -1 207C CALL F01ADF(nx,COM1,nx+1,IFAIL) 208 CALL DPOTRF('L',nx,COM1(1:nx,1:nx),nx,IFAIL) ! COM1 = L*L' 209 CALL DPOTRI('L',nx,COM1(1:nx,1:nx),nx,IFAIL) ! COM1 = DD^-1 210 211 DO 135 I=1,nx 212 DO 135 J=1,I 213 DI(I,J) = COM1(I,j) 214135 DI(J,I) = DI(I,J) 215 216C OMCDIC = I - OM(+1)*CS*DI*CS' (nx x nx) 217 DO 140 I=1,nx 218 DO 140 J=1,nx 219140 COM1(I,J) = SUM(OMC(I,:)*DI(:,J)) ! OM(+1)*CS*DI (nx x nx) 220 221 OMCDIC(:,:) = ZERO 222 DO 145 I=1,nx 223 OMCDIC(I,I) = ONE 224 DO 145 J=1,nx 225145 OMCDIC(I,J) = OMCDIC(I,J)-SUM(COM1(I,:)*CS(J,:)) 226 227C AOMCDIC = AA'*(I - OM(+1)*CS*DINV CS') (nx x nx) 228 DO 150 I=1,nx 229 DO 150 J=1,nx 230150 AOMCDIC(I,J) = SUM(AA(:,I)*OMCDIC(:,J)) 231 232C AOMCDICOM = AA'*(I - OM(+1)*CS*DINV*CS')*OM(+1) (nx x nx) 233 DO 160 I=1,nx 234 DO 160 J=1,nx 235160 AOMCDICOM(I,J) = SUM(AOMCDIC(I,:)*OM(IT+1,:,J)) 236 237C VV*H*F (ny x nx) 238 DO 170 I=1,iny 239 DO 170 J=1,nx 240170 VVHF(I,J) = SUM(VV(I,1:iny)*HF(1:iny,J)) 241 242C OM = AA*(I - OM(+1)*C*DI*C')*OM(+1)*AA' + 243C + F'*H'*VV*H*F 244 DO 180 I=1,nx 245 DO 180 J=1,nx 246180 OM(IT,I,J) = SUM(AOMCDICOM(I,:)*AA(:,J)) 247 + + SUM(HF(1:iny,I)*VVHF(1:iny,J)) 248 249C MU = AA'*(I - OM(+1)*C*DI* C')*MU(+1) + 250C - AA'*(I - OM C DINV C')*OM(+1)*LAM 251C + F'*H'*VV*(y(+1) - H*a - c*z) 252C LAM = a - B*H*a + B*[y(+1)-c*z] (nx x 1) 253 COM(1:iny,1) = 0.D0 254 DO 185 I=1,iny 255185 COM(I,1) = SUM(c(IYK(IT+1,I),1:nz,IS(1))*yk(IT+1,ny+1:ny+nz)) 256 257 DO 190 I=1,nx 258190 LAM(I) = a(I,IS(4)) - SUM(BH(I,1:nx)*a(1:nx,IS(4))) 259 + + SUM(B(I,1:iny)*(yk(IT+1,IYK(IT+1,1:iny)) 260 + - COM(1:iny,1))) 261 DO 200 I=1,nx 262200 MU(IT,I) = SUM(AOMCDIC(I,:)*MU(IT+1,:)) 263 + - SUM(AOMCDICOM(I,:)*LAM(:)) 264 + + SUM(VVHF(1:iny,I)*(yk(IT+1,IYK(IT+1,1:iny)) 265 # - Ha(1:iny)-COM(1:iny,1))) 266 267250 CONTINUE 268 269C --------------- 270C START SAMPLING 271C --------------- 272 ALLOCATE(DLL(nstot),PROB(nstot),PROBL(nstot),XT(0:nstot,nx), 273 1 PT(0:nstot,nx,nx),IRANK(nstot),PMUL(nstot,nstot,2)) 274 PMUL(:,:,1) = PMAT(:,:) ! one-step ahead 275 PMUL(:,:,2) = 0.D0 ! two-step ahead 276 DO 260 I = 1,nstot 277 DO 260 J = 1,nstot 278 DO 260 K = 1,nstot 279260 PMUL(I,J,2) = PMUL(I,J,2) + PMAT(I,K)*PMAT(K,J) 280 281C FEASIBLE Z-STATES 282 NIFS = 0 283 IFS(:) = 0 284 DO 265 K =1,nstot 285 IF (PE(K).GT.0.D0) THEN 286 NIFS = NIFS + 1 287 IFS(NIFS) = K 288265 ENDIF 289 dc(1:2) = 0 290 DO 2000 IT = 1,nobs 291 DO 1000 KKK = 1,NIFS 292 K = IFS(KKK) 293 CALL INT2SEQ(K,nv,INFOS,SEQ,IS(:)) 294 IF ((IT.LE.d(1)).AND.(d(1).GT.0)) THEN 295 DO 300 I = 1,d(1) 296300 CALL INT2SEQ(Z(I),nv,INFOS,SEQ,IND(I,:)) 297 IND(IT,:)= IS(:) 298 CALL IKF(d,ny,nz,nx,nu,ns,IND(1:d(1),:),yk(1:d(1),:), 299 1 IYK(1:d(1),:),c,H,G,a,F,R,Xdd(1:d(1),:), 300 2 Pdd(1:d(1),:,:),DLL(1:max(1,d(1)))) 301 XT(K,:) = Xdd(IT,:) ! xi(t|t) 302 PT(K,:,:) = Pdd(IT,:,:) ! P(t|t) 303 DLL(:) = ZERO ! log likelihood 304 ELSEIF ((IT.GT.d(1)).AND.(d(1).GT.0)) THEN 305C Input XT(0) PT(0), Output XT(K),PT(K),DLL(K) 306 Xdd(0,:) = XT(0,:) 307 Pdd(0,:,:) = PT(0,:,:) 308 CALL KF(1,dc,ny,nz,nx,nu,ns,IS,yk(IT,:),IYK(IT,:),c,H,G,a,F,R, 309 1 Xdd(0:1,:),Pdd(0:1,:,:),DLL(K)) 310 XT(K,:) = Xdd(1,:) 311 PT(K,:,:) = Pdd(1,:,:) 312 ELSEIF ((IT.EQ.1).AND.(d(1).EQ.0)) THEN 313 CALL IKF(d,ny,nz,nx,nu,ns,IS(:),yk(1,:),IYK(1,:),c,H,G,a,F,R, 314 1 Xdd(0,:),Pdd(0,:,:),DLL(K)) 315 CALL KF(1,dc,ny,nz,nx,nu,ns,IS(:),yk(1,:),IYK(1,:),c,H,G,a, 316 1 F,R,Xdd(0:1,:),Pdd(0:1,:,:),DLL(K)) ! log likelihood 317 XT(K,:) = Xdd(IT,:) ! xi(t|t) 318 PT(K,:,:) = Pdd(IT,:,:) ! P(t|t) 319 ELSEIF ((IT.GT.1).AND.(d(1).EQ.0)) THEN 320C Input XT(0) PT(0), Output XT(K),PT(K),DLL(K) 321 Xdd(0,:) = XT(0,:) 322 Pdd(0,:,:) = PT(0,:,:) 323 CALL KF(1,dc,ny,nz,nx,nu,ns,IS(:),yk(IT,:),IYK(IT,:),c,H,G,a, 324 1 F,R,Xdd(0:1,:),Pdd(0:1,:,:),DLL(K)) 325 XT(K,:) = Xdd(1,:) 326 PT(K,:,:) = Pdd(1,:,:) 327 ENDIF 328 329 SH(1) = Z(max(1,IT-1)) 330 SH(2) = K 331 SH(3) = Z(min(nobs,IT+1)) 332 PROBL(K) = DLL(K) 333 + + LEMMA4(OM(IT,:,:),MU(IT,:),XT(K,:),PT(K,:,:),nx) 334 + + MARKOVP(PMUL,PE,nstot,1,IT,nobs,SH) 335 3361000 CONTINUE 337 338C --------------------------------------------- 339C SAMPLING Z(t:t+h-1) using PROB 340C ISEQ is the sampled sequence - out of nstot 341C --------------------------------------------- 342C To prevent exp overflow 343 PROB(:) = 0.D0 344 IMAX = MAXLOC(PROBL(IFS(1:NIFS))) 345 KKK = IFS(IMAX(1)) 346 PROBL(IFS(1:NIFS)) = PROBL(IFS(1:NIFS))-PROBL(KKK) 347 PROB(IFS(1:NIFS)) = DEXP(PROBL(IFS(1:NIFS))) 348 # / SUM(DEXP(PROBL(IFS(1:NIFS)))) 349 350C U = G05CAF(U) ! Sampling from U(0,1) 351 U = genunf(0.D0,1.D0) 352 ISEQ = 1 353 AUX = PROB(1) 354 DO 310 WHILE (AUX.LT.U) 355 ISEQ = ISEQ + 1 356310 AUX = AUX + PROB(ISEQ) 357 358 Z(IT) = ISEQ 359 360 XT(0,:) = XT(ISEQ,:) 361 PT(0,:,:) = PT(ISEQ,:,:) 362 3632000 CONTINUE 364 365 DO I=1,nobs 366 CALL INT2SEQ(Z(I),nv,INFOS,SEQ,S(I,:)) 367 ENDDO 368 369 DEALLOCATE(DLL,PROB,PROBL,XT,PT,IRANK,PMUL) 370 RETURN 371 END 372 373C **** da butta ****************** 374c DO IT = 1,nobs-1 375c CALL INT2SEQ(Z(IT),nv,INFOS,SEQ,SS(IT,:)) 376c ENDDO 377 378c SS(nobs,:) = 1 379c CALL IKF(d,ny,nz,nx,nu,ns,SS(1:max(d(1),1),1:6), 380c 1 yk(1:max(d(1),1),1:ny+nz),IYK(1:max(d(1),1),1:ny+1), 381c 2 c,H,G,a,F,R,Xdd,Pdd,LIKE(1:max(d(1),1))) 382c XX(d(1),1:nx) = Xdd(max(d(1),1),1:nx) 383c PP(d(1),1:nx,1:nx) = Pdd(max(d(1),1),1:nx,1:nx) 384c CALL KF(nobs,d,ny,nz,nx,nu,ns,SS,yk,IYK,c,H,G,a,F,R,XX,PP,LIKE) 385 386c R1 = DEXP(SUM(LIKE))*P1(SS(nobs,4),SS(nobs-1,4)) 387 388c SS(nobs,4) = 2 389c CALL IKF(d,ny,nz,nx,nu,ns,SS(1:max(d(1),1),1:6), 390c 1 yk(1:max(d(1),1),1:ny+nz),IYK(1:max(d(1),1),1:ny+1), 391c 2 c,H,G,a,F,R,Xdd,Pdd,LIKE(1:max(d(1),1))) 392c XX(d(1),1:nx) = Xdd(max(d(1),1),1:nx) 393c PP(d(1),1:nx,1:nx) = Pdd(max(d(1),1),1:nx,1:nx) 394c CALL KF(nobs,d,ny,nz,nx,nu,ns,SS,yk,IYK,c,H,G,a,F,R,XX,PP,LIKE) 395c R2 = DEXP(SUM(LIKE))*P1(SS(nobs,4),SS(nobs-1,4)) 396c R1 = R1/(R1+R2) 397