1C --------------------------------------------------------------------
2C KS IMPLEMENTS THE KALMAN SMMOOTHER RECURSIONS in
3C Koopman (1997), JASA, 92, 440, 1630-38
4C Developed by A.Rossi, C.Planas and G.Fiorentini
5C
6C XS = E[x(t)|y(1),...,y(nobs)]
7C PS = V[x(t)|y(1),...,y(nobs)], t = 1,2,...,nobs
8C
9C State-space format:   y(t) = c(t)z(t) + H(t)x(t)   + G(t)u(t)
10C                       x(t) = a(t)     + F(t)x(t-1) + R(t)u(t)
11C
12C y(t) (ny x 1)          ny  = # of endogenous series
13C z(t) (nz x 1)          nz  = # of exogenous series
14C x(t) (nx x 1)          nx  = # of continous states
15C u(t) (nu x 1)          nu  = # of shocks
16C c(t) (ny x nz x ns1)   ns1 = # of states for c(t)
17C H(t) (ny x nx x ns2)   ns2 = # of states for S2(t)
18C G(t) (ny x nu x ns3)   ns3 = # of states for S3(t)
19C a(t) (nx x ns4)        ns4 = # of states for S4(t)
20C F(t) (nx x nx x ns5)   ns5 = # of states for S5(t)
21C R(t) (nx x nu x ns6)   ns6 = # of states for S6(t)
22C
23C d(1): order of integration of the system
24C d(2): number of non-stationary elements
25C
26C Copyright (C) 2010-2014 European Commission
27C
28C This file is part of Program DMM
29C
30C DMM is free software developed at the Joint Research Centre of the
31C European Commission: you can redistribute it and/or modify it under
32C the terms of the GNU General Public License as published by
33C the Free Software Foundation, either version 3 of the License, or
34C (at your option) any later version.
35C
36C DMM is distributed in the hope that it will be useful,
37C but WITHOUT ANY WARRANTY; without even the implied warranty of
38C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
39C GNU General Public License for more details.
40C
41C You should have received a copy of the GNU General Public License
42C along with DMM.  If not, see <http://www.gnu.org/licenses/>.
43C --------------------------------------------------------------------
44      SUBROUTINE KS(nobs,d,ny,nz,nx,nu,ns,S,yk,IYK,c,H,G,a,F,R,XS,PS)
45C INPUT
46	INTEGER nobs,d(2),ny,nz,nx,nu,ns(6),S(nobs,6),IYK(nobs,ny+1)
47	DOUBLE PRECISION yk(nobs,ny+nz),c(ny,max(nz,1),ns(1)),
48	1 H(ny,nx,ns(2)),G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)),
49     2 R(nx,nu,ns(6))
50
51C OUTPUT
52	DOUBLE PRECISION XS(nobs,nx),PS(nobs,nx,nx)
53
54C LOCALS
55	INTEGER imain,iny,I,J,IFAIL,FiRANK,ITIME
56	INTEGER IPIV(nx)
57	DOUBLE PRECISION TOL,SUMW1
58	DOUBLE PRECISION,ALLOCATABLE:: Pi(:,:,:),HPs(:,:),HPi(:,:),
59	1 Fi(:,:),Fs(:,:),Fim(:,:),Fsm(:,:),PHFs(:,:),PHFi(:,:),FFF(:,:),
60     2 Mi(:,:),Ms(:,:),Ci(:,:),Kst(:,:,:),Kit(:,:,:),W1(:),WORK(:),
61     3 WORK1(:),PFP(:,:),APPO(:,:),APPO1(:,:),COM(:,:),RG(:,:),
62	4 FP(:,:),HP1(:,:),V(:,:),CC(:,:),HPV(:,:),X0(:),P0(:,:),
63     5 RECR(:),RECRI(:),RECN(:,:),XT(:,:),PT(:,:,:),INN(:,:),
64	1 Vinv(:,:,:),Vis(:,:,:),Vii(:,:,:)
65
66C EXTERNAL SUBROUTINES
67      EXTERNAL INVF,INVFBIS,LYAP,DSYEV,DPOTRF,DPOTRI,DGETRF,DGETRI
68
69      ALLOCATE(Pi(d(1),nx,nx),HPs(ny,nx),HPi(ny,nx),
70	1 Fi(ny,ny),Fs(ny,ny),Fim(ny,ny),Fsm(ny,ny),
71     2 PHFs(nx,ny),PHFi(nx,ny),FFF(ny,ny),Mi(nx,ny),Ms(nx,ny),Ci(nx,nx),
72	3 Kst(d(1),nx,ny),Kit(d(1),nx,ny))
73
74	ALLOCATE(W1(ny),WORK(64*nx),WORK1(64*ny),
75	1 PFP(nx,nx),APPO(nx,nx),APPO1(nx,ny),COM(ny+1,ny),RG(nx,ny))
76
77	ALLOCATE(FP(nx,nx),HP1(ny,nx),V(ny,ny),CC(nx,nx),
78	1 HPV(nx,ny),X0(nx),P0(nx,nx),RECR(nx),RECRI(nx),RECN(nx,nx))
79
80	ALLOCATE(XT(nobs,nx),PT(nobs,nx,nx),INN(nobs,ny),
81	1 Vinv(nobs,ny,ny),Vis(d(1),ny,ny),Vii(d(1),ny,ny))
82
83	TOL = 1.D-3
84C Unconditional mean and variance
85      IF (d(1).EQ.0) THEN ! stationary models
86	 IF(SUM(ABS(a(:,S(1,4)))).EQ.0.D0) THEN
87	  XT(1,:) = 0.D0  ! X(1|0)
88	 ELSE
89	  APPO = -F(:,:,S(1,5))
90	  DO 1 I = 1,nx
911	  APPO(I,I) = 1.D0+APPO(I,I)
92C	  CALL F07ADF(nx,nx,APPO,nx,IPIV,IFAIL)
93C	  CALL F07AJF(nx,APPO,nx,IPIV,WORK,64*nx,IFAIL)
94	  CALL DGETRF(nx,nx,APPO,nx,IPIV,IFAIL)
95	  CALL DGETRI(nx,APPO,nx,IPIV,WORK,64*nx,IFAIL)
96
97        DO 2 I =1,nx
982	  XT(1,I) = SUM(APPO(I,:)*a(:,S(1,4))) ! inv(I-F)*a
99       ENDIF
100
101C P(1|0) - F*P(1|0)*F' = R*R'
102	 CALL LYAP(nx,nu,TOL,F(:,:,S(1,5)),R(:,:,S(1,6)),PT(1,:,:))
103	ELSE
104C -----------------------------------------------------------
105C Non-stationary models
106C Define X(1) = aa + A*eta + B*delta (A*B' = 0)
107C               eta~N(0,I), delta~N(0,k*I) k -> +inf
108C               X(1)~N(aa,P), P=Ps+k*Pi, Ps=AA', Pi=BB'.
109C CARE!! aa (uncond. mean),Ps, and Pi to be filled by users
110C -----------------------------------------------------------
111	 XT(1,1:nx)      = 0.D0  ! X(1|0)
112	 PT(1,1:nx,1:nx) = 0.D0  ! P(1|0)
113	 IF (d(2).LT.nx) THEN
114	  IF(SUM(ABS(a(d(2)+1:nx,S(1,4)))).NE.0.D0) THEN
115	   APPO(d(2)+1:nx,d(2)+1:nx) = -F(d(2)+1:nx,d(2)+1:nx,S(1,5))
116	   DO 3 I = d(2)+1,nx
1173	   APPO(I,I) = 1.D0+APPO(I,I)
118C	   CALL F07ADF(nx-d(2),nx-d(2),APPO(d(2)+1:nx,d(2)+1:nx),nx-d(2),
119C	1               IPIV(d(2)+1:nx),IFAIL)
120C	   CALL F07AJF(nx-d(2),APPO(d(2)+1:nx,d(2)+1:nx),nx-d(2),
121C	1               IPIV(d(2)+1:nx),WORK,64*nx,IFAIL)
122	   CALL DGETRF(nx-d(2),nx-d(2),APPO(d(2)+1:nx,d(2)+1:nx),nx-d(2),
123	1               IPIV(d(2)+1:nx),IFAIL)
124	   CALL DGETRI(nx-d(2),APPO(d(2)+1:nx,d(2)+1:nx),nx-d(2),
125	1               IPIV(d(2)+1:nx),WORK,64*nx,IFAIL)
126
127	   DO 4 I = d(2)+1,nx
1284	   XT(1,I) = SUM(APPO(I,d(2)+1:nx)*a(d(2)+1:nx,S(1,4))) ! inv(I-F)*a
129        ENDIF
130
131C Lyapunov eqn
132	  CALL LYAP(nx-d(2),nu,TOL,F(d(2)+1:nx,d(2)+1:nx,S(1,5)),
133	1            R(d(2)+1:nx,1:nu,S(1,6)),PT(1,d(2)+1:nx,d(2)+1:nx))
134       ENDIF
135
136	 Pi(:,:,:) = 0.D0
137	 DO 5 I = 1,d(2)
1385	 Pi(1,I,I) = 1.D0
139
140	 DO 200 imain = 1,d(1)
141	  iny = IYK(imain,ny+1)
142	  DO 30 I=1,iny
143        DO 30 J=1,nx
14430      HPs(I,J) = SUM(H(IYK(imain,I),:,S(imain,2))*PT(imain,:,J))
145
146	  DO 40 I=1,iny
147	   Fs(I,I) = SUM(HPs(I,:)*H(IYK(imain,I),:,S(imain,2)))
148     +   +SUM(G(IYK(imain,I),:,S(imain,3))*G(IYK(imain,I),:,S(imain,3)))
149         DO 40 J=1,I-1
150         Fs(I,J) = SUM(HPs(I,:)*H(IYK(imain,J),:,S(imain,2)))
151     +   +SUM(G(IYK(imain,I),:,S(imain,3))*G(IYK(imain,J),:,S(imain,3)))
15240       Fs(J,I) = Fs(I,J)
153
154  	  DO 50 I=1,iny
155        DO 50 J=1,nx
15650      HPi(I,J) = SUM(H(IYK(imain,I),:,S(imain,2))*Pi(imain,:,J))
157
158	  DO 60 I=1,iny
159         Fi(I,I) = SUM(HPi(I,:)*H(IYK(imain,I),:,S(imain,2)))
160	   DO 60 J=1,I-1
161         Fi(I,J) = SUM(HPi(I,:)*H(IYK(imain,J),:,S(imain,2)))
16260       Fi(J,I) = Fi(I,J)
163
164C --------------------------------------------------------------------------
165C Computes inverse of the innovation variance matrix
166C Cases: ny = 1, Fi is scalar >0 (or 0 not considered)
167C        ny > 1, Fi is full rank or singular (or 0 matrix not considered)
168C --------------------------------------------------------------------------
169        IF (iny.EQ.1) THEN
170	   Fsm = 0.D0
171	   Fim = 1.D0/Fi
172	   FFF = Fim*Fs*Fim
173	  ELSE
174
175	  IFAIL = -1
176	  COM(1:iny,1:iny) = Fi(1:iny,1:iny)
177C	  CALL F02FAF('N','U',iny,COM(1:iny,1:iny),iny,W1(1:iny),
178C	1              WORK1,64*iny,IFAIL)
179    	  CALL DSYEV('N','U',iny,COM(1:iny,1:iny),iny,W1(1:iny),
180	1              WORK1,64*iny,IFAIL)
181
182	  FiRANK = 0
183	  SUMW1  = SUM(ABS(W1(1:iny)))
184	  DO 70 I=1,iny
185	  W1(I) = W1(I)/SUMW1
18670      IF (W1(I).GT.1.D-10) FiRANK=FiRANK+1
187	  FiRANK = min(FiRANK,d(2))
188
189        IF(FiRANK.EQ.iny) THEN
190	    Fsm = 0.D0
191	    COM(1:iny,1:iny) = Fi(1:iny,1:iny)
192		IFAIL = -1
193C		CALL F01ADF(iny,COM(1:iny+1,1:iny),iny+1,IFAIL)
194          CALL DPOTRF('L',iny,COM(1:iny,1:iny),iny,IFAIL) ! COM = L*L'
195          CALL DPOTRI('L',iny,COM(1:iny,1:iny),iny,IFAIL) ! COM = VV^-1
196
197		DO 80 I=1,iny
198		Fim(I,I) = COM(I,I)
199		DO 80 J=1,I-1
200		Fim(I,J) = COM(I,J)
20180		Fim(J,I) = Fim(I,J)
202
203		DO 81 I=1,iny
204		DO 81 J=1,iny
20581		COM(I,J) = SUM(Fim(I,1:iny)*Fs(1:iny,J))  ! Fim x Fs
206
207		DO 82 I=1,iny
208		 FFF(I,I) = SUM(COM(I,1:iny)*Fim(1:iny,I))
209		 DO 82 J=1,I-1
210	     FFF(I,J) = SUM(COM(I,1:iny)*Fim(1:iny,J))  ! Fim x Fs x Fim
21182		 FFF(J,I) = FFF(I,J)
212	   ELSE
213	    SUMW1=0.D0
214		DO I=Firank+1,iny
215		 SUMW1 = SUMW1 + Fi(I,I)
216          ENDDO
217          IF (SUMW1.GT.0.D0) THEN
218		 CALL INVFBIS(Fs(1:iny,1:iny),Fi(1:iny,1:iny),iny,FiRANK,
219	1          Fsm(1:iny,1:iny),Fim(1:iny,1:iny),FFF(1:iny,1:iny))
220		ELSE
221		 CALL INVF(Fs(1:iny,1:iny),Fi(1:iny,1:iny),iny,FiRANK,
222	1          Fsm(1:iny,1:iny),Fim(1:iny,1:iny),FFF(1:iny,1:iny))
223          ENDIF
224	   ENDIF
225	  ENDIF
226	  Vis(imain,1:iny,1:iny) = Fsm(1:iny,1:iny)
227	  Vii(imain,1:iny,1:iny) = Fim(1:iny,1:iny)
228
229C ------------------------------------------------------------------
230C X(d|d) = X(d|d-1)+((Ps*H'+R*G')*Fsm+Pi*H'*Fim)*(Y(d)-H*X(d|d-1)-c)
231C ------------------------------------------------------------------
232	  DO 85 I = 1,nx
233	  DO 85 J = 1,iny
234 	  RG(I,J) =
235     #	 SUM(R(I,1:nu,S(imain,6))*G(IYK(imain,J),1:nu,S(imain,3)))
23685      HPs(J,I) = HPs(J,I) + RG(I,J) ! HPs = (Ps*H'+R*G')'
237
238	  DO 90 I = 1,nx
239	  DO 90 J = 1,iny
240 	  PHFs(I,J) = SUM(HPs(1:iny,I)*Fsm(1:iny,J))
24190	  PHFi(I,J) = SUM(HPi(1:iny,I)*Fim(1:iny,J))
242
243C Innovations
244	  DO 100 I=1,iny
245100	  INN(imain,I) = yk(imain,IYK(imain,I))
246     +   - SUM(H(IYK(imain,I),1:nx,S(imain,2))*XT(imain,1:nx))
247     +   - SUM(c(IYK(imain,I),1:nz,S(imain,1))*yk(imain,ny+1:ny+nz))
248
249	  DO 110 I=1,nx
250110     X0(I) = XT(imain,I)
251     +        + SUM((PHFs(I,1:iny)+PHFi(I,1:iny))*INN(imain,1:iny))
252
253C P(d|d) = P(d|d-1) + Pi*H'*Fim*Fs*Fim*H*Pi - Ps*H'*Fsm*H*Ps - Ps*H'*Fim*H*Pi - (Ps*H'*Fim*H*Pi)'
254C - Ps*H'*Fsm*H*Ps
255	  DO 120 I = 1,nx
256	   APPO(I,I) = -SUM(PHFs(I,1:iny)*HPs(1:iny,I))
257	   DO 120 J = 1,I-1
258	   APPO(I,J) = -SUM(PHFs(I,1:iny)*HPs(1:iny,J))
259120	   APPO(J,I) = APPO(I,J)
260
261C - Ps*H'*Fim*H*Pi - (Ps*H'*Fim*H*Pi)'
262	  DO 130 I = 1,nx
263	   APPO(I,I) = APPO(I,I) - SUM(HPs(1:iny,I)*PHFi(I,1:iny))
264     +             - SUM(PHFi(I,1:iny)*HPs(1:iny,I))
265	   DO 130 J = 1,I-1
266	   APPO(I,J) = APPO(I,J) - SUM(HPs(1:iny,I)*PHFi(J,1:iny))
267     +             - SUM(PHFi(I,1:iny)*HPs(1:iny,J))
268130	   APPO(J,I) = APPO(I,J)
269
270C Pi*H'*Fim*Fs*Fim*H*Pi
271	  DO 140 I = 1,nx
272	  DO 140 J = 1,iny
273140	  APPO1(I,J) = SUM(HPi(1:iny,I)*FFF(1:iny,J))
274
275	  DO 150 I = 1,nx
276	   PFP(I,I) = SUM(APPO1(I,1:iny)*HPi(1:iny,I))
277	   DO 150 J = 1,I-1
278	   PFP(I,J) = SUM(APPO1(I,1:iny)*HPi(1:iny,J))
279150      PFP(J,I) = PFP(I,J)
280
281	  P0(:,:) = PT(imain,:,:) + PFP(:,:) + APPO(:,:)
282
283C Mi = F*Pi*H'
284	  DO 151 I = 1,nx
285	  DO 151 J = 1,nx
286151	  PFP(I,J) = SUM(F(I,1:nx,S(imain,5))*Pi(imain,1:nx,J)) ! F*Pi
287
288  	  DO 152 I = 1,nx
289	  DO 152 J = 1,iny
290152	  Mi(I,J) = SUM(PFP(I,1:nx)*H(IYK(imain,J),1:nx,S(imain,2)))
291
292C Ms = F*Ps*H' + R*G'
293	  DO 153 I = 1,nx
294	  DO 153 J = 1,nx
295153	  FP(I,J) = SUM(F(I,:,S(imain,5))*PT(imain,:,J)) ! F*Ps
296
297	  DO 154 I = 1,nx
298	  DO 154 J = 1,iny
299154	  Ms(I,J)=RG(I,J)+SUM(FP(I,1:nx)*H(IYK(imain,J),1:nx,S(imain,2)))
300
301C Kit = Ms*Fim - Mi*FFF
302C Kst = Ms*Fsm + Mi*Fim
303  	  DO 155 I = 1,nx
304	  DO 155 J = 1,iny
305	  Kit(imain,I,J) = SUM(Ms(I,1:iny)*Fim(1:iny,J))
306     +                 - SUM(Mi(I,1:iny)*FFF(1:iny,J))
307155	  Kst(imain,I,J) = SUM(Ms(I,1:iny)*Fsm(1:iny,J))
308     +                 + SUM(Mi(I,1:iny)*Fim(1:iny,J))
309
310C ----------------------------------
311C Predictions X(t+1|t) and P(t+1|t)
312C ----------------------------------
313	  IF (imain.LT.d(1)) THEN
314C Pi+1 = F*Pi*F'-Ci
315
316C Ci = Mi*Fim*Mi'
317	   DO 164 I = 1,nx
318	   DO 164 J = 1,iny
319164	   RG(I,J) = SUM(Mi(I,1:iny)*Fim(1:iny,J))  ! Mi*Fim
320
321	   DO 166 I = 1,nx
322	   Ci(I,I) = SUM(RG(I,1:iny)*Mi(I,1:iny))
323	   DO 166 J = 1,I-1
324166	   Ci(I,J) = SUM(RG(I,1:iny)*Mi(J,1:iny))
325
326	   DO 168 I = 1,nx
327	   Pi(imain+1,I,I)=SUM(PFP(I,1:nx)*F(I,1:nx,S(imain+1,5)))-Ci(I,I)
328	   DO 168 J = 1,I-1
329	   Pi(imain+1,I,J)=SUM(PFP(I,1:nx)*F(J,1:nx,S(imain+1,5)))-Ci(I,J)
330168      Pi(imain+1,J,I) = Pi(imain+1,I,J)
331
332        ENDIF
333
334C X(t+1|t) = a + F*X(t|t)
335	  DO 170 I=1,nx
336170	  XT(imain+1,I) = a(I,S(imain+1,4))
337     +                + SUM(F(I,1:nx,S(imain+1,5))*X0(1:nx))
338
339C P(t+1|t) = F*PddF' + R*R'
340	  DO 172 I = 1,nx
341	  DO 172 J = 1,nx
342172	  APPO(I,J) = SUM(F(I,:,S(imain+1,5))*P0(:,J)) ! F*Pdd
343
344	  DO 180 I = 1,nx
345	   PT(imain+1,I,I) = SUM(APPO(I,1:nx)*F(I,1:nx,S(imain+1,5)))
346     +           + SUM(R(I,1:nu,S(imain+1,6))*R(I,1:nu,S(imain+1,6)))
347	   DO 180 J = 1,I-1
348	   PT(imain+1,I,J) = SUM(APPO(I,1:nx)*F(J,1:nx,S(imain+1,5)))
349     +           + SUM(R(I,1:nu,S(imain+1,6))*R(J,1:nu,S(imain+1,6)))
350180	   PT(imain+1,J,I) = PT(imain+1,I,J)
351
352200    CONTINUE
353	ENDIF
354
355	DO 400 imain = d(1)+1,nobs
356	 iny = IYK(imain,ny+1)
357C -------------------------------
358C Innovations: INN = yk-H*X1-c*z
359C -------------------------------
360       DO 210 I=1,iny
361210    INN(imain,I) = yk(imain,IYK(imain,I))
362     +       - SUM(H(IYK(imain,I),1:nx,S(imain,2))*XT(imain,:))
363     +       - SUM(c(IYK(imain,I),1:nz,S(imain,1))*yk(imain,ny+1:ny+nz))
364
365C ----------------------------------------------------------
366C Innovation variance V = H*P1*H' + G*G' + H*R*G' + G*R'*H'
367C ----------------------------------------------------------
368       DO 220 I=1,iny
369       DO 220 J=1,nx
370220    HP1(I,J) = SUM(H(IYK(imain,I),1:nx,S(imain,2))*PT(imain,1:nx,J))
371
372	 DO 221 I=1,nx
373	 DO 221 J=1,iny
374221	 RG(I,J)=SUM(R(I,1:nu,S(imain,6))*G(IYK(imain,J),1:nu,S(imain,3)))
375
376	 DO 222 I=1,iny
377	 DO 222 J=1,iny
378222	 COM(I,J)=SUM(H(IYK(imain,I),1:nx,S(imain,2))*RG(1:nx,J)) ! H*R*G'
379
380	 DO 230 I=1,iny
381        V(I,I) = SUM(HP1(I,1:nx)*H(IYK(imain,I),1:nx,S(imain,2)))
382     #         + SUM(G(IYK(imain,I),1:nu,S(imain,3))*
383     #           G(IYK(imain,I),1:nu,S(imain,3))) + 2.*COM(I,I)
384        DO 230 J=1,I-1
385        V(I,J) = SUM(HP1(I,1:nx)*H(IYK(imain,J),1:nx,S(imain,2)))+
386     #           SUM(G(IYK(imain,I),1:nu,S(imain,3))*
387     #           G(IYK(imain,J),1:nu,S(imain,3)))+COM(I,J)+COM(J,I)
388230     V(J,I) = V(I,J)
389
390C -------------------------------------------------------------------
391C  Updating equations:
392C  x0 = x1 + (P1*H'+R*G')*Vinv*INN
393C  p0 = p1 - (P1*H'+R*G')*Vinv*(P1*H'+R*G')'
394C -------------------------------------------------------------------
395       IF (iny.GT.0) THEN
396	  COM(1:iny,1:iny) = V(1:iny,1:iny)
397	  IFAIL = -1
398C	  CALL F01ADF(iny,COM(1:iny+1,1:iny),iny+1,IFAIL)
399        CALL DPOTRF('L',iny,COM(1:iny,1:iny),iny,IFAIL) ! COM = L*L'
400        CALL DPOTRI('L',iny,COM(1:iny,1:iny),iny,IFAIL) ! COM = VV^-1
401
402	  DO 240 I=1,iny
403	  Vinv(imain,I,I) = COM(I,I)
404	  DO 240 J=1,I-1
405	  Vinv(imain,I,J) = COM(I,J)
406240	  Vinv(imain,J,I) = Vinv(imain,I,J)
407
408	  DO 260 I=1,nx
409	  DO 260 J=1,iny
410260	  HPV(I,J) = SUM((HP1(1:iny,I)+RG(I,1:iny))*Vinv(imain,1:iny,J))
411
412     	  DO 270 I=1,nx
413270	  X0(I) = XT(imain,I)+SUM(HPV(I,1:iny)*INN(imain,1:iny))
414
415	  DO 280 I=1,nx
416	  P0(I,I) = PT(imain,I,I)
417     +          - SUM(HPV(I,1:iny)*(HP1(1:iny,I)+RG(I,1:iny)))
418	  DO 280 J=1,I-1
419	  P0(I,J) = PT(imain,I,J)
420     +	     - SUM(HPV(I,1:iny)*(HP1(1:iny,J)+RG(J,1:iny)))
421280	  P0(J,I) = P0(I,J)
422
423	 ELSE
424
425	  X0(1:nx)      = XT(imain,1:nx)
426	  P0(1:nx,1:nx) = PT(imain,1:nx,1:nx)
427
428	 ENDIF
429
430C ------------------------------------
431C Prediction       x1 = c+F*x0
432C Prediction var.  P1 = F*p0*F'+ R*R'
433C ------------------------------------
434	 IF (imain.LT.nobs) THEN
435	  DO 290 I=1,nx
436290     XT(imain+1,I) = a(I,S(imain+1,4))
437     +                + SUM(F(I,1:nx,S(imain+1,5))*X0(1:nx))
438
439	  DO 300 I=1,nx
440        DO 300 J=1,nx
441300     FP(I,J) = SUM(F(I,1:nx,S(imain+1,5))*P0(1:nx,J))
442
443        DO 310 I=1,nx
444        PT(imain+1,I,I) = SUM(FP(I,:)*F(I,:,S(imain+1,5)))
445     +           + SUM(R(I,1:nu,S(imain+1,6))*R(I,1:nu,S(imain+1,6)))
446	  DO 310 J=1,I-1
447        PT(imain+1,I,J) = SUM(FP(I,:)*F(J,:,S(imain+1,5)))
448     +          + SUM(R(I,1:nu,S(imain+1,6))*R(J,1:nu,S(imain+1,6)))
449310     PT(imain+1,J,I) = PT(imain+1,I,J)
450	 ENDIF
451400   CONTINUE
452
453C **** SMOOTHING BAKWARD RECURSIONS ****
454	RECR(1:nx)      = 0.D0
455	RECN(1:nx,1:nx) = 0.D0
456	DO 600 ITIME = nobs,d(1)+1,-1
457	 iny = IYK(ITIME,ny+1)
458
459C R*G' and H'*Vinv
460	 DO 420 J=1,iny
461	 DO 420 I=1,nx
462	 RG(I,J)=SUM(R(I,1:nu,S(ITIME,6))*G(IYK(ITIME,J),1:nu,S(ITIME,3)))
463420    PHFs(I,J) =
464     #  SUM(H(IYK(ITIME,1:iny),I,S(ITIME,2))*Vinv(ITIME,1:iny,J))
465
466C F(t+1)*P(t|t-1)
467	 DO 430 I=1,nx
468       DO 430 J=1,nx
469430    FP(I,J) = SUM(F(I,1:nx,S(min(nobs,ITIME+1),5))*PT(ITIME,1:nx,J))
470
471C H'*Vinv*H
472     	 DO 440 I=1,nx
473       APPO(I,I) = SUM(PHFs(I,1:iny)*H(IYK(ITIME,1:iny),I,S(ITIME,2)))
474	 DO 440 J=1,I-1
475       APPO(I,J) = SUM(PHFs(I,1:iny)*H(IYK(ITIME,1:iny),J,S(ITIME,2)))
476440    APPO(J,I) = APPO(I,J)
477
478C L(t) =  F(t+1)-(F(t+1)*P(t|t-1)*H(t)'+R(t)*G(t)')*Vinv(t)*H(t)
479    	 DO 450 I=1,nx
480       DO 450 J=1,nx
481450    PFP(I,J) = F(I,J,S(min(nobs,ITIME+1),5))
482     +          - SUM(FP(I,1:nx)*APPO(1:nx,J))
483     +          - SUM(RG(I,1:iny)*PHFs(J,1:iny))
484
485C r(t-1) = H(t)'*Vinv(t)*INN(t)+L(t)'*r(t)
486	 DO 460 I=1,nx
487460    WORK(I) = SUM(PFP(1:nx,I)*RECR(1:nx))
488
489	 DO 470 I=1,nx
490470	 RECR(I) = WORK(I) + SUM(PHFs(I,1:iny)*INN(ITIME,1:iny))
491
492C N(t-1) = H(t)'*Vinv(t)*H(t)+L(t)'*N(t)*L(t)
493     	 DO 480 I=1,nx
494       DO 480 J=1,nx
495480    CC(I,J) = SUM(PFP(1:nx,I)*RECN(1:nx,J)) ! L(t)'*N(t)
496
497	 DO 490 I=1,nx
498	  RECN(I,I) = APPO(I,I) + SUM(CC(I,1:nx)*PFP(1:nx,I))
499        DO 490 J=1,I-1
500	  RECN(I,J) = APPO(I,J) + SUM(CC(I,1:nx)*PFP(1:nx,J))
501490     RECN(J,I) = RECN(I,J)
502
503C X(t|T) = X(t|t-1) + P(t|t-1)*r(t-1)
504	 DO 500 I = 1,nx
505500	 XS(ITIME,I) = XT(ITIME,I) + SUM(PT(ITIME,I,1:nx)*RECR(1:nx))
506
507C P(t|T) = P(t|t-1) - P(t|t-1)*N(t-1)*P(t|t-1)
508    	 DO 510 I=1,nx
509       DO 510 J=1,nx
510510    CC(I,J) = SUM(PT(ITIME,I,1:nx)*RECN(1:nx,J)) ! P(t|t-1)*N(t-1)
511
512	 DO 520 I=1,nx
513	 PS(ITIME,I,I) = PT(ITIME,I,I) - SUM(CC(I,1:nx)*PT(ITIME,1:nx,I))
514       DO 520 J=1,I-1
515	 PS(ITIME,I,J) = PT(ITIME,I,J) - SUM(CC(I,1:nx)*PT(ITIME,1:nx,J))
516520    PS(ITIME,J,I) = PS(ITIME,I,J)
517
518600   CONTINUE
519
520C INITIAL KALMAN SAMOOTING
521	RECRI(1:nx) = 0.D0
522	DO 800 ITIME = d(1),1,-1
523	 iny = IYK(ITIME,ny+1)
524C L(t) =  F(t)-Kst(t)*H(t)
525    	 DO 610 I=1,nx
526       DO 610 J=1,nx
527610    PFP(I,J) = F(I,J,S(ITIME,5))
528     +   - SUM(Kst(ITIME,I,1:iny)*H(IYK(ITIME,1:iny),J,S(ITIME,2)))
529
530C r(t-1)  = H(t)'*Fsm*INN(t) + L(t)'*r(t)
531C ri(t-1) = H(t)'*Fim*INN(t) + L(t)'*ri(t) + Li(t)'*r(t)
532	 DO 620 I=1,nx
533       WORK(I)    = SUM(PFP(1:nx,I)*RECR(1:nx))  ! L(t)'*r(t)
534620    WORK(nx+I) = SUM(PFP(1:nx,I)*RECRI(1:nx)) ! L(t)'*ri(t)
535
536C Li(t) = -Kit(t)*H(t)
537   	 DO 621 I=1,nx
538       DO 621 J=1,nx
539621    PFP(I,J) =
540     #  - SUM(Kit(ITIME,I,1:iny)*H(IYK(ITIME,1:iny),J,S(ITIME,2)))
541
542C  L(t)'*ri(t) + Li(t)'*r(t)
543	 DO 622 I=1,nx
544622    WORK(nx+I) = WORK(nx+I)+SUM(PFP(1:nx,I)*RECR(1:nx))
545
546C H'*Fsm
547	 DO 625 I=1,nx
548       DO 625 J=1,iny
549625    PHFs(I,J) =
550     #  SUM(H(IYK(ITIME,1:iny),I,S(ITIME,2))*Vis(ITIME,1:iny,J))
551
552	 DO 630 I=1,nx
553630	 RECR(I) = WORK(I) + SUM(PHFs(I,1:iny)*INN(ITIME,1:iny))
554
555C H'*Fim
556	 DO 631 I=1,nx
557       DO 631 J=1,iny
558631    PHFs(I,J) =
559     #  SUM(H(IYK(ITIME,1:iny),I,S(ITIME,2))*Vii(ITIME,1:iny,J))
560
561	 DO 632 I=1,nx
562632	 RECRI(I) = WORK(nx+I) + SUM(PHFs(I,1:iny)*INN(ITIME,1:iny))
563
564C X(d|T) = X(d|d-1) + Psd*r(t-1) + Pid*ri(t-1)
565	 DO 660 I = 1,nx
566660	 XS(ITIME,I) = XT(ITIME,I) + SUM(PT(ITIME,I,1:nx)*RECR(1:nx))
567     +             + SUM(Pi(ITIME,I,1:nx)*RECRI(1:nx))
568
569800   CONTINUE
570
571	DEALLOCATE(Pi,HPs,HPi,Fi,Fs,Fim,Fsm,PHFs,PHFi,FFF,Mi,Ms,Ci,
572	1 Kst,Kit,W1,WORK,WORK1,PFP,APPO,APPO1,COM,RG,FP,HP1,V,CC,HPV,
573     2 X0,P0,RECR,RECRI,RECN,XT,PT,INN,Vinv,Vis,Vii)
574
575	RETURN
576	END
577
578C This is the variance and must be completed!!
579C N(t-1) = H(t)'*Vinv(t)*H(t)+L(t)'*N(t)*L(t)
580c     	 DO 640 I=1,nx
581c       DO 640 J=1,nx
582c640    CC(I,J) = SUM(PFP(1:nx,I)*RECN(1:nx,J)) ! L(t)'*N(t)
583c	 DO 650 I=1,nx
584c	  RECN(I,I) = APPO(I,I) + SUM(CC(I,1:nx)*PFP(1:nx,I))
585c        DO 650 J=1,I-1
586c	  RECN(I,J) = APPO(I,J) + SUM(CC(I,1:nx)*PFP(1:nx,J))
587c650     RECN(J,I) = RECN(I,J)
588C P(t|T) = P(t|t-1) - P(t|t-1)*N(t-1)*P(t|t-1)
589c    	 DO 670 I=1,nx
590c       DO 670 J=1,nx
591c670    CC(I,J) = SUM(PT(ITIME,I,1:nx)*RECN(1:nx,J)) ! P(t|t-1)*N(t-1)
592c
593c	 DO 680 I=1,nx
594c	  PS(ITIME,I,I) = PT(ITIME,I,I) - SUM(CC(I,1:nx)*PT(ITIME,1:nx,I))
595c       DO 680 J=1,I-1
596c	  PS(ITIME,I,J) = PT(ITIME,I,J) - SUM(CC(I,1:nx)*PT(ITIME,1:nx,J))
597c680    PS(ITIME,J,I) = PS(ITIME,I,J)
598
599