stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

fixing extended-path

time-shift
Michel Juillard 2015-05-25 09:43:21 +02:00
parent b3047c9742
commit a36fd53eff
5 changed files with 12 additions and 303 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
matlab/ep/extended_path.m
View File

@ -118,10 +118,12 @@ replic_nbr = ep.replic_nbr;
switch ep.innovation_distribution
case 'gaussian'
shocks = transpose(transpose(covariance_matrix_upper_cholesky)* ...
randn(effective_number_of_shocks,sample_size*replic_nbr));
randn(effective_number_of_shocks,sample_size* ...
replic_nbr));
shocks(:,positive_var_indx) = shocks;
case 'calibrated'
replic_nbr = 1;
shocks = zeros(sample_size,effective_number_of_shocks);
shocks = zeros(sample_size,M_.exo_nbr);
for i = 1:length(M_.unanticipated_det_shocks)
k = M_.unanticipated_det_shocks(i).periods;
ivar = M_.unanticipated_det_shocks(i).exo_id;
@ -132,7 +134,6 @@ switch ep.innovation_distribution
socks(k,ivar) = shocks(k,ivar) * v;
end
end
shocks = shocks(:,positive_var_indx);
otherwise
error(['extended_path:: ' ep.innovation_distribution ' distribution for the structural innovations is not (yet) implemented!'])
end
@ -203,7 +204,7 @@ while (t <= sample_size)
parfor k = 1:replic_nbr
exo_simul = repmat(oo_.exo_steady_state',periods+2,1);
% exo_simul(1:sample_size+3-t,:) = exo_simul_(t:end,:);
exo_simul(2,positive_var_indx) = exo_simul_(M_.maximum_lag+t,positive_var_indx) + ...
exo_simul(2,:) = exo_simul_(M_.maximum_lag+t,:) + ...
shocks((t-2)*replic_nbr+k,:);
initial_conditions = results{k}(:,t-1);
results{k}(:,t) = extended_path_core(ep.periods,endo_nbr,exo_nbr,positive_var_indx, ...
@ -219,8 +220,8 @@ while (t <= sample_size)
maximum_lead,1);
% exo_simul(1:sample_size+maximum_lag+maximum_lead-t+1,:) = ...
% exo_simul_(t:end,:);
exo_simul(maximum_lag+1,positive_var_indx) = ...
exo_simul_(maximum_lag+t,positive_var_indx) + shocks((t-2)*replic_nbr+k,:);
exo_simul(maximum_lag+1,:) = ...
exo_simul_(maximum_lag+t,:) + shocks((t-2)*replic_nbr+k,:);
initial_conditions = results{k}(:,t-1);
results{k}(:,t) = extended_path_core(ep.periods,endo_nbr,exo_nbr,positive_var_indx, ...
exo_simul,ep.init,initial_conditions,...

4
matlab/ep/setup_stochastic_perfect_foresight_model_solver.m
View File

@ -20,7 +20,9 @@ function pfm = setup_stochastic_perfect_foresight_model_solver(DynareModel,Dynar
pfm.lead_lag_incidence = DynareModel.lead_lag_incidence;
pfm.ny = DynareModel.endo_nbr;
pfm.Sigma = DynareModel.Sigma_e;
pfm.Omega = chol(pfm.Sigma,'upper'); % Sigma = Omega'*Omega
if det(pfm.Sigma) > 0
pfm.Omega = chol(pfm.Sigma,'upper'); % Sigma = Omega'*Omega
end
pfm.number_of_shocks = length(pfm.Sigma);
pfm.stochastic_order = DynareOptions.ep.stochastic.order;
pfm.max_lag = DynareModel.maximum_endo_lag;

4
matlab/perfect-foresight-models/sim1.m
View File

@ -93,8 +93,8 @@ o_periods = periods;
ZERO = zeros(length(i_upd),1);
h1 = clock ;
iA = zeros(periods*M_.NNZDerivatives(1),3);
for iter = 1:options_.simul.maxit
iA = zeros(periods*M.NNZDerivatives(1),3);
for iter = 1:options.simul.maxit
h2 = clock ;
i_rows = (1:ny)';

160
matlab/perfect_foresight_solver.m
View File

@ -1,160 +0,0 @@
function perfect_foresight_solver()
% Computes deterministic simulations
%
% INPUTS
% None
%
% OUTPUTS
% none
%
% ALGORITHM
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 1996-2014 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 <http://www.gnu.org/licenses/>.
global M_ options_ oo_
if options_.stack_solve_algo < 0 || options_.stack_solve_algo > 7
error('PERFECT_FORESIGHT_SOLVER: stack_solve_algo must be between 0 and 7')
end
if ~options_.block && ~options_.bytecode && options_.stack_solve_algo ~= 0 ...
&& options_.stack_solve_algo ~= 6 && options_.stack_solve_algo ~= 7
error('PERFECT_FORESIGHT_SOLVER: you must use stack_solve_algo=0 or stack_solve_algo=6 when not using block nor bytecode option')
end
if options_.block && ~options_.bytecode && options_.stack_solve_algo == 5
error('PERFECT_FORESIGHT_SOLVER: you can''t use stack_solve_algo = 5 without bytecode option')
end
if (options_.block || options_.bytecode) && options_.stack_solve_algo == 6
error('PERFECT_FORESIGHT_SOLVER: you can''t use stack_solve_algo = 6 with block or bytecode option')
end
if isoctave && options_.stack_solve_algo == 2
error('PERFECT_FORESIGHT_SOLVER: you can''t use stack_solve_algo = 2 under Octave')
end
if isempty(oo_.endo_simul) || any(size(oo_.endo_simul) ~= [ M_.endo_nbr, M_.maximum_lag+options_.periods+M_.maximum_lead ])
error('PERFECT_FORESIGHT_SOLVER: ''oo_.endo_simul'' has wrong size. Did you run ''perfect_foresight_setup'' ?')
end
if isempty(oo_.exo_simul) || any(size(oo_.exo_simul) ~= [ M_.maximum_lag+options_.periods+M_.maximum_lead, M_.exo_nbr ])
error('PERFECT_FORESIGHT_SOLVER: ''oo_.exo_simul'' has wrong size. Did you run ''perfect_foresight_setup'' ?')
end
if isempty(options_.scalv) || options_.scalv == 0
options_.scalv = oo_.steady_state;
end
options_.scalv= 1;
if options_.debug
model_static = str2func([M_.fname,'_static']);
for ii=1:size(oo_.exo_simul,1)
[residual(:,ii)] = model_static(oo_.steady_state, oo_.exo_simul(ii,:),M_.params);
end
problematic_periods=find(any(isinf(residual)) | any(isnan(residual)))-M_.maximum_endo_lag;
if ~isempty(problematic_periods)
period_string=num2str(problematic_periods(1));
for ii=2:length(problematic_periods)
period_string=[period_string, ', ', num2str(problematic_periods(ii))];
end
fprintf('\n\nWARNING: Value for the exogenous variable(s) in period(s) %s inconsistent with the static model.\n',period_string);
fprintf('WARNING: Check for division by 0.\n')
end
end
% Effectively compute simulation, possibly with homotopy
if options_.no_homotopy
[oo_.endo_simul,oo_.deterministic_simulation.status] = perfect_foresight_solver_core(M_,oo_,options_);
else
exosim = oo_.exo_simul;
exoinit = repmat(oo_.exo_steady_state',M_.maximum_lag+options_.periods+M_.maximum_lead,1);
endosim = oo_.endo_simul;
endoinit = repmat(oo_.steady_state, 1,M_.maximum_lag+options_.periods+M_.maximum_lead);
current_weight = 0; % Current weight of target point in convex combination
step = 1;
success_counter = 0;
while (step > options_.dynatol.x)
new_weight = current_weight + step; % Try this weight, and see if it succeeds
if new_weight >= 1
new_weight = 1; % Don't go beyond target point
step = new_weight - current_weight;
end
% Compute convex combination for exo path and initial/terminal endo conditions
% But take care of not overwriting the computed part of oo_.endo_simul
oo_.exo_simul = exosim*new_weight + exoinit*(1-new_weight);
endocombi = endosim*new_weight + endoinit*(1-new_weight);
oo_.endo_simul(:,1:M_.maximum_endo_lag) = endocombi(:,1:M_.maximum_endo_lag);
oo_.endo_simul(:,(end-M_.maximum_endo_lead):end) = endocombi(:,(end-M_.maximum_endo_lead):end);
saved_endo_simul = oo_.endo_simul;
[oo_.endo_simul,oo_.deterministic_simulation.status] = perfect_foresight_solver_core(M_,oo_,options_);
if oo_.deterministic_simulation.status == 1
current_weight = new_weight;
if current_weight >= 1
break
end
success_counter = success_counter + 1;
if success_counter >= 3
success_counter = 0;
step = step * 2;
disp([ 'Homotopy step succeeded, doubling step size (completed ' sprintf('%.1f', current_weight*100) '%, step size ' sprintf('%.3g', step) ')' ])
else
disp([ 'Homotopy step succeeded (completed ' sprintf('%.1f', current_weight*100) '%, step size ' sprintf('%.3g', step) ')' ])
end
else
oo_.endo_simul = saved_endo_simul;
success_counter = 0;
step = step / 2;
disp([ 'Homotopy step failed, halving step size (completed ' sprintf('%.1f', current_weight*100) '%, step size ' sprintf('%.3g', step) ')' ])
end
end
end
if oo_.deterministic_simulation.status == 1
disp('Perfect foresight solution found.')
else
warning('Failed to solve perfect foresight model')
end
dyn2vec;
if isnan(options_.initial_period)
initial_period = dates(1,1);
else
initial_period = options_.initial_period;
end
ts = dseries(transpose(oo_.endo_simul),initial_period,cellstr(M_.endo_names));
assignin('base', 'Simulated_time_series', ts);
end

134
matlab/perfect_foresight_solver_core.m
View File

@ -1,134 +0,0 @@
function [endo_simul, status] = perfect_foresight_solver_core(M,oo,options)
% Core function to compute deterministic simulations
%
% INPUTS
% M: model structure
% oo: output structure
% options: options structure
%
% OUTPUTS
% endo_simul: matrix endogenous variables
% deterministic_simulation: simulation status
%
% ALGORITHM
%
% various
%
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 1996-2015 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 <http://www.gnu.org/licenses/>.
endo_simul = oo.endo_simul;
status = 0;
deterministic_simulation = struct();
if(options.block)
if(options.bytecode)
[info, endo_simul] = bytecode('dynamic');
if info == 1
status = 0;
else
status = 1;
end
mexErrCheck('bytecode', info);
else
eval([M.fname '_dynamic']);
end
else
if(options.bytecode)
[info, endo_simul]=bytecode('dynamic');
if info == 1
status = 0;
else
status = 1;
end;
mexErrCheck('bytecode', info);
else
if M.maximum_endo_lead == 0
% Purely backward model
global oo_ options_
oo_ = oo;
options_ = options;
sim1_purely_backward;
endo_simul = oo_.endo_simul;
if oo_.deterministic_simulation.status == 1
status = 0;
end
elseif M.maximum_endo_lag == 0
% Purely forward model
global oo_ options_
oo_ = oo;
options_ = options;
sim1_purely_forward;
endo_simul = oo_.endo_simul;
if oo_.deterministic_simulation.status == 1
status = 0;
end
else
% General case
if options.stack_solve_algo == 0
oo = sim1(M,options,oo);
endo_simul = oo.endo_simul;
if oo.deterministic_simulation.status == 1
status = 0;
end
elseif options.stack_solve_algo == 6
global oo_ options_
oo_ = oo;
options_ = options;
sim1_lbj;
endo_simul = oo_.endo_simul;
if oo_.deterministic_simulation.status == 1
status = 0;
end
elseif options.stack_solve_algo == 7
periods = options.periods;
if ~isfield(options.lmmcp,'lb')
[lb,ub,pfm.eq_index] = get_complementarity_conditions(M);
options.lmmcp.lb = repmat(lb,periods,1);
options.lmmcp.ub = repmat(ub,periods,1);
end
y = endo_simul;
y0 = y(:,1);
yT = y(:,periods+2);
z = y(:,2:periods+1);
illi = M.lead_lag_incidence';
[i_cols,~,i_cols_j] = find(illi(:));
illi = illi(:,2:3);
[i_cols_J1,~,i_cols_1] = find(illi(:));
i_cols_T = nonzeros(M.lead_lag_incidence(1:2,:)');
[y,info] = dynare_solve(@perfect_foresight_problem,z(:),1, ...
str2func([M.fname '_dynamic']),y0,yT, ...
oo.exo_simul,M.params,oo.steady_state, ...
options.periods,M.endo_nbr,i_cols, ...
i_cols_J1, i_cols_1, i_cols_T, i_cols_j, ...
M.NNZDerivatives(1));
endo_simul = [y0 reshape(y,M.endo_nbr,periods) yT];
if info == 1
status = 1;
else
status = 0;
end;
end
end
end
end