Re: [DynareDev] [Dynare Git] Dynare branch, master, updated. a83c0e68c70df4481ab269a1de493eabb804df71
I wonder whether ierating the product of the transition matrix in forecast is not less accurate than taking the power of the matrix. But I agree that it must be faster for large matrices. Maybe we should decompose the matrix before the loop and only take the power of the eigenvalues.
Best
Michel
Stephane Adjemian wrote:
This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "Dynare".
The branch, master has been updated via a83c0e68c70df4481ab269a1de493eabb804df71 (commit) from a828c3d2c1edf1bc4465561a425a9591910c0dc3 (commit)
Those revisions listed above that are new to this repository have not appeared on any other notification email; so we list those revisions in full, below.
- Log -----------------------------------------------------------------
commit a83c0e68c70df4481ab269a1de493eabb804df71 Author: Stéphane Adjemian (Charybdis) stephane.adjemian@ens.fr Date: Fri Feb 12 12:31:32 2010 +0100
Bug fix and efficiency changes in smoothing routines.diff --git a/matlab/DiffuseKalmanSmoother1.m b/matlab/DiffuseKalmanSmoother1.m index 8e624ad..ef64185 100644 --- a/matlab/DiffuseKalmanSmoother1.m +++ b/matlab/DiffuseKalmanSmoother1.m @@ -59,7 +59,7 @@ spstar = size(Pstar1); v = zeros(pp,smpl); a = zeros(mm,smpl+1); atilde = zeros(mm,smpl); -aK = zeros(nk,mm,smpl+1); +aK = zeros(nk,mm,smpl+nk); iF = zeros(pp,pp,smpl); Fstar = zeros(pp,pp,smpl); iFinf = zeros(pp,pp,smpl); @@ -71,11 +71,11 @@ P = zeros(mm,mm,smpl+1); Pstar = zeros(spstar(1),spstar(2),smpl+1); Pstar(:,:,1) = Pstar1; Pinf = zeros(spinf(1),spinf(2),smpl+1); Pinf(:,:,1) = Pinf1; crit = options_.kalman_tol; -crit1 = 1.e-8; +crit1 = 1.e-8; steady = smpl; rr = size(Q,1); QQ = R*Q*transpose(R); -QRt = Q*transpose(R); +QRt = Q*transpose(R); alphahat = zeros(mm,smpl); etahat = zeros(rr,smpl); r = zeros(mm,smpl+1); @@ -99,7 +99,7 @@ while rank(Pinf(:,:,t+1),crit1) & t<smpl a(:,t+1) = T*atilde(:,t); aK(1,:,t+1) = a(:,t+1); for jnk=2:nk,
aK(jnk,:,t+jnk) = T^(jnk-1)*a(:,t+1);
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));@@ -129,9 +129,9 @@ while notsteady & t<smpl K(:,:,t) = T*PZI; L(:,:,t) = T-K(:,:,t)*Z; a(:,t+1) = T*atilde(:,t);
- aK(1,:,t+1) = a(:,t+1);
- aK(1,:,t+1) = a(:,t+1); for jnk=2:nk,
aK(jnk,:,t+jnk) = T^(jnk-1)*a(:,t+1);
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));@@ -154,7 +154,7 @@ while t<smpl a(:,t+1) = T*a(:,t) + K_s*v(:,t); aK(1,:,t+1) = a(:,t+1); for jnk=2:nk,
aK(jnk,:,t+jnk) = T^(jnk-1)*a(:,t+1);
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));end t = smpl+1; @@ -174,4 +174,4 @@ if d alphahat(:,t) = a(:,t) + Pstar(:,:,t)*r0(:,t) + Pinf(:,:,t)*r1(:,t); etahat(:,t) = QRt*r0(:,t); end -end +end \ No newline at end of file diff --git a/matlab/DiffuseKalmanSmoother1_Z.m b/matlab/DiffuseKalmanSmoother1_Z.m index f5614af..debf47c 100644 --- a/matlab/DiffuseKalmanSmoother1_Z.m +++ b/matlab/DiffuseKalmanSmoother1_Z.m @@ -105,8 +105,8 @@ while rank(Pinf(:,:,t+1),crit1) & t<smpl a(:,t+1) = T*atilde(:,t); aK(1,:,t+1) = a(:,t+1); % isn't a meaningless as long as we are in the diffuse part? MJ
- for jnk=2:nk,
aK(jnk,:,t+jnk) = T^(jnk-1)*a(:,t+1);
- for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));@@ -138,10 +138,13 @@ while notsteady & t<smpl L(:,:,t) = T-K(:,:,t)*Z; a(:,t+1) = T*atilde(:,t); Pf = P(:,:,t);
- aK(1,:,t+1) = a(:,t+1); for jnk=1:nk, Pf = T*Pf*T' + QQ;
aK(jnk,:,t+jnk) = T^jnk*atilde(:,t); PK(jnk,:,:,t+jnk) = Pf;
if jnk>1aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));end@@ -162,10 +165,13 @@ while t<smpl atilde(:,t) = a(:,t) + PZI*v(:,t); a(:,t+1) = T*atilde(:,t); Pf = P(:,:,t);
- for jnk=1:nk,
- aK(1,:,t+1) = a(:,t+1);
- for jnk=1:nk Pf = T*Pf*T' + QQ;
aK(jnk,:,t+jnk) = T^jnk*atilde(:,t); PK(jnk,:,:,t+jnk) = Pf;
if jnk>1aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));endend t = smpl+1; diff --git a/matlab/DiffuseKalmanSmoother3.m b/matlab/DiffuseKalmanSmoother3.m index 1127cf6..514890c 100644 --- a/matlab/DiffuseKalmanSmoother3.m +++ b/matlab/DiffuseKalmanSmoother3.m @@ -165,13 +165,14 @@ while newRank & t < smpl Pstar(:,:,t)=tril(Pstar(:,:,t))+transpose(tril(Pstar(:,:,t),-1)); end end
- a1(:,t+1) = T*a(:,t);
- for jnk=1:nk,
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
- a1(:,t+1) = T*a(:,t);
- aK(1,:,t+1) = a1(:,t+1);
- for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
- P0=Pinf(:,:,t+1);
- P0 = Pinf(:,:,t+1); if newRank, %newRank = any(diag(P0(mf,mf))>crit); newRank = rank(P0,crit1);
@@ -207,12 +208,15 @@ while notsteady & t<smpl P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1)); end end
- a1(:,t+1) = T*a(:,t);
- a1(:,t+1) = T*a(:,t); Pf = P(:,:,t);
- aK(1,:,t+1) = a1(:,t+1); for jnk=1:nk, Pf = T*Pf*T' + QQ;
aK(jnk,:,t+jnk) = T^jnk*a(:,t); PK(jnk,:,:,t+jnk) = Pf;
if jnk>1aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));end@@ -239,10 +243,13 @@ while t<smpl end a1(:,t+1) = T*a(:,t); Pf = P(:,:,t);
- for jnk=1:nk,
- aK(1,:,t+1) = a1(:,t+1);
- for jnk=1:nk Pf = T*Pf*T' + QQ;
aK(jnk,:,t+jnk) = T^jnk*a(:,t); PK(jnk,:,:,t+jnk) = Pf;
if jnk>1aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));endend ri=zeros(mm,1); diff --git a/matlab/DiffuseKalmanSmoother3_Z.m b/matlab/DiffuseKalmanSmoother3_Z.m index b30e3b3..c6fa2a2 100644 --- a/matlab/DiffuseKalmanSmoother3_Z.m +++ b/matlab/DiffuseKalmanSmoother3_Z.m @@ -91,7 +91,6 @@ Ki = zeros(mm,pp,smpl); Kstar = zeros(mm,pp,smpl_diff); P = zeros(mm,mm,smpl+1); P1 = P; -aK = zeros(nk,mm,smpl+nk); PK = zeros(nk,mm,mm,smpl+nk); Pstar = zeros(spstar(1),spstar(2),smpl_diff+1); Pstar(:,:,1) = Pstar1; Pinf = zeros(spinf(1),spinf(2),smpl_diff+1); Pinf(:,:,1) = Pinf1; @@ -162,9 +161,10 @@ while newRank & t < smpl Pstar(:,:,t)=tril(Pstar(:,:,t))+tril(Pstar(:,:,t),-1)'; end end
- a1(:,t+1) = T*a(:,t);
- for jnk=1:nk,
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
- a1(:,t+1) = T*a(:,t);
- aK(1,:,t+1) = a1(:,t+1);
- for jnk=2:nk
aK(jnk,:,t+jnk) = T*squueze(aK(jnk-1,:,t+jnk-1));@@ -205,11 +205,14 @@ while notsteady & t<smpl end end a1(:,t+1) = T*a(:,t);
- Pf = P(:,:,t);
- Pf = P(:,:,t);
- aK(1,:,t+1) = a1(:,t+1) for jnk=1:nk, Pf = T*Pf*T' + QQ;
aK(jnk,:,t+jnk) = T^jnk*a(:,t); PK(jnk,:,:,t+jnk) = Pf;
if jnk>1aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));enddiff --git a/matlab/DiffuseKalmanSmootherH1.m b/matlab/DiffuseKalmanSmootherH1.m index 362a7a6..8e16023 100644 --- a/matlab/DiffuseKalmanSmootherH1.m +++ b/matlab/DiffuseKalmanSmootherH1.m @@ -98,8 +98,9 @@ while rank(Pinf(:,:,t+1),crit1) & t<smpl atilde(:,t) = a(:,t) + PZI*v(:,t); Kinf(:,:,t) = T*PZI; a(:,t+1) = T*atilde(:,t);
- for jnk=1:nk,
aK(jnk,:,t+jnk) = T^(jnk-1)*a(:,t+1);
- aK(1,:,t+1) = a(:,t+1);
- for jnk=1:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));@@ -118,25 +119,26 @@ Pinf = Pinf(:,:,1:d); notsteady = 1; while notsteady & t<smpl t = t+1;
- v(:,t) = Y(:,t) - a(mf,t) - trend(:,t);
- v(:,t) = Y(:,t) - a(mf,t) - trend(:,t); P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1)); if rcond(P(mf,mf,t) + H) < crit
return- end
return- end iF(:,:,t) = inv(P(mf,mf,t) + H); PZI = P(:,mf,t)*iF(:,:,t); atilde(:,t) = a(:,t) + PZI*v(:,t); K(:,:,t) = T*PZI; L(:,:,t) = T-K(:,:,t)*Z; a(:,t+1) = T*atilde(:,t);
- for jnk=1:nk,
aK(jnk,:,t+jnk) = T^(jnk-1)*a(:,t+1);
- aK(1,:,t+1) = a(:,t+1);
- for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));end if t<smpl
- PZI_s = PZI;
- PZI_s = PZI; K_s = K(:,:,t); iF_s = iF(:,:,t); P_s = P(:,:,t+1);
@@ -151,8 +153,9 @@ while t<smpl v(:,t) = Y(:,t) - a(mf,t) - trend(:,t); atilde(:,t) = a(:,t) + PZI_s*v(:,t); a(:,t+1) = T*atilde(:,t);
- for jnk=1:nk,
aK(jnk,:,t+jnk) = T^(jnk-1)*a(:,t+1);
- aK(1,:,t+1) = a(:,t+1);
- for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));end t = smpl+1; diff --git a/matlab/DiffuseKalmanSmootherH1_Z.m b/matlab/DiffuseKalmanSmootherH1_Z.m index bbcdfd9..65fcd44 100644 --- a/matlab/DiffuseKalmanSmootherH1_Z.m +++ b/matlab/DiffuseKalmanSmootherH1_Z.m @@ -108,8 +108,8 @@ while rank(Pinf(:,:,t+1),crit1) & t<smpl a(:,t+1) = T*atilde(:,t); aK(1,:,t+1) = a(:,t+1); % isn't a meaningless as long as we are in the diffuse part? MJ
- for jnk=2:nk,
aK(jnk,:,t+jnk) = T^(jnk-1)*a(:,t+1);
- for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));@@ -141,10 +141,13 @@ while notsteady & t<smpl L(:,:,t) = T-K(:,:,t)*Z; a(:,t+1) = T*atilde(:,t); Pf = P(:,:,t);
- for jnk=1:nk,
- aK(1,:,t+1) = a(:,t+1);
- for jnk=1:nk Pf = T*Pf*T' + QQ;
aK(jnk,:,t+jnk) = T^jnk*atilde(:,t); PK(jnk,:,:,t+jnk) = Pf;
if jnk>1aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));end@@ -165,10 +168,13 @@ while t<smpl atilde(:,t) = a(:,t) + PZI*v(:,t); a(:,t+1) = T*atilde(:,t); Pf = P(:,:,t);
- for jnk=1:nk,
- aK(1,:,t+1) = a(:,t+1);
- for jnk=1:nk Pf = T*Pf*T' + QQ;
aK(jnk,:,t+jnk) = T^jnk*atilde(:,t); PK(jnk,:,:,t+jnk) = Pf;
if jnk>1aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));endend t = smpl+1; diff --git a/matlab/DiffuseKalmanSmootherH3.m b/matlab/DiffuseKalmanSmootherH3.m index c7dc47f..afa98b5 100644 --- a/matlab/DiffuseKalmanSmootherH3.m +++ b/matlab/DiffuseKalmanSmootherH3.m @@ -163,9 +163,10 @@ while newRank & t < smpl Pstar(:,:,t)=tril(Pstar(:,:,t))+transpose(tril(Pstar(:,:,t),-1)); end end
- a1(:,t+1) = T*a(:,t);
- for jnk=1:nk,
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
- a1(:,t+1) = T*a(:,t);
- aK(1,:,t+1) = a1(:,t+1)
- for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));@@ -206,8 +207,9 @@ while notsteady & t<smpl end end a1(:,t+1) = T*a(:,t);
- for jnk=1:nk,
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
- aK(1,:,t+1) = a1(:,t+1);
- for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));@@ -233,8 +235,9 @@ while t<smpl end end a1(:,t+1) = T*a(:,t);
- for jnk=1:nk,
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
- aK(1,:,t+1) = a1(:,t+1);
- for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));end ri=zeros(mm,1); diff --git a/matlab/DiffuseKalmanSmootherH3_Z.m b/matlab/DiffuseKalmanSmootherH3_Z.m index c01b833..5eaf948 100644 --- a/matlab/DiffuseKalmanSmootherH3_Z.m +++ b/matlab/DiffuseKalmanSmootherH3_Z.m @@ -91,7 +91,6 @@ L0 = zeros(mm,mm,pp,smpl_diff); Kstar = zeros(mm,pp,smpl_diff); P = zeros(mm,mm,smpl+1); P1 = P; -aK = zeros(nk,mm,smpl+nk); PK = zeros(nk,mm,mm,smpl+nk); Pstar = zeros(spstar(1),spstar(2),smpl_diff+1); Pstar(:,:,1) = Pstar1; Pinf = zeros(spinf(1),spinf(2),smpl_diff+1); Pinf(:,:,1) = Pinf1; @@ -163,9 +162,10 @@ while newRank & t < smpl Pstar(:,:,t)=tril(Pstar(:,:,t))+tril(Pstar(:,:,t),-1)'; end end
- a1(:,t+1) = T*a(:,t);
- for jnk=1:nk,
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
- a1(:,t+1) = T*a(:,t);
- aK(1,:,t+1) = a1(:,t+1);
- for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));@@ -207,10 +207,13 @@ while notsteady & t<smpl end a1(:,t+1) = T*a(:,t); Pf = P(:,:,t);
- aK(1,:,t+1) = a1(:,t+1); for jnk=1:nk, Pf = T*Pf*T' + QQ;
aK(jnk,:,t+jnk) = T^jnk*a(:,t); PK(jnk,:,:,t+jnk) = Pf;
if jnk>1aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1)) ;end@@ -238,11 +241,14 @@ while t<smpl end end a1(:,t+1) = T*a(:,t);
- Pf = P(:,:,t);
- Pf = P(:,:,t);
- aK(1,:,t+1) = a1(:,t+1); for jnk=1:nk, Pf = T*Pf*T' + QQ;
aK(jnk,:,t+jnk) = T^jnk*a(:,t); PK(jnk,:,:,t+jnk) = Pf;
if jnk>1aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));endend ri=zeros(mm,1); diff --git a/matlab/DiffuseKalmanSmootherH3corr.m b/matlab/DiffuseKalmanSmootherH3corr.m index 800362c..cfbd826 100644 --- a/matlab/DiffuseKalmanSmootherH3corr.m +++ b/matlab/DiffuseKalmanSmootherH3corr.m @@ -122,8 +122,9 @@ while newRank & t < smpl end end a1(:,t+1) = T*a(:,t);
- for jnk=1:nk,
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
- aK(1,:,t+1) = a1(:,t+1)
- for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));@@ -159,8 +160,9 @@ while notsteady & t<smpl end end a1(:,t+1) = T*a(:,t);
- for jnk=1:nk,
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
- aK(1,:,t+1) = a1(:,t+1);
- for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));@@ -186,8 +188,9 @@ while t<smpl end end a1(:,t+1) = T*a(:,t);
- for jnk=1:nk,
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
- aK(1,:,t+1) = a1(:,t+1);
- for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));end
diff --git a/matlab/missing_DiffuseKalmanSmoother1.m b/matlab/missing_DiffuseKalmanSmoother1.m index da7d840..5d76505 100644 --- a/matlab/missing_DiffuseKalmanSmoother1.m +++ b/matlab/missing_DiffuseKalmanSmoother1.m @@ -60,7 +60,7 @@ spstar = size(Pstar1); v = zeros(pp,smpl); a = zeros(mm,smpl+1); atilde = zeros(mm,smpl); -aK = zeros(nk,mm,smpl+1); +aK = zeros(nk,mm,smpl+nk); iF = zeros(pp,pp,smpl); Fstar = zeros(pp,pp,smpl); iFinf = zeros(pp,pp,smpl); @@ -113,9 +113,9 @@ while rank(Pinf(:,:,t+1),crit1) & t<smpl Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T)-T*Pinf(:,mf1,t)* ... transpose(Kinf(:,di,t)); end
- aK(1,:,t+1) = a(:,t+1);
- for jnk=2:nk,
aK(jnk,:,t+jnk) = T^(jnk-1)*a(:,t+1);
- aK(1,:,t+1) = a(:,t+1);
- for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));end d = t; @@ -150,7 +150,7 @@ while notsteady & t<smpl end aK(1,:,t+1) = a(:,t+1); for jnk=2:nk,
aK(jnk,:,t+jnk) = T^(jnk-1)*a(:,t+1);
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));diff --git a/matlab/missing_DiffuseKalmanSmoother1_Z.m b/matlab/missing_DiffuseKalmanSmoother1_Z.m index 097f1de..e640fa5 100644 --- a/matlab/missing_DiffuseKalmanSmoother1_Z.m +++ b/matlab/missing_DiffuseKalmanSmoother1_Z.m @@ -121,7 +121,7 @@ while rank(Pinf(:,:,t+1),crit1) & t<smpl aK(1,:,t+1) = a(:,t+1); % isn't a meaningless as long as we are in the diffuse part? MJ for jnk=2:nk,
aK(jnk,:,t+jnk) = T^(jnk-1)*a(:,t+1);
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));end d = t; @@ -157,10 +157,13 @@ while notsteady & t<smpl end a(:,t+1) = T*atilde(:,t); Pf = P(:,:,t);
- for jnk=1:nk,
- aK(1,:,t+1) = a(:,t+1);
- for jnk=1:nk Pf = T*Pf*T' + QQ;
aK(jnk,:,t+jnk) = T^jnk*atilde(:,t); PK(jnk,:,:,t+jnk) = Pf;
if jnk>1aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));endend diff --git a/matlab/missing_DiffuseKalmanSmoother3.m b/matlab/missing_DiffuseKalmanSmoother3.m index 8f2a041..c4c6df3 100644 --- a/matlab/missing_DiffuseKalmanSmoother3.m +++ b/matlab/missing_DiffuseKalmanSmoother3.m @@ -91,11 +91,11 @@ Pinf = zeros(spinf(1),spinf(2),smpl_diff+1); Pinf(:,:,1) = Pinf1; Pstar1 = Pstar; Pinf1 = Pinf; crit = options_.kalman_tol; -crit1 = 1.e-6; +crit1 = 1.e-6; steady = smpl; rr = size(Q,1); QQ = R*Q*transpose(R); -QRt = Q*transpose(R); +QRt = Q*transpose(R); alphahat = zeros(mm,smpl); etahat = zeros(rr,smpl); r = zeros(mm,smpl+1); @@ -167,16 +167,17 @@ while newRank & t < smpl Pstar(:,:,t)=tril(Pstar(:,:,t))+transpose(tril(Pstar(:,:,t),-1)); end end
- a1(:,t+1) = T*a(:,t);
- for jnk=1:nk,
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
- a1(:,t+1) = T*a(:,t);
- aK(1,:,t+1) = a1(:,t+1);
- for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));
newRank = rank(P0,crit1);
newRank = rank(P0,crit1);end
@@ -212,10 +213,13 @@ while notsteady & t<smpl end a1(:,t+1) = T*a(:,t); Pf = P(:,:,t);
- aK(1,:,t+1) = a1(:,t+1); for jnk=1:nk, Pf = T*Pf*T' + QQ;
aK(jnk,:,t+jnk) = T^jnk*a1(:,t); PK(jnk,:,:,t+jnk) = Pf;
if jnk>1aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));enddiff --git a/matlab/missing_DiffuseKalmanSmoother3_Z.m b/matlab/missing_DiffuseKalmanSmoother3_Z.m index 88cde2f..8a7313f 100644 --- a/matlab/missing_DiffuseKalmanSmoother3_Z.m +++ b/matlab/missing_DiffuseKalmanSmoother3_Z.m @@ -92,7 +92,6 @@ Ki = zeros(mm,pp,smpl); Kstar = zeros(mm,pp,smpl_diff); P = zeros(mm,mm,smpl+1); P1 = P; -aK = zeros(nk,mm,smpl+nk); PK = zeros(nk,mm,mm,smpl+nk); Pstar = zeros(spstar(1),spstar(2),smpl_diff+1); Pstar(:,:,1) = Pstar1; Pinf = zeros(spinf(1),spinf(2),smpl_diff+1); Pinf(:,:,1) = Pinf1; @@ -164,9 +163,10 @@ while newRank & t < smpl Pstar(:,:,t)=tril(Pstar(:,:,t))+tril(Pstar(:,:,t),-1)'; end end
- a1(:,t+1) = T*a(:,t);
- for jnk=1:nk,
aK(jnk,:,t+jnk) = T^jnk*a(:,t);
- a1(:,t+1) = T*a(:,t);
- aK(1,:,t+1) = a1(:,t+1);
- for jnk=2:nk
aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));@@ -209,10 +209,13 @@ while notsteady & t<smpl end a1(:,t+1) = T*a(:,t); Pf = P(:,:,t);
- for jnk=1:nk,
- aK(1,:,t+1) = a1(:,t+1);
- for jnk=1:nk Pf = T*Pf*T' + QQ;
aK(jnk,:,t+jnk) = T^jnk*a(:,t); PK(jnk,:,:,t+jnk) = Pf;
if jnk>1aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1));end
Summary of changes: matlab/DiffuseKalmanSmoother1.m | 16 ++++++++-------- matlab/DiffuseKalmanSmoother1_Z.m | 16 +++++++++++----- matlab/DiffuseKalmanSmoother3.m | 23 +++++++++++++++-------- matlab/DiffuseKalmanSmoother3_Z.m | 15 +++++++++------ matlab/DiffuseKalmanSmootherH1.m | 23 +++++++++++++---------- matlab/DiffuseKalmanSmootherH1_Z.m | 18 ++++++++++++------ matlab/DiffuseKalmanSmootherH3.m | 17 ++++++++++------- matlab/DiffuseKalmanSmootherH3_Z.m | 20 +++++++++++++------- matlab/DiffuseKalmanSmootherH3corr.m | 15 +++++++++------ matlab/missing_DiffuseKalmanSmoother1.m | 10 +++++----- matlab/missing_DiffuseKalmanSmoother1_Z.m | 9 ++++++--- matlab/missing_DiffuseKalmanSmoother3.m | 18 +++++++++++------- matlab/missing_DiffuseKalmanSmoother3_Z.m | 15 +++++++++------ 13 files changed, 131 insertions(+), 84 deletions(-)
hooks/post-receive
-
Michel Juillard