1function oo_ = realtime_shock_decomposition(M_,oo_,options_,varlist,bayestopt_,estim_params_) 2% function oo_ = realtime_shock_decomposition(M_,oo_,options_,varlist,bayestopt_,estim_params_) 3% Computes shocks contribution to a simulated trajectory. The fields set are 4% oo_.realtime_shock_decomposition, oo_.conditional_shock_decomposition and oo_.realtime_forecast_shock_decomposition. 5% Subfields are arrays n_var by nshock+2 by nperiods. The 6% first nshock columns store the respective shock contributions, column n+1 7% stores the role of the initial conditions, while column n+2 stores the 8% value of the smoothed variables. Both the variables and shocks are stored 9% in the order of declaration, i.e. M_.endo_names and M_.exo_names, respectively. 10% 11% INPUTS 12% M_: [structure] Definition of the model 13% oo_: [structure] Storage of results 14% options_: [structure] Options 15% varlist: [char] List of variables 16% bayestopt_: [structure] describing the priors 17% estim_params_: [structure] characterizing parameters to be estimated 18% 19% OUTPUTS 20% oo_: [structure] Storage of results 21% 22% SPECIAL REQUIREMENTS 23% none 24 25% Copyright (C) 2009-2020 Dynare Team 26% 27% This file is part of Dynare. 28% 29% Dynare is free software: you can redistribute it and/or modify 30% it under the terms of the GNU General Public License as published by 31% the Free Software Foundation, either version 3 of the License, or 32% (at your option) any later version. 33% 34% Dynare is distributed in the hope that it will be useful, 35% but WITHOUT ANY WARRANTY; without even the implied warranty of 36% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 37% GNU General Public License for more details. 38% 39% You should have received a copy of the GNU General Public License 40% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 41 42 43if isfield(oo_,'shock_decomposition_info') && isfield(oo_.shock_decomposition_info,'i_var') 44 if isfield (oo_,'realtime_conditional_shock_decomposition') ... 45 || isfield (oo_,'realtime_forecast_shock_decomposition') ... 46 || isfield (oo_,'realtime_shock_decomposition') ... 47 || isfield (oo_,'shock_decomposition') ... 48 || isfield (oo_,'conditional_shock_decomposition') ... 49 || isfield (oo_,'initval_decomposition') 50 error('realtime_shock_decomposition::squeezed shock decompositions are already stored in oo_') 51 end 52end 53 54with_epilogue = options_.shock_decomp.with_epilogue; 55 56% indices of endogenous variables 57if isempty(varlist) 58 varlist = M_.endo_names(1:M_.orig_endo_nbr); 59end 60 61[~, ~, index_uniques] = varlist_indices(varlist,M_.endo_names); 62varlist = varlist(index_uniques); 63 64% number of variables 65endo_nbr = M_.endo_nbr; 66 67% number of shocks 68nshocks = M_.exo_nbr; 69 70% parameter set 71parameter_set = options_.parameter_set; 72if isempty(parameter_set) 73 if isfield(oo_,'posterior_mean') 74 parameter_set = 'posterior_mean'; 75 elseif isfield(oo_,'mle_mode') 76 parameter_set = 'mle_mode'; 77 elseif isfield(oo_,'posterior') 78 parameter_set = 'posterior_mode'; 79 else 80 error(['realtime_shock_decomposition: option parameter_set is not specified ' ... 81 'and posterior mode is not available']) 82 end 83end 84 85presample = max(1,options_.presample-1); 86if isfield(options_.shock_decomp,'presample') 87 my_presample = max(1,options_.shock_decomp.presample); 88 presample = min(presample,my_presample); 89end 90% forecast_=0; 91forecast_ = options_.shock_decomp.forecast; 92forecast_params=0; 93if forecast_ && isfield(options_.shock_decomp,'forecast_params') 94 forecast_params = options_.shock_decomp.forecast_params; 95end 96fast_realtime = 0; 97if isfield(options_.shock_decomp,'fast_realtime') 98 fast_realtime = options_.shock_decomp.fast_realtime; 99end 100 101% save_realtime=0; 102save_realtime = options_.shock_decomp.save_realtime; 103% array of time points in the range options_.presample+1:options_.nobs 104 105zreal = zeros(endo_nbr+length(M_.epilogue_names)*with_epilogue,nshocks+2,options_.nobs+forecast_); 106zcond = zeros(endo_nbr+length(M_.epilogue_names)*with_epilogue,nshocks+2,options_.nobs); 107 108options_.selected_variables_only = 0; %make sure all variables are stored 109options_.plot_priors=0; 110init=1; 111nobs = options_.nobs; 112 113if forecast_ && any(forecast_params) 114 M1=M_; 115 M1.params = forecast_params; 116 [~,~,~,~,~,~,oo1] = dynare_resolve(M1,options_,oo_); 117end 118 119if fast_realtime 120 skipline() 121 skipline() 122 running_text = 'Fast realtime shock decomposition '; 123 newString=sprintf(running_text); 124 fprintf(['%s'],newString); 125 options_.nobs=fast_realtime; 126 [oo0,M_,~,~,Smoothed_Variables_deviation_from_mean0] = evaluate_smoother(parameter_set,varlist,M_,oo_,options_,bayestopt_,estim_params_); 127 gend0 = size(oo0.SmoothedShocks.(M_.exo_names{1}),1); 128 prctdone=0.5; 129 if isoctave 130 printf([running_text,' %3.f%% done\r'], prctdone*100); 131 else 132 s0=repmat('\b',1,length(newString)+1); 133 newString=sprintf([running_text,' %3.1f%% done'], prctdone*100); 134 fprintf([s0,'%s'],newString); 135 end 136 options_.nobs=nobs; 137 [oo2,M_,~,~,Smoothed_Variables_deviation_from_mean2] = evaluate_smoother(parameter_set,varlist,M_,oo_,options_,bayestopt_,estim_params_); 138 prctdone=1; 139 if isoctave 140 printf([running_text,' %3.f%% done\r'], prctdone*100); 141 else 142 s0=repmat('\b',1,length(newString)+1); 143 newString=sprintf([running_text,' %3.1f%% done'], prctdone*100); 144 fprintf([s0,'%s'],newString); 145 end 146end 147 148skipline() 149skipline() 150running_text = 'Realtime shock decomposition '; 151newString=sprintf(running_text); 152fprintf(['%s'],newString); 153 154for j=presample+1:nobs 155 % evalin('base',['options_.nobs=' int2str(j) ';']) 156 options_.nobs=j; 157 if ~fast_realtime 158 [oo,M_,~,~,Smoothed_Variables_deviation_from_mean] = evaluate_smoother(parameter_set,varlist,M_,oo_,options_,bayestopt_,estim_params_); 159 gend = size(oo.SmoothedShocks.(M_.exo_names{1}),1); 160 else 161 gend = gend0+j-fast_realtime; 162 if j>fast_realtime 163 oo=oo2; 164 Smoothed_Variables_deviation_from_mean = Smoothed_Variables_deviation_from_mean2(:,1:gend); 165 else 166 oo=oo0; 167 Smoothed_Variables_deviation_from_mean = Smoothed_Variables_deviation_from_mean0(:,1:gend); 168 end 169 end 170 % reduced form 171 dr = oo.dr; 172 173 % data reordering 174 order_var = dr.order_var; 175 inv_order_var = dr.inv_order_var; 176 177 178 % coefficients 179 A = dr.ghx; 180 B = dr.ghu; 181 182 if forecast_ 183 if any(forecast_params) 184 Af = oo1.dr.ghx; 185 Bf = oo1.dr.ghu; 186 else 187 Af = A; 188 Bf = B; 189 end 190 end 191 192 % initialization 193 epsilon=NaN(nshocks,gend); 194 for i = 1:nshocks 195 epsilon(i,:) = oo.SmoothedShocks.(M_.exo_names{i})(1:gend); 196 end 197 epsilon=[epsilon zeros(nshocks,forecast_)]; 198 199 z = zeros(endo_nbr,nshocks+2,gend+forecast_); 200 201 z(:,end,1:gend) = Smoothed_Variables_deviation_from_mean; 202 203 maximum_lag = M_.maximum_lag; 204 205 k2 = dr.kstate(find(dr.kstate(:,2) <= maximum_lag+1),[1 2]); 206 i_state = order_var(k2(:,1))+(min(i,maximum_lag)+1-k2(:,2))*M_.endo_nbr; 207 for i=1:gend+forecast_ 208 if i > 1 && i <= maximum_lag+1 209 lags = min(i-1,maximum_lag):-1:1; 210 end 211 212 if i > 1 213 tempx = permute(z(:,1:nshocks,lags),[1 3 2]); 214 m = min(i-1,maximum_lag); 215 tempx = [reshape(tempx,endo_nbr*m,nshocks); zeros(endo_nbr*(maximum_lag-i+1),nshocks)]; 216 if i > gend 217 z(:,nshocks+2,i) = Af(inv_order_var,:)*z(i_state,nshocks+2,lags); 218 % z(:,nshocks+2,i) = A(inv_order_var,:)*permute(z(i_state,nshocks+2,lags),[1 3 2]); 219 z(:,1:nshocks,i) = Af(inv_order_var,:)*tempx(i_state,:); 220 else 221 z(:,1:nshocks,i) = A(inv_order_var,:)*tempx(i_state,:); 222 end 223 lags = lags+1; 224 z(:,1:nshocks,i) = z(:,1:nshocks,i) + B(inv_order_var,:).*repmat(epsilon(:,i)',endo_nbr,1); 225 end 226 227 % z(:,1:nshocks,i) = z(:,1:nshocks,i) + B(inv_order_var,:).*repmat(epsilon(:,i)',endo_nbr,1); 228 z(:,nshocks+1,i) = z(:,nshocks+2,i) - sum(z(:,1:nshocks,i),2); 229 end 230 231 if with_epilogue 232 [z, epilogue_steady_state] = epilogue_shock_decomposition(z, M_, oo_); 233 if ~isfield(oo_,'shock_decomposition_info') || ~isfield(oo_.shock_decomposition_info,'epilogue_steady_state') 234 oo_.shock_decomposition_info.epilogue_steady_state = epilogue_steady_state; 235 end 236 end 237 %% conditional shock decomp 1 step ahead 238 z1 = zeros(endo_nbr,nshocks+2); 239 z1(:,end) = Smoothed_Variables_deviation_from_mean(:,gend); 240 for i=gend 241 242 z1(:,1:nshocks) = z1(:,1:nshocks) + B(inv_order_var,:).*repmat(epsilon(:,i)',endo_nbr,1); 243 z1(:,nshocks+1) = z1(:,nshocks+2) - sum(z1(:,1:nshocks),2); 244 end 245 if with_epilogue 246 clear ztmp0 247 ztmp0(:,1,:) = Smoothed_Variables_deviation_from_mean(:,1:gend-1); 248 ztmp0(:,2,:) = Smoothed_Variables_deviation_from_mean(:,1:gend-1); 249 ztmp = cat(3,cat(2,zeros(endo_nbr,nshocks,gend-1),ztmp0),z1); 250% ztmp = cat(3,zeros(endo_nbr,nshocks+2,40),ztmp); % pad with zeros in presample 251 z1 = epilogue_shock_decomposition(ztmp, M_, oo_); 252 z1=squeeze(z1(:,:,end)); 253 end 254 %% 255 256 %% conditional shock decomp k step ahead 257 if forecast_ && forecast_<j 258 zn = zeros(endo_nbr,nshocks+2,forecast_+1); 259 zn(:,end,1:forecast_+1) = Smoothed_Variables_deviation_from_mean(:,gend-forecast_:gend); 260 for i=1:forecast_+1 261 if i > 1 && i <= maximum_lag+1 262 lags = min(i-1,maximum_lag):-1:1; 263 end 264 265 if i > 1 266 tempx = permute(zn(:,1:nshocks,lags),[1 3 2]); 267 m = min(i-1,maximum_lag); 268 tempx = [reshape(tempx,endo_nbr*m,nshocks); zeros(endo_nbr*(maximum_lag-i+1-1),nshocks)]; 269 zn(:,1:nshocks,i) = A(inv_order_var,:)*tempx(i_state,:); 270 lags = lags+1; 271 zn(:,1:nshocks,i) = zn(:,1:nshocks,i) + B(inv_order_var,:).*repmat(epsilon(:,i+gend-forecast_-1)',endo_nbr,1); 272 end 273 274 % zn(:,1:nshocks,i) = zn(:,1:nshocks,i) + B(inv_order_var,:).*repmat(epsilon(:,i+gend-forecast_-1)',endo_nbr,1); 275 zn(:,nshocks+1,i) = zn(:,nshocks+2,i) - sum(zn(:,1:nshocks,i),2); 276 end 277 if with_epilogue 278 clear ztmp0 279 ztmp0(:,1,:) = Smoothed_Variables_deviation_from_mean(:,1:gend-forecast_-1); 280 ztmp0(:,2,:) = Smoothed_Variables_deviation_from_mean(:,1:gend-forecast_-1); 281 ztmp = cat(3,cat(2,zeros(endo_nbr,nshocks,gend-forecast_-1),ztmp0),zn); 282% ztmp = cat(3,zeros(endo_nbr,nshocks+2,40),ztmp); % pad with zeros (st state) in presample 283 zn = epilogue_shock_decomposition(ztmp, M_, oo_); 284 zn=squeeze(zn(:,:,end-forecast_:end)); 285 end 286 if ismember(j-forecast_,save_realtime) 287 oo_.conditional_shock_decomposition.(['time_' int2str(j-forecast_)])=zn; 288 end 289 end 290 %% 291 292 if init 293 zreal(:,:,1:j) = z(:,:,1:j); 294 elseif j<nobs 295 zreal(:,:,j) = z(:,:,gend); 296 else 297 zreal(:,:,j:end) = z(:,:,gend:end); 298 end 299 zcond(:,:,j) = z1; 300 if ismember(j,save_realtime) 301 oo_.realtime_shock_decomposition.(['time_' int2str(j)])=z; 302 end 303 304 if forecast_ 305 zfrcst(:,:,j+1) = z(:,:,gend+1); 306 ootmp.realtime_forecast_shock_decomposition.(['time_' int2str(j)])=z(:,:,gend:end); 307 if ismember(j,save_realtime) 308 oo_.realtime_forecast_shock_decomposition.(['time_' int2str(j)]) = ... 309 ootmp.realtime_forecast_shock_decomposition.(['time_' int2str(j)]); 310 end 311 if j>forecast_+presample 312 %% realtime conditional shock decomp k step ahead 313 ootmp.realtime_conditional_shock_decomposition.(['time_' int2str(j-forecast_)]) = ... 314 zreal(:,:,j-forecast_:j) - ... 315 ootmp.realtime_forecast_shock_decomposition.(['time_' int2str(j-forecast_)]); 316 ootmp.realtime_conditional_shock_decomposition.(['time_' int2str(j-forecast_)])(:,end-1,:) = ... 317 ootmp.realtime_conditional_shock_decomposition.(['time_' int2str(j-forecast_)])(:,end-1,:) + ... 318 ootmp.realtime_forecast_shock_decomposition.(['time_' int2str(j-forecast_)])(:,end,:); 319 ootmp.realtime_conditional_shock_decomposition.(['time_' int2str(j-forecast_)])(:,end,:) = ... 320 zreal(:,end,j-forecast_:j); 321 if ismember(j-forecast_,save_realtime) 322 oo_.realtime_conditional_shock_decomposition.(['time_' int2str(j-forecast_)]) = ... 323 ootmp.realtime_conditional_shock_decomposition.(['time_' int2str(j-forecast_)]); 324 end 325 if j==nobs 326 for my_forecast_=(forecast_-1):-1:1 327 ootmp.realtime_conditional_shock_decomposition.(['time_' int2str(j-my_forecast_)]) = ... 328 zreal(:,:,j-my_forecast_:j) - ... 329 ootmp.realtime_forecast_shock_decomposition.(['time_' int2str(j-my_forecast_)])(:,:,1:my_forecast_+1); 330 ootmp.realtime_conditional_shock_decomposition.(['time_' int2str(j-my_forecast_)])(:,end-1,:) = ... 331 ootmp.realtime_conditional_shock_decomposition.(['time_' int2str(j-my_forecast_)])(:,end-1,:) + ... 332 ootmp.realtime_forecast_shock_decomposition.(['time_' int2str(j-my_forecast_)])(:,end,1:my_forecast_+1); 333 ootmp.realtime_conditional_shock_decomposition.(['time_' int2str(j-my_forecast_)])(:,end,:) = ... 334 zreal(:,end,j-my_forecast_:j); 335 if ismember(j-my_forecast_,save_realtime) 336 oo_.realtime_conditional_shock_decomposition.(['time_' int2str(j-my_forecast_)]) = ... 337 ootmp.realtime_conditional_shock_decomposition.(['time_' int2str(j-my_forecast_)]); 338 end 339 end 340 end 341 342 end 343 end 344 345 prctdone=(j-presample)/(nobs-presample); 346 if isoctave 347 printf([running_text,' %3.f%% done\r'], prctdone*100); 348 else 349 s0=repmat('\b',1,length(newString)); 350 newString=sprintf([running_text,' %3.1f%% done'], prctdone*100); 351 fprintf([s0,'%s'],newString); 352 end 353 init=0; 354end 355oo_.realtime_shock_decomposition.pool = zreal; 356oo_.conditional_shock_decomposition.pool = zcond; 357if forecast_ 358 oo_.realtime_forecast_shock_decomposition.pool = zfrcst; 359end 360oo_.gui.ran_realtime_shock_decomposition = true; 361 362skipline() 363