1! Importance sampling filtering of non-gaussian model 2 3subroutine isamplefilter(yt, ymiss, timevar, zt, tt, rtv, qt, a1, p1,p1inf, u, dist, & 4p, n, m, r, theta, maxiter,rankp,convtol, nnd,nsim,epsplus,etaplus,& 5aplus1,c,tol,info,antithetics,w,sim,simwhat,simdim, expected, htol) 6 7 implicit none 8 9 integer, intent(in) :: p,m, r, n,nnd,antithetics,nsim& 10 ,simwhat,simdim,rankp, expected 11 integer, intent(in), dimension(p) :: dist 12 integer, intent(in), dimension(n,p) :: ymiss 13 integer, intent(in), dimension(5) :: timevar 14 integer, intent(inout) ::info, maxiter 15 integer :: t, j,i,k,maxiter2,maxitermax,info2 16 double precision, intent(in) :: convtol,tol 17 double precision, intent(inout) :: htol 18 double precision, intent(in), dimension(n,p) :: u 19 double precision, intent(in), dimension(n,p) :: yt 20 double precision, intent(in), dimension(p,m,(n-1)*timevar(1)+1) :: zt 21 double precision, intent(in), dimension(m,m,(n-1)*timevar(3)+1) :: tt 22 double precision, intent(in), dimension(m,r,(n-1)*timevar(4)+1) :: rtv 23 double precision, intent(in), dimension(r,r,(n-1)*timevar(5)+1) :: qt 24 double precision, intent(in), dimension(m) :: a1 25 double precision, intent(in), dimension(m,m) :: p1,p1inf 26 double precision, intent(in),dimension(nsim) :: c 27 double precision, intent(inout), dimension(p,n,nsim) :: epsplus 28 double precision, intent(inout), dimension(r,n,nsim) :: etaplus 29 double precision, intent(inout), dimension(m,nsim) :: aplus1 30 double precision, intent(inout), dimension(n,p) :: theta 31 double precision, dimension(p,p,n) :: ht 32 double precision, intent(inout), dimension(simdim,n,3 * nsim * antithetics + nsim) :: sim 33 double precision, dimension(p,(3 * nsim * antithetics + nsim)*(5-simwhat)) :: tsim 34 double precision, dimension(n,p) :: ytilde 35 double precision, dimension(n) :: tmp 36 double precision, dimension(n,3 * nsim * antithetics + nsim) :: w 37 double precision, external :: ddot 38 39 integer, dimension(n,p) :: ymiss2 40 double precision, dimension(p,n,nsim) :: epsplus2 41 double precision, dimension(r,n,nsim) :: etaplus2 42 double precision, dimension(m,nsim) :: aplus12 43 double precision, dimension(simdim,n,3 * nsim * antithetics + nsim) :: sim2 44 double precision :: diff 45 double precision :: lik 46 47 external approx, simgaussian 48 49 ht=0.0d0 50 ytilde=0.0d0 51 w=1.0d0 52 53 epsplus2 = epsplus 54 etaplus2 = etaplus 55 aplus12 = aplus1 56 ymiss2 = ymiss 57 ymiss2(1,:) = 1 58 59 call simgaussian(ymiss2(1,:),timevar, ytilde(1,:), zt(:,:,1), & 60 ht(:,:,1), tt(:,:,1), rtv(:,:,1), & 61 qt(:,:,1), a1, p1, p1inf, nnd,nsim, epsplus2(:,1,:), etaplus2(:,1,:), aplus12(:,:), & 62 p, 1, m, r, info,rankp,tol,sim(:,1,:),c,simwhat,simdim,antithetics) 63 64 65 if(info /= 0) then 66 return 67 end if 68 maxitermax = 0 69 70 do i = 1, (n-1) ! increase time 71 72 ht = 0.0d0 73 ytilde = 0.0d0 74 75 maxiter2 = maxiter 76 info2 = 0 77 ! approximate 78 call approx(yt(1:i,:), ymiss(1:i,:), timevar, zt(:,:,1:((i-1)*timevar(1)+1)), & 79 tt(:,:,1:((i-1)*timevar(3)+1)), rtv(:,:,1:((i-1)*timevar(4)+1)), ht(:,:,1:i),& 80 qt(:,:,1:((i-1)*timevar(5)+1)), a1, p1,p1inf, p,i,m,r,& 81 theta(1:i,:), u(1:i,:), ytilde(1:i,:), dist,maxiter2,tol,rankp,convtol,diff,lik,& 82 info2, expected, htol) 83 84 if(info2 .ne. 0 .and. info2 .ne. 3) then !check for errors in approximating algorithm 85 info = info2 86 return 87 end if 88 89 if(maxiter2.GT.maxitermax) then 90 maxitermax = maxiter2 91 end if 92 epsplus2 = epsplus 93 etaplus2 = etaplus 94 aplus12 = aplus1 95 ymiss2 = ymiss 96 ymiss2(i+1,:) = 1 97 sim2=0.0d0 98 info2 = 0 99 ! simulate signals 100 call simgaussian(ymiss2(1:(i+1),:),timevar, ytilde(1:(i+1),:), zt(:,:,1:(i*timevar(1)+1)), & 101 ht(:,:,1:(i+1)), tt(:,:,1:(i*timevar(3)+1)), rtv(:,:,1:(i*timevar(4)+1)), & 102 qt(:,:,1:(i*timevar(5)+1)), a1, p1, p1inf, nnd,nsim, epsplus2(:,1:(i+1),:), & 103 etaplus2(:,1:(i+1),:), aplus12(:,:),p, i+1, m, r, info2,rankp,tol,& 104 sim2(:,1:(i+1),:),c,simwhat,simdim,antithetics) 105 106 if(info2 /= 0) then 107 info = info2 108 return 109 end if 110 111 ! compute importance weights 112 113 if(simwhat.EQ.5) then 114 do j=1,p 115 select case(dist(j)) 116 case(2) !poisson 117 tmp(1:i) = exp(theta(1:i,j)) 118 do t=1,i 119 if(ymiss2(t,j) .EQ. 0) then 120 w(i+1,:) = w(i+1,:)*exp(yt(t,j)*(sim2(j,t,:)-theta(t,j))-& 121 u(t,j)*(exp(sim2(j,t,:))-tmp(t)))/& 122 exp(-0.5d0/ht(j,j,t)*( (ytilde(t,j)-sim2(j,t,:))**2 - (ytilde(t,j)-theta(t,j))**2)) 123 end if 124 end do 125 case(3) !binomial 126 tmp(1:i) = log(1.0d0+exp(theta(1:i,j))) 127 do t=1,i 128 if(ymiss2(t,j) .EQ. 0) then 129 130 w(i+1,:) = w(i+1,:)*exp( yt(t,j)*(sim2(j,t,:)-theta(t,j))-& 131 u(t,j)*(log(1.0d0+exp(sim2(j,t,:)))-tmp(t)))/& 132 exp(-0.5d0/ht(j,j,t)*( (ytilde(t,j)-sim2(j,t,:))**2 -(ytilde(t,j)-theta(t,j))**2)) 133 134 end if 135 end do 136 case(4) ! gamma 137 tmp(1:i) = exp(-theta(1:i,j)) 138 do t=1,i 139 if(ymiss2(t,j) .EQ. 0) then 140 w(i+1,:) = w(i+1,:)*exp( u(t,j)*(yt(t,j)*(tmp(t)-exp(-sim2(j,t,:)))& 141 +theta(t,j)-sim2(j,t,:)))/& 142 exp(-0.5d0/ht(j,j,t)*( (ytilde(t,j)-sim2(j,t,:))**2 - (ytilde(t,j)-theta(t,j))**2)) 143 end if 144 end do 145 case(5) !negbin 146 tmp(1:i) = exp(theta(1:i,j)) 147 do t=1,i 148 if(ymiss2(t,j) .EQ. 0) then 149 w(i+1,:) = w(i+1,:)*exp(yt(t,j)*(sim2(j,t,:)-theta(t,j)) +& 150 (yt(t,j)+u(t,j))*log((u(t,j)+tmp(t))/(u(t,j)+exp(sim2(j,t,:)))))/& 151 exp(-0.5d0/ht(j,j,t)*( (ytilde(t,j)-sim2(j,t,:))**2 - (ytilde(t,j)-theta(t,j))**2)) 152 end if 153 end do 154 end select 155 end do 156 157 else 158 do j=1,p 159 select case(dist(j)) 160 case(2) !poisson 161 tmp(1:i) = exp(theta(1:i,j)) 162 do t=1,i 163 if(ymiss2(t,j) .EQ. 0) then 164 do k=1,3 * nsim * antithetics + nsim 165 tsim(j,k) = ddot(m,zt(j,:,(t-1)*timevar(1)+1),1,sim2(:,t,k),1) 166 end do 167 w(i+1,:) = w(i+1,:)*exp(yt(t,j)*(tsim(j,:)-theta(t,j))-& 168 u(t,j)*(exp(tsim(j,:))-tmp(t)))/& 169 exp(-0.5d0/ht(j,j,t)*( (ytilde(t,j)-tsim(j,:))**2 - (ytilde(t,j)-theta(t,j))**2)) 170 171 end if 172 end do 173 case(3) !binomial 174 tmp(1:i) = log(1.0d0+exp(theta(1:i,j))) 175 do t=1,i 176 if(ymiss2(t,j) .EQ. 0) then 177 do k=1,3 * nsim * antithetics + nsim 178 tsim(j,k) = ddot(m,zt(j,:,(t-1)*timevar(1)+1),1,sim2(:,t,k),1) 179 end do 180 w(i+1,:) = w(i+1,:)*exp( yt(t,j)*(tsim(j,:)-theta(t,j))-& 181 u(t,j)*(log(1.0d0+exp(tsim(j,:)))-tmp(t)))/& 182 exp(-0.5d0/ht(j,j,t)*( (ytilde(t,j)-tsim(j,:))**2 - (ytilde(t,j)-theta(t,j))**2)) 183 184 end if 185 end do 186 case(4) ! gamma 187 tmp(1:i) = exp(-theta(1:i,j)) 188 do t=1,i 189 if(ymiss2(t,j) .EQ. 0) then 190 do k=1,3 * nsim * antithetics + nsim 191 tsim(j,k) = ddot(m,zt(j,:,(t-1)*timevar(1)+1),1,sim2(:,t,k),1) 192 end do 193 w(i+1,:) = w(i+1,:)*exp( u(t,j)*(yt(t,j)*(tmp(t)-exp(-tsim(j,:)))+theta(t,j)-tsim(j,:)))/& 194 exp(-0.5d0/ht(j,j,t)*( (ytilde(t,j)-tsim(j,:))**2 - (ytilde(t,j)-theta(t,j))**2)) 195 end if 196 end do 197 case(5) !negbin 198 tmp(1:i) = exp(theta(1:i,j)) 199 do t=1,i 200 if(ymiss2(t,j) .EQ. 0) then 201 do k=1,3 * nsim * antithetics + nsim 202 tsim(j,k) = ddot(m,zt(j,:,(t-1)*timevar(1)+1),1,sim2(:,t,k),1) 203 end do 204 w(i+1,:) = w(i+1,:)*exp(yt(t,j)*(tsim(j,:)-theta(t,j)) +& 205 (yt(t,j)+u(t,j))*log((u(t,j)+tmp(t))/(u(t,j)+exp(tsim(j,:)))))/& 206 exp(-0.5d0/ht(j,j,t)*( (ytilde(t,j)-tsim(j,:))**2 - (ytilde(t,j)-theta(t,j))**2)) 207 end if 208 end do 209 end select 210 end do 211 end if 212 sim(:,i+1,:) = sim2(:,i+1,:) 213 end do 214 maxiter=maxitermax 215 216end subroutine isamplefilter 217