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