diff --git a/matlab/ep/extended_path.m b/matlab/ep/extended_path.m
index 806880e87390b7101cb87bfc602e0ab6995be409..4335e9492d140fdf31efb93ab3be31e581218cdc 100644
--- a/matlab/ep/extended_path.m
+++ b/matlab/ep/extended_path.m
@@ -1,17 +1,17 @@
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;
diff --git a/matlab/ep/homotopic_steps.m b/matlab/ep/homotopic_steps.m
index 696b88014b97375f8c9e5a6f6387b8c905e03437..3801e003e417c5b50aa6d760d5f6ea3dd1838b5b 100644
--- a/matlab/ep/homotopic_steps.m
+++ b/matlab/ep/homotopic_steps.m
@@ -1,20 +1,36 @@
-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