From f441b9d7ae6fa9bb881303c9e105f8b2b8f8faef Mon Sep 17 00:00:00 2001
From: Marco Ratto <marco.ratto@jrc.ec.europa.eu>
Date: Thu, 16 Nov 2017 16:55:58 +0100
Subject: [PATCH] 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
e97e5c340757098a2904b24aaeb5f8011812cdf3 that led to
2f9dc092855bbdb0b3fdd970b6ab551842c5080a It is tested and completely fixes
the issue highlighted in #1312.
(cherry picked from commit cf8213f7a0d71c7520e7e0f0904f826cfc8026f5)
---
matlab/DsgeSmoother.m | 2 +-
matlab/compute_Pinf_Pstar.m | 83 ++++++++++++++++++++++++++++++++++---
matlab/dsge_likelihood.m | 2 +-
3 files changed, 80 insertions(+), 7 deletions(-)
diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m
index 6466ad2f1..d83aa07d7 100644
--- a/matlab/DsgeSmoother.m
+++ b/matlab/DsgeSmoother.m
@@ -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);
diff --git a/matlab/compute_Pinf_Pstar.m b/matlab/compute_Pinf_Pstar.m
index d847ffacc..4eb5fbc87 100644
--- a/matlab/compute_Pinf_Pstar.m
+++ b/matlab/compute_Pinf_Pstar.m
@@ -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
diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m
index 5be6aa2ff..16cc82117 100644
--- a/matlab/dsge_likelihood.m
+++ b/matlab/dsge_likelihood.m
@@ -358,7 +358,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;
--
GitLab