1function [out, steady_state] = plot_shock_decomposition(M_,oo_,options_,varlist) 2% function plot_shock_decomposition(M_,oo_,options_,varlist) 3% Plots the results of shock_decomposition 4% 5% INPUTS 6% M_: [structure] Definition of the model 7% oo_: [structure] Storage of results 8% options_: [structure] Options 9% varlist: [char] List of variables 10% 11% SPECIAL REQUIREMENTS 12% none 13 14% Copyright (C) 2016-2019 Dynare Team 15% 16% This file is part of Dynare. 17% 18% Dynare is free software: you can redistribute it and/or modify 19% it under the terms of the GNU General Public License as published by 20% the Free Software Foundation, either version 3 of the License, or 21% (at your option) any later version. 22% 23% Dynare is distributed in the hope that it will be useful, 24% but WITHOUT ANY WARRANTY; without even the implied warranty of 25% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 26% GNU General Public License for more details. 27% 28% You should have received a copy of the GNU General Public License 29% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 30 31options_.nodisplay = options_.plot_shock_decomp.nodisplay; 32options_.graph_format = options_.plot_shock_decomp.graph_format; 33 34if ~isfield(oo_,'shock_decomposition_info') 35 oo_.shock_decomposition_info = struct(); 36end 37if ~isfield(oo_,'plot_shock_decomposition_info') 38 oo_.plot_shock_decomposition_info = struct(); 39end 40 41out=oo_; 42% indices of endogenous variables 43exist_varlist = 1; 44if size(varlist,1) == 0 45 exist_varlist = 0; 46 if size( M_.endo_names,1) >= M_.orig_endo_nbr 47 varlist = M_.endo_names(1:M_.orig_endo_nbr); 48 else 49 varlist = M_.endo_names; 50 end 51end 52 53if isfield(options_.plot_shock_decomp,'init2shocks') % private trap for uimenu calls 54 init2shocks=options_.plot_shock_decomp.init2shocks; 55else 56 init2shocks=[]; 57end 58if ~isempty(init2shocks) 59 init2shocks=M_.init2shocks.(init2shocks); 60end 61 62 63epilogue_decomp=false; 64if exist_varlist && any(ismember(varlist,M_.epilogue_names)) 65 epilogue_decomp=true; 66 M_.endo_names = [M_.endo_names;M_.epilogue_names]; 67 M_.endo_names_tex = [M_.endo_names_tex;M_.epilogue_names]; 68 M_.endo_nbr = length( M_.endo_names ); 69end 70if isfield(oo_.shock_decomposition_info,'i_var') && (M_.endo_nbr>=M_.orig_endo_nbr) 71 if max(oo_.shock_decomposition_info.i_var)>M_.orig_endo_nbr 72 epilogue_decomp=true; 73 M_.endo_names = [M_.endo_names;M_.epilogue_names]; 74 M_.endo_names_tex = [M_.endo_names_tex;M_.epilogue_names]; 75 M_.endo_nbr = length( M_.endo_names ); 76 end 77 M_.endo_names = M_.endo_names(oo_.shock_decomposition_info.i_var,:); 78 M_.endo_names_tex = M_.endo_names_tex(oo_.shock_decomposition_info.i_var,:); 79 M_.endo_nbr = length( oo_.shock_decomposition_info.i_var ); 80end 81 82try 83 [i_var,nvar,index_uniques] = varlist_indices(varlist,M_.endo_names); 84catch ME 85 if isfield(oo_.shock_decomposition_info,'i_var') 86 warning('shock decomp results for some input variable was not stored: I recompute all decompositions') 87 M_ = evalin('base','M_'); 88 bayestopt_ = evalin('base','bayestopt_'); 89 estim_params_ = evalin('base','estim_params_'); 90 options_.no_graph.shock_decomposition=1; % force nograph in computing decompositions! 91 oo_.shock_decomposition_info = rmfield(oo_.shock_decomposition_info,'i_var'); 92 var_list_ = char(); 93 disp('recomputing shock decomposition ...') 94 [oo_,M_]= shock_decomposition(M_,oo_,options_,var_list_,bayestopt_,estim_params_); 95 if isfield(oo_,'realtime_shock_decomposition') || options_.plot_shock_decomp.realtime 96 disp('recomputing realtime shock decomposition ...') 97 oo_ = realtime_shock_decomposition(M_,oo_,options_,var_list_,bayestopt_,estim_params_); 98 end 99 if isfield(oo_,'initval_decomposition') 100 disp('recomputing initval shock decomposition ...') 101 oo_ = initial_condition_decomposition(M_,oo_,options_,0,bayestopt_,estim_params_); 102 end 103 [i_var,nvar,index_uniques] = varlist_indices(varlist,M_.endo_names); 104 out = oo_; 105 else 106 rethrow(ME) 107 end 108end 109 110varlist = varlist(index_uniques); 111 112if ~isfield(out.shock_decomposition_info,'i_var') && exist_varlist 113 if ~isfield(out.plot_shock_decomposition_info,'i_var') 114 out.plot_shock_decomposition_info.i_var = i_var; 115 else 116 out.plot_shock_decomposition_info.i_var = unique([i_var(:); out.plot_shock_decomposition_info.i_var(:)]); 117 end 118end 119 120type=options_.plot_shock_decomp.type; 121if isequal(type, 'aoa') && isfield(options_.plot_shock_decomp,'q2a') && isstruct(options_.plot_shock_decomp.q2a) 122 q2avec=options_.plot_shock_decomp.q2a; 123 if nvar>1 124 for jv = 1:nvar 125 my_varlist = varlist(jv); 126 indv = strcmp(my_varlist,{q2avec.qname}); 127 options_.plot_shock_decomp.q2a = q2avec(indv); 128 plot_shock_decomposition(M_,oo_,options_,my_varlist); 129 end 130 return 131 else 132 indv = strcmp(varlist,{q2avec.qname}); 133 options_.plot_shock_decomp.q2a = q2avec(indv); 134 end 135end 136 137% number of shocks 138nshocks = M_.exo_nbr; 139fig_name=''; 140 141if isfield(options_.plot_shock_decomp,'diff') % private trap for uimenu calls 142 differentiate_decomp=options_.plot_shock_decomp.diff; 143else 144 differentiate_decomp=0; 145end 146if isfield(options_.plot_shock_decomp,'flip') % private trap for uimenu calls 147 flip_decomp=options_.plot_shock_decomp.flip; 148else 149 flip_decomp=0; 150end 151if isfield(options_.plot_shock_decomp,'expand') % private trap for uimenu calls 152 expand=options_.plot_shock_decomp.expand; 153else 154 expand=0; 155 options_.plot_shock_decomp.expand=0; 156end 157if ~isfield(options_.plot_shock_decomp,'init_cond_decomp') 158 options_.plot_shock_decomp.init_cond_decomp=0; 159end 160 161options_.plot_shock_decomp.initval=0; 162if ~isempty(options_.plot_shock_decomp.fig_name) 163 fig_name=[' ' options_.plot_shock_decomp.fig_name]; 164 if length(fig_name)>=8 && strcmp(fig_name(end-6:end),'initval') 165 options_.plot_shock_decomp.initval=1; 166 end 167end 168 169detail_plot=options_.plot_shock_decomp.detail_plot; 170realtime_= options_.plot_shock_decomp.realtime; 171vintage_ = options_.plot_shock_decomp.vintage; 172forecast_ = options_.shock_decomp.forecast; 173steadystate = options_.plot_shock_decomp.steadystate; 174write_xls = options_.plot_shock_decomp.write_xls; 175 176if vintage_ 177 forecast_ = min(forecast_,options_.nobs-vintage_); 178end 179 180initial_date = options_.initial_date; 181if isempty(initial_date) 182 if isempty(type) 183 % we assume annual model 184 initial_date = dates('1Y'); 185 else 186 % we assume the sample starts in Q1 187 initial_date = dates('1Q1'); 188 end 189end 190 191if isfield(options_.plot_shock_decomp,'q2a') % private trap for aoa calls 192 q2a=options_.plot_shock_decomp.q2a; 193 if isstruct(q2a) && isempty(fieldnames(q2a)) 194 q2a=0; 195 end 196else 197 q2a=0; 198end 199 200switch realtime_ 201 202 case 0 203 if ~expand 204 z = oo_.shock_decomposition; 205 end 206 fig_name1=fig_name; 207 208 case 1 % realtime 209 if vintage_ 210 if ~expand 211 z = oo_.realtime_shock_decomposition.(['time_' int2str(vintage_)]); 212 end 213 fig_name1=[fig_name ' realtime (vintage ' char(initial_date+vintage_-1) ')']; 214 else 215 if ~expand 216 z = oo_.realtime_shock_decomposition.pool; 217 end 218 fig_name1=[fig_name ' realtime (rolling)']; 219 end 220 221 case 2 % conditional 222 if vintage_ 223 if ~expand 224 z = oo_.realtime_conditional_shock_decomposition.(['time_' int2str(vintage_)]); 225 end 226 initial_date = initial_date+vintage_-1; 227 fig_name1=[fig_name ' ' int2str(forecast_) '-step ahead conditional forecast (given ' char(initial_date) ')']; 228 else 229 if ~expand 230 z = oo_.conditional_shock_decomposition.pool; 231 end 232 fig_name1=[fig_name ' 1-step ahead conditional forecast (rolling)']; 233 end 234 235 case 3 % forecast 236 if vintage_ 237 if ~expand 238 z = oo_.realtime_forecast_shock_decomposition.(['time_' int2str(vintage_)]); 239 end 240 initial_date = initial_date+vintage_-1; 241 fig_name1=[fig_name ' ' int2str(forecast_) '-step ahead forecast (given ' char(initial_date) ')']; 242 else 243 if ~expand 244 z = oo_.realtime_forecast_shock_decomposition.pool; 245 end 246 fig_name1=[fig_name ' 1-step ahead forecast (rolling)']; 247 end 248end 249 250 251if ~isempty(init2shocks) && ~expand 252 n=size(init2shocks,1); 253 M_.exo_names_init=M_.exo_names; 254 for i=1:n 255 j=strmatch(init2shocks{i}{1},M_.endo_names,'exact'); 256 if ~isempty(init2shocks{i}{2}) 257 jj=strmatch(init2shocks{i}{2},M_.exo_names,'exact'); 258 M_.exo_names_init{jj}=[M_.exo_names_init{jj} ' + ' M_.endo_names{j}]; 259 z(:,jj,:)= z(:,jj,:) + oo_.initval_decomposition (:,j,:); 260 else 261 z(:,end,:)= z(:,end,:) - oo_.initval_decomposition (:,j,:); 262 end 263 z(:,end-1,:)= z(:,end-1,:) - oo_.initval_decomposition (:,j,:); 264 265 end 266end 267 268if ~epilogue_decomp 269 if isfield(oo_.dr,'ys') 270 steady_state = oo_.dr.ys; 271 else 272 steady_state = oo_.steady_state; 273 end 274else 275 steady_state = oo_.shock_decomposition_info.epilogue_steady_state; 276end 277 278if isequal(type,'aoa') && isstruct(q2a) 279 if expand 280 za = options_.plot_shock_decomp.zfull; 281 genda = size(za,3); 282 endo_nbra = size(za,1); 283 [za, shock_names, M_] = make_the_groups(za,genda,endo_nbra,nshocks,M_, options_); 284 end 285 if realtime_ 286 % take all dates where realtime is saved 287 qqq=options_.initial_date+options_.shock_decomp.save_realtime(:)-1; 288 % take the first Q4 of saved realtime 289 t0=min(options_.shock_decomp.save_realtime(qqq.time(:,2)==4)); 290 if isempty(t0) 291 error('the realtime decompositions are not stored in Q4! Please check your dates and settings.') 292 end 293 if ~isfield(q2a,'type') % private trap for aoa calls 294 q2a.type=1; 295 end 296 if ~isfield(q2a,'islog') % private trap for aoa calls 297 q2a.islog=0; 298 end 299 if ~isfield(q2a,'GYTREND0') % private trap for aoa calls 300 q2a.GYTREND0=0; 301 end 302 if ~isfield(q2a,'aux') % private trap for aoa calls 303 q2a.aux=0; 304 end 305 if ~isfield(q2a,'cumfix') % private trap for aoa calls 306 q2a.cumfix=1; 307 end 308 if ~isfield(q2a,'plot') % private trap for aoa calls 309 q2a.plot=1; % growth rate 310 end 311 312 if ~expand 313 if options_.plot_shock_decomp.interactive && ~isempty(options_.plot_shock_decomp.use_shock_groups) 314 mygroup = options_.plot_shock_decomp.use_shock_groups; 315 options_.plot_shock_decomp.use_shock_groups=''; 316 zafull = ... 317 annualized_shock_decomposition(oo_,M_, options_, i_var, t0, options_.nobs, realtime_, vintage_, steady_state,q2a); 318 options_.plot_shock_decomp.use_shock_groups = mygroup; 319 end 320 [za, endo_names, endo_names_tex, steady_state, i_var, oo_] = ... 321 annualized_shock_decomposition(oo_,M_, options_, i_var, t0, options_.nobs, realtime_, vintage_, steady_state,q2a); 322 end 323 end 324end 325 326if ~expand 327 fig_name = fig_name1; 328end 329 330if options_.plot_shock_decomp.use_shock_groups 331 fig_name=[fig_name ' group ' options_.plot_shock_decomp.use_shock_groups]; 332 if expand 333 z = options_.plot_shock_decomp.zfull; 334 endo_names = options_.plot_shock_decomp.endo_names; 335 endo_names_tex = options_.plot_shock_decomp.endo_names_tex; 336 steady_state=steady_state(i_var); 337 i_var = 1; 338 if ~(isequal(type,'aoa') && isstruct(q2a)) 339 gend = size(z,3); 340 endo_nbr = size(z,1); 341 [z, shock_names, M_] = make_the_groups(z,gend,endo_nbr,nshocks,M_,options_); 342 M_.endo_names = endo_names; 343 M_.endo_names_tex = endo_names_tex; 344 else 345 % here we know we only have one variable to handle 346 if isstruct(q2a.aux) && ischar(q2a.aux.y) 347 steady_state_aux = get_mean(q2a.aux.y); 348 q2a.aux.y=repmat(steady_state_aux,16,1); 349 q2a.aux.yss=steady_state_aux; 350 end 351 [~, yssa, ~, gyssa] = ... 352 quarterly2annual(repmat(steady_state,16,1),steady_state,q2a.GYTREND0,q2a.type,q2a.islog,q2a.aux); 353 if q2a.plot==1 354 steady_state = gyssa; 355 else 356 steady_state = yssa; 357 end 358 end 359 else 360 gend = size(z,3); 361 zfull = z; 362 endo_nbr = size(z,1); 363 [z, shock_names, M_] = make_the_groups(z,gend,endo_nbr,nshocks,M_,options_); 364 end 365 if ~isempty(init2shocks) && ~expand 366 M_.exo_names=M_.exo_names_init; 367 end 368else 369 if ~isempty(init2shocks) && ~expand 370 M_.exo_names=M_.exo_names_init; 371 end 372 shock_names = M_.exo_names; 373end 374 375if ~expand 376 if flip_decomp 377 fig_name=[fig_name ' flip']; 378 end 379 if differentiate_decomp 380 fig_name=[fig_name ' diff']; 381 end 382 if ~isempty(init2shocks) 383 fig_name=[fig_name ' init2shocks']; 384 end 385end 386 387func = @(x) colorspace('RGB->Lab',x); 388MAP = distinguishable_colors(size(z,2)-1,'w',func); 389% MAP = [MAP; MAP(end,:)]; 390MAP(end,:) = [0.7 0.7 0.7]; 391% MAP = [MAP; [0.7 0.7 0.7]; [0.3 0.3 0.3]]; 392 393if isempty(options_.plot_shock_decomp.colormap) 394 options_.plot_shock_decomp.colormap = MAP; 395end 396 397if differentiate_decomp 398 z(:,:,2:end) = z(:,:,2:end)-z(:,:,1:end-1); 399 z(:,:,1) = nan; 400 steady_state = steady_state*0; 401end 402 403switch type 404 405 case '' % default 406 407 case 'qoq' 408 409 case 'yoy' 410 z=z(:,:,1:end-3)+z(:,:,2:end-2)+z(:,:,3:end-1)+z(:,:,4:end); 411 if ~isempty(initial_date) 412 initial_date = initial_date+3; 413 else 414 initial_date = dates('1Q4'); 415 end 416 steady_state = 4*steady_state; 417 418 case 'aoa' 419 420 if isempty(initial_date) 421 t0=1; % we assume the sample starts Q1 of 1st year 422 initial_date = dates('1Y'); 423 else 424 initial_date0 = dates([int2str(initial_date.time(1)) 'Y']); 425 if initial_date.time(2)==1 % the first year is full 426 t0=1; 427 initial_date1=initial_date0; 428 else 429 t0=(4-initial_date.time(2)+2); % 1st period of the 1st full year in sample 430 initial_date1=initial_date0+1; 431 end 432 end 433 if realtime_ == 0 434 t0=t0+4-1; % we start in Q4 of the first full year 435 end 436 if isempty(options_.plot_shock_decomp.plot_init_date) && realtime_ == 0 437 options_.plot_shock_decomp.plot_init_date=initial_date+t0; 438 end 439 if isstruct(q2a) 440 if realtime_ == 0 441 if ~isfield(q2a,'type') % private trap for aoa calls 442 q2a.type=1; 443 end 444 if ~isfield(q2a,'islog') % private trap for aoa calls 445 q2a.islog=0; 446 end 447 if ~isfield(q2a,'GYTREND0') % private trap for aoa calls 448 q2a.GYTREND0=0; 449 end 450 if ~isfield(q2a,'aux') % private trap for aoa calls 451 q2a.aux=0; 452 end 453 if ~isfield(q2a,'cumfix') % private trap for aoa calls 454 q2a.cumfix=1; 455 end 456 if ~isfield(q2a,'plot') % private trap for aoa calls 457 q2a.plot=1; % growth rate 458 end 459 460 if ~expand 461 if isstruct(q2a.aux) && ischar(q2a.aux.y) 462 opts=options_; 463 opts.plot_shock_decomp.type='qoq'; 464 opts.plot_shock_decomp.use_shock_groups=[]; 465 [y_aux, steady_state_aux] = plot_shock_decomposition(M_,oo_,opts,q2a.aux.y); 466 q2a.aux.y=y_aux; 467 q2a.aux.yss=steady_state_aux; 468 end 469 i_var0 = i_var; 470 [za, endo_names, endo_names_tex, steady_state, i_var, oo_] = ... 471 annualized_shock_decomposition(z,M_, options_, i_var, t0, options_.nobs, realtime_, vintage_, steady_state,q2a); 472 if options_.plot_shock_decomp.interactive && ~isempty(options_.plot_shock_decomp.use_shock_groups) 473 mygroup = options_.plot_shock_decomp.use_shock_groups; 474 options_.plot_shock_decomp.use_shock_groups=''; 475 zafull = ... 476 annualized_shock_decomposition(zfull(i_var0,:,:),M_, options_, i_var, t0, options_.nobs, realtime_, vintage_, steady_state,q2a); 477 options_.plot_shock_decomp.use_shock_groups = mygroup; 478 end 479 end 480 end 481 z = za; 482 if options_.plot_shock_decomp.interactive && ~isempty(options_.plot_shock_decomp.use_shock_groups) 483 zfull = zafull; 484 end 485 M_.endo_names = endo_names; 486 M_.endo_names_tex = endo_names_tex; 487 % endo_nbr = size(z,1); 488 if realtime_<2 || vintage_ == 0 489 initial_date = initial_date1; 490 else 491 initial_date = initial_date0; 492 end 493 else 494 % this is for quarterly-annualized variables already defined in model, so we can just take Q4 495 t0=4-initial_date.time(2)+1; 496 initial_date = initial_date0; 497 z=z(:,:,t0:4:end); 498 end 499 500 if ~isempty(options_.plot_shock_decomp.plot_init_date) 501 options_.plot_shock_decomp.plot_init_date = dates([int2str(options_.plot_shock_decomp.plot_init_date.time(1)) 'Y']); 502 end 503 if ~isempty(options_.plot_shock_decomp.plot_end_date) 504 options_.plot_shock_decomp.plot_end_date = dates([int2str(options_.plot_shock_decomp.plot_end_date.time(1)) 'Y']); 505 end 506 507 508 otherwise 509 510 error('plot_shock_decomposition:: Wrong type') 511 512end 513 514if flip_decomp 515 z = -z; 516 steady_state = - steady_state; 517end 518if steadystate 519 options_.plot_shock_decomp.steady_state=steady_state; 520end 521 522if nargout == 2 523 out=z(i_var,:,:); 524 steady_state = steady_state(i_var); 525 return 526end 527 528% here we crop data if needed 529my_initial_date = initial_date; 530a = 1; 531b = size(z,3); 532if ~isempty(options_.plot_shock_decomp.plot_init_date) 533 my_initial_date = max(initial_date,options_.plot_shock_decomp.plot_init_date); 534 a = find((initial_date:initial_date+b-1)==options_.plot_shock_decomp.plot_init_date); 535end 536if ~isempty(options_.plot_shock_decomp.plot_end_date) 537 if options_.plot_shock_decomp.plot_end_date<=(max(initial_date:initial_date+b-1)) 538 b = find((initial_date:initial_date+b-1)==options_.plot_shock_decomp.plot_end_date); 539 else 540 warning('You set plot_end_date larger than smoother size!!'); 541 end 542end 543z = z(:,:,a:b); 544% end crop data 545 546options_.plot_shock_decomp.fig_name=fig_name; 547options_.plot_shock_decomp.orig_varlist = varlist; 548if options_.plot_shock_decomp.interactive && ~isempty(options_.plot_shock_decomp.use_shock_groups) 549 options_.plot_shock_decomp.zfull = zfull; 550end 551 552if ~options_.no_graph.plot_shock_decomposition 553 if detail_plot 554 graph_decomp_detail(z, shock_names, M_.endo_names, i_var, my_initial_date, M_, options_); 555 else 556 graph_decomp(z, shock_names, M_.endo_names, i_var, my_initial_date, M_, options_); 557 end 558end 559 560if write_xls 561 WriteShockDecomp2Excel(z,shock_names,M_.endo_names,i_var,my_initial_date,M_,options_,options_.plot_shock_decomp); 562end 563 564end 565 566function [z, shock_names, M_] = make_the_groups(z,gend,endo_nbr,nshocks,M_,options_) 567 568shock_groups = M_.shock_groups.(options_.plot_shock_decomp.use_shock_groups); 569shock_ind = fieldnames(shock_groups); 570ngroups = length(shock_ind); 571shock_names = shock_ind; 572shock_varexo = shock_ind; 573for i=1:ngroups 574 shock_names{i} = (shock_groups.(shock_ind{i}).label); 575 if isfield(M_,'exo_names_init') 576 shock_varexo{i} = (shock_groups.(shock_ind{i}).shocks); 577 end 578end 579zz = zeros(endo_nbr,ngroups+2,gend); 580kcum=[]; 581for i=1:ngroups 582 indx=0; 583 for j = shock_groups.(shock_ind{i}).shocks 584 k = find(strcmp(j,cellstr(M_.exo_names))); 585 if isfield(M_,'exo_names_init') 586 indx=indx+1; 587 shock_varexo{i}{indx} = M_.exo_names_init{k}; 588 end 589 zz(:,i,:) = zz(:,i,:) + z(:,k,:); 590 z(:,k,:) = 0; 591 kcum = [kcum k]; 592 end 593 if isfield(M_,'exo_names_init') 594 shock_groups.(shock_ind{i}).shocks = shock_varexo{i}; 595 end 596end 597zothers = sum(z(:,1:nshocks,:),2); 598shock_groups.(['group' int2str(ngroups+1)]).label = 'Others'; 599shock_groups.(['group' int2str(ngroups+1)]).shocks = cellstr(M_.exo_names(find(~ismember([1:M_.exo_nbr],kcum)),:))'; 600M_.shock_groups.(options_.plot_shock_decomp.use_shock_groups)=shock_groups; 601if any(any(zothers)) 602 shock_names = [shock_names; {'Others + Initial Values'}]; 603end 604zz(:,ngroups+1,:) = sum(z(:,1:nshocks+1,:),2); 605zz(:,ngroups+2,:) = z(:,nshocks+2,:); 606z = zz; 607 608end 609