From 1ea763193093309ca414b2e11886541949ddbea1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Fri, 27 Dec 2013 18:35:53 +0100 Subject: [PATCH] Added new option for perfect foresight simulations (sim1 routine, available only with stack_solve_algo==0). Try to reduce the size of the nonlinear system of equations by skipping the (last) periods for wich the residuals are already (almost) zero. The number of periods is not constant during the Newton, the effective number of periods for each iteration of the Newton is available in oo_.deterministic_simulation.vperiods. --- matlab/global_initialization.m | 1 + matlab/sim1.m | 31 +++++-- preprocessor/DynareBison.yy | 4 +- preprocessor/DynareFlex.ll | 1 + tests/Makefile.am | 1 + tests/deterministic_simulations/rbc_det6.mod | 85 ++++++++++++++++++++ 6 files changed, 115 insertions(+), 8 deletions(-) create mode 100644 tests/deterministic_simulations/rbc_det6.mod diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index 04b2242a2..f00cb2615 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -273,6 +273,7 @@ options_.linear = 0; options_.stack_solve_algo = 0; options_.markowitz = 0.5; options_.minimal_solving_periods = 1; +options_.endogenous_terminal_period = 0; % Solution options_.order = 2; diff --git a/matlab/sim1.m b/matlab/sim1.m index 677bf07d9..fcf445a4d 100644 --- a/matlab/sim1.m +++ b/matlab/sim1.m @@ -32,6 +32,9 @@ function sim1() global M_ options_ oo_ +endogenous_terminal_period = options_.endogenous_terminal_period; +vperiods = options_.periods*ones(1,options_.simul.maxit); + lead_lag_incidence = M_.lead_lag_incidence; ny = M_.endo_nbr; @@ -79,7 +82,6 @@ z = Y(find(lead_lag_incidence')); A = sparse([],[],[],periods*ny,periods*ny,periods*nnz(jacobian)); res = zeros(periods*ny,1); - h1 = clock ; for iter = 1:options_.simul.maxit @@ -88,11 +90,10 @@ for iter = 1:options_.simul.maxit i_rows = 1:ny; i_cols = find(lead_lag_incidence'); i_cols_A = i_cols; - + for it = (M_.maximum_lag+1):(M_.maximum_lag+periods) - [d1,jacobian] = model_dynamic(Y(i_cols),exo_simul, params, ... - steady_state,it); + [d1,jacobian] = model_dynamic(Y(i_cols), exo_simul, params, steady_state,it); if it == M_.maximum_lag+periods && it == M_.maximum_lag+1 A(i_rows,i_cols_A0) = jacobian(:,i_cols_0); elseif it == M_.maximum_lag+periods @@ -105,8 +106,17 @@ for iter = 1:options_.simul.maxit res(i_rows) = d1; + if endogenous_terminal_period && iter>1 + dr = max(abs(d1)); + if 10000*dr M_.maximum_lag+1 i_cols_A = i_cols_A + ny; end @@ -125,7 +135,12 @@ for iter = 1:options_.simul.maxit break end - dy = -A\res; + if endogenous_terminal_period && iter>1 + dy = zeros(length(i_upd),1); + dy(1:i_rows(end)) = -A(1:i_rows(end),1:i_rows(end))\res(1:i_rows(end)); + else + dy = -A\res; + end Y(i_upd) = Y(i_upd) + dy; @@ -137,6 +152,7 @@ if stop oo_.deterministic_simulation.status = 0;% NaN or Inf occurred oo_.deterministic_simulation.error = err; oo_.deterministic_simulation.iterations = iter; + oo_.deterministic_simulation.periods = vperiods(1:iter); oo_.endo_simul = reshape(Y,ny,periods+2); skipline(); fprintf('\nSimulation terminated after %d iterations.\n',iter); @@ -150,6 +166,7 @@ if stop oo_.deterministic_simulation.status = 1;% Convergency obtained. oo_.deterministic_simulation.error = err; oo_.deterministic_simulation.iterations = iter; + oo_.deterministic_simulation.periods = vperiods(1:iter); oo_.endo_simul = reshape(Y,ny,periods+2); end elseif ~stop @@ -159,9 +176,9 @@ elseif ~stop fprintf('WARNING : maximum number of iterations is reached (modify options_.simul.maxit).\n') ; oo_.deterministic_simulation.status = 0;% more iterations are needed. oo_.deterministic_simulation.error = err; + oo_.deterministic_simulation.periods = vperiods(1:iter); %oo_.deterministic_simulation.errors = c/abs(err) oo_.deterministic_simulation.iterations = options_.simul.maxit; end disp (['-----------------------------------------------------']) ; -skipline(); - +skipline(); \ No newline at end of file diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy index 3e95184c9..8e806bd90 100644 --- a/preprocessor/DynareBison.yy +++ b/preprocessor/DynareBison.yy @@ -122,7 +122,7 @@ class ParsingDriver; %token QUOTED_STRING %token QZ_CRITERIUM QZ_ZERO_THRESHOLD FULL DSGE_VAR DSGE_VARLAG DSGE_PRIOR_WEIGHT TRUNCATE %token RELATIVE_IRF REPLIC SIMUL_REPLIC RPLOT SAVE_PARAMS_AND_STEADY_STATE PARAMETER_UNCERTAINTY -%token SHOCKS SHOCK_DECOMPOSITION SIGMA_E SIMUL SIMUL_ALGO SIMUL_SEED SIMULATION_TYPE +%token SHOCKS SHOCK_DECOMPOSITION SIGMA_E SIMUL SIMUL_ALGO SIMUL_SEED SIMULATION_TYPE ENDOGENOUS_TERMINAL_PERIOD %token SMOOTHER SQUARE_ROOT_SOLVER STACK_SOLVE_ALGO STEADY_STATE_MODEL STOCHASTIC SOLVE_ALGO SOLVER_PERIODS %token STDERR STEADY STOCH_SIMUL SURPRISE SYLVESTER SYLVESTER_FIXED_POINT_TOL REGIMES REGIME %token TEX RAMSEY_POLICY PLANNER_DISCOUNT DISCRETIONARY_POLICY DISCRETIONARY_TOL @@ -986,6 +986,7 @@ simul_options : o_periods | o_markowitz | o_minimal_solving_periods | o_simul_maxit + | o_endogenous_terminal_period ; external_function : EXTERNAL_FUNCTION '(' external_function_options_list ')' ';' @@ -2329,6 +2330,7 @@ o_simul_algo : SIMUL_ALGO EQUAL INT_NUMBER { driver.error("simul_algo=1 option is no longer supported"); }; o_stack_solve_algo : STACK_SOLVE_ALGO EQUAL INT_NUMBER { driver.option_num("stack_solve_algo", $3); }; +o_endogenous_terminal_period : ENDOGENOUS_TERMINAL_PERIOD { driver.option_num("endogenous_terminal_period", "1"); }; o_linear : LINEAR { driver.linear(); }; o_order : ORDER EQUAL INT_NUMBER { driver.option_num("order", $3); }; o_replic : REPLIC EQUAL INT_NUMBER { driver.option_num("replic", $3); }; diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll index 9b6b5ed46..98dda2663 100644 --- a/preprocessor/DynareFlex.ll +++ b/preprocessor/DynareFlex.ll @@ -295,6 +295,7 @@ DATE -?[0-9]+([YyAa]|[Mm]([1-9]|1[0-2])|[Qq][1-4]|[Ww]([1-9]{1}|[1-4][0-9]|5[0-2 nocorr {return token::NOCORR;} optim {return token::OPTIM;} periods {return token::PERIODS;} +endogenous_terminal_period {return token::ENDOGENOUS_TERMINAL_PERIOD;} sub_draws {return token::SUB_DRAWS;} minimal_solving_periods {return token::MINIMAL_SOLVING_PERIODS;} markowitz {return token::MARKOWITZ;} diff --git a/tests/Makefile.am b/tests/Makefile.am index 31a42eddc..ec51bad78 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -150,6 +150,7 @@ MODFILES = \ deterministic_simulations/rbc_det3.mod \ deterministic_simulations/rbc_det4.mod \ deterministic_simulations/rbc_det5.mod \ + deterministic_simulations/rbc_det6.mod \ walsh.mod \ measurement_errors/fs2000_corr_me_ml_mcmc/fs2000_corr_ME.mod \ trend_var/fs2000_nonstationary.mod \ diff --git a/tests/deterministic_simulations/rbc_det6.mod b/tests/deterministic_simulations/rbc_det6.mod new file mode 100644 index 000000000..2d46248da --- /dev/null +++ b/tests/deterministic_simulations/rbc_det6.mod @@ -0,0 +1,85 @@ +var Capital, Output, Labour, Consumption, Efficiency, efficiency, ExpectedTerm; + +varexo EfficiencyInnovation; + +parameters beta, theta, tau, alpha, psi, delta, rho, effstar, sigma2; + +beta = 0.9900; +theta = 0.3570; +tau = 2.0000; +alpha = 0.4500; +psi = -0.1000; +delta = 0.0200; +rho = 0.8000; +effstar = 1.0000; +sigma2 = 0; + +model(use_dll); + + // Eq. n°1: + efficiency = rho*efficiency(-1) + EfficiencyInnovation; + + // Eq. n°2: + Efficiency = effstar*exp(efficiency); + + // Eq. n°3: + Output = Efficiency*(alpha*(Capital(-1)^psi)+(1-alpha)*(Labour^psi))^(1/psi); + + // Eq. n°4: + Capital = Output-Consumption + (1-delta)*Capital(-1); + + // Eq. n°5: + ((1-theta)/theta)*(Consumption/(1-Labour)) - (1-alpha)*(Output/Labour)^(1-psi); + + // Eq. n°6: + (((Consumption^theta)*((1-Labour)^(1-theta)))^(1-tau))/Consumption = ExpectedTerm(1); + + // Eq. n°7: + ExpectedTerm = beta*((((Consumption^theta)*((1-Labour)^(1-theta)))^(1-tau))/Consumption)*(alpha*((Output/Capital(-1))^(1-psi))+(1-delta)); + +end; + +steady_state_model; +efficiency = EfficiencyInnovation/(1-rho); +Efficiency = effstar*exp(efficiency); +Output_per_unit_of_Capital=((1/beta-1+delta)/alpha)^(1/(1-psi)); +Consumption_per_unit_of_Capital=Output_per_unit_of_Capital-delta; +Labour_per_unit_of_Capital=(((Output_per_unit_of_Capital/Efficiency)^psi-alpha)/(1-alpha))^(1/psi); +Output_per_unit_of_Labour=Output_per_unit_of_Capital/Labour_per_unit_of_Capital; +Consumption_per_unit_of_Labour=Consumption_per_unit_of_Capital/Labour_per_unit_of_Capital; + +% Compute steady state share of capital. +ShareOfCapital=alpha/(alpha+(1-alpha)*Labour_per_unit_of_Capital^psi); + +% Compute steady state of the endogenous variables. +Labour=1/(1+Consumption_per_unit_of_Labour/((1-alpha)*theta/(1-theta)*Output_per_unit_of_Labour^(1-psi))); +Consumption=Consumption_per_unit_of_Labour*Labour; +Capital=Labour/Labour_per_unit_of_Capital; +Output=Output_per_unit_of_Capital*Capital; +ExpectedTerm=beta*((((Consumption^theta)*((1-Labour)^(1-theta)))^(1-tau))/Consumption) + *(alpha*((Output/Capital)^(1-psi))+1-delta); +end; + +steady; + +ik = varlist_indices('Capital',M_.endo_names); +CapitalSS = oo_.steady_state(ik); + +histval; +Capital(0) = CapitalSS/2; +end; + +//options_.endogenous_terminal_period = 0; +simul(periods=500); +fff = oo_.endo_simul; + +oo_.deterministic_simulation + +//options_.endogenous_terminal_period = 1; +simul(periods=500, endogenous_terminal_period); +ggg = oo_.endo_simul; + +oo_.deterministic_simulation + +t1 = fff-ggg; +t2 = max(max(t1));