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

* Changed name of ct_ --> options_.terminal_condition. The default value is zero 

(the terminal condition is y_{T} = y^{\star}), other possible values
  are  1  (terminal  condition   is  y_{T+1}=y_{T})  and  2  (terminal
  condition is  y_{T+1}=TransitionMatrix*y_{T}, where TransitionMatrix
  is given by the first order approximation of the reduced form model).
* Added mode options_.terminal_condition=2 in perfect_foresight_simulation.m.



git-svn-id: https://www.dynare.org/svn/dynare/trunk@3176 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
stepan 2009-11-30 16:24:57 +00:00
parent 622a2d7ec2
commit 9a23cee31c
4 changed files with 97 additions and 73 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
matlab/global_initialization.m
View File

@ -30,8 +30,7 @@ function global_initialization()
global oo_ M_ options_ ct_ global oo_ M_ options_ ct_
ct_=0; options_.terminal_condition = 0;
options_.rplottype = 0; options_.rplottype = 0;
options_.smpl = 0; options_.smpl = 0;
options_.dynatol = 0.00001; options_.dynatol = 0.00001;

155
matlab/perfect_foresight_simulation.m
View File

@ -1,9 +1,11 @@
function info = perfect_foresight_simulation(init) 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 % Performs deterministic simulations with lead or lag on one period
% %
% INPUTS % 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 % OUTPUTS
% none % none
% %
@ -33,40 +35,44 @@ function info = perfect_foresight_simulation(init)
% You should have received a copy of the GNU General Public License % You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>. % along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global M_ options_ oo_ global M_ options_ it_
global ct_ it_
persistent flag_init persistent flag_init
persistent lead_lag_incidence dynamic_model ny nyp nyf nrs nrc iyf iyp isp is isf isf1 iz icf persistent lead_lag_incidence dynamic_model ny nyp nyf nrs nrc iyf iyp isp is isf isf1 iz icf
persistent ghx
if nargin==1
flag_init = [];
end
if isempty(flag_init) if isempty(flag_init)
lead_lag_incidence = M_.lead_lag_incidence; lead_lag_incidence = M_.lead_lag_incidence;
dynamic_model = [M_.fname '_dynamic']; dynamic_model = [M_.fname '_dynamic'];
ny = size(oo_.endo_simul,1); ny = size(endo_simul,1);
nyp = nnz(lead_lag_incidence(1,:)); nyp = nnz(lead_lag_incidence(1,:));% number of lagged variables.
nyf = nnz(lead_lag_incidence(3,:)); nyf = nnz(lead_lag_incidence(3,:));% number of leaded variables.
nrs = ny+nyp+nyf+1; nrs = ny+nyp+nyf+1;
nrc = nyf+1; nrc = nyf+1;
iyf = find(lead_lag_incidence(3,:)>0); iyf = find(lead_lag_incidence(3,:)>0);% indices for leaded variables.
iyp = find(lead_lag_incidence(1,:)>0); iyp = find(lead_lag_incidence(1,:)>0);% indices for lagged variables.
isp = 1:nyp; isp = 1:nyp;
is = (nyp+1):(ny+nyp); is = (nyp+1):(nyp+ny); % Indices for contemporaneaous variables.
isf = iyf+nyp; isf = iyf+nyp;
isf1 = (nyp+ny+1):(nyf+nyp+ny+1); isf1 = (nyp+ny+1):(nyf+nyp+ny+1);
iz = 1:(ny+nyp+nyf); iz = 1:(ny+nyp+nyf);
icf = 1:size(iyf,2); icf = 1:size(iyf,2);
flag_init = 1; flag_init = 1;
if nargin==1 end
return
end
end
endo_simul = oo_.endo_simul; if nargin<3
periods = options_.periods; 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 ; stop = 0 ;
it_init = M_.maximum_lag+1; it_init = M_.maximum_lag+1;
@ -80,63 +86,80 @@ function info = perfect_foresight_simulation(init)
last_line = options_.maxit_; last_line = options_.maxit_;
error_growth = 0; error_growth = 0;
h1 = clock; h1 = clock;
for iter = 1:options_.maxit_ for iter = 1:options_.maxit_
h2 = clock; h2 = clock;
if ct_ == 0 if options_.terminal_condition
c = zeros(ny*periods,nrc);
else
c = zeros(ny*(periods+1),nrc); c = zeros(ny*(periods+1),nrc);
else
c = zeros(ny*periods,nrc);
end end
it_ = it_init ; it_ = it_init;
z = [ endo_simul(iyp,it_-1) ; endo_simul(:,it_) ; endo_simul(iyf,it_+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 = [jacobian(:,iz) , -d1];
ic = 1:ny; ic = 1:ny;
icp = iyp; icp = iyp;
c(ic,:) = jacobian(:,is)\jacobian(:,isf1) ; 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)]; 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 = [jacobian(:,iz) , -d1];
jacobian(:,[isf nrs]) = jacobian(:,[isf nrs])-jacobian(:,isp)*c(icp,:); 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; 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 = bksup0(c,ny,nrc,iyf,icf,periods);
c = reshape(c,ny,periods+1); c = reshape(c,ny,periods+1);
endo_simul(:,it_init+(0:periods)) = endo_simul(:,it_init+(0:periods))+options_.slowc*c; endo_simul(:,it_init+(0:periods)) = endo_simul(:,it_init+(0:periods))+options_.slowc*c;
else else% Terminal condition is Y_{T}=Y^{\star}
c = bksup0(c,ny,nrc,iyf,icf,periods); c = bksup0(c,ny,nrc,iyf,icf,periods);
c = reshape(c,ny,periods); c = reshape(c,ny,periods);
endo_simul(:,it_init+(0:periods-1)) = endo_simul(:,it_init+(0:periods-1))+options_.slowc*c; endo_simul(:,it_init+(0:periods-1)) = endo_simul(:,it_init+(0:periods-1))+options_.slowc*c;
end end
err = max(max(abs(c))); err = max(max(abs(c)));
info.iterations.time(iter) = etime(clock,h2); info.iterations.time(iter) = etime(clock,h2);
info.iterations.error(iter) = err; info.iterations.error(iter) = err;
if iter>1 % if iter>1
error_growth = error_growth + (info.iterations.error(iter)>info.iterations.error(iter-1)); % error_growth = error_growth + (info.iterations.error(iter)>info.iterations.error(iter-1));
end % end
if isnan(err) || error_growth>3 % if isnan(err) || error_growth>3
last_line = iter; % 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 break
end end
if err < options_.dynatol end
stop = 1;
info.time = etime(clock,h1); if options_.terminal_condition==2
info.error = err; distance_to_steady_state = abs(((endo_simul(:,end-1)-endo_simul(:,end))./endo_simul(:,end)))*100;
info.iterations.time = info.iterations.time(1:iter); disp('Distance to steady state at the end is (in percentage):')
info.iterations.error = info.iterations.error(1:iter); distance_to_steady_state
oo_.endo_simul = endo_simul; end
break
end
end
if ~stop if ~stop
info.time = etime(clock,h1); info.time = etime(clock,h1);
@ -144,4 +167,6 @@ function info = perfect_foresight_simulation(init)
info.convergence = 0; info.convergence = 0;
info.iterations.time = info.iterations.time(1:last_line); info.iterations.time = info.iterations.time(1:last_line);
info.iterations.error = info.iterations.error(1:last_line); info.iterations.error = info.iterations.error(1:last_line);
info.iterations.error
endo_simul = [ ];
end end

