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