1function data_set = det_cond_forecast(varargin) 2% Computes conditional forecasts using the extended path method. 3% 4% INPUTS 5% o plan [structure] A structure describing the different shocks and the endogenous varibales, the date of the shocks and the path of the shock. 6% The plan structure is created by the functions init_plan, basic_plan and flip_plan 7% o [dataset] [dseries] A dserie containing the initial values of the shocks and the endogenous variables (usually the dseries generated by the smoother). 8% o [starting_date] [dates] The first date of the forecast. 9% 10% 11% OUTPUTS 12% dataset [dseries] Returns a dseries containing the forecasted endgenous variables and shocks 13% 14 15% Copyright (C) 2013-2018 Dynare Team 16% 17% This file is part of Dynare. 18% 19% Dynare is free software: you can redistribute it and/or modify 20% it under the terms of the GNU General Public License as published by 21% the Free Software Foundation, either version 3 of the License, or 22% (at your option) any later version. 23% 24% Dynare is distributed in the hope that it will be useful, 25% but WITHOUT ANY WARRANTY; without even the implied warranty of 26% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 27% GNU General Public License for more details. 28% 29% You should have received a copy of the GNU General Public License 30% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 31 32global options_ oo_ M_ 33pp = 2; 34initial_conditions = oo_.steady_state; 35verbosity = options_.verbosity; 36if options_.periods == 0 37 options_.periods = 25; 38end 39if isempty(options_.qz_criterium) 40 options_.qz_criterium = 1+1e-6; 41end 42%We have to get an initial guess for the conditional forecast 43% and terminal conditions for the non-stationary variables, we 44% use the first order approximation of the rational expectation solution. 45if ~isfield(oo_,'dr') || ~isfield(oo_.dr,'ghx') 46 fprintf('computing the first order solution of the model as initial guess...'); 47 dr = struct(); 48 oo_.dr=set_state_space(dr,M_,options_); 49 options_.order = 1; 50 [dr,Info,M_,options_,oo_] = resol(0,M_,options_,oo_); 51 fprintf('done\n'); 52end 53b_surprise = 0; 54b_pf = 0; 55surprise = 0; 56pf = 0; 57is_shock = []; 58is_constraint = []; 59if length(varargin) > 3 60 % regular way to call 61 constrained_paths = varargin{1}; 62 max_periods_simulation = size(constrained_paths, 2); 63 constrained_vars = varargin{2}; 64 options_cond_fcst = varargin{3}; 65 constrained_perfect_foresight = varargin{4}; 66 constraint_index = cell(max_periods_simulation,1); 67 nvars = length(constrained_vars); 68 for i = 1:max_periods_simulation 69 constraint_index{i} = 1:nvars; 70 end 71 direct_mode = 0; 72 shocks_present = 0; 73 controlled_varexo = options_cond_fcst.controlled_varexo; 74 nvarexo = size(controlled_varexo, 1); 75 options_cond_fcst.controlled_varexo = zeros(nvarexo,1); 76 exo_names = M_.exo_names; 77 for i = 1:nvarexo 78 j = find(strcmp(controlled_varexo(i,:), exo_names)); 79 if ~isempty(j) 80 options_cond_fcst.controlled_varexo(i) = j; 81 else 82 error(['Unknown exogenous variable ' controlled_varexo(i,:)]); 83 end 84 end 85 86else 87 % alternative way to call: plan, dset, dates_of_frcst 88 plan = varargin{1}; 89 if length(varargin) >= 2 90 dset = varargin{2}; 91 if ~isa(dset,'dseries') 92 error('the second argmuent should be a dseries'); 93 end 94 if length(varargin) >= 3 95 range = varargin{3}; 96 if ~isa(range,'dates') 97 error('the third argmuent should be a dates'); 98 end 99 %if (range(range.ndat) > dset.time(dset.nobs) ) 100 if (range(range.ndat) > dset.dates(dset.nobs)+1 ) 101 s1 = strings(dset.dates(dset.nobs)); 102 s2 = strings(range(range.ndat)); 103 error(['the dseries ' inputname(2) ' finish at time ' s1{1} ' before the last period of forecast ' s2{1}]); 104 end 105 106 sym_dset = dset(dates(-range(1)):dates(range(range.ndat))); 107 periods = options_.periods + M_.maximum_lag + M_.maximum_lead; 108 total_periods = periods + range.ndat; 109 if isfield(oo_, 'exo_simul') 110 if size(oo_.exo_simul, 1) ~= total_periods 111 oo_.exo_simul = repmat(oo_.exo_steady_state',total_periods,1); 112 end 113 else 114 oo_.exo_simul = repmat(oo_.exo_steady_state',total_periods,1); 115 end 116 117 oo_.endo_simul = repmat(oo_.steady_state, 1, total_periods); 118 119 for i = 1:sym_dset.vobs 120 iy = find(strcmp(strtrim(sym_dset.name{i}), strtrim(plan.endo_names))); 121 if ~isempty(iy) 122 oo_.endo_simul(iy,1:sym_dset.nobs) = sym_dset.data(:, i); 123 initial_conditions(iy) = sym_dset.data(1, i); 124 else 125 ix = find(strcmp(strtrim(sym_dset.name{i}), strtrim(plan.exo_names))); 126 if ~isempty(ix) 127 oo_.exo_simul(1, ix) = sym_dset.data(1, i)'; 128 else 129 %warning(['The variable ' sym_dset.name{i} ' in the dataset ' inputname(2) ' is not a endogenous neither an exogenous variable!!']); 130 end 131 end 132 end 133 for i = 1:length(M_.aux_vars) 134 if M_.aux_vars(i).type == 1 %lag variable 135 iy = find(strcmp(M_.endo_names{M_.aux_vars(i).orig_index}, sym_dset.name)); 136 if ~isempty(iy) 137 oo_.endo_simul(M_.aux_vars(i).endo_index, 1:sym_dset.nobs) = dset(dates(range(1) + (M_.aux_vars(i).orig_lead_lag - 1))).data(:,iy); 138 initial_conditions(M_.aux_vars(i).endo_index) = dset(dates(range(1) + (M_.aux_vars(i).orig_lead_lag - 1))).data(:,iy); 139 else 140 warning(['The variable auxiliary ' M_.endo_names{M_.aux_vars(i).endo_index} ' associated to the variable ' M_.endo_names{M_.aux_vars(i).orig_index} ' do not appear in the dataset']); 141 end 142 else 143 oo_.endo_simul(M_.aux_vars(i).endo_index, 1:sym_dset.nobs) = repmat(oo_.steady_state(M_.aux_vars(i).endo_index), 1, range.ndat + 1); 144 end 145 end 146 %Compute the initial path using the the steady-state 147 % steady-state 148 %for jj = 2 : (options_.periods + 2) 149 for jj = 2 : (range.ndat + 2) 150 oo_.endo_simul(:, jj) = oo_.steady_state; 151 end 152 missings = isnan(oo_.endo_simul(:,1)); 153 if any(missings) 154 for jj = 1:M_.endo_nbr 155 if missings(jj) 156 oo_.endo_simul(jj,1) = oo_.steady_state(jj,1); 157 end 158 end 159 end 160 161 if options_.bytecode 162 save_options_dynatol_f = options_.dynatol.f; 163 options_.dynatol.f = 1e-7; 164 [Info, endo, exo] = bytecode('extended_path', plan, oo_.endo_simul, oo_.exo_simul, M_.params, oo_.steady_state, options_.periods); 165 options_.dynatol.f = save_options_dynatol_f; 166 167 if Info == 0 168 oo_.endo_simul = endo; 169 oo_.exo_simul = exo; 170 end 171 endo = endo'; 172 endo_l = size(endo(1+M_.maximum_lag:end,:),1); 173 jrng = dates(plan.date(1)):dates(plan.date(1)+endo_l); 174 data_set = dseries(nan(endo_l, dset.vobs), plan.date(1), dset.name); 175 for i = 1:length(dset.name) 176 pos = find(strcmp(dset.name{i},plan.endo_names)); 177 if ~isempty(pos) 178 data_set.(dset.name{i}) = dseries(endo(1+M_.maximum_lag:end,pos), plan.date(1), dset.name{i}); 179 else 180 pos = find(strcmp(dset.name{i},plan.exo_names)); 181 if ~isempty(pos) 182 data_set{dset.name{i}} = dseries(exo(1+M_.maximum_lag:end,pos), plan.date(1),dset.name{i}); 183 end 184 end 185 end 186 data_set = [dset(dset.dates(1):(plan.date(1)-1)) ; data_set]; 187 for i=1:M_.exo_nbr 188 pos = find(strcmp(M_.exo_names{i}, dset.name)); 189 if isempty(pos) 190 data_set{M_.exo_names{i}} = dseries(exo(1+M_.maximum_lag:end,i), plan.date(1), M_.exo_names{i}); 191 else 192 data_set{M_.exo_names{i}}(plan.date(1):plan.date(1)+ (size(exo, 1) - M_.maximum_lag)) = exo(1+M_.maximum_lag:end,i); 193 end 194 end 195 data_set = merge(dset(dset.dates(1):(plan.date(1)-1)), data_set); 196 return 197 union_names = union(data_set.name, dset.name); 198 dif = setdiff(union_names, data_set.name); 199 data_set_nobs = data_set.nobs; 200 for i = 1:length(dif) 201 data_set{dif{i}} = dseries(nan(data_set_nobs,1),plan.date(1), dif(i), dif(i)); 202 end 203 dif = setdiff(union_names, dset.name); 204 dset_nobs = dset.nobs; 205 for i = 1:length(dif) 206 dset{dif{i}} = dseries(nan(dset_nobs,1),dset.dates(1), dif(i), dif(i)); 207 end 208 data_set = [dset(dset.dates(1):(plan.date(1)-1)) ; data_set]; 209 return 210 end 211 else 212 error('impossible case'); 213 end 214 215 else 216 oo_.exo_simul = repmat(oo_.exo_steady_state',options_.periods+2,1); 217 oo_.endo_simul = repmat(oo_.steady_state, 1, options_.periods+2); 218 end 219 220 direct_mode = 1; 221 constrained_paths = plan.constrained_paths_; 222 constrained_vars = plan.constrained_vars_; 223 options_cond_fcst = plan.options_cond_fcst_; 224 constrained_perfect_foresight = plan.constrained_perfect_foresight_; 225 constrained_periods = plan.constrained_date_; 226 if ~isempty(plan.shock_paths_) 227 shock_paths = plan.shock_paths_; 228 shock_vars = plan.shock_vars_; 229 shock_perfect_foresight = plan.shock_perfect_foresight_; 230 shock_periods = plan.shock_date_; 231 shocks_present = 1; 232 else 233 shocks_present = 0; 234 end 235 236 total_periods = plan.date; 237 238end 239 240if ~isfield(options_cond_fcst,'periods') || isempty(options_cond_fcst.periods) 241 options_cond_fcst.periods = 100; 242end 243 244options_.periods = 10; 245 246if direct_mode == 1 247 n_periods = length(constrained_periods); 248 is_constraint = zeros(length(total_periods), n_periods); 249 constrained_paths_cell = constrained_paths; 250 clear constrained_paths; 251 constrained_paths = zeros(n_periods, length(total_periods)); 252 max_periods_simulation = 0; 253 for i = 1:n_periods 254 period_i = constrained_periods{i}; 255 %period_i 256 tp = total_periods(1); 257 if size(period_i) > 1 258 init_periods = period_i(1); 259 tp_end = period_i(end); 260 else 261 init_periods = period_i; 262 tp_end = period_i; 263 end 264 tp0 = tp; 265 while tp < init_periods 266 tp = tp + 1; 267 end 268 j = 0; 269 while tp <= tp_end 270 is_constraint(tp - tp0 + 1, i) = 1; 271 constrained_paths(i, tp - tp0 + 1) = constrained_paths_cell{i}(j + 1); 272 tp = tp + 1; 273 j = j + 1; 274 end 275 if tp - tp0 > max_periods_simulation 276 max_periods_simulation = tp - tp0; 277 end 278 end 279 n_nnz = length(sum(is_constraint,2)); 280 if n_nnz > 0 281 constraint_index = cell(n_nnz,1); 282 for i= 1:n_nnz 283 constraint_index{i} = find(is_constraint(i,:)); 284 end 285 end 286 if shocks_present 287 n_periods = length(shock_periods); 288 shock_paths_cell = shock_paths; 289 clear shock_paths; 290 shock_paths = zeros(n_periods, length(total_periods)); 291 is_shock = zeros(length(total_periods), n_periods); 292 for i = 1:n_periods 293 period_i = shock_periods{i}; 294 %period_i 295 tp = total_periods(1); 296 if size(period_i) > 1 297 init_periods = period_i(1); 298 tp_end = period_i(end); 299 else 300 init_periods = period_i; 301 tp_end = period_i; 302 end 303 tp0 = tp; 304 while tp < init_periods 305 tp = tp + 1; 306 end 307 j = 0; 308 while tp <= tp_end 309 is_shock(tp - tp0 + 1, i) = 1; 310 shock_paths(i, tp - tp0 + 1) = shock_paths_cell{i}(j + 1); 311 tp = tp + 1; 312 j = j + 1; 313 end 314 if tp - tp0 > max_periods_simulation 315 max_periods_simulation = tp - tp0; 316 end 317 end 318 n_nnz = length(sum(is_shock,2)); 319 if n_nnz > 0 320 shock_index = cell(n_nnz, 1); 321 for i= 1:n_nnz 322 shock_index{i} = find(is_shock(i,:)); 323 end 324 end 325 end 326else 327 is_constraint = ones(size(constrained_paths)); 328end 329 330 331 332maximum_lag = M_.maximum_lag; 333 334ys = oo_.steady_state; 335ny = size(ys,1); 336xs = [oo_.exo_steady_state ; oo_.exo_det_steady_state]; 337nx = size(xs,1); 338 339constrained_periods = max_periods_simulation; 340n_endo_constrained = size(constrained_vars,1); 341if isfield(options_cond_fcst,'controlled_varexo') 342 n_control_exo = size(options_cond_fcst.controlled_varexo, 1); 343 if n_control_exo ~= n_endo_constrained 344 error(['det_cond_forecast: the number of exogenous controlled variables (' int2str(n_control_exo) ') has to be equal to the number of constrained endogenous variabes (' int2str(n_endo_constrained) ')']) 345 end 346else 347 error('det_cond_forecast: to run a deterministic conditional forecast you have to specified the exogenous variables controlled using the option controlled_varexo in forecast command'); 348end 349 350% if n_endo_constrained == 0 351% options_.ep.use_bytecode = options_.bytecode; 352% data_set = extended_path(initial_conditions, max_periods_simulation); 353% end 354 355if length(varargin) >= 1 356 controlled_varexo = options_cond_fcst.controlled_varexo; 357else 358 exo_names = M_.exo_names; 359 controlled_varexo = zeros(1,n_control_exo); 360 for i = 1:nx 361 for j=1:n_control_exo 362 if strcmp(exo_names{i}, deblank(options_cond_fcst.controlled_varexo(j,:))) 363 controlled_varexo(j) = i; 364 end 365 end 366 end 367end 368 369%todo check if zero => error message 370 371save_options_initval_file = options_.initval_file; 372options_.initval_file = '__'; 373 374[pos_constrained_pf, ~] = find(constrained_perfect_foresight); 375indx_endo_solve_pf = constrained_vars(pos_constrained_pf); 376if isempty(indx_endo_solve_pf) 377 pf = 0; 378else 379 pf = length(indx_endo_solve_pf); 380end 381indx_endo_solve_surprise = setdiff(constrained_vars, indx_endo_solve_pf); 382 383if isempty(indx_endo_solve_surprise) 384 surprise = 0; 385else 386 surprise = length(indx_endo_solve_surprise); 387end 388 389 390eps = options_.solve_tolf; 391maxit = options_.simul.maxit; 392 393 394past_val = 0; 395save_options_periods = options_.periods; 396options_.periods = options_cond_fcst.periods; 397save_options_dynatol_f = options_.dynatol.f; 398options_.dynatol.f = 1e-8; 399eps1 = 1e-7;%1e-4; 400exo = zeros(maximum_lag + options_cond_fcst.periods, nx); 401endo = zeros(maximum_lag + options_cond_fcst.periods, ny); 402endo(1,:) = oo_.steady_state'; 403 404 405% if all the endogenous paths are perfectly anticipated we do not need to 406% implement the extended path 407if pf && ~surprise 408 time_index_constraint = maximum_lag + 1:maximum_lag + constrained_periods; 409 second_system_size = pf * constrained_periods; 410 J = zeros(second_system_size,second_system_size); 411 r = zeros(second_system_size,1); 412 indx_endo = zeros(second_system_size,1); 413 col_count = 1; 414 for j = 1:length(constrained_vars) 415 indx_endo(col_count : col_count + constrained_periods - 1) = constrained_vars(j) + (time_index_constraint - 1) * ny; 416 col_count = col_count + constrained_periods; 417 end 418 oo_=make_ex_(M_,options_,oo_); 419 oo_=make_y_(M_,options_,oo_); 420 it = 1; 421 convg = 0; 422 normra = 1e+50; 423 while ~convg && it <= maxit 424 disp('---------------------------------------------------------------------------------------------'); 425 disp(['iteration ' int2str(it)]); 426 not_achieved = 1; 427 alpha = 1; 428 while not_achieved 429 simul(); 430 result = sum(sum(isfinite(oo_.endo_simul(:,time_index_constraint)))) == ny * constrained_periods; 431 if result 432 y = oo_.endo_simul(constrained_vars, time_index_constraint); 433 ys = oo_.endo_simul(indx_endo); 434 col_count = 1; 435 for j = 1:length(constrained_vars) 436 y = oo_.endo_simul(constrained_vars(j), time_index_constraint); 437 r(col_count:col_count+constrained_periods - 1) = (y - constrained_paths(j,1:constrained_periods))'; 438 col_count = col_count + constrained_periods; 439 end 440 normr = norm(r, 1); 441 end 442 if (~oo_.deterministic_simulation.status || ~result || normr > normra) && it > 1 443 not_achieved = 1; 444 alpha = alpha / 2; 445 col_count = 1; 446 for j = controlled_varexo' 447 oo_.exo_simul(time_index_constraint,j) = (old_exo(:,j) + alpha * D_exo(col_count: (col_count + constrained_periods - 1))); 448 col_count = col_count + constrained_periods; 449 end 450 disp(['Divergence in Newton: reducing the path length alpha=' num2str(alpha)]); 451 oo_.endo_simul = repmat(oo_.steady_state, 1, options_cond_fcst.periods + 2); 452 else 453 not_achieved = 0; 454 end 455 end 456 457 per = 30; 458 z = oo_.endo_simul(:, 1 : per + 2 ); 459 zx = oo_.exo_simul(1: per + 2,:); 460 g1 = spalloc(M_.endo_nbr * (per ), M_.endo_nbr * (per ), 3* M_.endo_nbr * per ); 461 g1_x = spalloc(M_.endo_nbr * (per ), M_.exo_nbr, M_.endo_nbr * (per )* M_.exo_nbr ); 462 lag_indx = find(M_.lead_lag_incidence(1,:)); 463 cur_indx = M_.endo_nbr + find(M_.lead_lag_incidence(2,:)); 464 lead_indx = 2 * M_.endo_nbr + find(M_.lead_lag_incidence(3,:)); 465 cum_l1 = 0; 466 cum_index_d_y_x = []; 467 indx_x = []; 468 for k = 1 : per 469 if k == 1 470 if (isfield(M_,'block_structure')) 471 data1 = M_.block_structure.block; 472 Size = length(M_.block_structure.block); 473 else 474 data1 = M_; 475 Size = 1; 476 end 477 data1 = M_; 478 if (options_.bytecode) 479 [chck, zz, data1]= bytecode('dynamic','evaluate', z, zx, M_.params, oo_.steady_state, k, data1); 480 else 481 [zz, g1b] = feval([M_.fname '.dynamic'], z', zx, M_.params, oo_.steady_state, k); 482 data1.g1_x = g1b(:,end - M_.exo_nbr + 1:end); 483 data1.g1 = g1b(:,1 : end - M_.exo_nbr); 484 chck = 0; 485 end 486 mexErrCheck('bytecode', chck); 487 end 488 if k == 1 489 g1(1:M_.endo_nbr,-M_.endo_nbr + [cur_indx lead_indx]) = data1.g1(:,M_.nspred + 1:end); 490 elseif k == per 491 g1(M_.endo_nbr * (k - 1) + 1 :M_.endo_nbr * k,M_.endo_nbr * (k -2) + [lag_indx cur_indx]) = data1.g1(:,1:M_.nspred + M_.endo_nbr); 492 else 493 g1(M_.endo_nbr * (k - 1) + 1 :M_.endo_nbr * k, M_.endo_nbr * (k -2) + [lag_indx cur_indx lead_indx]) = data1.g1; 494 end 495 l2 = 1; 496 pf_c = 1; 497 if k <= constrained_periods 498 for l = constraint_index{k} 499 l1 = controlled_varexo(l); 500 g1_x(M_.endo_nbr * (k - 1) + 1:M_.endo_nbr * k,1 + cum_l1) = data1.g1_x(:,l1); 501 if k == 1 502 indx_x(l2) = l ; 503 l2 = l2 + 1; 504 for ii = 2:constrained_periods 505 indx_x(l2) = length(controlled_varexo) + pf * (ii - 2) + constraint_index{k + ii - 1}(pf_c); 506 l2 = l2 + 1; 507 end 508 pf_c = pf_c + 1; 509 cum_index_d_y_x = [cum_index_d_y_x; constrained_vars(l)]; 510 else 511 cum_index_d_y_x = [cum_index_d_y_x; constrained_vars(l) + (k - 1) * M_.endo_nbr]; 512 end 513 cum_l1 = cum_l1 + length(l1); 514 end 515 end 516 end 517 518 d_y_x = - g1 \ g1_x; 519 520 cum_l1 = 0; 521 count_col = 1; 522 cum_index_J = 1:length(cum_index_d_y_x(indx_x)); 523 J= zeros(length(cum_index_J)); 524 for j1 = 1:length(controlled_varexo) 525 cum_l1 = 0; 526 for k = 1:(constrained_periods) 527 l1 = constraint_index{k}; 528 l1 = find(constrained_perfect_foresight(l1) | (k == 1)); 529 if constraint_index{k}( j1) 530 J(cum_index_J,count_col) = d_y_x(cum_index_d_y_x(indx_x),indx_x(count_col)); 531 count_col = count_col + 1; 532 end 533 cum_l1 = cum_l1 + length(l1); 534 end 535 cum_l1 = cum_l1 + length(constrained_vars(j1)); 536 end 537 538 539 % col_count = 1; 540 % for j = controlled_varexo' 541 % for time = time_index_constraint 542 % saved = oo_.exo_simul(time,j); 543 % oo_.exo_simul(time,j) = oo_.exo_simul(time,j) + eps1; 544 % simul(); 545 % J1(:,col_count) = (oo_.endo_simul(indx_endo) - ys) / eps1; 546 % oo_.exo_simul(time,j) = saved; 547 % col_count = col_count + 1; 548 % end 549 % end 550 % J1 551 % sdfmlksdf; 552 553 disp(['iteration ' int2str(it) ' error = ' num2str(normr)]); 554 555 if normr <= eps 556 convg = 1; 557 disp('convergence achieved'); 558 else 559 % Newton update on exogenous shocks 560 old_exo = oo_.exo_simul(time_index_constraint,:); 561 D_exo = - J \ r; 562 col_count = 1; 563 %constrained_periods 564 for j = controlled_varexo' 565 oo_.exo_simul(time_index_constraint,j) = oo_.exo_simul(time_index_constraint,j) + D_exo(col_count: (col_count + constrained_periods - 1)); 566 col_count = col_count + constrained_periods - 1; 567 end 568 end 569 it = it + 1; 570 normra = normr; 571 end 572 endo = oo_.endo_simul'; 573 exo = oo_.exo_simul; 574else 575 for t = 1:constrained_periods 576 577 if direct_mode && ~isempty(is_constraint) 578 [pos_constrained_pf, ~] = find(constrained_perfect_foresight .* is_constraint(t, :)'); 579 indx_endo_solve_pf = constrained_vars(pos_constrained_pf); 580 if isempty(indx_endo_solve_pf) 581 pf = 0; 582 else 583 pf = length(indx_endo_solve_pf); 584 end 585 586 [pos_constrained_surprise, ~] = find((1-constrained_perfect_foresight) .* is_constraint(t, :)'); 587 indx_endo_solve_surprise = constrained_vars(pos_constrained_surprise); 588 589 if isempty(indx_endo_solve_surprise) 590 surprise = 0; 591 else 592 surprise = length(indx_endo_solve_surprise); 593 end 594 end 595 596 if direct_mode && ~isempty(is_shock) 597 [pos_shock_pf, ~] = find(shock_perfect_foresight .* is_shock(t, :)'); 598 indx_endo_solve_pf = shock_vars(pos_shock_pf); 599 if isempty(indx_endo_solve_pf) 600 b_pf = 0; 601 else 602 b_pf = length(indx_endo_solve_pf); 603 end 604 605 [pos_shock_surprise, ~] = find((1-shock_perfect_foresight) .* is_shock(t, :)'); 606 indx_endo_solve_surprise = shock_vars(pos_shock_surprise); 607 608 if isempty(indx_endo_solve_surprise) 609 b_surprise = 0; 610 else 611 b_surprise = length(indx_endo_solve_surprise); 612 end 613 end 614 615 disp('==============================================================================================='); 616 disp(['t=' int2str(t) ' conditional (surprise=' int2str(surprise) ' perfect foresight=' int2str(pf) ') unconditional (surprise=' int2str(b_surprise) ' perfect foresight=' int2str(b_pf) ')']); 617 disp('==============================================================================================='); 618 if t == 1 619 oo_=make_ex_(M_,options_,oo_); 620 if maximum_lag > 0 621 exo_init = oo_.exo_simul; 622 else 623 exo_init = zeros(size(oo_.exo_simul)); 624 end 625 oo_=make_y_(M_,options_,oo_); 626 end 627 %exo_init 628 oo_.exo_simul = exo_init; 629 oo_.endo_simul(:,1) = initial_conditions; 630 631 time_index_constraint = maximum_lag + 1:maximum_lag + constrained_periods - t + 1; 632 633 if direct_mode 634 nb_shocks = length(plan.shock_vars_); 635 if nb_shocks > 0 636 shock_index_t = shock_index{t}; 637 shock_vars_t = shock_vars(shock_index_t); 638 for shock_indx = shock_index_t 639 if shock_perfect_foresight(shock_indx) 640 oo_.exo_simul(time_index_constraint,shock_vars_t(shock_indx)) = shock_paths(shock_indx,t:constrained_periods); 641 else 642 oo_.exo_simul(maximum_lag + 1,shock_vars_t(shock_indx)) = shock_paths(shock_indx,t); 643 end 644 end 645 end 646 end 647 %Compute the initial path using the first order solution around the 648 % steady-state 649 oo_.endo_simul(:, 1) = oo_.endo_simul(:, 1) - oo_.steady_state; 650 for jj = 2 : (options_.periods + 2) 651 oo_.endo_simul(:,jj) = oo_.dr.ghx(oo_.dr.inv_order_var,:) * oo_.endo_simul(oo_.dr.state_var, jj-1) + oo_.dr.ghu * oo_.exo_simul(jj,:)'; 652 end 653 for jj = 1 : (options_.periods + 2) 654 oo_.endo_simul(:,jj) = oo_.steady_state + oo_.endo_simul(:,jj); 655 end 656 657 constraint_index_t = constraint_index{t}; 658 controlled_varexo_t = controlled_varexo(constraint_index_t); 659 constrained_vars_t = constrained_vars(constraint_index_t); 660 661 second_system_size = surprise + pf * (constrained_periods - t + 1); 662 indx_endo = zeros(second_system_size,1); 663 col_count = 1; 664 for j = 1:length(constrained_vars_t) 665 if constrained_perfect_foresight(j) 666 indx_endo(col_count : col_count + constrained_periods - t) = constrained_vars(j) + (time_index_constraint - 1) * ny; 667 col_count = col_count + constrained_periods - t + 1; 668 else 669 indx_endo(col_count) = constrained_vars(j) + maximum_lag * ny; 670 col_count = col_count + 1; 671 end 672 end 673 674 r = zeros(second_system_size,1); 675 676 convg = 0; 677 it = 1; 678 while ~convg && it <= maxit 679 disp('-----------------------------------------------------------------------------------------------'); 680 disp(['iteration ' int2str(it) ' time ' int2str(t)]); 681 not_achieved = 1; 682 alpha = 1; 683 while not_achieved 684 perfect_foresight_setup; 685 perfect_foresight_solver; 686 result = sum(sum(isfinite(oo_.endo_simul(:,time_index_constraint)))) == ny * length(time_index_constraint); 687 if (~oo_.deterministic_simulation.status || ~result) && it > 1 688 not_achieved = 1; 689 alpha = alpha / 2; 690 col_count = 1; 691 for j1 = constraint_index_t 692 j = controlled_varexo(j1); 693 if constrained_perfect_foresight(j1) 694 oo_.exo_simul(time_index_constraint,j) = (old_exo(time_index_constraint,j) + alpha * D_exo(col_count: col_count + constrained_periods - t)); 695 col_count = col_count + constrained_periods - t + 1; 696 else 697 oo_.exo_simul(maximum_lag + 1,j) = old_exo(maximum_lag + 1,j) + alpha * D_exo(col_count); 698 col_count = col_count + 1; 699 end 700 end 701 disp(['Divergence in Newton: reducing the path length alpha=' num2str(alpha) ' result=' num2str(result) ' sum(sum)=' num2str(sum(sum(isfinite(oo_.endo_simul(:,time_index_constraint))))) ' ny * length(time_index_constraint)=' num2str(ny * length(time_index_constraint)) ' oo_.deterministic_simulation.status=' num2str(oo_.deterministic_simulation.status)]); 702 oo_.endo_simul = [initial_conditions repmat(oo_.steady_state, 1, options_cond_fcst.periods + 1)]; 703 else 704 not_achieved = 0; 705 end 706 end 707 if t==constrained_periods 708 ycc = oo_.endo_simul; 709 end 710 yc = oo_.endo_simul(:,maximum_lag + 1); 711 ys = oo_.endo_simul(indx_endo); 712 713 col_count = 1; 714 for j = constraint_index_t 715 if constrained_perfect_foresight(j) 716 y = oo_.endo_simul(constrained_vars(j), time_index_constraint); 717 r(col_count:col_count+constrained_periods - t) = (y - constrained_paths(j,t:constrained_periods))'; 718 col_count = col_count + constrained_periods - t + 1; 719 else 720 y = yc(constrained_vars(j)); 721 r(col_count) = y - constrained_paths(j,t); 722 col_count = col_count + 1; 723 end 724 end 725 726 disp('computation of derivatives w.r. to exogenous shocks'); 727 728 per = 30; 729 z = oo_.endo_simul(:, 1 : per + 2 ); 730 zx = oo_.exo_simul(1: per + 2,:); 731 g1 = spalloc(M_.endo_nbr * (per ), M_.endo_nbr * (per ), 3* M_.endo_nbr * per ); 732 g1_x = spalloc(M_.endo_nbr * (per ), M_.exo_nbr, M_.endo_nbr * (per )* M_.exo_nbr ); 733 lag_indx = find(M_.lead_lag_incidence(1,:)); 734 cur_indx = M_.endo_nbr + find(M_.lead_lag_incidence(2,:)); 735 lead_indx = 2 * M_.endo_nbr + find(M_.lead_lag_incidence(3,:)); 736 cum_l1 = 0; 737 %indx_x = zeros(length(constraint_index_t)* constrained_periods, 1); 738 cum_index_d_y_x = []; 739 indx_x = []; 740 for k = 1 : per 741 if k == 1 742 if (isfield(M_,'block_structure')) 743 data1 = M_.block_structure.block; 744 Size = length(M_.block_structure.block); 745 else 746 data1 = M_; 747 Size = 1; 748 end 749 data1 = M_; 750 if (options_.bytecode) 751 [chck, zz, data1]= bytecode('dynamic','evaluate', z, zx, M_.params, oo_.steady_state, k, data1); 752 else 753 [zz, g1b] = feval([M_.fname '.dynamic'], z', zx, M_.params, oo_.steady_state, k); 754 data1.g1_x = g1b(:,end - M_.exo_nbr + 1:end); 755 data1.g1 = g1b(:,1 : end - M_.exo_nbr); 756 chck = 0; 757 end 758 mexErrCheck('bytecode', chck); 759 end 760 if k == 1 761 g1(1:M_.endo_nbr,-M_.endo_nbr + [cur_indx lead_indx]) = data1.g1(:,M_.nspred + 1:end); 762 elseif k == per 763 g1(M_.endo_nbr * (k - 1) + 1 :M_.endo_nbr * k,M_.endo_nbr * (k -2) + [lag_indx cur_indx]) = data1.g1(:,1:M_.nspred + M_.endo_nbr); 764 else 765 g1(M_.endo_nbr * (k - 1) + 1 :M_.endo_nbr * k, M_.endo_nbr * (k -2) + [lag_indx cur_indx lead_indx]) = data1.g1; 766 end 767 768 l2 = 1; 769 pf_c = 1; 770 if t + k - 1 <= constrained_periods 771 for l = constraint_index{t + k - 1} 772 l1 = controlled_varexo(l); 773 if (k == 1) || ((k > 1) && constrained_perfect_foresight(l)) 774 g1_x(M_.endo_nbr * (k - 1) + 1:M_.endo_nbr * k,1 + cum_l1) = data1.g1_x(:,l1); 775 if k == 1 776 if constrained_perfect_foresight(l) 777 indx_x(l2) = l ; 778 l2 = l2 + 1; 779 for ii = 2:constrained_periods - t + 1 780 indx_x(l2) = length(controlled_varexo) + pf * (ii - 2) + constraint_index{k + ii - 1}(pf_c); 781 l2 = l2 + 1; 782 end 783 pf_c = pf_c + 1; 784 else 785 indx_x(l2) = l; 786 l2 = l2 + 1; 787 end 788 cum_index_d_y_x = [cum_index_d_y_x; constrained_vars_t(l)]; 789 else 790 cum_index_d_y_x = [cum_index_d_y_x; constrained_vars_t(l) + (k - 1) * M_.endo_nbr]; 791 end 792 cum_l1 = cum_l1 + length(l1); 793 end 794 end 795 end 796 end 797 798 d_y_x = - g1 \ g1_x; 799 800 801 cum_l1 = 0; 802 count_col = 1; 803 cum_index_J = 1:length(cum_index_d_y_x(indx_x)); 804 J= zeros(length(cum_index_J)); 805 for j1 = constraint_index_t 806 if constrained_perfect_foresight(j1) 807 cum_l1 = 0; 808 for k = 1:(constrained_periods - t + 1) 809 l1 = constraint_index{k}; 810 l1 = find(constrained_perfect_foresight(l1) | (k == 1)); 811 if constraint_index{k}( j1) 812 J(cum_index_J,count_col) = d_y_x(cum_index_d_y_x(indx_x),indx_x(count_col)); 813 count_col = count_col + 1; 814 end 815 cum_l1 = cum_l1 + length(l1); 816 end 817 else 818 J(cum_index_J,count_col) = d_y_x(cum_index_d_y_x(indx_x),indx_x(count_col)); 819 count_col = count_col + 1; 820 end 821 cum_l1 = cum_l1 + length(constrained_vars_t(j1)); 822 end 823 824 825 % % Numerical computation of the derivatives in the second systme 826 % J1 = []; 827 % col_count = 1; 828 % for j = constraint_index_t 829 % j_pos = controlled_varexo(j); 830 % if constrained_perfect_foresight(j) 831 % for time = time_index_constraint 832 % saved = oo_.exo_simul(time,j_pos); 833 % oo_.exo_simul(time,j_pos) = oo_.exo_simul(time,j_pos) + eps1; 834 % simul(); 835 % J1(:,col_count) = (oo_.endo_simul(indx_endo) - ys) / eps1; 836 % oo_.exo_simul(time,j_pos) = saved; 837 % col_count = col_count + 1; 838 % end 839 % else 840 % saved = oo_.exo_simul(maximum_lag+1,j_pos); 841 % oo_.exo_simul(maximum_lag+1,j_pos) = oo_.exo_simul(maximum_lag+1,j_pos) + eps1; 842 % simul(); 843 % % indx_endo 844 % J1(:,col_count) = (oo_.endo_simul(indx_endo) - ys) / eps1; 845 % % J(:,col_count) = (oo_.endo_simul((pp - 1) * M_.endo_nbr + 1: pp * M_.endo_nbr) - ys) / eps1; 846 % oo_.exo_simul(maximum_lag+1,j_pos) = saved; 847 % col_count = col_count + 1; 848 % end 849 % end 850 % disp('J1'); 851 % disp(full(J1)); 852 % sdfmlk; 853 854 855 normr = norm(r, 1); 856 857 disp(['iteration ' int2str(it) ' error = ' num2str(normr) ' at time ' int2str(t)]); 858 859 if normr <= eps 860 convg = 1; 861 disp('convergence achieved'); 862 else 863 % Newton update on exogenous shocks 864 try 865 D_exo = - J \ r; 866 catch 867 [V, D] = eig(full(J)); 868 ev = diag(D); 869 [ev abs(ev)] 870 z_root = find(abs(ev) < 1e-4); 871 z_root 872 disp(V(:,z_root)); 873 end 874 old_exo = oo_.exo_simul; 875 col_count = 1; 876 for j = constraint_index_t 877 j_pos=controlled_varexo(j); 878 if constrained_perfect_foresight(j) 879 oo_.exo_simul(time_index_constraint,j_pos) = (oo_.exo_simul(time_index_constraint,j_pos) + D_exo(col_count:(col_count + constrained_periods - t) )); 880 col_count = col_count + constrained_periods - t + 1; 881 else 882 oo_.exo_simul(maximum_lag + 1,j_pos) = oo_.exo_simul(maximum_lag + 1,j_pos) + D_exo(col_count); 883 col_count = col_count + 1; 884 end 885 end 886 end 887 it = it + 1; 888 end 889 if ~convg 890 error(['convergence not achived at time ' int2str(t) ' after ' int2str(it) ' iterations']); 891 end 892 for j = constraint_index_t 893 j_pos=controlled_varexo(j); 894 if constrained_perfect_foresight(j) 895 % in case of mixed surprise and perfect foresight on the 896 % endogenous path, at each date all the exogenous paths have to be 897 % stored. The paths are stacked in exo. 898 for time = time_index_constraint; 899 exo(past_val + time,j_pos) = oo_.exo_simul(time,j_pos); 900 end 901 else 902 exo(maximum_lag + t,j_pos) = oo_.exo_simul(maximum_lag + 1,j_pos); 903 end 904 end 905 past_val = past_val + length(time_index_constraint); 906 if t < constrained_periods 907 endo(maximum_lag + t,:) = yc; 908 else 909 endo(maximum_lag + t :maximum_lag + options_cond_fcst.periods ,:) = ycc(:,maximum_lag + 1:maximum_lag + options_cond_fcst.periods - constrained_periods + 1)'; 910 end 911 initial_conditions = yc; 912 if maximum_lag > 0 913 exo_init(1,:) = exo(maximum_lag + t,:); 914 end 915 end 916end 917options_.periods = save_options_periods; 918options_.dynatol.f = save_options_dynatol_f; 919options_.initval_file = save_options_initval_file; 920options_.verbosity = verbosity; 921oo_.endo_simul = endo'; 922oo_.exo_simul = exo; 923if direct_mode 924 data_set = dseries([endo(2:end,1:M_.orig_endo_nbr) exo(2:end,:)], total_periods(1), {plan.endo_names{:} plan.exo_names{:}}, {plan.endo_names{:} plan.exo_names{:}}); 925end 926