6
matlab/sim1.m
View File

@ -33,7 +33,7 @@ function sim1
% along with Dynare. If not, see <http://www.gnu.org/licenses/>. % along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global M_ options_ oo_ global M_ options_ oo_
global iyp iyf ct_ M_ it_ c global iyp iyf M_ it_ c
lead_lag_incidence = M_.lead_lag_incidence; lead_lag_incidence = M_.lead_lag_incidence;
@ -61,7 +61,7 @@ h1 = clock ;
for iter = 1:options_.maxit_ for iter = 1:options_.maxit_
h2 = clock ; h2 = clock ;
if ct_ == 0 if options_.terminal_condition == 0
c = zeros(ny*options_.periods,nrc) ; c = zeros(ny*options_.periods,nrc) ;
else else
c = zeros(ny*(options_.periods+1),nrc) ; c = zeros(ny*(options_.periods+1),nrc) ;
@ -84,7 +84,7 @@ for iter = 1:options_.maxit_
c (ic,:) = jacobian(:,is)\jacobian(:,isf1) ; c (ic,:) = jacobian(:,is)\jacobian(:,isf1) ;
end end
if ct_ == 1 if options_.terminal_condition == 1
s = eye(ny) ; s = eye(ny) ;
s(:,isf) = s(:,isf)+c(ic,1:nyf) ; s(:,isf) = s(:,isf)+c(ic,1:nyf) ;
ic = ic + ny ; ic = ic + ny ;

6
matlab/simk.m
View File

@ -35,7 +35,7 @@ function simk
% along with Dynare. If not, see <http://www.gnu.org/licenses/>. % along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global M_ options_ oo_ global M_ options_ oo_
global it_ iyr0 ct_ broyden_ global it_ iyr0 broyden_
%func_name = [M_.fname '_static']; %func_name = [M_.fname '_static'];
nk = M_.maximum_endo_lag + M_.maximum_endo_lead + 1 ; nk = M_.maximum_endo_lag + M_.maximum_endo_lead + 1 ;
@ -283,7 +283,7 @@ for iter = 1:options_.maxit_
iyr = iyr + ny ; iyr = iyr + ny ;
icr0 = icr0 + ny ; icr0 = icr0 + ny ;
end end
if ct_ == 1 if options_.terminal_condition == 1
ofs = (((it_-M_.maximum_lag-2)*ny+1)-1)*ncc*8 ; ofs = (((it_-M_.maximum_lag-2)*ny+1)-1)*ncc*8 ;
junk = fseek(fid,ofs,-1) ; junk = fseek(fid,ofs,-1) ;
@ -303,7 +303,7 @@ for iter = 1:options_.maxit_
end end
end end
oo_.endo_simul = reshape(oo_.endo_simul,ny,options_.periods+M_.maximum_lag+M_.maximum_endo_lead) ; 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 ; hbacsup = clock ;
c = bksupk(ny,fid,ncc,icc1) ; c = bksupk(ny,fid,ncc,icc1) ;
hbacsup = etime(clock,hbacsup) ; hbacsup = etime(clock,hbacsup) ;