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