diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index 814af6c69..7155eb0e1 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -30,8 +30,7 @@ function global_initialization()
global oo_ M_ options_ ct_
- ct_=0;
-
+ options_.terminal_condition = 0;
options_.rplottype = 0;
options_.smpl = 0;
options_.dynatol = 0.00001;
diff --git a/matlab/perfect_foresight_simulation.m b/matlab/perfect_foresight_simulation.m
index b776efd1b..047c91b59 100644
--- a/matlab/perfect_foresight_simulation.m
+++ b/matlab/perfect_foresight_simulation.m
@@ -1,9 +1,11 @@
-function info = perfect_foresight_simulation(init)
- % performs deterministic simulations with lead or lag on one period
+function [info,endo_simul] = perfect_foresight_simulation(endo_simul,exo_simul,compute_linear_solution,steady_state)
+ % Performs deterministic simulations with lead or lag on one period
%
% INPUTS
- % none
- %
+ % endo_simul [double] n*T matrix, where n is the number of endogenous variables.
+ % exo_simul [double] q*T matrix, where q is the number of shocks.
+ % compute_linear_solution [integer] scalar equal to zero or one.
+ %
% OUTPUTS
% none
%
@@ -33,40 +35,44 @@ function info = perfect_foresight_simulation(init)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see .
- global M_ options_ oo_
- global ct_ it_
+ global M_ options_ it_
persistent flag_init
persistent lead_lag_incidence dynamic_model ny nyp nyf nrs nrc iyf iyp isp is isf isf1 iz icf
-
- if nargin==1
- flag_init = [];
- end
+ persistent ghx
if isempty(flag_init)
lead_lag_incidence = M_.lead_lag_incidence;
- dynamic_model = [M_.fname '_dynamic'];
- ny = size(oo_.endo_simul,1);
- nyp = nnz(lead_lag_incidence(1,:));
- nyf = nnz(lead_lag_incidence(3,:));
+ dynamic_model = [M_.fname '_dynamic'];
+ ny = size(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;
nrc = nyf+1;
- iyf = find(lead_lag_incidence(3,:)>0);
- iyp = find(lead_lag_incidence(1,:)>0);
- isp = 1:nyp;
- is = (nyp+1):(ny+nyp);
+ iyf = find(lead_lag_incidence(3,:)>0);% indices for leaded variables.
+ iyp = find(lead_lag_incidence(1,:)>0);% indices for lagged variables.
+ isp = 1:nyp;
+ is = (nyp+1):(nyp+ny); % Indices for contemporaneaous variables.
isf = iyf+nyp;
isf1 = (nyp+ny+1):(nyf+nyp+ny+1);
- iz = 1:(ny+nyp+nyf);
- icf = 1:size(iyf,2);
- flag_init = 1;
- if nargin==1
- return
- end
- end
+ iz = 1:(ny+nyp+nyf);
+ icf = 1:size(iyf,2);
+ flag_init = 1;
+ end
- endo_simul = oo_.endo_simul;
- periods = options_.periods;
+ if nargin<3
+ compute_linear_solution = 0;
+ if nargin<4
+ error('The steady state (fourth input argument) is missing!');
+ end
+ end
+
+ if compute_linear_solution
+ [dr,info]=resol(steady_state,0);
+ ghx = dr.ghx(end-dr.nfwrd+1:end,:);
+ end
+
+ periods = options_.periods;
stop = 0 ;
it_init = M_.maximum_lag+1;
@@ -80,63 +86,80 @@ function info = perfect_foresight_simulation(init)
last_line = options_.maxit_;
error_growth = 0;
- h1 = clock;
+ h1 = clock;
for iter = 1:options_.maxit_
- h2 = clock;
- if ct_ == 0
- c = zeros(ny*periods,nrc);
- else
+ h2 = clock;
+ if options_.terminal_condition
c = zeros(ny*(periods+1),nrc);
+ else
+ c = zeros(ny*periods,nrc);
end
- it_ = it_init ;
+ it_ = it_init;
z = [ endo_simul(iyp,it_-1) ; endo_simul(:,it_) ; endo_simul(iyf,it_+1) ];
- [d1,jacobian] = feval(dynamic_model,z,oo_.exo_simul, M_.params, it_);
+ [d1,jacobian] = feval(dynamic_model,z,exo_simul, M_.params, it_);
jacobian = [jacobian(:,iz) , -d1];
- ic = 1:ny;
- icp = iyp;
+ ic = 1:ny;
+ icp = iyp;
c(ic,:) = jacobian(:,is)\jacobian(:,isf1) ;
- for it_ = it_init+(1:periods-1)
+ 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,oo_.exo_simul, M_.params, it_);
+ [d1,jacobian] = feval(dynamic_model,z,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;
- c(ic,:) = jacobian(:,is)\jacobian(:,isf1);
- end
- if ct_ == 1
- s = eye(ny);
- s(:,isf) = s(:,isf)+c(ic,1:nyf);
ic = ic + ny;
- c(ic,nrc) = s\c(:,nrc);
+ icp = icp + ny;
+ c(ic,:) = jacobian(:,is)\jacobian(:,isf1);
+ end
+ if options_.terminal_condition
+ if options_.terminal_condition==1% Terminal condition is Y_{T} = Y_{T+1}
+ s = eye(ny);
+ 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_);
+ jacobian = [jacobian(:,iz) -d1];
+ jacobian(:,[isf nrs]) = jacobian(:,[isf nrs])-jacobian(:,isp)*c(icp,:) ;
+ ic = ic + ny;
+ icp = icp + ny;
+ s = jacobian(:,is);
+ s(:,iyp-nyp) = s(:,iyp-nyp)+jacobian(:,isf)*ghx;
+ c (ic,:) = s\jacobian(:,isf1);
+ end
c = bksup0(c,ny,nrc,iyf,icf,periods);
c = reshape(c,ny,periods+1);
endo_simul(:,it_init+(0:periods)) = endo_simul(:,it_init+(0:periods))+options_.slowc*c;
- else
- c = bksup0(c,ny,nrc,iyf,icf,periods);
- c = reshape(c,ny,periods);
+ 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;
- end
+ 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;
+ % 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);
+ info.error = err;
+ info.iterations.time = info.iterations.time(1:iter);
+ info.iterations.error = info.iterations.error(1:iter);
break
end
- if err < options_.dynatol
- stop = 1;
- info.time = etime(clock,h1);
- info.error = err;
- info.iterations.time = info.iterations.time(1:iter);
- info.iterations.error = info.iterations.error(1:iter);
- oo_.endo_simul = endo_simul;
- break
- end
- end
+ end
+
+ if options_.terminal_condition==2
+ distance_to_steady_state = abs(((endo_simul(:,end-1)-endo_simul(:,end))./endo_simul(:,end)))*100;
+ disp('Distance to steady state at the end is (in percentage):')
+ distance_to_steady_state
+ end
if ~stop
info.time = etime(clock,h1);
@@ -144,4 +167,6 @@ function info = perfect_foresight_simulation(init)
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 = [ ];
end
\ No newline at end of file
diff --git a/matlab/sim1.m b/matlab/sim1.m
index 3190f3137..70f4583ee 100644
--- a/matlab/sim1.m
+++ b/matlab/sim1.m
@@ -33,7 +33,7 @@ function sim1
% along with Dynare. If not, see .
global M_ options_ oo_
-global iyp iyf ct_ M_ it_ c
+global iyp iyf M_ it_ c
lead_lag_incidence = M_.lead_lag_incidence;
@@ -61,7 +61,7 @@ h1 = clock ;
for iter = 1:options_.maxit_
h2 = clock ;
- if ct_ == 0
+ if options_.terminal_condition == 0
c = zeros(ny*options_.periods,nrc) ;
else
c = zeros(ny*(options_.periods+1),nrc) ;
@@ -84,7 +84,7 @@ for iter = 1:options_.maxit_
c (ic,:) = jacobian(:,is)\jacobian(:,isf1) ;
end
- if ct_ == 1
+ if options_.terminal_condition == 1
s = eye(ny) ;
s(:,isf) = s(:,isf)+c(ic,1:nyf) ;
ic = ic + ny ;
diff --git a/matlab/simk.m b/matlab/simk.m
index 0982e9a8f..1767aaa1a 100644
--- a/matlab/simk.m
+++ b/matlab/simk.m
@@ -35,7 +35,7 @@ function simk
% along with Dynare. If not, see .
global M_ options_ oo_
-global it_ iyr0 ct_ broyden_
+global it_ iyr0 broyden_
%func_name = [M_.fname '_static'];
nk = M_.maximum_endo_lag + M_.maximum_endo_lead + 1 ;
@@ -283,7 +283,7 @@ for iter = 1:options_.maxit_
iyr = iyr + ny ;
icr0 = icr0 + ny ;
end
- if ct_ == 1
+ if options_.terminal_condition == 1
ofs = (((it_-M_.maximum_lag-2)*ny+1)-1)*ncc*8 ;
junk = fseek(fid,ofs,-1) ;
@@ -303,7 +303,7 @@ for iter = 1:options_.maxit_
end
end
oo_.endo_simul = reshape(oo_.endo_simul,ny,options_.periods+M_.maximum_lag+M_.maximum_endo_lead) ;
- if ct_ == 1
+ if options_.terminal_condition == 1
hbacsup = clock ;
c = bksupk(ny,fid,ncc,icc1) ;
hbacsup = etime(clock,hbacsup) ;