schur_statespace_transformation

 function [Z,ST,QT,R1,Pstar,Pinf] = schur_statespace(mf,T,R,Q,qz_criterium)
 Kitagawa transformation of state space system with a quasi-triangular
 transition matrix with unit roots at the top. Computation of Pstar and
 Pinf for Durbin and Koopman Diffuse filter
 
 INPUTS 
   mf           [integer]    vector of indices of observed variables in
                             state vector
   T            [double]     matrix of transition
   R            [double]     matrix of structural shock effects
   Q            [double]     matrix of covariance of structural shocks
   qz_criterium [double]     numerical criterium for unit roots   
  
 OUTPUTS 
   Z            [double]     transformed matrix of measurement equation
   ST           [double]     tranformed matrix of transition
   R1           [double]     tranformed matrix of structural shock effects
   QT           [double]     matrix of Schur vectors
   Pstar        [double]     matrix of covariance of stationary part
   Pinf         [double]     matrix of covariance initialization for
                             nonstationary part    
    
 ALGORITHM 
   Real Schur transformation of transition equation

 SPECIAL REQUIREMENTS
   None

0001 function [Z,ST,R1,QT,Pstar,Pinf] = schur_statespace_transformation(mf,T,R,Q,qz_criterium)
0002 % function [Z,ST,QT,R1,Pstar,Pinf] = schur_statespace(mf,T,R,Q,qz_criterium)
0003 % Kitagawa transformation of state space system with a quasi-triangular
0004 % transition matrix with unit roots at the top. Computation of Pstar and
0005 % Pinf for Durbin and Koopman Diffuse filter
0006 %
0007 % INPUTS
0008 %   mf           [integer]    vector of indices of observed variables in
0009 %                             state vector
0010 %   T            [double]     matrix of transition
0011 %   R            [double]     matrix of structural shock effects
0012 %   Q            [double]     matrix of covariance of structural shocks
0013 %   qz_criterium [double]     numerical criterium for unit roots
0014 %
0015 % OUTPUTS
0016 %   Z            [double]     transformed matrix of measurement equation
0017 %   ST           [double]     tranformed matrix of transition
0018 %   R1           [double]     tranformed matrix of structural shock effects
0019 %   QT           [double]     matrix of Schur vectors
0020 %   Pstar        [double]     matrix of covariance of stationary part
0021 %   Pinf         [double]     matrix of covariance initialization for
0022 %                             nonstationary part
0023 %
0024 % ALGORITHM
0025 %   Real Schur transformation of transition equation
0026 %
0027 % SPECIAL REQUIREMENTS
0028 %   None
0029 
0030 % Copyright (C) 2006-2011 Dynare Team
0031 %
0032 % This file is part of Dynare.
0033 %
0034 % Dynare is free software: you can redistribute it and/or modify
0035 % it under the terms of the GNU General Public License as published by
0036 % the Free Software Foundation, either version 3 of the License, or
0037 % (at your option) any later version.
0038 %
0039 % Dynare is distributed in the hope that it will be useful,
0040 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0041 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0042 % GNU General Public License for more details.
0043 %
0044 % You should have received a copy of the GNU General Public License
0045 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0046 
0047 np = size(T,1);
0048 [QT,ST] = schur(T);
0049 e1 = abs(ordeig(ST)) > 2-qz_criterium;
0050 [QT,ST] = ordschur(QT,ST,e1);
0051 k = find(abs(ordeig(ST)) > 2-qz_criterium);
0052 nk = length(k);
0053 nk1 = nk+1;
0054 Pstar = zeros(np,np);
0055 B = QT'*R*Q*R'*QT;
0056 i = np;
0057 while i >= nk+2
0058     if ST(i,i-1) == 0
0059         if i == np
0060             c = zeros(np-nk,1);
0061         else
0062             c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+...
0063                 ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i);
0064         end
0065         q = eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i);
0066         Pstar(nk1:i,i) = q\(B(nk1:i,i)+c);
0067         Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)';
0068         i = i - 1;
0069     else
0070         if i == np
0071             c = zeros(np-nk,1);
0072             c1 = zeros(np-nk,1);
0073         else
0074             c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+...
0075                 ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i)+...
0076                 ST(i,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1);
0077             c1 = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i-1,i+1:end)')+...
0078                  ST(i-1,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1)+...
0079                  ST(i-1,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i);
0080         end
0081         q = [eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i) -ST(nk1:i,nk1:i)*ST(i,i-1);...
0082              -ST(nk1:i,nk1:i)*ST(i-1,i) eye(i-nk)-ST(nk1:i,nk1:i)*ST(i-1,i-1)];
0083         z =  q\[B(nk1:i,i)+c;B(nk1:i,i-1)+c1];
0084         Pstar(nk1:i,i) = z(1:(i-nk));
0085         Pstar(nk1:i,i-1) = z(i-nk+1:end);
0086         Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)';
0087         Pstar(i-1,nk1:i-2) = Pstar(nk1:i-2,i-1)';
0088         i = i - 2;
0089     end
0090 end
0091 if i == nk+1
0092     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);
0093     Pstar(nk1,nk1)=(B(nk1,nk1)+c)/(1-ST(nk1,nk1)*ST(nk1,nk1));
0094 end
0095 
0096 Z = QT(mf,:);
0097 R1 = QT'*R;
0098 
0099 % stochastic trends with no influence on observed variables are
0100 % arbitrarily initialized to zero
0101 Pinf = zeros(np,np);
0102 Pinf(1:nk,1:nk) = eye(nk);
0103 [QQ,RR,EE] = qr(Z*ST(:,1:nk),0);
0104 k = find(abs(diag([RR; zeros(nk-size(Z,1),size(RR,2))])) < 1e-8);
0105 if length(k) > 0
0106     k1 = EE(:,k);
0107     dd =ones(nk,1);
0108     dd(k1) = zeros(length(k1),1);
0109     Pinf(1:nk,1:nk) = diag(dd);
0110 end