Commit cf8213f7 authored by Marco Ratto's avatar Marco Ratto Committed by Stéphane Adjemian (Charybdis)

fix to the Kitigawa transformation that allows to reduce the computing time of...

fix to the Kitigawa transformation that allows to reduce the computing time of the likelihood in large models, with a lot of static variables, by 30-50%.
This fixes the bug in e97e5c34 that led to 2f9dc092
It is tested and completely fixes the issue highlighted in #1312.
parent 60efdd49
......@@ -183,7 +183,7 @@ elseif options_.lik_init == 3 % Diffuse Kalman filter
Z = [Z, eye(vobs)];
end
end
[Pstar,Pinf] = compute_Pinf_Pstar(mf,T,R,Q,options_.qz_criterium,oo_.dr.restrict_var_list);
[Pstar,Pinf] = compute_Pinf_Pstar(mf,T,R,Q,options_.qz_criterium);
elseif options_.lik_init == 4 % Start from the solution of the Riccati equation.
[err, Pstar] = kalman_steady_state(transpose(T),R*Q*transpose(R),transpose(build_selection_matrix(mf,np,vobs)),H);
mexErrCheck('kalman_steady_state',err);
......
......@@ -48,8 +48,24 @@ function [Pstar,Pinf] = compute_Pinf_Pstar(mf,T,R,Q,qz_criterium, restrict_colum
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
np = size(T,1);
if nargin == 6
indx = restrict_columns;
indx0=find(~ismember([1:np],indx));
else
% indx=(find(max(abs(T))>0));
% indx0=(find(max(abs(T))==0));
indx=(find(max(abs(T))>=1.e-10));
indx0=(find(max(abs(T))<1.e-10));
end
np0=length(indx0);
Tbkp = T;
T0=T(indx0,indx); %static variables vs. dynamic ones
R0=R(indx0,:); % matrix of shocks for static variables
% perform Kitagawa transformation
% perform Kitagawa transformation only for non-zero columns of T
T=T(indx,indx);
R=R(indx,:);
np = size(T,1);
[QT,ST] = schur(T);
e1 = abs(ordeig(ST)) > 2-qz_criterium;
[QT,ST] = ordschur(QT,ST,e1);
......@@ -59,7 +75,6 @@ nk1 = nk+1;
Pstar = zeros(np,np);
R1 = QT'*R;
B = R1*Q*R1';
% computes variance of stationary block (lower right)
i = np;
while i >= nk+2
if ST(i,i-1) == 0
......@@ -100,13 +115,71 @@ if i == nk+1
Pstar(nk1,nk1)=(B(nk1,nk1)+c)/(1-ST(nk1,nk1)*ST(nk1,nk1));
end
if np0
ST1=ST;
% now I recover stationarized static variables
% using
% ss = s-A*z and
% z-z(-1) (growth rates of unit roots) only depends on stationary variables
Pstar = blkdiag(zeros(np0),Pstar);
ST = [zeros(length(Pstar),length(indx0)) [T0*QT ;ST]];
R1 = [R0; R1];
% build the matrix for stationarized variables
STinf = ST(np0+1:np0+nk,np0+1:np0+nk);
iSTinf = inv(STinf);
ST0=ST;
ST0(:,1:np0+nk)=0; % stationarized static + 1st difference only respond to lagged stationary states
ST00 = ST(1:np0,np0+1:np0+nk);
% A\B is the matrix division of A into B, which is roughly the
% same as INV(A)*B
ST0(1:np0,np0+nk+1:end) = ST(1:np0,np0+nk+1:end)-ST00*(iSTinf*ST(np0+1:np0+nk,np0+nk+1:end)); % snip non-stationary part
R10 = R1;
R10(1:np0,:) = R1(1:np0,:)-ST00*(iSTinf*R1(np0+1:np0+nk,:)); % snip non-stationary part
% kill non-stationary part before projecting Pstar
ST0(np0+1:np0+nk,:)=0;
R10(np0+1:np0+nk,:)=0; % is this questionable???? IT HAS TO in order to match Michel's version!!!
% project Pstar onto static x
Pstar = ...
ST0*Pstar*ST0'+ ...
R10*Q*R10';
% QT(1:np0,np0+1:np0+nk) = QT(1:np0,np0+1:np0+nk)+ST(1:np0,np0+1:np0+nk); %%% is this questionable ????
% reorder QT entries
else
STinf = ST(np0+1:np0+nk,np0+1:np0+nk);
end
% stochastic trends with no influence on observed variables are
% arbitrarily initialized to zero
Pinf = zeros(np,np);
Pinf(1:nk,1:nk) = eye(nk);
for k = 1:nk
if norm(QT(mf,:)*ST(:,k)) < 1e-8
Pinf(k,k) = 0;
if np0
STtriu = STinf-eye(nk);
% A\B is the matrix division of A into B, which is roughly the
% same as INV(A)*B
STinf0 = ST00*(eye(nk)-iSTinf*STtriu);
Pinf = blkdiag(zeros(np0),Pinf);
QT = blkdiag(eye(np0),QT);
QTinf = QT;
QTinf(1:np0,np0+1:np0+nk) = STinf0;
QTinf([indx0(:); indx(:)],:) = QTinf;
STinf1 = [zeros(np0+np,np0) [STinf0; eye(nk); zeros(np-nk,nk)] zeros(np0+np,np-nk)];
for k = 1:nk
if norm(QTinf(mf,:)*ST([indx0(:); indx(:)],k+np0)) < 1e-8
Pinf(k+np0,k+np0) = 0;
end
end
Pinf = STinf1*Pinf*STinf1';
QT([indx0(:); indx(:)],:) = QT;
else
for k = 1:nk
if norm(QT(mf,:)*ST(:,k)) < 1e-8
Pinf(k+np0,k+np0) = 0;
end
end
end
......
......@@ -372,7 +372,7 @@ switch DynareOptions.lik_init
error(['The model requires Diffuse filter, but you specified a different Kalman filter. You must set options_.kalman_algo ' ...
'to 0 (default), 3 or 4'])
end
[Pstar,Pinf] = compute_Pinf_Pstar(Z,T,R,Q,DynareOptions.qz_criterium,[1:length(T)]);
[Pstar,Pinf] = compute_Pinf_Pstar(Z,T,R,Q,DynareOptions.qz_criterium);
Z =zeros(length(BayesInfo.mf),size(T,1));
for i = 1:length(BayesInfo.mf)
Z(i,BayesInfo.mf(i))=1;
......
Markdown is supported
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