1!simulate gaussian state space model
2subroutine simgaussian(ymiss,timevar, yt, zt, ht, tt, rtv, qt, a1, p1, &
3p1inf, nnd,nsim, epsplus, etaplus, aplus1, p, n, m, r, info,rankp,&
4tol,sim,c,simwhat,simdim,antithetics)
5
6    implicit none
7
8    integer, intent(in) :: p, m, r, n, nsim,nnd,simdim,simwhat,antithetics,rankp
9    integer, intent(in), dimension(n,p) :: ymiss
10    integer, intent(in), dimension(5) :: timevar
11    integer, intent(inout) :: info
12    integer ::  t, i, d, j,k,tv,l,rankp2
13    double precision, intent(in) :: tol
14    double precision, intent(in), dimension(n,p) :: yt
15    double precision, intent(in), dimension(p,m,(n-1)*timevar(1)+1) :: zt
16    double precision, intent(in), dimension(p,p,(n-1)*timevar(2)+1) :: ht
17    double precision, intent(in), dimension(m,m,(n-1)*timevar(3)+1) :: tt
18    double precision, intent(in), dimension(m,r,(n-1)*timevar(4)+1) :: rtv
19    double precision, intent(in), dimension(r,r,(n-1)*timevar(5)+1) :: qt
20    double precision, intent(in), dimension(m) :: a1
21    double precision, intent(in), dimension(m,m) ::  p1,p1inf
22    double precision, intent(in),dimension(nsim) :: c
23    double precision, intent(inout), dimension(simdim,n,3 * nsim * antithetics + nsim) :: sim
24    double precision, intent(inout), dimension(r,n,nsim) :: etaplus
25    double precision, intent(inout), dimension(p,n,nsim) :: epsplus
26    double precision, intent(inout), dimension(m,nsim) :: aplus1
27
28    double precision, dimension(n,p) :: yplus
29    double precision, dimension(m,n+1) :: aplus
30    double precision, dimension(m) :: ahat
31    double precision, dimension(m) :: aplushat
32    double precision, dimension(p,n) :: ft,finf
33    double precision, dimension(m,p,n) :: kt,kinf
34    double precision, dimension(r,r,(n-1)*timevar(5)+1) :: cholqt
35    double precision, dimension(m,m) :: cholp1
36    double precision, dimension(r,r) :: rcholtmp
37    double precision, dimension(m) :: rt0,rt1
38    double precision, dimension(m,r) :: mr
39    double precision, dimension(m,m,(n-1)*max(timevar(4),timevar(5))+1) :: rqr
40    double precision, dimension(r,n) :: etahat
41    double precision, dimension(p,n) :: epshat
42    double precision, dimension(r,n) :: etaplushat
43    double precision, dimension(p,n) :: epsplushat
44    double precision, dimension(r,n-1,4) :: etatmp
45    double precision, dimension(m,n,4) :: alphatmp
46    double precision, external :: ddot
47    logical needeps
48
49    external dsymm, dgemm, smoothsim, dsymv, ldl, dtrmv,smoothsimfast, dgemv
50
51
52
53    needeps = simwhat.EQ.1 .or. simwhat.EQ.3
54    if(needeps) then
55        epshat = 0.0d0
56    end if
57
58    !compute rqr
59    tv= max(timevar(4),timevar(5))
60    do t=1, (n-1)*tv+1
61        call dsymm('r','u',m,r,1.0d0,qt(:,:,(t-1)*timevar(5)+1),r,rtv(:,:,(t-1)*timevar(4)+1),m,0.0d0,mr,m)
62        call dgemm('n','t',m,m,r,1.0d0,mr,m,rtv(:,:,(t-1)*timevar(4)+1),m,0.0d0,rqr(:,:,t),m)
63    end do
64
65
66
67    rankp2 = rankp
68    call smoothsim(yt, ymiss, timevar, zt, ht,tt, rtv,qt,rqr, a1, p1, p1inf, &
69    d, j, p, m, n, r,tol,rankp2,ft,finf,kt,kinf,epshat,etahat,rt0,rt1,needeps)
70        !simwhat = 1: epsilon, 2: eta, 3: both, 4: state, 5: signal, 6: observations
71    if(simwhat .GT. 3) then
72        ahat = a1
73        call dsymv('l',m,1.0d0,p1,m,rt0,1,1.0d0,ahat,1)
74        if(d .GT. 0) then
75            call dsymv('l',m,1.0d0,p1inf,m,rt1,1,1.0d0,ahat,1)
76        end if
77    end if
78
79
80
81    do t = 1, (n-1)*timevar(5)+1
82        if(r.EQ.1) then
83            cholqt(1,1,t)=sqrt(qt(1,1,t))
84        else
85            rcholtmp = qt(:,:,t)
86            call ldl(rcholtmp,r,tol,info)
87            if(info .NE. 0) then
88                info = -2
89                return
90            end if
91            do i=1,r
92                cholqt(i,i,t)=sqrt(rcholtmp(i,i))
93            end do
94            do i=1,r-1
95                cholqt((i+1):r,i,t) = rcholtmp((i+1):r,i)*cholqt(i,i,t)
96            end do
97        end if
98    end do
99
100
101    if(nnd.GT.0) then
102        if(m.EQ.1) then
103            cholp1(1,1)=sqrt(p1(1,1))
104        else
105            cholp1 = p1
106            call ldl(cholp1,m,tol,info)
107            if(info .NE. 0) then
108                info=-3
109                return
110            end if
111            do i=1,m
112                cholp1(i,i)=sqrt(cholp1(i,i))
113            end do
114            do i=1,m-1
115                cholp1((i+1):m,i) = cholp1((i+1):m,i)*cholp1(i,i)
116            end do
117        end if
118    end if
119
120    do i = 1, nsim
121        aplus = 0.0d0
122        aplus(:,1) = a1
123        if(nnd.GT.0) then
124            call dtrmv('l','n','n',m,cholp1,m,aplus1(:,i),1)
125            aplus(:,1) = aplus(:,1)+aplus1(:,i)
126        end if
127
128        do t = 1, n
129            do k = 1, p
130                epsplus(k,t,i) = epsplus(k,t,i)*sqrt(ht(k,k,(t-1)*timevar(2)+1))
131                if(ymiss(t,k).EQ.0) then
132                    yplus(t,k) = epsplus(k,t,i) + ddot(m,zt(k,:,(t-1)*timevar(1)+1),1,aplus(:,t),1)
133                end if
134            end do
135            call dtrmv('l','n','n',r,cholqt(:,:,(t-1)*timevar(5)+1),r,etaplus(:,t,i),1)
136            call dgemv('n',m,m,1.0d0,tt(:,:,(t-1)*timevar(3)+1),m,aplus(:,t),1,0.0d0,aplus(:,t+1),1)
137            call dgemv('n',m,r,1.0d0,rtv(:,:,(t-1)*timevar(4)+1),m,etaplus(:,t,i),1,1.0d0,aplus(:,t+1),1)
138        end do
139
140        call smoothsimfast(yplus, ymiss, timevar, zt, ht,tt, rtv,qt,a1, ft,kt,&
141        finf, kinf, d, j, p, m, n,r,epsplushat,etaplushat,rt0,rt1,needeps)
142
143        !simwhat = 1: epsilon, 2: eta, 3: both, 4: state, 5: signal, 6: observations
144        select case(simwhat)
145            case(1)
146                sim(:,:,i) = epshat - epsplushat + epsplus(:,:,i)
147                if(antithetics .EQ. 1) then
148                    sim(:,:,i+nsim) = epshat + epsplushat - epsplus(:,:,i)
149                    sim(:,:,i+2*nsim) = epshat + c(i)*(sim(:,:,i)-epshat)
150                    sim(:,:,i+3*nsim) = epshat + c(i)*(sim(:,:,i+nsim)-epshat)
151                end if
152            case(2)
153                sim(:,:,i) = etahat - etaplushat + etaplus(:,:,i)
154                if(antithetics .EQ. 1) then
155                    sim(:,:,i+nsim) = etahat + etaplushat - etaplus(:,:,i)
156                    sim(:,:,i+2*nsim) = etahat + c(i)*(sim(:,:,i)-etahat)
157                    sim(:,:,i+3*nsim) = etahat + c(i)*(sim(:,:,i+nsim)-etahat)
158                end if
159            case(3)
160                sim(1:p,:,i) = epshat - epsplushat + epsplus(:,:,i)
161                sim((p+1):,:,i) = etahat - etaplushat + etaplus(:,:,i)
162                if(antithetics .EQ. 1) then
163                    sim(1:p,:,i+nsim) = epshat + epsplushat - epsplus(:,:,i)
164                    sim(1:p,:,i+2*nsim) = epshat + c(i)*(sim(1:p,:,i)-epshat)
165                    sim(1:p,:,i+3*nsim) = epshat + c(i)*(sim(1:p,:,i+nsim)-epshat)
166                    sim((p+1):,:,i+nsim) = etahat + etaplushat - etaplus(:,:,i)
167                    sim((p+1):,:,i+2*nsim) = etahat + c(i)*(sim((p+1):,:,i)-etahat)
168                    sim((p+1):,:,i+3*nsim) = etahat + c(i)*(sim((p+1):,:,i+nsim)-etahat)
169                end if
170            case(4)
171                aplushat = a1
172
173                call dsymv('l',m,1.0d0,p1,m,rt0,1,1.0d0,aplushat,1)
174                if(d .GT. 0) then
175                    call dsymv('l',m,1.0d0,p1inf,m,rt1,1,1.0d0,aplushat,1)
176                end if
177
178                sim(:,1,i) = ahat - aplushat + aplus(:,1)
179                etatmp(:,:,1) = etahat(:,1:(n-1)) - etaplushat(:,1:(n-1)) + etaplus(:,1:(n-1),i)
180                if(antithetics .EQ. 1) then
181                    sim(:,1,i+nsim) = ahat + aplushat - aplus(:,1)
182                    sim(:,1,i+2*nsim) = ahat+ c(i)*(sim(:,1,i)-ahat)
183                    sim(:,1,i+3*nsim) = ahat+ c(i)*(sim(:,1,i+nsim)-ahat)
184
185                    etatmp(:,:,2) = etahat(:,1:(n-1)) + etaplushat(:,1:(n-1)) - etaplus(:,1:(n-1),i)
186                    etatmp(:,:,3) = etahat(:,1:(n-1)) + c(i)*(etatmp(:,:,1)-etahat(:,1:(n-1)))
187                    etatmp(:,:,4) = etahat(:,1:(n-1)) + c(i)*(etatmp(:,:,2)-etahat(:,1:(n-1)))
188                end if
189                do k = 1, 3*antithetics+1
190                    do t = 2, n
191                        call dgemv('n',m,m,1.0d0,tt(:,:,(t-2)*timevar(3)+1),m,sim(:,t-1,i+(k-1)*nsim),&
192                        1,0.0d0,sim(:,t,i+(k-1)*nsim),1)
193                        call dgemv('n',m,r,1.0d0,rtv(:,:,(t-2)*timevar(4)+1),m,etatmp(:,t-1,k),1,&
194                        1.0d0,sim(:,t,i+(k-1)*nsim),1)
195                    end do
196                end do
197            case(5)
198                aplushat = a1
199
200                call dsymv('l',m,1.0d0,p1,m,rt0,1,1.0d0,aplushat,1)
201                if(d .GT. 0) then
202                    call dsymv('l',m,1.0d0,p1inf,m,rt1,1,1.0d0,aplushat,1)
203                end if
204
205                alphatmp(:,1,1) = ahat - aplushat + aplus(:,1)
206                etatmp(:,:,1) = etahat(:,1:(n-1)) - etaplushat(:,1:(n-1)) + etaplus(:,1:(n-1),i)
207                call dgemv('n',p,m,1.0d0,zt(:,:,1),p,alphatmp(:,1,1),1,0.0d0,sim(:,1,i),1)
208
209                if(antithetics .EQ. 1) then
210                    alphatmp(:,1,2) = ahat + aplushat - aplus(:,1)
211                    alphatmp(:,1,3) = ahat + c(i)*(alphatmp(:,1,1)-ahat)
212                    alphatmp(:,1,4) = ahat + c(i)*(alphatmp(:,1,2)-ahat)
213                    call dgemv('n',p,m,1.0d0,zt(:,:,1),p,alphatmp(:,1,2),1,0.0d0,sim(:,1,i+nsim),1)
214                    call dgemv('n',p,m,1.0d0,zt(:,:,1),p,alphatmp(:,1,3),1,0.0d0,sim(:,1,i+2*nsim),1)
215                    call dgemv('n',p,m,1.0d0,zt(:,:,1),p,alphatmp(:,1,4),1,0.0d0,sim(:,1,i+3*nsim),1)
216
217                    etatmp(:,:,2) = etahat(:,1:(n-1)) + etaplushat(:,1:(n-1)) - etaplus(:,1:(n-1),i)
218                    etatmp(:,:,3) = etahat(:,1:(n-1)) + c(i)*(etatmp(:,:,1)-etahat(:,1:(n-1)))
219                    etatmp(:,:,4) = etahat(:,1:(n-1)) + c(i)*(etatmp(:,:,2)-etahat(:,1:(n-1)))
220
221                end if
222                do k = 1, 3*antithetics+1
223                    do t = 2, n
224                        call dgemv('n',m,m,1.0d0,tt(:,:,(t-2)*timevar(3)+1),m,alphatmp(:,t-1,k),&
225                        1,0.0d0,alphatmp(:,t,k),1)
226                        call dgemv('n',m,r,1.0d0,rtv(:,:,(t-2)*timevar(4)+1),m,etatmp(:,t-1,k),1,&
227                        1.0d0,alphatmp(:,t,k),1)
228                        call dgemv('n',p,m,1.0d0,zt(:,:,(t-1)*timevar(1)+1),p,&
229                        alphatmp(:,t,k),1,0.0d0,sim(:,t,i+(k-1)*nsim),1)
230                    end do
231                end do
232            case(6)
233                aplushat = a1
234
235                call dsymv('l',m,1.0d0,p1,m,rt0,1,1.0d0,aplushat,1)
236                if(d .GT. 0) then
237                    call dsymv('l',m,1.0d0,p1inf,m,rt1,1,1.0d0,aplushat,1)
238                end if
239
240                alphatmp(:,1,1) = ahat - aplushat + aplus(:,1)
241                etatmp(:,:,1) = etahat(:,1:(n-1)) - etaplushat(:,1:(n-1)) + etaplus(:,1:(n-1),i)
242
243                if(antithetics .EQ. 1) then
244                    alphatmp(:,1,2) = ahat + aplushat - aplus(:,1)
245                    alphatmp(:,1,3) = ahat+ c(i)*(alphatmp(:,1,1)-ahat)
246                    alphatmp(:,1,4) = ahat+ c(i)*(alphatmp(:,1,2)-ahat)
247
248                    etatmp(:,:,2) = etahat(:,1:(n-1)) + etaplushat(:,1:(n-1)) - etaplus(:,1:(n-1),i)
249                    etatmp(:,:,3) = etahat(:,1:(n-1)) + c(i)*(etatmp(:,:,1)-etahat(:,1:(n-1)))
250                    etatmp(:,:,4) = etahat(:,1:(n-1)) + c(i)*(etatmp(:,:,2)-etahat(:,1:(n-1)))
251
252                end if
253                do k = 1, 3*antithetics+1
254                    do l = 1, p
255                        if(ymiss(1,l).GT.0) then
256                            sim(l,1,i+(k-1)*nsim) = ddot(m,zt(l,:,1),1,alphatmp(:,1,k),1)
257                        end if
258                    end do
259                    do t = 2, n
260                        call dgemv('n',m,m,1.0d0,tt(:,:,(t-2)*timevar(3)+1),m,alphatmp(:,t-1,k),&
261                        1,0.0d0,alphatmp(:,t,k),1)
262                        call dgemv('n',m,r,1.0d0,rtv(:,:,(t-2)*timevar(4)+1),m,etatmp(:,t-1,k),1,&
263                        1.0d0,alphatmp(:,t,k),1)
264                        do l = 1, p
265                            if(ymiss(t,l).GT.0) then
266                                sim(l,t,i+(k-1)*nsim) = ddot(m,zt(l,:,(t-1)*timevar(1)+1),1,alphatmp(:,t,k),1)
267                            end if
268                        end do
269                    end do
270                end do
271        end select
272    end do
273
274end subroutine simgaussian
275
276