1function [CheckCO,minns,minSYS] = get_minimal_state_representation(SYS, derivs_flag) 2% Derives and checks the minimal state representation 3% Let x = A*x(-1) + B*u and y = C*x(-1) + D*u be a linear state space 4% system, then this function computes the following representation 5% xmin = minA*xmin(-1) + minB*u and and y=minC*xmin(-1) + minD*u 6% 7% ------------------------------------------------------------------------- 8% INPUTS 9% SYS [structure] 10% with the following necessary fields: 11% A: [nspred by nspred] in DR order 12% Transition matrix for all state variables 13% B: [nspred by exo_nbr] in DR order 14% Transition matrix mapping shocks today to states today 15% C: [varobs_nbr by nspred] in DR order 16% Measurement matrix linking control/observable variables to states 17% D: [varobs_nbr by exo_nbr] in DR order 18% Measurement matrix mapping shocks today to controls/observables today 19% and optional fields: 20% dA: [nspred by nspred by totparam_nbr] in DR order 21% Jacobian (wrt to all parameters) of transition matrix A 22% dB: [nspred by exo_nbr by totparam_nbr] in DR order 23% Jacobian (wrt to all parameters) of transition matrix B 24% dC: [varobs_nbr by nspred by totparam_nbr] in DR order 25% Jacobian (wrt to all parameters) of measurement matrix C 26% dD: [varobs_nbr by exo_nbr by totparam_nbr] in DR order 27% Jacobian (wrt to all parameters) of measurement matrix D 28% derivs_flag [scalar] 29% (optional) indicator whether to output parameter derivatives 30% ------------------------------------------------------------------------- 31% OUTPUTS 32% CheckCO: [scalar] 33% equals to 1 if minimal state representation is found 34% minns: [scalar] 35% length of minimal state vector 36% SYS [structure] 37% with the following fields: 38% minA: [minns by minns] in DR-order 39% transition matrix A for evolution of minimal state vector 40% minB: [minns by exo_nbr] in DR-order 41% transition matrix B for evolution of minimal state vector 42% minC: [varobs_nbr by minns] in DR-order 43% measurement matrix C for evolution of controls, depending on minimal state vector only 44% minD: [varobs_nbr by minns] in DR-order 45% measurement matrix D for evolution of controls, depending on minimal state vector only 46% dminA: [minns by minns by totparam_nbr] in DR order 47% Jacobian (wrt to all parameters) of transition matrix minA 48% dminB: [minns by exo_nbr by totparam_nbr] in DR order 49% Jacobian (wrt to all parameters) of transition matrix minB 50% dminC: [varobs_nbr by minns by totparam_nbr] in DR order 51% Jacobian (wrt to all parameters) of measurement matrix minC 52% dminD: [varobs_nbr by u_nbr by totparam_nbr] in DR order 53% Jacobian (wrt to all parameters) of measurement matrix minD 54% ------------------------------------------------------------------------- 55% This function is called by 56% * get_identification_jacobians.m (previously getJJ.m) 57% ------------------------------------------------------------------------- 58% This function calls 59% * check_minimality (embedded) 60% * minrealold (embedded) 61% ========================================================================= 62% Copyright (C) 2019-2020 Dynare Team 63% 64% This file is part of Dynare. 65% 66% Dynare is free software: you can redistribute it and/or modify 67% it under the terms of the GNU General Public License as published by 68% the Free Software Foundation, either version 3 of the License, or 69% (at your option) any later version. 70% 71% Dynare is distributed in the hope that it will be useful, 72% but WITHOUT ANY WARRANTY; without even the implied warranty of 73% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 74% GNU General Public License for more details. 75% 76% You should have received a copy of the GNU General Public License 77% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 78% ========================================================================= 79if nargin == 1 80 derivs_flag = 0; 81end 82realsmall = 1e-7; 83[nspred,exo_nbr] = size(SYS.B); 84varobs_nbr = size(SYS.C,1); 85 86% Check controllability and observability conditions for full state vector 87CheckCO = check_minimality(SYS.A,SYS.B,SYS.C); 88if CheckCO == 1 % If model is already minimal, we are finished 89 minns = nspred; 90 minSYS = SYS; 91else 92 %Model is not minimal 93 try 94 minreal_flag = 1; 95 % In future we will use SLICOT TB01PD.f mex file [to do @wmutschl], currently use workaround 96 [mSYS,U] = minrealold(SYS,realsmall); 97 minns = size(mSYS.A,1); 98 CheckCO = check_minimality(mSYS.A,mSYS.B,mSYS.C); 99 if CheckCO 100 minSYS.A = mSYS.A; 101 minSYS.B = mSYS.B; 102 minSYS.C = mSYS.C; 103 minSYS.D = mSYS.D; 104 if derivs_flag 105 totparam_nbr = size(SYS.dA,3); 106 minSYS.dA = zeros(minns,minns,totparam_nbr); 107 minSYS.dB = zeros(minns,exo_nbr,totparam_nbr); 108 minSYS.dC = zeros(varobs_nbr,minns,totparam_nbr); 109 % Note that orthogonal matrix U is such that (U*dA*U',U*dB,dC*U') is a Kalman decomposition of (dA,dB,dC) % 110 for jp=1:totparam_nbr 111 dA_tmp = U*SYS.dA(:,:,jp)*U'; 112 dB_tmp = U*SYS.dB(:,:,jp); 113 dC_tmp = SYS.dC(:,:,jp)*U'; 114 minSYS.dA(:,:,jp) = dA_tmp(1:minns,1:minns); 115 minSYS.dB(:,:,jp) = dB_tmp(1:minns,:); 116 minSYS.dC(:,:,jp) = dC_tmp(:,1:minns); 117 end 118 minSYS.dD = SYS.dD; 119 end 120 end 121 catch 122 minreal_flag = 0; % if something went wrong use below procedure 123 end 124 125 if minreal_flag == 0 126 fprintf('Use naive brute-force approach to find minimal state space system.\n These computations may be inaccurate/wrong as ''minreal'' is not available, please raise an issue on GitLab or the forum\n') 127 % create indices for unnecessary states 128 exogstateindex = find(abs(sum(SYS.A,1))>realsmall); 129 minns = length(exogstateindex); 130 % remove unnecessary states from solution matrices 131 A_2 = SYS.A(exogstateindex,exogstateindex); 132 B_2 = SYS.B(exogstateindex,:); 133 C_2 = SYS.C(:,exogstateindex); 134 D = SYS.D; 135 ind_A2 = exogstateindex; 136 % minimal representation 137 minSYS.A = A_2; 138 minSYS.B = B_2; 139 minSYS.C = C_2; 140 minSYS.D = D; 141 142 % Check controllability and observability conditions 143 CheckCO = check_minimality(minSYS.A,minSYS.B,minSYS.C); 144 145 if CheckCO ~=1 146 % do brute-force search 147 j=1; 148 while (CheckCO==0 && j<minns) 149 combis = nchoosek(1:minns,j); 150 i=1; 151 while i <= size(combis,1) 152 ind_A2 = exogstateindex; 153 ind_A2(combis(j,:)) = []; 154 % remove unnecessary states from solution matrices 155 A_2 = SYS.A(ind_A2,ind_A2); 156 B_2 = SYS.B(ind_A2,:); 157 C_2 = SYS.C(:,ind_A2); 158 D = SYS.D; 159 % minimal representation 160 minSYS.A = A_2; 161 minSYS.B = B_2; 162 minSYS.C = C_2; 163 minSYS.D = D; 164 % Check controllability and observability conditions 165 CheckCO = check_minimality(minSYS.A,minSYS.B,minSYS.C); 166 if CheckCO == 1 167 minns = length(ind_A2); 168 break; 169 end 170 i=i+1; 171 end 172 j=j+1; 173 end 174 end 175 if derivs_flag 176 minSYS.dA = SYS.dA(ind_A2,ind_A2,:); 177 minSYS.dB = SYS.dB(ind_A2,:,:); 178 minSYS.dC = SYS.dC(:,ind_A2,:); 179 minSYS.dD = SYS.dD; 180 end 181 end 182end 183 184function CheckCO = check_minimality(a,b,c) 185%% This function computes the controllability and the observability matrices of the ABCD system and checks if the system is minimal 186% 187% Inputs: Solution matrices A,B,C of ABCD representation of a DSGE model 188% Outputs: CheckCO: equals 1, if both rank conditions for observability and controllability are fullfilled 189n = size(a,1); 190CC = b; % Initialize controllability matrix 191OO = c; % Initialize observability matrix 192if n >= 2 193 for indn = 1:1:n-1 194 CC = [CC, (a^indn)*b]; % Set up controllability matrix 195 OO = [OO; c*(a^indn)]; % Set up observability matrix 196 end 197end 198CheckC = (rank(full(CC))==n); % Check rank of controllability matrix 199CheckO = (rank(full(OO))==n); % Check rank of observability matrix 200CheckCO = CheckC&CheckO; % equals 1 if minimal state 201end % check_minimality end 202 203function [mSYS,U] = minrealold(SYS,tol) 204 % This is a temporary replacement for minreal, will be replaced by a mex file from SLICOT TB01PD.f soon 205 a = SYS.A; 206 b = SYS.B; 207 c = SYS.C; 208 [ns,nu] = size(b); 209 [am,bm,cm,U,k] = ControllabilityStaircaseRosenbrock(a,b,c,tol); 210 kk = sum(k); 211 nu = ns - kk; 212 nn = nu; 213 am = am(nu+1:ns,nu+1:ns); 214 bm = bm(nu+1:ns,:); 215 cm = cm(:,nu+1:ns); 216 ns = ns - nu; 217 if ns 218 [am,bm,cm,t,k] = ObservabilityStaircaseRosenbrock(am,bm,cm,tol); 219 kk = sum(k); 220 nu = ns - kk; 221 nn = nn + nu; 222 am = am(nu+1:ns,nu+1:ns); 223 bm = bm(nu+1:ns,:); 224 cm = cm(:,nu+1:ns); 225 end 226 mSYS.A = am; 227 mSYS.B = bm; 228 mSYS.C = cm; 229 mSYS.D = SYS.D; 230end 231 232 233function [abar,bbar,cbar,t,k] = ObservabilityStaircaseRosenbrock(a,b,c,tol) 234 %Observability staircase form 235 [aa,bb,cc,t,k] = ControllabilityStaircaseRosenbrock(a',c',b',tol); 236 abar = aa'; bbar = cc'; cbar = bb'; 237end 238 239function [abar,bbar,cbar,t,k] = ControllabilityStaircaseRosenbrock(a, b, c, tol) 240 % Controllability staircase algorithm of Rosenbrock, 1968 241 [ra,ca] = size(a); 242 [rb,cb] = size(b); 243 ptjn1 = eye(ra); 244 ajn1 = a; 245 bjn1 = b; 246 rojn1 = cb; 247 deltajn1 = 0; 248 sigmajn1 = ra; 249 k = zeros(1,ra); 250 if nargin == 3 251 tol = ra*norm(a,1)*eps; 252 end 253 for jj = 1 : ra 254 [uj,sj,vj] = svd(bjn1); 255 [rsj,csj] = size(sj); 256 %p = flip(eye(rsj),2); 257 p = eye(rsj); 258 p = p(:,end:-1:1); 259 p = permute(p,[2 1 3:ndims(eye(rsj))]); 260 uj = uj*p; 261 bb = uj'*bjn1; 262 roj = rank(bb,tol); 263 [rbb,cbb] = size(bb); 264 sigmaj = rbb - roj; 265 sigmajn1 = sigmaj; 266 k(jj) = roj; 267 if roj == 0, break, end 268 if sigmaj == 0, break, end 269 abxy = uj' * ajn1 * uj; 270 aj = abxy(1:sigmaj,1:sigmaj); 271 bj = abxy(1:sigmaj,sigmaj+1:sigmaj+roj); 272 ajn1 = aj; 273 bjn1 = bj; 274 [ruj,cuj] = size(uj); 275 ptj = ptjn1 * ... 276 [uj zeros(ruj,deltajn1); ... 277 zeros(deltajn1,cuj) eye(deltajn1)]; 278 ptjn1 = ptj; 279 deltaj = deltajn1 + roj; 280 deltajn1 = deltaj; 281 end 282 t = ptjn1'; 283 abar = t * a * t'; 284 bbar = t * b; 285 cbar = c * t'; 286end 287 288end % Main function end 289