1!simulate gaussian state space model 2subroutine simgaussian(ymiss,timevar, yt, zt, ht, tt, rtv, qt, a1, p1, & 3p1inf, nnd,nsim, epsplus, etaplus, aplus1, p, n, m, r, info,rankp,& 4tol,sim,c,simwhat,simdim,antithetics) 5 6 implicit none 7 8 integer, intent(in) :: p, m, r, n, nsim,nnd,simdim,simwhat,antithetics,rankp 9 integer, intent(in), dimension(n,p) :: ymiss 10 integer, intent(in), dimension(5) :: timevar 11 integer, intent(inout) :: info 12 integer :: t, i, d, j,k,tv,l,rankp2 13 double precision, intent(in) :: tol 14 double precision, intent(in), dimension(n,p) :: yt 15 double precision, intent(in), dimension(p,m,(n-1)*timevar(1)+1) :: zt 16 double precision, intent(in), dimension(p,p,(n-1)*timevar(2)+1) :: ht 17 double precision, intent(in), dimension(m,m,(n-1)*timevar(3)+1) :: tt 18 double precision, intent(in), dimension(m,r,(n-1)*timevar(4)+1) :: rtv 19 double precision, intent(in), dimension(r,r,(n-1)*timevar(5)+1) :: qt 20 double precision, intent(in), dimension(m) :: a1 21 double precision, intent(in), dimension(m,m) :: p1,p1inf 22 double precision, intent(in),dimension(nsim) :: c 23 double precision, intent(inout), dimension(simdim,n,3 * nsim * antithetics + nsim) :: sim 24 double precision, intent(inout), dimension(r,n,nsim) :: etaplus 25 double precision, intent(inout), dimension(p,n,nsim) :: epsplus 26 double precision, intent(inout), dimension(m,nsim) :: aplus1 27 28 double precision, dimension(n,p) :: yplus 29 double precision, dimension(m,n+1) :: aplus 30 double precision, dimension(m) :: ahat 31 double precision, dimension(m) :: aplushat 32 double precision, dimension(p,n) :: ft,finf 33 double precision, dimension(m,p,n) :: kt,kinf 34 double precision, dimension(r,r,(n-1)*timevar(5)+1) :: cholqt 35 double precision, dimension(m,m) :: cholp1 36 double precision, dimension(r,r) :: rcholtmp 37 double precision, dimension(m) :: rt0,rt1 38 double precision, dimension(m,r) :: mr 39 double precision, dimension(m,m,(n-1)*max(timevar(4),timevar(5))+1) :: rqr 40 double precision, dimension(r,n) :: etahat 41 double precision, dimension(p,n) :: epshat 42 double precision, dimension(r,n) :: etaplushat 43 double precision, dimension(p,n) :: epsplushat 44 double precision, dimension(r,n-1,4) :: etatmp 45 double precision, dimension(m,n,4) :: alphatmp 46 double precision, external :: ddot 47 logical needeps 48 49 external dsymm, dgemm, smoothsim, dsymv, ldl, dtrmv,smoothsimfast, dgemv 50 51 52 53 needeps = simwhat.EQ.1 .or. simwhat.EQ.3 54 if(needeps) then 55 epshat = 0.0d0 56 end if 57 58 !compute rqr 59 tv= max(timevar(4),timevar(5)) 60 do t=1, (n-1)*tv+1 61 call dsymm('r','u',m,r,1.0d0,qt(:,:,(t-1)*timevar(5)+1),r,rtv(:,:,(t-1)*timevar(4)+1),m,0.0d0,mr,m) 62 call dgemm('n','t',m,m,r,1.0d0,mr,m,rtv(:,:,(t-1)*timevar(4)+1),m,0.0d0,rqr(:,:,t),m) 63 end do 64 65 66 67 rankp2 = rankp 68 call smoothsim(yt, ymiss, timevar, zt, ht,tt, rtv,qt,rqr, a1, p1, p1inf, & 69 d, j, p, m, n, r,tol,rankp2,ft,finf,kt,kinf,epshat,etahat,rt0,rt1,needeps) 70 !simwhat = 1: epsilon, 2: eta, 3: both, 4: state, 5: signal, 6: observations 71 if(simwhat .GT. 3) then 72 ahat = a1 73 call dsymv('l',m,1.0d0,p1,m,rt0,1,1.0d0,ahat,1) 74 if(d .GT. 0) then 75 call dsymv('l',m,1.0d0,p1inf,m,rt1,1,1.0d0,ahat,1) 76 end if 77 end if 78 79 80 81 do t = 1, (n-1)*timevar(5)+1 82 if(r.EQ.1) then 83 cholqt(1,1,t)=sqrt(qt(1,1,t)) 84 else 85 rcholtmp = qt(:,:,t) 86 call ldl(rcholtmp,r,tol,info) 87 if(info .NE. 0) then 88 info = -2 89 return 90 end if 91 do i=1,r 92 cholqt(i,i,t)=sqrt(rcholtmp(i,i)) 93 end do 94 do i=1,r-1 95 cholqt((i+1):r,i,t) = rcholtmp((i+1):r,i)*cholqt(i,i,t) 96 end do 97 end if 98 end do 99 100 101 if(nnd.GT.0) then 102 if(m.EQ.1) then 103 cholp1(1,1)=sqrt(p1(1,1)) 104 else 105 cholp1 = p1 106 call ldl(cholp1,m,tol,info) 107 if(info .NE. 0) then 108 info=-3 109 return 110 end if 111 do i=1,m 112 cholp1(i,i)=sqrt(cholp1(i,i)) 113 end do 114 do i=1,m-1 115 cholp1((i+1):m,i) = cholp1((i+1):m,i)*cholp1(i,i) 116 end do 117 end if 118 end if 119 120 do i = 1, nsim 121 aplus = 0.0d0 122 aplus(:,1) = a1 123 if(nnd.GT.0) then 124 call dtrmv('l','n','n',m,cholp1,m,aplus1(:,i),1) 125 aplus(:,1) = aplus(:,1)+aplus1(:,i) 126 end if 127 128 do t = 1, n 129 do k = 1, p 130 epsplus(k,t,i) = epsplus(k,t,i)*sqrt(ht(k,k,(t-1)*timevar(2)+1)) 131 if(ymiss(t,k).EQ.0) then 132 yplus(t,k) = epsplus(k,t,i) + ddot(m,zt(k,:,(t-1)*timevar(1)+1),1,aplus(:,t),1) 133 end if 134 end do 135 call dtrmv('l','n','n',r,cholqt(:,:,(t-1)*timevar(5)+1),r,etaplus(:,t,i),1) 136 call dgemv('n',m,m,1.0d0,tt(:,:,(t-1)*timevar(3)+1),m,aplus(:,t),1,0.0d0,aplus(:,t+1),1) 137 call dgemv('n',m,r,1.0d0,rtv(:,:,(t-1)*timevar(4)+1),m,etaplus(:,t,i),1,1.0d0,aplus(:,t+1),1) 138 end do 139 140 call smoothsimfast(yplus, ymiss, timevar, zt, ht,tt, rtv,qt,a1, ft,kt,& 141 finf, kinf, d, j, p, m, n,r,epsplushat,etaplushat,rt0,rt1,needeps) 142 143 !simwhat = 1: epsilon, 2: eta, 3: both, 4: state, 5: signal, 6: observations 144 select case(simwhat) 145 case(1) 146 sim(:,:,i) = epshat - epsplushat + epsplus(:,:,i) 147 if(antithetics .EQ. 1) then 148 sim(:,:,i+nsim) = epshat + epsplushat - epsplus(:,:,i) 149 sim(:,:,i+2*nsim) = epshat + c(i)*(sim(:,:,i)-epshat) 150 sim(:,:,i+3*nsim) = epshat + c(i)*(sim(:,:,i+nsim)-epshat) 151 end if 152 case(2) 153 sim(:,:,i) = etahat - etaplushat + etaplus(:,:,i) 154 if(antithetics .EQ. 1) then 155 sim(:,:,i+nsim) = etahat + etaplushat - etaplus(:,:,i) 156 sim(:,:,i+2*nsim) = etahat + c(i)*(sim(:,:,i)-etahat) 157 sim(:,:,i+3*nsim) = etahat + c(i)*(sim(:,:,i+nsim)-etahat) 158 end if 159 case(3) 160 sim(1:p,:,i) = epshat - epsplushat + epsplus(:,:,i) 161 sim((p+1):,:,i) = etahat - etaplushat + etaplus(:,:,i) 162 if(antithetics .EQ. 1) then 163 sim(1:p,:,i+nsim) = epshat + epsplushat - epsplus(:,:,i) 164 sim(1:p,:,i+2*nsim) = epshat + c(i)*(sim(1:p,:,i)-epshat) 165 sim(1:p,:,i+3*nsim) = epshat + c(i)*(sim(1:p,:,i+nsim)-epshat) 166 sim((p+1):,:,i+nsim) = etahat + etaplushat - etaplus(:,:,i) 167 sim((p+1):,:,i+2*nsim) = etahat + c(i)*(sim((p+1):,:,i)-etahat) 168 sim((p+1):,:,i+3*nsim) = etahat + c(i)*(sim((p+1):,:,i+nsim)-etahat) 169 end if 170 case(4) 171 aplushat = a1 172 173 call dsymv('l',m,1.0d0,p1,m,rt0,1,1.0d0,aplushat,1) 174 if(d .GT. 0) then 175 call dsymv('l',m,1.0d0,p1inf,m,rt1,1,1.0d0,aplushat,1) 176 end if 177 178 sim(:,1,i) = ahat - aplushat + aplus(:,1) 179 etatmp(:,:,1) = etahat(:,1:(n-1)) - etaplushat(:,1:(n-1)) + etaplus(:,1:(n-1),i) 180 if(antithetics .EQ. 1) then 181 sim(:,1,i+nsim) = ahat + aplushat - aplus(:,1) 182 sim(:,1,i+2*nsim) = ahat+ c(i)*(sim(:,1,i)-ahat) 183 sim(:,1,i+3*nsim) = ahat+ c(i)*(sim(:,1,i+nsim)-ahat) 184 185 etatmp(:,:,2) = etahat(:,1:(n-1)) + etaplushat(:,1:(n-1)) - etaplus(:,1:(n-1),i) 186 etatmp(:,:,3) = etahat(:,1:(n-1)) + c(i)*(etatmp(:,:,1)-etahat(:,1:(n-1))) 187 etatmp(:,:,4) = etahat(:,1:(n-1)) + c(i)*(etatmp(:,:,2)-etahat(:,1:(n-1))) 188 end if 189 do k = 1, 3*antithetics+1 190 do t = 2, n 191 call dgemv('n',m,m,1.0d0,tt(:,:,(t-2)*timevar(3)+1),m,sim(:,t-1,i+(k-1)*nsim),& 192 1,0.0d0,sim(:,t,i+(k-1)*nsim),1) 193 call dgemv('n',m,r,1.0d0,rtv(:,:,(t-2)*timevar(4)+1),m,etatmp(:,t-1,k),1,& 194 1.0d0,sim(:,t,i+(k-1)*nsim),1) 195 end do 196 end do 197 case(5) 198 aplushat = a1 199 200 call dsymv('l',m,1.0d0,p1,m,rt0,1,1.0d0,aplushat,1) 201 if(d .GT. 0) then 202 call dsymv('l',m,1.0d0,p1inf,m,rt1,1,1.0d0,aplushat,1) 203 end if 204 205 alphatmp(:,1,1) = ahat - aplushat + aplus(:,1) 206 etatmp(:,:,1) = etahat(:,1:(n-1)) - etaplushat(:,1:(n-1)) + etaplus(:,1:(n-1),i) 207 call dgemv('n',p,m,1.0d0,zt(:,:,1),p,alphatmp(:,1,1),1,0.0d0,sim(:,1,i),1) 208 209 if(antithetics .EQ. 1) then 210 alphatmp(:,1,2) = ahat + aplushat - aplus(:,1) 211 alphatmp(:,1,3) = ahat + c(i)*(alphatmp(:,1,1)-ahat) 212 alphatmp(:,1,4) = ahat + c(i)*(alphatmp(:,1,2)-ahat) 213 call dgemv('n',p,m,1.0d0,zt(:,:,1),p,alphatmp(:,1,2),1,0.0d0,sim(:,1,i+nsim),1) 214 call dgemv('n',p,m,1.0d0,zt(:,:,1),p,alphatmp(:,1,3),1,0.0d0,sim(:,1,i+2*nsim),1) 215 call dgemv('n',p,m,1.0d0,zt(:,:,1),p,alphatmp(:,1,4),1,0.0d0,sim(:,1,i+3*nsim),1) 216 217 etatmp(:,:,2) = etahat(:,1:(n-1)) + etaplushat(:,1:(n-1)) - etaplus(:,1:(n-1),i) 218 etatmp(:,:,3) = etahat(:,1:(n-1)) + c(i)*(etatmp(:,:,1)-etahat(:,1:(n-1))) 219 etatmp(:,:,4) = etahat(:,1:(n-1)) + c(i)*(etatmp(:,:,2)-etahat(:,1:(n-1))) 220 221 end if 222 do k = 1, 3*antithetics+1 223 do t = 2, n 224 call dgemv('n',m,m,1.0d0,tt(:,:,(t-2)*timevar(3)+1),m,alphatmp(:,t-1,k),& 225 1,0.0d0,alphatmp(:,t,k),1) 226 call dgemv('n',m,r,1.0d0,rtv(:,:,(t-2)*timevar(4)+1),m,etatmp(:,t-1,k),1,& 227 1.0d0,alphatmp(:,t,k),1) 228 call dgemv('n',p,m,1.0d0,zt(:,:,(t-1)*timevar(1)+1),p,& 229 alphatmp(:,t,k),1,0.0d0,sim(:,t,i+(k-1)*nsim),1) 230 end do 231 end do 232 case(6) 233 aplushat = a1 234 235 call dsymv('l',m,1.0d0,p1,m,rt0,1,1.0d0,aplushat,1) 236 if(d .GT. 0) then 237 call dsymv('l',m,1.0d0,p1inf,m,rt1,1,1.0d0,aplushat,1) 238 end if 239 240 alphatmp(:,1,1) = ahat - aplushat + aplus(:,1) 241 etatmp(:,:,1) = etahat(:,1:(n-1)) - etaplushat(:,1:(n-1)) + etaplus(:,1:(n-1),i) 242 243 if(antithetics .EQ. 1) then 244 alphatmp(:,1,2) = ahat + aplushat - aplus(:,1) 245 alphatmp(:,1,3) = ahat+ c(i)*(alphatmp(:,1,1)-ahat) 246 alphatmp(:,1,4) = ahat+ c(i)*(alphatmp(:,1,2)-ahat) 247 248 etatmp(:,:,2) = etahat(:,1:(n-1)) + etaplushat(:,1:(n-1)) - etaplus(:,1:(n-1),i) 249 etatmp(:,:,3) = etahat(:,1:(n-1)) + c(i)*(etatmp(:,:,1)-etahat(:,1:(n-1))) 250 etatmp(:,:,4) = etahat(:,1:(n-1)) + c(i)*(etatmp(:,:,2)-etahat(:,1:(n-1))) 251 252 end if 253 do k = 1, 3*antithetics+1 254 do l = 1, p 255 if(ymiss(1,l).GT.0) then 256 sim(l,1,i+(k-1)*nsim) = ddot(m,zt(l,:,1),1,alphatmp(:,1,k),1) 257 end if 258 end do 259 do t = 2, n 260 call dgemv('n',m,m,1.0d0,tt(:,:,(t-2)*timevar(3)+1),m,alphatmp(:,t-1,k),& 261 1,0.0d0,alphatmp(:,t,k),1) 262 call dgemv('n',m,r,1.0d0,rtv(:,:,(t-2)*timevar(4)+1),m,etatmp(:,t-1,k),1,& 263 1.0d0,alphatmp(:,t,k),1) 264 do l = 1, p 265 if(ymiss(t,l).GT.0) then 266 sim(l,t,i+(k-1)*nsim) = ddot(m,zt(l,:,(t-1)*timevar(1)+1),1,alphatmp(:,t,k),1) 267 end if 268 end do 269 end do 270 end do 271 end select 272 end do 273 274end subroutine simgaussian 275 276