diff --git a/doc/dynare.texi b/doc/dynare.texi index 77debc657..5506670f0 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -2795,7 +2795,7 @@ Algorithm used for computing the solution. Possible values are: @item 0 Newton method to solve simultaneously all the equations for every -period, see @cite{Juillard (1996)} (Default). +period, using sparse matrices (Default). @item 1 Use a Newton algorithm with a sparse LU solver at each iteration @@ -2821,6 +2821,13 @@ declaration}). Use a Newton algorithm with a sparse Gaussian elimination (SPE) solver at each iteration (requires @code{bytecode} option, @pxref{Model declaration}). + +@item 6 +Use the historical algorithm proposed in @cite{Juillard (1996)}: it is +slower than @code{stack_solve_algo=0}, but may be less memory consuming +on big models (not available with @code{bytecode} and/or @code{block} +options). + @end table @item markowitz = @var{DOUBLE} diff --git a/matlab/sim1.m b/matlab/sim1.m index b8e27bc4e..d88a66341 100644 --- a/matlab/sim1.m +++ b/matlab/sim1.m @@ -1,21 +1,19 @@ function sim1 % function sim1 -% performs deterministic simulations with lead or lag on one period +% Performs deterministic simulations with lead or lag on one period. +% Uses sparse matrices. % % INPUTS % ... % OUTPUTS % ... % ALGORITHM -% Laffargue, Boucekkine, Juillard (LBJ) -% see Juillard (1996) Dynare: A program for the resolution and -% simulation of dynamic models with forward variables through the use -% of a relaxation algorithm. CEPREMAP. Couverture Orange. 9602. +% ... % % SPECIAL REQUIREMENTS % None. -% Copyright (C) 1996-2010 Dynare Team +% Copyright (C) 1996-2012 Dynare Team % % This file is part of Dynare. % @@ -36,13 +34,19 @@ global M_ options_ oo_ lead_lag_incidence = M_.lead_lag_incidence; -ny = size(oo_.endo_simul,1) ; +ny = M_.endo_nbr; + +max_lag = M_.maximum_endo_lag; + nyp = nnz(lead_lag_incidence(1,:)) ; -nyf = nnz(lead_lag_incidence(3,:)) ; -nrs = ny+nyp+nyf+1 ; -nrc = nyf+1 ; -iyf = find(lead_lag_incidence(3,:)>0) ; iyp = find(lead_lag_incidence(1,:)>0) ; +ny0 = nnz(lead_lag_incidence(2,:)) ; +iy0 = find(lead_lag_incidence(2,:)>0) ; +nyf = nnz(lead_lag_incidence(3,:)) ; +iyf = find(lead_lag_incidence(3,:)>0) ; + +nd = nyp+ny0+nyf; +nrc = nyf+1 ; isp = [1:nyp] ; is = [nyp+1:ny+nyp] ; isf = iyf+nyp ; @@ -50,57 +54,63 @@ isf1 = [nyp+ny+1:nyf+nyp+ny+1] ; stop = 0 ; iz = [1:ny+nyp+nyf]; +periods = options_.periods +steady_state = oo_.steady_state; +params = M_.params; +endo_simul = oo_.endo_simul; +exo_simul = oo_.exo_simul; +i_cols_1 = nonzeros(lead_lag_incidence(2:3,:)'); +i_cols_A1 = find(lead_lag_incidence(2:3,:)'); +i_cols_T = nonzeros(lead_lag_incidence(1:2,:)'); +i_cols_j = 1:nd; +i_upd = ny+(1:periods*ny); + +Y = endo_simul(:); + disp (['-----------------------------------------------------']) ; disp (['MODEL SIMULATION :']) ; fprintf('\n') ; -it_init = M_.maximum_lag+1 ; +model_dynamic = str2func([M_.fname,'_dynamic']); +z = Y(find(lead_lag_incidence')); +[d1,jacobian] = model_dynamic(z,oo_.exo_simul, params, ... + steady_state,2); + +A = sparse([],[],[],periods*ny,periods*ny,periods*nnz(jacobian)); +res = zeros(periods*ny,1); + + h1 = clock ; for iter = 1:options_.maxit_ h2 = clock ; - if options_.terminal_condition == 0 - c = zeros(ny*options_.periods,nrc) ; - else - c = zeros(ny*(options_.periods+1),nrc) ; - end + i_rows = 1:ny; + i_cols = find(lead_lag_incidence'); + i_cols_A = i_cols; - it_ = it_init ; - z = [oo_.endo_simul(iyp,it_-1) ; oo_.endo_simul(:,it_) ; oo_.endo_simul(iyf,it_+1)] ; - [d1,jacobian] = feval([M_.fname '_dynamic'],z,oo_.exo_simul, M_.params, oo_.steady_state,it_); - jacobian = [jacobian(:,iz) -d1] ; - ic = [1:ny] ; - icp = iyp ; - c (ic,:) = jacobian(:,is)\jacobian(:,isf1) ; - for it_ = it_init+(1:options_.periods-1) - z = [oo_.endo_simul(iyp,it_-1) ; oo_.endo_simul(:,it_) ; oo_.endo_simul(iyf,it_+1)] ; - [d1,jacobian] = feval([M_.fname '_dynamic'],z,oo_.exo_simul, ... - M_.params, oo_.steady_state, 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) ; + for it = 2:(periods+1) + + [d1,jacobian] = model_dynamic(Y(i_cols),exo_simul, params, ... + steady_state,it); + if it == 2 + A(i_rows,i_cols_A1) = jacobian(:,i_cols_1); + elseif it == periods+1 + A(i_rows,i_cols_A(i_cols_T)) = jacobian(:,i_cols_T); + else + A(i_rows,i_cols_A) = jacobian(:,i_cols_j); + end + + res(i_rows) = d1; + + i_rows = i_rows + ny; + i_cols = i_cols + ny; + if it > 2 + i_cols_A = i_cols_A + ny; + end end - - if options_.terminal_condition == 1 - s = eye(ny) ; - s(:,isf) = s(:,isf)+c(ic,1:nyf) ; - ic = ic + ny ; - c(ic,nrc) = s\c(ic,nrc) ; - c = bksup1(c,ny,nrc,iyf,options_.periods) ; - c = reshape(c,ny,options_.periods+1) ; - oo_.endo_simul(:,it_init+(0:options_.periods)) = oo_.endo_simul(:,it_init+(0:options_.periods))+options_.slowc*c ; - else - c = bksup1(c,ny,nrc,iyf,options_.periods) ; - c = reshape(c,ny,options_.periods) ; - oo_.endo_simul(:,it_init+(0:options_.periods-1)) = oo_.endo_simul(:,it_init+(0:options_.periods-1))+options_.slowc*c ; - end - - err = max(max(abs(c./options_.scalv'))); - disp([num2str(iter) ' - err = ' num2str(err)]) ; - disp([' Time of iteration :' num2str(etime(clock,h2))]) ; + + err = max(abs(res)); if err < options_.dynatol.f stop = 1 ; @@ -112,10 +122,17 @@ for iter = 1:options_.maxit_ oo_.deterministic_simulation.status = 1;% Convergency obtained. oo_.deterministic_simulation.error = err; oo_.deterministic_simulation.iterations = iter; + oo_.endo_simul = reshape(Y,ny,periods+2); break end + + dy = -A\res; + + Y(i_upd) = Y(i_upd) + dy; + end + if ~stop fprintf('\n') ; disp([' Total time of simulation :' num2str(etime(clock,h1))]) ; diff --git a/matlab/sim1a.m b/matlab/sim1_lbj.m similarity index 50% rename from matlab/sim1a.m rename to matlab/sim1_lbj.m index d552a2a92..57184583a 100644 --- a/matlab/sim1a.m +++ b/matlab/sim1_lbj.m @@ -1,19 +1,22 @@ -function sim1a -% function sim1a -% Performs deterministic simulations with lead or lag on one period -% Alternative algorithm to the one implemented in sim1 +function sim1_lbj +% function sim1_lbj +% performs deterministic simulations with lead or lag on one period +% using the historical LBJ algorithm % % INPUTS % ... % OUTPUTS % ... % ALGORITHM -% ... +% Laffargue, Boucekkine, Juillard (LBJ) +% see Juillard (1996) Dynare: A program for the resolution and +% simulation of dynamic models with forward variables through the use +% of a relaxation algorithm. CEPREMAP. Couverture Orange. 9602. % % SPECIAL REQUIREMENTS % None. -% Copyright (C) 1996-2012 Dynare Team +% Copyright (C) 1996-2010 Dynare Team % % This file is part of Dynare. % @@ -34,19 +37,13 @@ global M_ options_ oo_ lead_lag_incidence = M_.lead_lag_incidence; -ny = M_.endo_nbr; - -max_lag = M_.maximum_endo_lag; - +ny = size(oo_.endo_simul,1) ; nyp = nnz(lead_lag_incidence(1,:)) ; -iyp = find(lead_lag_incidence(1,:)>0) ; -ny0 = nnz(lead_lag_incidence(2,:)) ; -iy0 = find(lead_lag_incidence(2,:)>0) ; nyf = nnz(lead_lag_incidence(3,:)) ; -iyf = find(lead_lag_incidence(3,:)>0) ; - -nd = nyp+ny0+nyf; +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] ; isf = iyf+nyp ; @@ -54,63 +51,57 @@ isf1 = [nyp+ny+1:nyf+nyp+ny+1] ; stop = 0 ; iz = [1:ny+nyp+nyf]; -periods = options_.periods -steady_state = oo_.steady_state; -params = M_.params; -endo_simul = oo_.endo_simul; -exo_simul = oo_.exo_simul; -i_cols_1 = nonzeros(lead_lag_incidence(2:3,:)'); -i_cols_A1 = find(lead_lag_incidence(2:3,:)'); -i_cols_T = nonzeros(lead_lag_incidence(1:2,:)'); -i_cols_j = 1:nd; -i_upd = ny+(1:periods*ny); - -Y = endo_simul(:); - disp (['-----------------------------------------------------']) ; disp (['MODEL SIMULATION :']) ; fprintf('\n') ; +it_init = M_.maximum_lag+1 ; -model_dynamic = str2func([M_.fname,'_dynamic']); -z = Y(find(lead_lag_incidence')); -[d1,jacobian] = model_dynamic(z,oo_.exo_simul, params, ... - steady_state,2); - -A = sparse([],[],[],periods*ny,periods*ny,periods*nnz(jacobian)); -res = zeros(periods*ny,1); - - h1 = clock ; for iter = 1:options_.maxit_ h2 = clock ; - i_rows = 1:ny; - i_cols = find(lead_lag_incidence'); - i_cols_A = i_cols; - - for it = 2:(periods+1) - - [d1,jacobian] = model_dynamic(Y(i_cols),exo_simul, params, ... - steady_state,it); - if it == 2 - A(i_rows,i_cols_A1) = jacobian(:,i_cols_1); - elseif it == periods+1 - A(i_rows,i_cols_A(i_cols_T)) = jacobian(:,i_cols_T); - else - A(i_rows,i_cols_A) = jacobian(:,i_cols_j); - end - - res(i_rows) = d1; - - i_rows = i_rows + ny; - i_cols = i_cols + ny; - if it > 2 - i_cols_A = i_cols_A + ny; - end + if options_.terminal_condition == 0 + c = zeros(ny*options_.periods,nrc) ; + else + c = zeros(ny*(options_.periods+1),nrc) ; end - - err = max(abs(res)); + + it_ = it_init ; + z = [oo_.endo_simul(iyp,it_-1) ; oo_.endo_simul(:,it_) ; oo_.endo_simul(iyf,it_+1)] ; + [d1,jacobian] = feval([M_.fname '_dynamic'],z,oo_.exo_simul, M_.params, oo_.steady_state,it_); + jacobian = [jacobian(:,iz) -d1] ; + ic = [1:ny] ; + icp = iyp ; + c (ic,:) = jacobian(:,is)\jacobian(:,isf1) ; + for it_ = it_init+(1:options_.periods-1) + z = [oo_.endo_simul(iyp,it_-1) ; oo_.endo_simul(:,it_) ; oo_.endo_simul(iyf,it_+1)] ; + [d1,jacobian] = feval([M_.fname '_dynamic'],z,oo_.exo_simul, ... + M_.params, oo_.steady_state, 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 options_.terminal_condition == 1 + s = eye(ny) ; + s(:,isf) = s(:,isf)+c(ic,1:nyf) ; + ic = ic + ny ; + c(ic,nrc) = s\c(ic,nrc) ; + c = bksup1(c,ny,nrc,iyf,options_.periods) ; + c = reshape(c,ny,options_.periods+1) ; + oo_.endo_simul(:,it_init+(0:options_.periods)) = oo_.endo_simul(:,it_init+(0:options_.periods))+options_.slowc*c ; + else + c = bksup1(c,ny,nrc,iyf,options_.periods) ; + c = reshape(c,ny,options_.periods) ; + oo_.endo_simul(:,it_init+(0:options_.periods-1)) = oo_.endo_simul(:,it_init+(0:options_.periods-1))+options_.slowc*c ; + end + + err = max(max(abs(c./options_.scalv'))); + disp([num2str(iter) ' - err = ' num2str(err)]) ; + disp([' Time of iteration :' num2str(etime(clock,h2))]) ; if err < options_.dynatol.f stop = 1 ; @@ -122,17 +113,10 @@ for iter = 1:options_.maxit_ oo_.deterministic_simulation.status = 1;% Convergency obtained. oo_.deterministic_simulation.error = err; oo_.deterministic_simulation.iterations = iter; - oo_.endo_simul = reshape(Y,ny,periods+2); break end - - dy = -A\res; - - Y(i_upd) = Y(i_upd) + dy; - end - if ~stop fprintf('\n') ; disp([' Total time of simulation :' num2str(etime(clock,h1))]) ; diff --git a/matlab/simul.m b/matlab/simul.m index a0acb8512..cf96981e7 100644 --- a/matlab/simul.m +++ b/matlab/simul.m @@ -33,18 +33,23 @@ global M_ options_ oo_ test_for_deep_parameters_calibration(M_); -if options_.stack_solve_algo < 0 || options_.stack_solve_algo > 5 - error('SIMUL: stack_solve_algo must be between 0 and 5') +if options_.stack_solve_algo < 0 || options_.stack_solve_algo > 6 + error('SIMUL: stack_solve_algo must be between 0 and 6') end -if ~options_.block && ~options_.bytecode && options_.stack_solve_algo ~= 0 - error('SIMUL: you must use stack_solve_algo=0 when not using block nor bytecode option') +if ~options_.block && ~options_.bytecode && options_.stack_solve_algo ~= 0 ... + && options_.stack_solve_algo ~= 6 + error('SIMUL: you must use stack_solve_algo=0 or stack_solve_algo=6 when not using block nor bytecode option') end if options_.block && ~options_.bytecode && options_.stack_solve_algo == 5 error('SIMUL: you can''t use stack_solve_algo = 5 without bytecode option') end +if (options_.block || options_.bytecode) && options_.stack_solve_algo == 6 + error('SIMUL: you can''t use stack_solve_algo = 6 with block or bytecode option') +end + if exist('OCTAVE_VERSION') && options_.stack_solve_algo == 2 error('SIMUL: you can''t use stack_solve_algo = 2 under Octave') end @@ -91,7 +96,11 @@ else elseif M_.maximum_endo_lag == 0 % Purely forward model sim1_purely_forward; else % General case - sim1; + if options_.stack_solve_algo == 0 + sim1; + else % stack_solve_algo = 6 + sim1_lbj; + end end end; end; diff --git a/tests/run_test_matlab.m b/tests/run_test_matlab.m index c5a8c0755..c861d7fdc 100644 --- a/tests/run_test_matlab.m +++ b/tests/run_test_matlab.m @@ -78,7 +78,7 @@ for blockFlag = 0:1 default_stack_solve_algo = 0; if ~blockFlag && ~bytecodeFlag solve_algos = 1:4; - stack_solve_algos = 0; + stack_solve_algos = [0 6]; elseif blockFlag && ~bytecodeFlag solve_algos = [1:4 6:8]; stack_solve_algos = 0:4; diff --git a/tests/run_test_octave.m b/tests/run_test_octave.m index 66f0c0cbd..26742a6a2 100644 --- a/tests/run_test_octave.m +++ b/tests/run_test_octave.m @@ -75,7 +75,7 @@ for blockFlag = 0:1 default_stack_solve_algo = 0; if !blockFlag && !bytecodeFlag solve_algos = 0:4; - stack_solve_algos = 0; + stack_solve_algos = [0 6]; elseif blockFlag && !bytecodeFlag solve_algos = [0:4 6 8]; stack_solve_algos = [0 1 3 4];