diff --git a/matlab/extended_path.m b/matlab/extended_path.m
index d9ee33f84..f217c952a 100644
--- a/matlab/extended_path.m
+++ b/matlab/extended_path.m
@@ -31,14 +31,14 @@ function time_series = extended_path(initial_conditions,sample_size,init)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see .
global M_ oo_ options_
-
+
% Set default initial conditions.
if isempty(initial_conditions)
initial_conditions = repmat(oo_.steady_state,1,M_.maximum_lag);
end
-
-% Copy sample_size to periods.
-options_.periods = sample_size;
+
+% Set the number of periods for the deterministic solver.
+options_.periods = 100;
% Initialize the exogenous variables.
make_ex_;
@@ -60,8 +60,8 @@ tdx = M_.maximum_lag+1;
norme = 0;
% Set verbose option
-verbose = 0;
-
+verbose = 1;
+
for t=1:sample_size
shocks = exp(randn(1,number_of_structural_innovations)*covariance_matrix_upper_cholesky-.5*variances(positive_var_indx)');
oo_.exo_simul(tdx,positive_var_indx) = shocks;
@@ -72,7 +72,7 @@ for t=1:sample_size
info
end
if ~info.convergence
- info = homotopic_steps(tdx,positive_var_indx,shocks,norme,.2);
+ info = homotopic_steps(tdx,positive_var_indx,shocks,norme,.5);
if verbose
norme
info
@@ -88,49 +88,50 @@ for t=1:sample_size
oo_.endo_simul(:,1:end-1) = oo_.endo_simul(:,2:end);
oo_.endo_simul(:,end) = oo_.steady_state;
end
-
+
+
function info = homotopic_steps(tdx,positive_var_indx,shocks,init_weight,step)
-global oo_
-weight = init_weight;
-verbose = 0;
-iter = 0;
-time = 0;
-reduce_step = 0;
-while iter<=100 && weight<=1
- iter = iter+1;
- old_weight = weight;
- weight = weight+step;
- oo_.exo_simul(tdx,positive_var_indx) = weight*shocks+(1-weight);
- info = perfect_foresight_simulation;
- time = time+info.time;
- if verbose
- [iter,step]
- [info.iterations.time,info.iterations.error]
- end
- if ~info.convergence
+ global oo_
+ weight = init_weight;
+ verbose = 1;
+ iter = 0;
+ time = 0;
+ reduce_step = 0;
+ while iter<=100 && weight<=1
+ iter = iter+1;
+ old_weight = weight;
+ weight = weight+step;
+ oo_.exo_simul(tdx,positive_var_indx) = weight*shocks+(1-weight);
+ info = perfect_foresight_simulation;
+ time = time+info.time;
if verbose
- disp('Reduce step size!')
+ [iter,step]
+ [info.iterations.time,info.iterations.error]
end
- reduce_step = 1;
- break
- else
- if length(info.iterations.error)<5
+ if ~info.convergence
if verbose
- disp('Increase step size!')
+ disp('Reduce step size!')
+ end
+ reduce_step = 1;
+ break
+ else
+ if length(info.iterations.error)<5
+ if verbose
+ disp('Increase step size!')
+ end
+ step = step*1.5;
end
- step = step*1.5;
end
end
-end
-if reduce_step
- step=step/1.5;
- info = homotopic_steps(tdx,positive_var_indx,shocks,old_weight,step);
- time = time+info.time;
-elseif weight<1 && iter<100
- oo_.exo_simul(tdx,positive_var_indx) = shocks;
- info = perfect_foresight_simulation;
- info.time = info.time+time;
-else
- info.time = time;
-end
\ No newline at end of file
+ if reduce_step
+ step=step/1.5;
+ info = homotopic_steps(tdx,positive_var_indx,shocks,old_weight,step);
+ time = time+info.time;
+ elseif weight<1 && iter<100
+ oo_.exo_simul(tdx,positive_var_indx) = shocks;
+ info = perfect_foresight_simulation;
+ info.time = info.time+time;
+ else
+ info.time = time;
+ end
\ No newline at end of file
diff --git a/matlab/perfect_foresight_simulation.m b/matlab/perfect_foresight_simulation.m
index b00ec2c99..214e742f0 100644
--- a/matlab/perfect_foresight_simulation.m
+++ b/matlab/perfect_foresight_simulation.m
@@ -1,4 +1,4 @@
-function [info,endo_simul] = perfect_foresight_simulation(endo_simul,exo_simul,compute_linear_solution,steady_state)
+function info = perfect_foresight_simulation(compute_linear_solution,steady_state)
% Performs deterministic simulations with lead or lag on one period
%
% INPUTS
@@ -18,7 +18,7 @@ function [info,endo_simul] = perfect_foresight_simulation(endo_simul,exo_simul,c
% SPECIAL REQUIREMENTS
% None.
-% Copyright (C) 1996-2009 Dynare Team
+% Copyright (C) 2009 Dynare Team
%
% This file is part of Dynare.
%
@@ -34,17 +34,19 @@ function [info,endo_simul] = perfect_foresight_simulation(endo_simul,exo_simul,c
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see .
-
-global M_ options_ it_
-
+
+global M_ options_ it_ oo_
+
persistent flag_init
persistent lead_lag_incidence dynamic_model ny nyp nyf nrs nrc iyf iyp isp is isf isf1 iz icf
persistent ghx
+initial_endo_simul = oo_.endo_simul;
+
if isempty(flag_init)
lead_lag_incidence = M_.lead_lag_incidence;
dynamic_model = [M_.fname '_dynamic'];
- ny = size(endo_simul,1);
+ ny = size(oo_.endo_simul,1);
nyp = nnz(lead_lag_incidence(1,:));% number of lagged variables.
nyf = nnz(lead_lag_incidence(3,:));% number of leaded variables.
nrs = ny+nyp+nyf+1;
@@ -60,10 +62,11 @@ if isempty(flag_init)
flag_init = 1;
end
-if nargin<3
+if nargin<1
compute_linear_solution = 0;
- if nargin<4
- error('The steady state (fourth input argument) is missing!');
+else
+ if nargin<2
+ error('The steady state (second input argument) is missing!');
end
end
@@ -95,16 +98,16 @@ for iter = 1:options_.maxit_
c = zeros(ny*periods,nrc);
end
it_ = it_init;
- z = [ endo_simul(iyp,it_-1) ; endo_simul(:,it_) ; endo_simul(iyf,it_+1) ];
- [d1,jacobian] = feval(dynamic_model,z,exo_simul, M_.params, it_);
+ z = [ oo_.endo_simul(iyp,it_-1) ; oo_.endo_simul(:,it_) ; oo_.endo_simul(iyf,it_+1) ];
+ [d1,jacobian] = feval(dynamic_model,z,oo_.exo_simul, M_.params, it_);
jacobian = [jacobian(:,iz) , -d1];
ic = 1:ny;
icp = iyp;
c(ic,:) = jacobian(:,is)\jacobian(:,isf1) ;
for it_ = it_init+(1:periods-1)
- z = [ endo_simul(iyp,it_-1) ; endo_simul(:,it_) ; endo_simul(iyf,it_+1)];
- [d1,jacobian] = feval(dynamic_model,z,exo_simul, M_.params, it_);
- jacobian = [jacobian(:,iz) , -d1];
+ z = [ oo_.endo_simul(iyp,it_-1) ; oo_.endo_simul(:,it_) ; oo_.endo_simul(iyf,it_+1)];
+ [d1,jacobian] = feval(dynamic_model,z,oo_.exo_simul, M_.params, it_);
+ jacobian = [jacobian(:,iz) , -d1];
jacobian(:,[isf nrs]) = jacobian(:,[isf nrs])-jacobian(:,isp)*c(icp,:);
ic = ic + ny;
icp = icp + ny;
@@ -116,14 +119,18 @@ for iter = 1:options_.maxit_
s(:,isf) = s(:,isf)+c(ic,1:nyf);
ic = ic + ny;
c(ic,nrc) = s\c(ic,nrc);
- else% Terminal condition is Y_{T}-Y^{\star} = TransitionMatrix*(Y_{T+1}-Y^{\star})
- z = [ endo_simul(iyp,it_-1) ; endo_simul(:,it_) ; endo_simul(iyf,it_+1) ] ;
- [d1,jacobian] = feval(dynamic_model,z,exo_simul, M_.params, it_);
+ else% Terminal condition is Y_{T+1}-Y^{\star} = TransitionMatrix*(Y_{T}-Y^{\star})
+ z = [ oo_.endo_simul(iyp,it_-1) ; oo_.endo_simul(:,it_) ; oo_.endo_simul(iyf,it_+1) ] ;
+ [d1,jacobian] = feval(dynamic_model,z,oo_.exo_simul, M_.params, it_);
jacobian = [jacobian(:,iz) -d1];
jacobian(:,[isf nrs]) = jacobian(:,[isf nrs])-jacobian(:,isp)*c(icp,:) ;
ic = ic + ny;
icp = icp + ny;
s = jacobian(:,is);
+ iyp
+ nyp
+ isf
+ iyp-nyp
s(:,iyp-nyp) = s(:,iyp-nyp)+jacobian(:,isf)*ghx;
c (ic,:) = s\jacobian(:,isf1);
end
@@ -133,18 +140,18 @@ for iter = 1:options_.maxit_
else% Terminal condition is Y_{T}=Y^{\star}
c = bksup0(c,ny,nrc,iyf,icf,periods);
c = reshape(c,ny,periods);
- endo_simul(:,it_init+(0:periods-1)) = endo_simul(:,it_init+(0:periods-1))+options_.slowc*c;
+ oo_.endo_simul(:,it_init+(0:periods-1)) = oo_.endo_simul(:,it_init+(0:periods-1))+options_.slowc*c;
end
err = max(max(abs(c)));
info.iterations.time(iter) = etime(clock,h2);
info.iterations.error(iter) = err;
- % if iter>1
- % error_growth = error_growth + (info.iterations.error(iter)>info.iterations.error(iter-1));
- % end
- % if isnan(err) || error_growth>3
- % last_line = iter;
- % break
- % end
+ if iter>1
+ error_growth = error_growth + (info.iterations.error(iter)>info.iterations.error(iter-1));
+ end
+ if isnan(err) || error_growth>3
+ last_line = iter;
+ break
+ end
if err < options_.dynatol
stop = 1;
info.time = etime(clock,h1);
@@ -156,7 +163,7 @@ for iter = 1:options_.maxit_
end
if options_.terminal_condition==2
- distance_to_steady_state = abs(((endo_simul(:,end-1)-endo_simul(:,end))./endo_simul(:,end)))*100;
+ distance_to_steady_state = abs(((oo_.endo_simul(:,end-1)-oo_.endo_simul(:,end))./oo_.endo_simul(:,end)))*100;
disp('Distance to steady state at the end is (in percentage):')
distance_to_steady_state
end
@@ -167,6 +174,5 @@ if ~stop
info.convergence = 0;
info.iterations.time = info.iterations.time(1:last_line);
info.iterations.error = info.iterations.error(1:last_line);
- info.iterations.error
- endo_simul = [ ];
+ oo_.endo_simul = initial_endo_simul;
end
\ No newline at end of file