File: compute_Pinf_Pstar.m

package info (click to toggle)
dynare 6.4-2
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 67,648 kB
  • sloc: cpp: 79,109; ansic: 28,917; objc: 12,430; yacc: 4,528; pascal: 1,993; lex: 1,441; sh: 1,129; python: 634; makefile: 626; lisp: 163; xml: 18
file content (187 lines) | stat: -rw-r--r-- 6,737 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
function [Pstar,Pinf] = compute_Pinf_Pstar(mf,T,R,Q,qz_criterium, restrict_columns)
% [Pstar,Pinf] = compute_Pinf_Pstar(mf,T,R,Q,qz_criterium, restrict_columns)
% Kitagawa transformation of state space system with a quasi-triangular
% transition matrix with unit roots at the top, but excluding zero columns of the transition matrix.
% Computation of Pstar and Pinf for Durbin and Koopman Diffuse filter
%
% The transformed state space is
%     y = [ss; z; x];
%     s = static variables (zero columns of T)
%     z = unit roots
%     x = stable roots
%     ss = s - z = stationarized static variables
%
% 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
%   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

% Copyright © 2006-2023 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.

np = size(T,1);
if iszero(Q)
    % this may happen if users set Q=0 and use heteroskedastic shocks to set
    % variances period by period
    % this in practice triggers a form of conditional filter where states
    % are initialized at st. state with zero variances
    Pstar=T*0;
    Pinf=T*0;
    return
end

if nargin == 6
    indx = restrict_columns;
    indx0=find(~ismember(1:np,indx));
else
    indx=(find(max(abs(T))>=1.e-10));
    indx0=(find(max(abs(T))<1.e-10));
end
np0=length(indx0);
T0=T(indx0,indx); % static variables vs. dynamic ones
R0=R(indx0,:);    % matrix of shocks for static variables

% 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);
k = find(abs(ordeig(ST)) > 2-qz_criterium);
nk = length(k);
nk1 = nk+1;
Pstar = zeros(np,np);
R1 = QT'*R;
B = R1*Q*R1';
i = np;
while i >= nk+2
    if ST(i,i-1) == 0
        if i == np
            c = zeros(np-nk,1);
        else
            c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+...
                ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i);
        end
        q = eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i);
        Pstar(nk1:i,i) = q\(B(nk1:i,i)+c);
        Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)';
        i = i - 1;
    else
        if i == np
            c = zeros(np-nk,1);
            c1 = zeros(np-nk,1);
        else
            c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+...
                ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i)+...
                ST(i,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1);
            c1 = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i-1,i+1:end)')+...
                 ST(i-1,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1)+...
                 ST(i-1,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i);
        end
        q = [eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i) -ST(nk1:i,nk1:i)*ST(i,i-1);...
             -ST(nk1:i,nk1:i)*ST(i-1,i) eye(i-nk)-ST(nk1:i,nk1:i)*ST(i-1,i-1)];
        z =  q\[B(nk1:i,i)+c;B(nk1:i,i-1)+c1];
        Pstar(nk1:i,i) = z(1:(i-nk));
        Pstar(nk1:i,i-1) = z(i-nk+1:end);
        Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)';
        Pstar(i-1,nk1:i-2) = Pstar(nk1:i-2,i-1)';
        i = i - 2;
    end
end
if i == nk+1
    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);
    Pstar(nk1,nk1)=(B(nk1,nk1)+c)/(1-ST(nk1,nk1)*ST(nk1,nk1));
end

if np0
    % 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);
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)];
    mf = ismember([indx0(:); indx(:)],mf);
    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

Pinf = QT*Pinf*QT';
Pstar = QT*Pstar*QT';