1function [Pstar,Pinf] = compute_Pinf_Pstar(mf,T,R,Q,qz_criterium, restrict_columns) 2% function [Z,ST,QT,R1,Pstar,Pinf] = schur_statespace_transformation(mf,T,R,Q,qz_criterium, restrict_columns) 3% Kitagawa transformation of state space system with a quasi-triangular 4% transition matrix with unit roots at the top, but excluding zero columns of the transition matrix. 5% Computation of Pstar and Pinf for Durbin and Koopman Diffuse filter 6% 7% The transformed state space is 8% y = [ss; z; x]; 9% s = static variables (zero columns of T) 10% z = unit roots 11% x = stable roots 12% ss = s - z = stationarized static variables 13% 14% INPUTS 15% mf [integer] vector of indices of observed variables in 16% state vector 17% T [double] matrix of transition 18% R [double] matrix of structural shock effects 19% Q [double] matrix of covariance of structural shocks 20% qz_criterium [double] numerical criterium for unit roots 21% 22% OUTPUTS 23% Pstar [double] matrix of covariance of stationary part 24% Pinf [double] matrix of covariance initialization for 25% nonstationary part 26% 27% ALGORITHM 28% Real Schur transformation of transition equation 29% 30% SPECIAL REQUIREMENTS 31% None 32 33% Copyright (C) 2006-2018 Dynare Team 34% 35% This file is part of Dynare. 36% 37% Dynare is free software: you can redistribute it and/or modify 38% it under the terms of the GNU General Public License as published by 39% the Free Software Foundation, either version 3 of the License, or 40% (at your option) any later version. 41% 42% Dynare is distributed in the hope that it will be useful, 43% but WITHOUT ANY WARRANTY; without even the implied warranty of 44% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 45% GNU General Public License for more details. 46% 47% You should have received a copy of the GNU General Public License 48% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 49 50np = size(T,1); 51if nargin == 6 52 indx = restrict_columns; 53 indx0=find(~ismember([1:np],indx)); 54else 55 indx=(find(max(abs(T))>=1.e-10)); 56 indx0=(find(max(abs(T))<1.e-10)); 57end 58np0=length(indx0); 59Tbkp = T; 60T0=T(indx0,indx); % static variables vs. dynamic ones 61R0=R(indx0,:); % matrix of shocks for static variables 62 63% Perform Kitagawa transformation only for non-zero columns of T 64T=T(indx,indx); 65R=R(indx,:); 66np = size(T,1); 67[QT,ST] = schur(T); 68e1 = abs(ordeig(ST)) > 2-qz_criterium; 69[QT,ST] = ordschur(QT,ST,e1); 70k = find(abs(ordeig(ST)) > 2-qz_criterium); 71nk = length(k); 72nk1 = nk+1; 73Pstar = zeros(np,np); 74R1 = QT'*R; 75B = R1*Q*R1'; 76i = np; 77while i >= nk+2 78 if ST(i,i-1) == 0 79 if i == np 80 c = zeros(np-nk,1); 81 else 82 c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+... 83 ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i); 84 end 85 q = eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i); 86 Pstar(nk1:i,i) = q\(B(nk1:i,i)+c); 87 Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)'; 88 i = i - 1; 89 else 90 if i == np 91 c = zeros(np-nk,1); 92 c1 = zeros(np-nk,1); 93 else 94 c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+... 95 ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i)+... 96 ST(i,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1); 97 c1 = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i-1,i+1:end)')+... 98 ST(i-1,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1)+... 99 ST(i-1,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i); 100 end 101 q = [eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i) -ST(nk1:i,nk1:i)*ST(i,i-1);... 102 -ST(nk1:i,nk1:i)*ST(i-1,i) eye(i-nk)-ST(nk1:i,nk1:i)*ST(i-1,i-1)]; 103 z = q\[B(nk1:i,i)+c;B(nk1:i,i-1)+c1]; 104 Pstar(nk1:i,i) = z(1:(i-nk)); 105 Pstar(nk1:i,i-1) = z(i-nk+1:end); 106 Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)'; 107 Pstar(i-1,nk1:i-2) = Pstar(nk1:i-2,i-1)'; 108 i = i - 2; 109 end 110end 111if i == nk+1 112 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); 113 Pstar(nk1,nk1)=(B(nk1,nk1)+c)/(1-ST(nk1,nk1)*ST(nk1,nk1)); 114end 115 116if np0 117 ST1=ST; 118 % Now I recover stationarized static variables using 119 % ss = s-A*z 120 % and 121 % z-z(-1) (growth rates of unit roots) only depends on stationary variables 122 Pstar = blkdiag(zeros(np0),Pstar); 123 ST = [zeros(length(Pstar),length(indx0)) [T0*QT ;ST]]; 124 R1 = [R0; R1]; 125 % Build the matrix for stationarized variables 126 STinf = ST(np0+1:np0+nk,np0+1:np0+nk); 127 iSTinf = inv(STinf); 128 ST0=ST; 129 ST0(:,1:np0+nk)=0; % stationarized static + 1st difference only respond to lagged stationary states 130 ST00 = ST(1:np0,np0+1:np0+nk); 131 % A\B is the matrix division of A into B, which is roughly the 132 % same as INV(A)*B 133 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 134 R10 = R1; 135 R10(1:np0,:) = R1(1:np0,:)-ST00*(iSTinf*R1(np0+1:np0+nk,:)); % snip non-stationary part 136 % Kill non-stationary part before projecting Pstar 137 ST0(np0+1:np0+nk,:)=0; 138 R10(np0+1:np0+nk,:)=0; % is this questionable???? IT HAS TO in order to match Michel's version!!! 139 % project Pstar onto static x 140 Pstar = ST0*Pstar*ST0'+R10*Q*R10'; 141 % QT(1:np0,np0+1:np0+nk) = QT(1:np0,np0+1:np0+nk)+ST(1:np0,np0+1:np0+nk); %%% is this questionable ???? 142 % reorder QT entries 143else 144 STinf = ST(np0+1:np0+nk,np0+1:np0+nk); 145end 146 147% stochastic trends with no influence on observed variables are 148% arbitrarily initialized to zero 149Pinf = zeros(np,np); 150Pinf(1:nk,1:nk) = eye(nk); 151if np0 152 STtriu = STinf-eye(nk); 153 % A\B is the matrix division of A into B, which is roughly the 154 % same as INV(A)*B 155 STinf0 = ST00*(eye(nk)-iSTinf*STtriu); 156 Pinf = blkdiag(zeros(np0),Pinf); 157 QT = blkdiag(eye(np0),QT); 158 QTinf = QT; 159 QTinf(1:np0,np0+1:np0+nk) = STinf0; 160 QTinf([indx0(:); indx(:)],:) = QTinf; 161 STinf1 = [zeros(np0+np,np0) [STinf0; eye(nk); zeros(np-nk,nk)] zeros(np0+np,np-nk)]; 162 for k = 1:nk 163 if norm(QTinf(mf,:)*ST([indx0(:); indx(:)],k+np0)) < 1e-8 164 Pinf(k+np0,k+np0) = 0; 165 end 166 end 167 Pinf = STinf1*Pinf*STinf1'; 168 QT([indx0(:); indx(:)],:) = QT; 169else 170 for k = 1:nk 171 if norm(QT(mf,:)*ST(:,k)) < 1e-8 172 Pinf(k+np0,k+np0) = 0; 173 end 174 end 175end 176 177Pinf = QT*Pinf*QT'; 178Pstar = QT*Pstar*QT'; 179