Commit a34ce9e3 authored by michel's avatar michel
Browse files

v4 bug corrections for new diffuse algorithm

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1684 ac1d8469-bf42-47a9-8791-bf33cf982152
parent d8d191bb
......@@ -83,7 +83,7 @@ while rank(Pinf(:,:,t+1),crit1) & t<smpl
Linf(:,:,t) = T - Kinf(:,:,t)*Z;
Fstar(:,:,t) = Z*Pstar(:,:,t)*Z';
Kstar(:,:,t) = (T*Pstar(:,:,t)*Z'-Kinf(:,:,t)*Fstar(:,:,t))*iFinf(:,:,t);
Pstar(:,:,t+1) = T*Pstar(:,:,t)*T'-T*Pstar(:,:,t)*Z'*Kinf(:,:,t)'-Kinf(:,:,t)*F*Kstar(:,:,t) + QQ;
Pstar(:,:,t+1) = T*Pstar(:,:,t)*T'-T*Pstar(:,:,t)*Z'*Kinf(:,:,t)'-Kinf(:,:,t)*F*Kstar(:,:,t)' + QQ;
Pinf(:,:,t+1) = T*Pinf(:,:,t)*T'-T*Pinf(:,:,t)*Z'*Kinf(:,:,t)';
end
d = t;
......@@ -99,7 +99,7 @@ while notsteady & t<smpl
t = t+1;
v(:,t) = Y(:,t) - Z*a(:,t);
P(:,:,t)=tril(P(:,:,t))+transpose(tril(P(:,:,t),-1));
F = Z*P(:,:,t)*Z'
F = Z*P(:,:,t)*Z';
if rcond(F) < crit
return
end
......@@ -150,7 +150,7 @@ if d
etahat(:,t) = QRt*r0(:,t);
end
r0_0 = Linf(:,:,1)'*r0(:,1);
r1_0 = Z'*(iFinf(:,:,1)*v(:,1)-Kstar(:,:,1)*r0(:,1)) + Linf(:,:,1)'*r1(:,1);
r1_0 = Z'*(iFinf(:,:,1)*v(:,1)-Kstar(:,:,1)'*r0(:,1)) + Linf(:,:,1)'*r1(:,1);
alphahat(:,1) = a(:,1) + Pstar(:,:,1)*r0_0 + Pinf(:,:,1)*r1_0;
etahat(:,1) = QRt*r0(:,1);
else
......
......@@ -253,7 +253,7 @@ if d
L0(:,:,i,1)'*r0_0 + Linf(:,:,i,1)'*r1_0;
r0_0 = Linf(:,:,i,1)'*r0_0;
elseif Fstar(i,1) > crit, % step needed when Finf=0
r0_0=transpose(Z(i,:))/Fstar(i,1)*v(i,1)+Li(:,:,i,1)'*r0_0;
r0_0=Z(i,:)'/Fstar(i,1)*v(i,1)+Li(:,:,i,1)'*r0_0;
end
end
alphahat(:,1) = a1(:,1) + Pstar1(:,:,1)*r0_0 + Pinf1(:,:,1)*r1_0;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment