Commit 6328a44f authored by Stéphane Adjemian's avatar Stéphane Adjemian
Browse files

Streamlined extended paths routines.

*  Removed the  necessity (for  the  user) to  run stoch_simul  bebore
executing the exetended path routine (when options_.init>0).

*  The  value  of  options_.ep.init  defines the  mix  (used  for  the
initialization of  the perfect foresight solver)  between the previous
perfect foresight  solution and  the path obtained  with an  order one
perturbation approach.

* Removed timing related statements.

*  Changed homotopy set-up  for stochastic  extended path:  add future
multivariate innovations one by one.

* Endogeneously increase step_length in the homotopy routine.

* Removed homotopy_2 related code.
parent d86daa01
function time_series = extended_path(initial_conditions,sample_size)
% Stochastic simulation of a non linear DSGE model using the Extended Path method (Fair and Taylor 1983). A time
% series of size T is obtained by solving T perfect foresight models.
%
% series of size T is obtained by solving T perfect foresight models.
%
% INPUTS
% o initial_conditions [double] m*nlags array, where m is the number of endogenous variables in the model and
% nlags is the maximum number of lags.
% o sample_size [integer] scalar, size of the sample to be simulated.
%
%
% OUTPUTS
% o time_series [double] m*sample_size array, the simulations.
%
%
% ALGORITHM
%
%
% SPECIAL REQUIREMENTS
% Copyright (C) 2009, 2010, 2011 Dynare Team
......@@ -58,16 +58,14 @@ options_.stack_solve_algo = options_.ep.stack_solve_algo;
%
% REMARK. It is assumed that the user did run the same mod file with stoch_simul(order=1) and save
% all the globals in a mat file called linear_reduced_form.mat;
if options_.ep.init
lrf = load('linear_reduced_form','oo_');
oo_.dr = lrf.oo_.dr; clear('lrf');
if options_.ep.init==2
lambda = .8;
end
options_.order = 1;
[dr,Info,M_,options_,oo_] = resol(1,M_,options_,oo_);
end
% Do not use a minimal number of perdiods for the perfect foresight solver (with bytecode and blocks)
options_.minimal_solving_period = options_.ep.periods;
options_.minimal_solving_period = 100;%options_.ep.periods;
% Get indices of variables with non zero steady state
idx = find(abs(oo_.steady_state)>1e-6);
......@@ -83,12 +81,12 @@ make_y_;
time_series = zeros(M_.endo_nbr,sample_size);
% Set the covariance matrix of the structural innovations.
variances = diag(M_.Sigma_e);
positive_var_indx = find(variances>0);
variances = diag(M_.Sigma_e);
positive_var_indx = find(variances>0);
effective_number_of_shocks = length(positive_var_indx);
stdd = sqrt(variances(positive_var_indx));
covariance_matrix = M_.Sigma_e(positive_var_indx,positive_var_indx);
covariance_matrix_upper_cholesky = chol(covariance_matrix);
covariance_matrix = M_.Sigma_e(positive_var_indx,positive_var_indx);
covariance_matrix_upper_cholesky = chol(covariance_matrix);
% Set seed.
if options_.ep.set_dynare_seed_to_default
......@@ -98,11 +96,11 @@ end
% Simulate shocks.
switch options_.ep.innovation_distribution
case 'gaussian'
oo_.ep.shocks = randn(sample_size,effective_number_of_shocks)*covariance_matrix_upper_cholesky;
oo_.ep.shocks = randn(sample_size,effective_number_of_shocks)*covariance_matrix_upper_cholesky;
otherwise
error(['extended_path:: ' options_.ep.innovation_distribution ' distribution for the structural innovations is not (yet) implemented!'])
end
% Set future shocks (Stochastic Extended Path approach)
if options_.ep.stochastic.status
switch options_.ep.stochastic.method
......@@ -186,13 +184,12 @@ while (t<sample_size)
oo_.exo_simul(2+options_.ep.stochastic.order+1:2+options_.ep.stochastic.order+options_.ep.stochastic.scramble,positive_var_indx) = ...
randn(options_.ep.stochastic.scramble,effective_number_of_shocks)*covariance_matrix_upper_cholesky;
end
if options_.ep.init% Compute first order solution...
initial_path = simult_(initial_conditions,oo_.dr,oo_.exo_simul(2:end,:),1);
if options_.ep.init==1
oo_.endo_simul(:,1:end-1) = initial_path(:,1:end-1);% Last column is the steady state.
elseif options_.ep.init==2
oo_.endo_simul(:,1:end-1) = initial_path(:,1:end-1)*lambda+oo_.endo_simul(:,1:end-1)*(1-lambda);
end
if options_.ep.init% Compute first order solution (Perturbation)...
ex = zeros(size(oo_.endo_simul,2),size(oo_.exo_simul,2));
ex(1:size(oo_.exo_simul,1),:) = oo_.exo_simul;
oo_.exo_simul = ex;
initial_path = simult_(initial_conditions,dr,oo_.exo_simul(2:end,:),1);
oo_.endo_simul(:,1:end-1) = initial_path(:,1:end-1)*options_.ep.init+oo_.endo_simul(:,1:end-1)*(1-options_.ep.init);
end
% Solve a perfect foresight model (using bytecoded version).
increase_periods = 0;
......@@ -200,10 +197,8 @@ while (t<sample_size)
while 1
if ~increase_periods
t0 = tic;
[flag,tmp] = bytecode('dynamic');
ctime = toc(t0);
[flag,tmp] = bytecode('dynamic');
info.convergence = ~flag;
info.time = ctime;
end
if verbosity
if info.convergence
......@@ -240,7 +235,7 @@ while (t<sample_size)
break
else
options_.periods = options_.periods + options_.ep.step;
options_.minimal_solving_period = options_.periods;
%options_.minimal_solving_period = 100;%options_.periods;
increase_periods = increase_periods + 1;
if verbosity
if t<10
......@@ -263,13 +258,11 @@ while (t<sample_size)
end
t0 = tic;
[flag,tmp] = bytecode('dynamic');
ctime = toc(t0);
info.time = info.time+ctime;
if info.convergence
maxdiff = max(max(abs(tmp(:,2:options_.ep.fp)-tmp_old(:,2:options_.ep.fp))));
if maxdiff<options_.dynatol.x
options_.periods = options_.ep.periods;
options_.minimal_solving_period = options_.periods;
%options_.minimal_solving_period = 100;%options_.periods;
oo_.exo_simul = oo_.exo_simul(1:(options_.periods+2),:);
break
end
......@@ -297,12 +290,19 @@ while (t<sample_size)
end
end
if ~info.convergence% If the previous step was unsuccesfull, use an homotopic approach
[INFO,tmp] = homotopic_steps(.5,.01,t);
% Cumulate time.
info.time = ctime+INFO.time;
[INFO,tmp] = homotopic_steps(.5,.01);
if (~isstruct(INFO) && isnan(INFO)) || ~INFO.convergence
disp('Homotopy:: No convergence of the perfect foresight model solver!')
error('I am not able to simulate this model!');
[INFO,tmp] = homotopic_steps(0,.01);
if ~INFO.convergence
disp('Homotopy:: No convergence of the perfect foresight model solver!')
error('I am not able to simulate this model!');
else
info.convergence = 1;
oo_.endo_simul = tmp;
if verbosity && info.convergence
disp('Homotopy:: Convergence of the perfect foresight model solver!')
end
end
else
info.convergence = 1;
oo_.endo_simul = tmp;
......@@ -316,14 +316,10 @@ while (t<sample_size)
% Save results of the perfect foresight model solver.
if options_.ep.memory
mArray1(:,:,s,t) = oo_.endo_simul(:,1:100);
mArrat2(:,:,s,t) = transpose(oo_.exo_simul(1:100,:));
mArrat2(:,:,s,t) = transpose(oo_.exo_simul(1:100,:));
end
time_series(:,t) = time_series(:,t)+ www(s)*oo_.endo_simul(:,2);
%save('simulated_paths.mat','time_series');
% Set initial condition for the nex round.
%initial_conditions = oo_.endo_simul(:,2);
end
%oo_.endo_simul = oo_.endo_simul(:,1:options_.periods+M_.maximum_endo_lag+M_.maximum_endo_lead);
oo_.endo_simul(:,1:end-1) = oo_.endo_simul(:,2:end);
oo_.endo_simul(:,1) = time_series(:,t);
oo_.endo_simul(:,end) = oo_.steady_state;
......
function [info,tmp] = homotopic_steps(initial_weight,step_length,time)
function [info,tmp] = homotopic_steps(initial_weight,step_length)
global oo_ options_ M_
% Set increase and decrease factors.
increase_factor = 5.0;
decrease_factor = 0.2;
% Save current state of oo_.endo_simul and oo_.exo_simul.
endo_simul = oo_.endo_simul;
exxo_simul = oo_.exo_simul;
% Reset exo_simul to zero.
oo_.exo_simul = zeros(size(oo_.exo_simul));
[idx,jdx] = find(abs(exxo_simul)>1e-12);
idx = unique(idx);
if idx(1)-2
% The first scalar in idx must be equal to two.
error('extended_path::homotopy: Something is wrong in oo_.exo_simul!')
end
stochastic_extended_path_depth = length(idx);
if stochastic_extended_path_depth>1
for i=2:length(idx)
if (idx(i)-idx(i-1)-1)
error('extended_path::homotopy: Something is wrong in oo_.exo_simul!')
end
end
end
initial_step_length = step_length;
max_iter = 1000/step_length;
weight = initial_weight;
verbose = 1;
iter = 0;
ctime = 0;
verbose = options_.ep.debug;
reduce_step_flag = 0;
......@@ -22,180 +38,100 @@ if verbose
format long
end
homotopy_1 = 1; % Only innovations are rescaled. Starting from weight equal to initial_weight.
homotopy_2 = 0; % Only innovations are rescaled. Starting from weight equal to zero.
disp(' ')
if homotopy_1
for d=1:stochastic_extended_path_depth
% (re)Set iter.
iter = 0;
% (re)Set iter.
jter = 0;
% (re)Set weight.
weight = initial_weight;
% (re)Set exo_simul to zero.
oo_.exo_simul = zeros(size(oo_.exo_simul));
while weight<1
iter = iter+1;
oo_.exo_simul = weight*exxo_simul;
t0 = tic;
oo_.exo_simul(idx(d),:) = weight*exxo_simul(idx(d),:);
if d>1
oo_.exo_simul(1:idx(d-1),:) = exxo_simul(1:idx(d-1),:);
end
[flag,tmp] = bytecode('dynamic');
TeaTime = toc(t0);
ctime = ctime+TeaTime;
info.convergence = ~flag;
info.convergence = ~flag;% Equal to one if the perfect foresight solver converged for the current value of weight.
if verbose
if ~info.convergence
disp(['Iteration n° ' int2str(iter) ', weight is ' num2str(weight,8) ', Convergence problem!' ])
else
if info.convergence
disp(['Iteration n° ' int2str(iter) ', weight is ' num2str(weight,8) ', Ok!' ])
else
disp(['Iteration n° ' int2str(iter) ', weight is ' num2str(weight,8) ', Convergence problem!' ])
end
end
if ~info.convergence
if abs(weight-initial_weight)<1e-12% First iterations.
if info.convergence
if d<stochastic_extended_path_depth
oo_.endo_simul = tmp;
end
jter = jter + 1;
if jter>3
if verbose
disp('I am increasing the step length!')
end
step_length=step_length*increase_factor;
jter = 0;
end
if abs(1-weight)<options_.dynatol.x;
break
end
weight = weight+step_length;
else% Perfect foresight solver failed for the current value of weight.
if initial_weight>0 && abs(weight-initial_weight)<1e-12% First iteration, the initial weight is too high.
if verbose
disp('I am reducing the initial weight!')
end
initial_weight = initial_weight/1.1;
initial_weight = initial_weight/2;
weight = initial_weight;
if weight<1e-12
homotopy_1 = 0;
homotopy_2 = 1;
break
oo_.endo_simul = endo_simul;
oo_.exo_simul = exxo_simul;
info.convergence = 0;
info.depth = d;
tmp = [];
return
end
continue
else% A good initial weight has been obtained. In case of convergence problem we have to reduce the step length.
else% Initial weight is OK, but the perfect foresight solver failed on some subsequent iteration.
if verbose
disp('I am reducing the step length!')
end
weight = weight-step_length;
step_length=step_length/10;
jter = 0;
if weight>0
weight = weight-step_length;
end
step_length=step_length*decrease_factor;
weight = weight+step_length;
if 10*step_length<options_.dynatol.x
homotopy_1 = 0;
homotopy_2 = 0;
if step_length<options_.dynatol.x
break
end
continue
end
else
oo_.endo_simul = tmp;
info.time = ctime;
if abs(1-weight)<=1e-12;
homotopy_1 = 0;
homotopy_2 = 0;
break
end
weight = weight+step_length;
end
if iter>max_iter
info = NaN;
return
end
end
if weight<1 && homotopy_1
if weight<1
oo_.exo_simul = exxo_simul;
t0 = tic;
[flag,tmp] = bytecode('dynamic');
TeaTime = toc(t0);
ctime = ctime+TeaTime;
info.convergence = ~flag;
info.time = ctime;
if info.convergence
oo_.endo_simul = tmp;
homotopy_1 = 0;
homotopy_2 = 0;
return
else
if step_length>1e-12
if verbose
disp('I am reducing step length!')
end
step_length=step_length/2;
else
weight = initial_weight;
step_length = initial_step_length;
info = NaN;
homotopy_2 = 1;
homotopy_1 = 0;
end
end
end
end
iter = 0;
weight = 0;
if homotopy_2
while weight<1
iter = iter+1;
oo_.exo_simul(2,:) = weight*exxo_simul(2,:);
if time==1
oo_.endo_simul = repmat(oo_.steady_state,1,size(oo_.endo_simul,2));
else
oo_.endo_simul = endo_simul;
end
t0 = tic;
[flag,tmp] = bytecode('dynamic');
TeaTime = toc(t0);
ctime = ctime+TeaTime;
old_weight = weight;
info.convergence = ~flag;
if verbose
if ~info.convergence
disp(['Iteration n° ' int2str(iter) ', weight is ' num2str(old_weight,8) ', Convergence problem!' ])
else
disp(['Iteration n° ' int2str(iter) ', weight is ' num2str(old_weight,8) ', Ok!' ])
end
end
if ~info.convergence
if iter==1
disp('I am not able to simulate this model!')
disp('There is something wrong with the initial condition of the homotopic')
disp('approach...')
error(' ')
else
if verbose
disp('I am reducing the step length!')
end
step_length=step_length/10;
if 10*step_length<options_.dynatol.x
homotopy_1 = 0;
homotopy_2 = 0;
break
end
weight = old_weight+step_length;
end
else
oo_.endo_simul = tmp;
info.time = ctime;
if abs(1-weight)<=1e-12;
homotopy_2 = 1;
break
end
weight = weight+step_length;
step_length = initial_step_length;
end
if iter>max_iter
info = NaN;
return
end
end
if weight<1 && homotopy_2
oo_.exo_simul(2,:) = exxo_simul(2,:);
t0 = tic;
[flag,tmp] = bytecode('dynamic');
TeaTime = toc(t0);
ctime = ctime+TeaTime;
info.convergence = ~flag;
info.time = ctime;
if info.convergence
oo_.endo_simul = tmp;
homotopy_1 = 0;
else
if step_length>1e-12
if verbose
disp('I am reducing step length!')
end
step_length=step_length/2;
if step_length>options_.dynatol.x
oo_.endo_simul = endo_simul;
oo_.exo_simul = exxo_simul;
info.convergence = 0;
info.depth = d;
tmp = [];
return
else
weight = initial_weight;
step_length = initial_step_length;
info = NaN;
homotopy_2 = 1;
homotopy_1 = 0;
error('extended_path::homotopy: Oups! I did my best, but I am not able to simulate this model...')
end
end
end
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment