1! Filtering of non-gaussian model
2
3subroutine ngfilter(yt, ymiss, timevar, zt, tt, rtv, qt, a1, p1,p1inf, u, theta,&
4dist, p,n, m, r, rankp, nnd,nsim,epsplus,etaplus,aplus1,c,tol,info,maxiter,&
5convtol,alphahat,alphavar,thetahat,thetavar,yhat,yvar,smootha,smooths,smoothy,&
6expected, htol)
7
8    implicit none
9
10    integer, intent(in) ::  p,m, r, n,nnd,nsim,rankp,smootha,smooths,smoothy&
11    , expected
12    integer, intent(in), dimension(p) :: dist
13    integer, intent(in), dimension(n,p) :: ymiss
14    integer, intent(in), dimension(5) :: timevar
15    integer, intent(inout) :: info,maxiter
16    integer ::  t, j
17    double precision, intent(in) :: tol,convtol
18    double precision, intent(inout) :: htol
19    double precision, intent(inout), dimension(n,p) :: theta
20    double precision, intent(in), dimension(n,p) :: u
21    double precision, intent(in), dimension(n,p) :: yt
22    double precision, intent(in), dimension(p,m,(n-1)*timevar(1)+1) :: zt
23    double precision, intent(in), dimension(m,m,(n-1)*timevar(3)+1) :: tt
24    double precision, intent(in), dimension(m,r,(n-1)*timevar(4)+1) :: rtv
25    double precision, intent(in), dimension(r,r,(n-1)*timevar(5)+1) :: qt
26    double precision, intent(in), dimension(m) :: a1
27    double precision, intent(in), dimension(m,m) ::  p1,p1inf
28    double precision, intent(in),dimension(nsim) :: c
29    double precision, intent(inout), dimension(p,n,nsim) :: epsplus
30    double precision, intent(inout), dimension(r,n,nsim) :: etaplus
31    double precision, intent(inout), dimension(m,nsim) :: aplus1
32    double precision, intent(inout), dimension((m-1)*smootha+1,(n-1)*smootha+1) :: alphahat
33    double precision, intent(inout), dimension((m-1)*smootha+1,(m-1)*smootha+1,(n-1)*smootha+1) :: alphavar
34    double precision, intent(inout), dimension((p-1)*smoothy+1,(n-1)*smoothy+1) :: yhat
35    double precision, intent(inout), dimension((p-1)*smoothy+1,(p-1)*smoothy+1,(n-1)*smoothy+1) :: yvar
36    double precision, intent(inout), dimension((p-1)*smooths+1,(n-1)*smooths+1) :: thetahat
37    double precision, intent(inout), dimension((p-1)*smooths+1,(p-1)*smooths+1,(n-1)*smooths+1) :: thetavar
38    double precision, dimension(smootha*m+(1-smootha)*p,n,4*nsim) :: sim
39    double precision, dimension(smootha*p+(1-smootha),n,4*nsim*smootha+(1-smootha)) :: osim
40    double precision, dimension(n,4*nsim) :: w
41    double precision, dimension(p,m) :: pm
42    double precision, external :: ddot
43
44    external isamplefilter, covmeanwprotect, dgemv, dsymm, dgemm, covmeanw
45
46    if(smootha.EQ.1) then
47
48        call isamplefilter(yt, ymiss, timevar, zt, tt, rtv, qt, a1, p1,p1inf, u, dist, &
49        p, n, m, r, theta, maxiter,rankp,convtol, nnd,nsim,epsplus,etaplus,&
50        aplus1,c,tol,info,1,w,sim,4,m, expected, htol)
51
52        if(info .ne. 0 .and. info .ne. 3) then
53            return
54        end if
55
56        do t=1,n
57            w(t,:) = w(t,:)/sum(w(t,:))
58            call covmeanwprotect(sim(:,t,:),w(t,:),m,1,4*nsim,alphahat(:,t),alphavar(:,:,t))
59        end do
60
61        if(smooths.EQ.1) then
62            do t = 1, n
63                call dgemv('n',p,m,1.0d0,zt(:,:,(t-1)*timevar(1)+1),p,alphahat(:,t),1,0.0d0,thetahat(:,t),1)
64                call dsymm('r','u',p,m,1.0d0,alphavar(:,:,t),m,zt(:,:,(t-1)*timevar(1)+1),p,0.0d0,pm,p)
65                call dgemm('n','t',p,p,m,1.0d0,pm,p,zt(:,:,(t-1)*timevar(1)+1),p,0.0d0,thetavar(:,:,t),p)
66            end do
67        end if
68
69        if(smoothy.EQ.1) then
70            do t = 1, n
71                do j = 1,p
72                    call dgemv('t',m,4*nsim,1.0d0,sim(:,t,:),m,zt(j,:,(t-1)*timevar(1)+1),1,0.0d0,osim(j,t,:),1)
73                end do
74            end do
75
76            do j= 1,p
77                select case(dist(j))
78                    case(1)
79
80                    case(2)
81                        do t=1, n
82                            osim(j,t,:) = exp(osim(j,t,:))*u(t,j)
83                        end do
84                    case(3)
85                        osim(j,:,:) = exp(osim(j,:,:))/(1.0d0+exp(osim(j,:,:)))
86                    case default
87                        osim(j,:,:) = exp(osim(j,:,:))
88                end select
89            end do
90            do t=1,n
91                call covmeanw(osim(:,t,:),w(t,:),p,1,4*nsim,yhat(:,t),yvar(:,:,t))
92            end do
93        end if
94    else
95        call isamplefilter(yt, ymiss, timevar, zt, tt, rtv, qt, a1, p1,p1inf, u, dist, &
96        p, n, m, r, theta, maxiter,rankp,convtol, nnd,nsim,epsplus,etaplus,&
97        aplus1,c,tol,info,1,w,sim,5,p, expected, htol)
98
99        if(info .ne. 0 .and. info .ne. 3) then
100            return
101        end if
102
103        do t=1,n
104            w(t,:) = w(t,:)/sum(w(t,:))
105        end do
106
107        if(smooths.EQ.1) then
108            do t=1,n
109                call covmeanwprotect(sim(:,t,:),w(t,:),p,1,4*nsim,thetahat(:,t),thetavar(:,:,t))
110            end do
111        end if
112
113        if(smoothy.EQ.1) then
114            do j= 1,p
115                select case(dist(j))
116                    case(1)
117
118                    case(2)
119                        do t=1, n
120                            sim(j,t,:) = exp(sim(j,t,:))*u(t,j)
121                        end do
122                    case(3)
123                        sim(j,:,:) = exp(sim(j,:,:))/(1.0d0+exp(sim(j,:,:)))
124                    case default
125                        sim(j,:,:) = exp(sim(j,:,:))
126                end select
127            end do
128            do t=1,n
129                call covmeanw(sim(:,t,:),w(t,:),p,1,4*nsim,yhat(:,t),yvar(:,:,t))
130            end do
131        end if
132
133    end if
134end subroutine ngfilter
135