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)); end Linf(:,:,t) = T - Kinf(:,:,t)*Z; Fstar(:,:,t) = Pstar(mf,mf,t); @@ -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)); end P(:,:,t+1) = T*P(:,:,t)*transpose(T)-T*P(:,mf,t)*transpose(K(:,:,t)) + QQ; notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit); @@ -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 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)); end Linf(:,:,t) = T - Kinf(:,:,t)*Z; Fstar(:,:,t) = Z*Pstar(:,:,t)*Z'; @@ -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>1 + aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1)); + end end P(:,:,t+1) = T*P(:,:,t)*T'-T*P(:,:,t)*Z'*K(:,:,t)' + QQ; notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit); @@ -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>1 + aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1)); + end end end 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)); end Pstar(:,:,t+1) = T*Pstar(:,:,t)*transpose(T)+ QQ; Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T); - 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>1 + aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1)); + end end P(:,:,t+1) = T*P(:,:,t)*transpose(T) + QQ; notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit); @@ -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>1 + aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1)); + end end end 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)); end Pstar(:,:,t+1) = T*Pstar(:,:,t)*T'+ QQ; Pinf(:,:,t+1) = T*Pinf(:,:,t)*T'; @@ -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>1 + aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1)); + end end P(:,:,t+1) = T*P(:,:,t)*T' + QQ; % notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit); diff --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)); end Linf(:,:,t) = T - Kinf(:,:,t)*Z; Fstar(:,:,t) = Pstar(mf,mf,t) + H; @@ -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 P(:,:,t+1) = T*P(:,:,t)*transpose(T)-T*P(:,mf,t)*transpose(K(:,:,t)) + QQ; notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit); 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 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)); end Linf(:,:,t) = T - Kinf(:,:,t)*Z; Fstar(:,:,t) = Z*Pstar(:,:,t)*Z' + H; @@ -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>1 + aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1)); + end end P(:,:,t+1) = T*P(:,:,t)*T'-T*P(:,:,t)*Z'*K(:,:,t)' + QQ; notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit); @@ -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>1 + aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1)); + end end end 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)); end Pstar(:,:,t+1) = T*Pstar(:,:,t)*transpose(T)+ QQ; Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T); @@ -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)); end P(:,:,t+1) = T*P(:,:,t)*transpose(T) + QQ; notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit); @@ -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 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)); end Pstar(:,:,t+1) = T*Pstar(:,:,t)*T'+ QQ; Pinf(:,:,t+1) = T*Pinf(:,:,t)*T'; @@ -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>1 + aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1)) ; + end end P(:,:,t+1) = T*P(:,:,t)*T' + QQ; notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit); @@ -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>1 + aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1)); + end end end 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)); end Pstar(:,:,t+1) = T*Pstar(:,:,t)*transpose(T)+ QQ; Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T); @@ -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)); end P(:,:,t+1) = T*P(:,:,t)*transpose(T) + QQ; notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit); @@ -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 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 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)); end P(:,:,t+1) = T*P(:,:,t)*transpose(T)-T*P(:,mf,t)*transpose(K(:,:,t)) + QQ; % notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit); 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 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>1 + aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1)); + end end % notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit); end 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)); end Pstar(:,:,t+1) = T*Pstar(:,:,t)*transpose(T)+ QQ; Pinf(:,:,t+1) = T*Pinf(:,:,t)*transpose(T); P0=Pinf(:,:,t+1); if newRank, %newRank = any(diag(P0(mf,mf))>crit); - newRank = rank(P0,crit1); + newRank = rank(P0,crit1); end 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>1 + aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1)); + end end P(:,:,t+1) = T*P(:,:,t)*transpose(T) + QQ; % notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit); diff --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)); end Pstar(:,:,t+1) = T*Pstar(:,:,t)*T'+ QQ; Pinf(:,:,t+1) = T*Pinf(:,:,t)*T'; @@ -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>1 + aK(jnk,:,t+jnk) = T*squeeze(aK(jnk-1,:,t+jnk-1)); + end end P(:,:,t+1) = T*P(:,:,t)*T' + QQ; % notsteady = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<crit);
-----------------------------------------------------------------------
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
participants (1)
-
Michel Juillard