v4 bug correction and correct smoother for new diffuse algorithm
git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1682 ac1d8469-bf42-47a9-8791-bf33cf982152time-shift
parent
d01ed3ca16
commit
4443d56136
|
@ -97,7 +97,7 @@ while newRank & t < smpl
|
||||||
v(i,t) = Y(i,t)-Zi*a(:,t);
|
v(i,t) = Y(i,t)-Zi*a(:,t);
|
||||||
Fstar(i,t) = Zi*Pstar(:,:,t)*Zi';
|
Fstar(i,t) = Zi*Pstar(:,:,t)*Zi';
|
||||||
Finf(i,t) = Zi*Pinf(:,:,t)*Zi';
|
Finf(i,t) = Zi*Pinf(:,:,t)*Zi';
|
||||||
Kstar(:,i,t) = Pstar(:,:,t)*Zi;
|
Kstar(:,i,t) = Pstar(:,:,t)*Zi';
|
||||||
if Finf(i,t) > crit & newRank
|
if Finf(i,t) > crit & newRank
|
||||||
icc=icc+1;
|
icc=icc+1;
|
||||||
Kinf(:,i,t) = Pinf(:,:,t)*Zi';
|
Kinf(:,i,t) = Pinf(:,:,t)*Zi';
|
||||||
|
@ -168,7 +168,7 @@ while notsteady & t<smpl
|
||||||
P(:,:,t)=tril(P(:,:,t))+tril(P(:,:,t),-1)';
|
P(:,:,t)=tril(P(:,:,t))+tril(P(:,:,t),-1)';
|
||||||
P1(:,:,t) = P(:,:,t);
|
P1(:,:,t) = P(:,:,t);
|
||||||
for i=1:pp
|
for i=1:pp
|
||||||
Zi = Z(i,:)'
|
Zi = Z(i,:);
|
||||||
v(i,t) = Y(i,t) - Zi*a(:,t);
|
v(i,t) = Y(i,t) - Zi*a(:,t);
|
||||||
Fi(i,t) = Zi*P(:,:,t)*Zi';
|
Fi(i,t) = Zi*P(:,:,t)*Zi';
|
||||||
Ki(:,i,t) = P(:,:,t)*Zi';
|
Ki(:,i,t) = P(:,:,t)*Zi';
|
||||||
|
@ -200,7 +200,7 @@ while t<smpl
|
||||||
t=t+1;
|
t=t+1;
|
||||||
a1(:,t) = a(:,t);
|
a1(:,t) = a(:,t);
|
||||||
for i=1:pp
|
for i=1:pp
|
||||||
Zi = Z(i,:)';
|
Zi = Z(i,:);
|
||||||
v(i,t) = Y(i,t) - Zi*a(:,t);
|
v(i,t) = Y(i,t) - Zi*a(:,t);
|
||||||
if Fi_s(i) > crit
|
if Fi_s(i) > crit
|
||||||
a(:,t) = a(:,t) + Ki_s(:,i)*v(i,t)/Fi_s(i);
|
a(:,t) = a(:,t) + Ki_s(:,i)*v(i,t)/Fi_s(i);
|
||||||
|
|
|
@ -149,7 +149,7 @@ while t < smpl
|
||||||
t = t+1;
|
t = t+1;
|
||||||
Pstar = oldP;
|
Pstar = oldP;
|
||||||
for i=1:pp
|
for i=1:pp
|
||||||
Zi = Z(i,i);
|
Zi = Z(i,:);
|
||||||
v(i) = Y(i,t) - Zi*a;
|
v(i) = Y(i,t) - Zi*a;
|
||||||
Fi = Zi*Pstar*Zi';
|
Fi = Zi*Pstar*Zi';
|
||||||
if Fi > crit
|
if Fi > crit
|
||||||
|
|
|
@ -82,18 +82,71 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R] = Dsge
|
||||||
Pstar = 10*eye(np);
|
Pstar = 10*eye(np);
|
||||||
Pinf = [];
|
Pinf = [];
|
||||||
elseif options_.lik_init == 3 % Diffuse Kalman filter
|
elseif options_.lik_init == 3 % Diffuse Kalman filter
|
||||||
Pstar = zeros(np,np);
|
if options_.kalman_algo < 4
|
||||||
ivs = bayestopt_.var_list_stationary;
|
Pstar = zeros(np,np);
|
||||||
Pstar(ivs,ivs) = lyapunov_symm(T(ivs,ivs),R(ivs,:)*Q* ...
|
ivs = bayestopt_.restrict_var_list_stationary;
|
||||||
transpose(R(ivs,:)));
|
R1 = R(ivs,:);
|
||||||
% Pinf = bayestopt_.Pinf;
|
Pstar(ivs,ivs) = lyapunov_symm(T(ivs,ivs),R1*Q*R1');
|
||||||
% by M. Ratto
|
% Pinf = bayestopt_.Pinf;
|
||||||
RR=T(:,find(~ismember([1:np],ivs)));
|
% by M. Ratto
|
||||||
i=find(abs(RR)>1.e-10);
|
RR=T(:,bayestopt_.restrict_var_list_nonstationary);
|
||||||
R0=zeros(size(RR));
|
i=find(abs(RR)>1.e-10);
|
||||||
R0(i)=sign(RR(i));
|
R0=zeros(size(RR));
|
||||||
Pinf=R0*R0';
|
R0(i)=sign(RR(i));
|
||||||
% by M. Ratto
|
Pinf=R0*R0';
|
||||||
|
% by M. Ratto
|
||||||
|
else
|
||||||
|
[QT,ST] = schur(T);
|
||||||
|
e1 = abs(ordeig(ST)) > 2-options_.qz_criterium;
|
||||||
|
[QT,ST] = ordschur(QT,ST,e1);
|
||||||
|
k = find(abs(ordeig(ST)) > 2-options_.qz_criterium);
|
||||||
|
nk = length(k);
|
||||||
|
nk1 = nk+1;
|
||||||
|
Pinf = zeros(np,np);
|
||||||
|
Pinf(1:nk,1:nk) = eye(nk);
|
||||||
|
Pstar = zeros(np,np);
|
||||||
|
B = QT'*R*Q*R'*QT;
|
||||||
|
for i=np:-1:nk+2
|
||||||
|
if ST(i,i-1) == 0
|
||||||
|
if i == np
|
||||||
|
c = zeros(np-nk,1);
|
||||||
|
else
|
||||||
|
c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+...
|
||||||
|
ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i);
|
||||||
|
end
|
||||||
|
q = eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i);
|
||||||
|
Pstar(nk1:i,i) = q\(B(nk1:i,i)+c);
|
||||||
|
Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)';
|
||||||
|
else
|
||||||
|
if i == np
|
||||||
|
c = zeros(np-nk,1);
|
||||||
|
c1 = zeros(np-nk,1);
|
||||||
|
else
|
||||||
|
c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+...
|
||||||
|
ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i)+...
|
||||||
|
ST(i,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1);
|
||||||
|
c1 = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i-1,i+1:end)')+...
|
||||||
|
ST(i-1,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1)+...
|
||||||
|
ST(i-1,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i);
|
||||||
|
end
|
||||||
|
q = [eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i) -ST(nk1:i,nk1:i)*ST(i,i-1);...
|
||||||
|
-ST(nk1:i,nk1:i)*ST(i-1,i) eye(i-nk)-ST(nk1:i,nk1:i)*ST(i-1,i-1)];
|
||||||
|
z = q\[B(nk1:i,i)+c;B(nk1:i,i-1)+c1];
|
||||||
|
Pstar(nk1:i,i) = z(1:(i-nk));
|
||||||
|
Pstar(nk1:i,i-1) = z(i-nk+1:end);
|
||||||
|
Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)';
|
||||||
|
Pstar(i-1,nk1:i-2) = Pstar(nk1:i-2,i-1)';
|
||||||
|
i = i - 1;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
if i == nk+2
|
||||||
|
c = ST(nk+1,:)*(Pstar(:,nk+2:end)*ST(nk1,nk+2:end)')+ST(nk1,nk1)*ST(nk1,nk+2:end)*Pstar(nk+2:end,nk1);
|
||||||
|
Pstar(nk1,nk1)=(B(nk1,nk1)+c)/(1-ST(nk1,nk1)*ST(nk1,nk1));
|
||||||
|
end
|
||||||
|
|
||||||
|
Z = QT(mf,:);
|
||||||
|
R1 = QT'*R;
|
||||||
|
end
|
||||||
end
|
end
|
||||||
% -----------------------------------------------------------------------------
|
% -----------------------------------------------------------------------------
|
||||||
% 4. Kalman smoother
|
% 4. Kalman smoother
|
||||||
|
@ -115,5 +168,15 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R] = Dsge
|
||||||
end
|
end
|
||||||
elseif options_.kalman_algo == 3
|
elseif options_.kalman_algo == 3
|
||||||
[alphahat,etahat,ahat,aK] = DiffuseKalmanSmoother3(T,R,Q,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
|
[alphahat,etahat,ahat,aK] = DiffuseKalmanSmoother3(T,R,Q,Pinf,Pstar,Y,trend,nobs,np,smpl,mf);
|
||||||
|
elseif options_.kalman_algo == 4
|
||||||
|
data1 = Y - trend;
|
||||||
|
[alphahat,etahat,ahat,aK] = DiffuseKalmanSmoother1_Z(ST,Z,R1,Q,Pinf,Pstar,data1,nobs,np,smpl);
|
||||||
|
alphahat = QT*alphahat;
|
||||||
|
ahat = QT*ahat;
|
||||||
|
elseif options_.kalman_algo == 5
|
||||||
|
data1 = Y - trend;
|
||||||
|
[alphahat,etahat,ahat,aK] = DiffuseKalmanSmoother3_Z(ST,Z,R1,Q,Pinf,Pstar,data1,nobs,np,smpl);
|
||||||
|
alphahat = QT*alphahat;
|
||||||
|
ahat = QT*ahat;
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
|
Loading…
Reference in New Issue