1!simulate gaussian state space model 2subroutine simgaussianuncond(timevar, zt, ht, tt, rtv, qt, a1, p1, & 3nnd,nsim, epsplus, etaplus, aplus1, p, n, m, r, info,& 4tol,sim,c,simwhat,simdim,antithetics) 5 6 implicit none 7 8 integer, intent(in) :: p, m, r, n, nsim,nnd,simdim,simwhat,antithetics 9 integer, intent(in), dimension(5) :: timevar 10 integer, intent(inout) :: info 11 integer :: t, i,k,l 12 double precision, intent(in) :: tol 13 double precision, intent(in), dimension(p,m,(n-1)*timevar(1)+1) :: zt 14 double precision, intent(in), dimension(p,p,(n-1)*timevar(2)+1) :: ht 15 double precision, intent(in), dimension(m,m,(n-1)*timevar(3)+1) :: tt 16 double precision, intent(in), dimension(m,r,(n-1)*timevar(4)+1) :: rtv 17 double precision, intent(in), dimension(r,r,(n-1)*timevar(5)+1) :: qt 18 double precision, intent(in), dimension(m) :: a1 19 double precision, intent(in), dimension(m,m) :: p1 20 double precision, intent(in),dimension(nsim) :: c 21 double precision, intent(inout), dimension(simdim,n,3 * nsim * antithetics + nsim) :: sim 22 double precision, intent(inout), dimension(r,n,nsim) :: etaplus 23 double precision, intent(inout), dimension(p,n,nsim) :: epsplus 24 double precision, intent(inout), dimension(m,nsim) :: aplus1 25 26 double precision, dimension(m,n+1) :: aplus 27 double precision, dimension(r,r,(n-1)*timevar(5)+1) :: cholqt 28 double precision, dimension(m,m) :: cholp1 29 double precision, dimension(r,r) :: rcholtmp 30 double precision, dimension(r,n-1,4) :: etatmp 31 double precision, dimension(m,n,4) :: alphatmp 32 double precision, external :: ddot 33 34 external dsymm, dgemm, dsymv, ldl, dtrmv, dgemv 35 36 37 38 do t = 1, (n-1)*timevar(5)+1 39 if(r.EQ.1) then 40 cholqt(1,1,t)=sqrt(qt(1,1,t)) 41 else 42 rcholtmp = qt(:,:,t) 43 call ldl(rcholtmp,r,tol,info) 44 if(info .NE. 0) then 45 info = -2 46 return 47 end if 48 do i=1,r 49 cholqt(i,i,t)=sqrt(rcholtmp(i,i)) 50 end do 51 do i=1,r-1 52 cholqt((i+1):r,i,t) = rcholtmp((i+1):r,i)*cholqt(i,i,t) 53 end do 54 end if 55 end do 56 57 58 if(nnd.GT.0) then 59 if(m.EQ.1) then 60 cholp1(1,1)=sqrt(p1(1,1)) 61 else 62 cholp1 = p1 63 call ldl(cholp1,m,tol,info) 64 if(info .NE. 0) then 65 info=-3 66 return 67 end if 68 do i=1,m 69 cholp1(i,i)=sqrt(cholp1(i,i)) 70 end do 71 do i=1,m-1 72 cholp1((i+1):m,i) = cholp1((i+1):m,i)*cholp1(i,i) 73 end do 74 end if 75 end if 76 77 78 do i = 1, nsim 79 aplus = 0.0d0 80 aplus(:,1) = a1 81 if(nnd.GT.0) then 82 call dtrmv('l','n','n',m,cholp1,m,aplus1(:,i),1) 83 aplus(:,1) = aplus(:,1)+aplus1(:,i) 84 end if 85 86 do t = 1, n 87 do k = 1, p 88 epsplus(k,t,i) = epsplus(k,t,i)*sqrt(ht(k,k,(t-1)*timevar(2)+1)) 89 end do 90 call dtrmv('l','n','n',r,cholqt(:,:,(t-1)*timevar(5)+1),r,etaplus(:,t,i),1) 91 call dgemv('n',m,m,1.0d0,tt(:,:,(t-1)*timevar(3)+1),m,aplus(:,t),1,0.0d0,aplus(:,t+1),1) 92 call dgemv('n',m,r,1.0d0,rtv(:,:,(t-1)*timevar(4)+1),m,etaplus(:,t,i),1,1.0d0,aplus(:,t+1),1) 93 end do 94 95 96 !simwhat = 1: epsilon, 2: eta, 3: both, 4: state, 5: signal, 6: observations 97 select case(simwhat) 98 case(1) 99 sim(:,:,i) = epsplus(:,:,i) 100 if(antithetics .EQ. 1) then 101 sim(:,:,i+nsim) = -epsplus(:,:,i) 102 sim(:,:,i+2*nsim) = c(i)*sim(:,:,i) 103 sim(:,:,i+3*nsim) = c(i)*sim(:,:,i+nsim) 104 end if 105 case(2) 106 sim(:,:,i) = etaplus(:,:,i) 107 if(antithetics .EQ. 1) then 108 sim(:,:,i+nsim) = -etaplus(:,:,i) 109 sim(:,:,i+2*nsim) = c(i)*sim(:,:,i) 110 sim(:,:,i+3*nsim) = c(i)*sim(:,:,i+nsim) 111 end if 112 case(3) 113 sim(1:p,:,i) = epsplus(:,:,i) 114 sim((p+1):,:,i) = etaplus(:,:,i) 115 if(antithetics .EQ. 1) then 116 sim(1:p,:,i+nsim) = -epsplus(:,:,i) 117 sim(1:p,:,i+2*nsim) = c(i)*sim(1:p,:,i) 118 sim(1:p,:,i+3*nsim) = c(i)*sim(1:p,:,i+nsim) 119 sim((p+1):,:,i+nsim) = -etaplus(:,:,i) 120 sim((p+1):,:,i+2*nsim) = c(i)*sim((p+1):,:,i) 121 sim((p+1):,:,i+3*nsim) = c(i)*sim((p+1):,:,i+nsim) 122 end if 123 case(4) 124 sim(:,1,i) = aplus(:,1) 125 etatmp(:,:,1) = etaplus(:,1:(n-1),i) 126 if(antithetics .EQ. 1) then 127 sim(:,1,i+nsim) = -aplus(:,1) 128 sim(:,1,i+2*nsim) = c(i)*sim(:,1,i) 129 sim(:,1,i+3*nsim) = c(i)*sim(:,1,i+nsim) 130 131 etatmp(:,:,2) = -etaplus(:,1:(n-1),i) 132 etatmp(:,:,3) = c(i)*etatmp(:,:,1) 133 etatmp(:,:,4) = c(i)*etatmp(:,:,2) 134 end if 135 do k = 1, 3*antithetics+1 136 do t = 2, n 137 call dgemv('n',m,m,1.0d0,tt(:,:,(t-2)*timevar(3)+1),m,sim(:,t-1,i+(k-1)*nsim),& 138 1,0.0d0,sim(:,t,i+(k-1)*nsim),1) 139 call dgemv('n',m,r,1.0d0,rtv(:,:,(t-2)*timevar(4)+1),m,etatmp(:,t-1,k),1,& 140 1.0d0,sim(:,t,i+(k-1)*nsim),1) 141 end do 142 end do 143 case(5) 144 145 alphatmp(:,1,1) = aplus(:,1) 146 etatmp(:,:,1) = etaplus(:,1:(n-1),i) 147 call dgemv('n',p,m,1.0d0,zt(:,:,1),p,alphatmp(:,1,1),1,0.0d0,sim(:,1,i),1) 148 149 if(antithetics .EQ. 1) then 150 alphatmp(:,1,2) = -aplus(:,1) 151 alphatmp(:,1,3) = c(i)*alphatmp(:,1,1) 152 alphatmp(:,1,4) = c(i)*alphatmp(:,1,2) 153 call dgemv('n',p,m,1.0d0,zt(:,:,1),p,alphatmp(:,1,2),1,0.0d0,sim(:,1,i+nsim),1) 154 call dgemv('n',p,m,1.0d0,zt(:,:,1),p,alphatmp(:,1,3),1,0.0d0,sim(:,1,i+2*nsim),1) 155 call dgemv('n',p,m,1.0d0,zt(:,:,1),p,alphatmp(:,1,4),1,0.0d0,sim(:,1,i+3*nsim),1) 156 157 etatmp(:,:,2) = -etaplus(:,1:(n-1),i) 158 etatmp(:,:,3) = c(i)*etatmp(:,:,1) 159 etatmp(:,:,4) = c(i)*etatmp(:,:,2) 160 161 end if 162 do k = 1, 3*antithetics+1 163 do t = 2, n 164 call dgemv('n',m,m,1.0d0,tt(:,:,(t-2)*timevar(3)+1),m,alphatmp(:,t-1,k),& 165 1,0.0d0,alphatmp(:,t,k),1) 166 call dgemv('n',m,r,1.0d0,rtv(:,:,(t-2)*timevar(4)+1),m,etatmp(:,t-1,k),1,& 167 1.0d0,alphatmp(:,t,k),1) 168 call dgemv('n',p,m,1.0d0,zt(:,:,(t-1)*timevar(1)+1),p,& 169 alphatmp(:,t,k),1,0.0d0,sim(:,t,i+(k-1)*nsim),1) 170 end do 171 end do 172 case(6) 173 174 175 alphatmp(:,1,1) = aplus(:,1) 176 etatmp(:,:,1) = etaplus(:,1:(n-1),i) 177 178 if(antithetics .EQ. 1) then 179 alphatmp(:,1,2) = -aplus(:,1) 180 alphatmp(:,1,3) = c(i)*alphatmp(:,1,1) 181 alphatmp(:,1,4) = c(i)*alphatmp(:,1,2) 182 183 etatmp(:,:,2) = -etaplus(:,1:(n-1),i) 184 etatmp(:,:,3) = c(i)*etatmp(:,:,1) 185 etatmp(:,:,4) = c(i)*etatmp(:,:,2) 186 187 end if 188 do k = 1, 3*antithetics+1 189 do l = 1, p 190 sim(l,1,i+(k-1)*nsim) = ddot(m,zt(l,:,1),1,alphatmp(:,1,k),1) 191 192 end do 193 do t = 2, n 194 call dgemv('n',m,m,1.0d0,tt(:,:,(t-2)*timevar(3)+1),m,alphatmp(:,t-1,k),& 195 1,0.0d0,alphatmp(:,t,k),1) 196 call dgemv('n',m,r,1.0d0,rtv(:,:,(t-2)*timevar(4)+1),m,etatmp(:,t-1,k),1,& 197 1.0d0,alphatmp(:,t,k),1) 198 do l = 1, p 199 200 sim(l,t,i+(k-1)*nsim) = ddot(m,zt(l,:,(t-1)*timevar(1)+1),1,alphatmp(:,t,k),1) 201 202 end do 203 end do 204 end do 205 end select 206 end do 207 208end subroutine simgaussianuncond 209 210