diff --git a/src/ComputingTasks.cc b/src/ComputingTasks.cc index 40d3c002..5dc5623f 100644 --- a/src/ComputingTasks.cc +++ b/src/ComputingTasks.cc @@ -683,6 +683,58 @@ ForecastStatement::writeJsonOutput(ostream &output) const output << "}"; } +DetCondForecast::DetCondForecast(const SymbolList &symbol_list_arg, + const OptionsList &options_list_arg, + const bool linear_decomposition_arg) : + options_list(options_list_arg), + symbol_list(symbol_list_arg), + linear_decomposition(linear_decomposition_arg) +{ + +} + +void +DetCondForecast::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const +{ + options_list.writeOutput(output); + if (linear_decomposition) + { + output << "first_order_solution_to_compute = 1;" << endl; + output << "if exist('oo_')" << endl; + output << " if isfield(oo_, 'dr')" << endl; + output << " if isfield(oo_.dr, 'ghx') && isfield(oo_.dr, 'ghu') && isfield(oo_.dr, 'state_var') && isfield(oo_.dr, 'order_var')" << endl; + output << " first_order_solution_to_compute = 0;" << endl; + output << " end;" << endl; + output << " end;" << endl; + output << "end;" << endl; + output << "if first_order_solution_to_compute" << endl; + output << " fprintf('%s','Computing the first order solution ...');" << endl; + output << " options_.nograph = 1;" << endl; + output << " options_.order = 1;" << endl; + output << " options_.noprint = 1;" << endl; + output << " options_.nocorr = 1;" << endl; + output << " options_.nomoments = 1;" << endl; + output << " options_.nodecomposition = 1;" << endl; + output << " options_.nofunctions = 1;" << endl; + output << " options_.irf = 0;" << endl; + output << " tmp_periods = options_.periods;" << endl; + output << " options_.periods = 0;" << endl; + output << " var_list_ = char();" << endl; + output << " info = stoch_simul(var_list_);" << endl; + output << " fprintf('%s\\n','done');" << endl; + output << " options_.periods = tmp_periods;" << endl; + output << "end;" << endl; + } + vector symbols = symbol_list.get_symbols(); + if (symbols.size() > 0) + output << symbols[1] << " = det_cond_forecast(" ; + for (unsigned int i = 0; i < symbols.size() - 1; i++) + output << symbols[i] << ", "; + if (symbols.size() > 0) + output << symbols[symbols.size() - 1]; + output << ");" << endl; +} + RamseyModelStatement::RamseyModelStatement(OptionsList options_list_arg) : options_list(move(options_list_arg)) { diff --git a/src/ComputingTasks.hh b/src/ComputingTasks.hh index 4985c2b7..6fb1264b 100644 --- a/src/ComputingTasks.hh +++ b/src/ComputingTasks.hh @@ -95,6 +95,21 @@ public: void writeJsonOutput(ostream &output) const override; }; + +class DetCondForecast : public Statement +{ +private: + const OptionsList options_list; + const SymbolList symbol_list; + const bool linear_decomposition; +public: + DetCondForecast(const SymbolList &symbol_list_arg, + const OptionsList &options_list_arg, + const bool linear_decompositiontion_arg); + //virtual void checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings); + virtual void writeOutput(ostream &output, const string &basename, bool minimal_workspace) const; +}; + class ModelInfoStatement : public Statement { private: diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc index 6484a26c..b7dbe15b 100644 --- a/src/DynamicModel.cc +++ b/src/DynamicModel.cc @@ -1043,7 +1043,7 @@ DynamicModel::writeModelEquationsCode(const string &basename, const map_idx_t &m } void -DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_idx_t &map_idx) const +DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_idx_t &map_idx, const bool linear_decomposition) const { struct Uff_l { @@ -1068,10 +1068,12 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id deriv_node_temp_terms_t tef_terms; vector feedback_variables; bool file_open = false; - + string main_name; boost::filesystem::create_directories(basename + "/model/bytecode"); - - string main_name = basename + "/model/bytecode/dynamic.cod"; + if (linear_decomposition) + main_name = basename + "/model/bytecode/non_linear.cod"; + else + main_name = basename + "/model/bytecode/dynamic.cod"; code_file.open(main_name, ios::out | ios::binary | ios::ate); if (!code_file.is_open()) { @@ -1104,7 +1106,7 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id || simulation_type == SOLVE_BACKWARD_COMPLETE || simulation_type == SOLVE_FORWARD_COMPLETE) { Write_Inf_To_Bin_File_Block(basename, block, u_count_int, file_open, - simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE); + simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE, linear_decomposition); file_open = true; } map>, expr_t> tmp_block_endo_derivative; @@ -1187,13 +1189,16 @@ DynamicModel::writeModelEquationsCode_Block(const string &basename, const map_id ); fbeginblock.write(code_file, instruction_number); + if (linear_decomposition) + compileTemporaryTerms(code_file, instruction_number, temporary_terms, map_idx, true, false); + // The equations for (i = 0; i < (int) block_size; i++) { //The Temporary terms temporary_terms_t tt2; tt2.clear(); - if (v_temporary_terms[block][i].size()) + if (v_temporary_terms[block][i].size() && !linear_decomposition) { for (auto it : v_temporary_terms[block][i]) { @@ -1715,11 +1720,17 @@ DynamicModel::setNonZeroHessianEquations(map &eqs) void DynamicModel::Write_Inf_To_Bin_File_Block(const string &basename, const int &num, - int &u_count_int, bool &file_open, bool is_two_boundaries) const + int &u_count_int, bool &file_open, bool is_two_boundaries, const bool linear_decomposition) const { int j; std::ofstream SaveCode; - string filename = basename + "/model/bytecode/dynamic.bin"; + string filename; + + if (!linear_decomposition) + filename = basename + "/model/bytecode/dynamic.bin"; + else + filename = basename + "/model/bytecode/non_linear.bin"; + if (file_open) SaveCode.open(filename, ios::out | ios::in | ios::binary | ios::ate); else @@ -2814,7 +2825,7 @@ DynamicModel::writeDynamicModel(const string &basename, ostream &DynamicOutput, } void -DynamicModel::writeOutput(ostream &output, const string &basename, bool block_decomposition, bool byte_code, bool use_dll, int order, bool estimation_present, bool compute_xrefs, bool julia) const +DynamicModel::writeOutput(ostream &output, const string &basename, bool block_decomposition, bool linear_decomposition, bool byte_code, bool use_dll, int order, bool estimation_present, bool compute_xrefs, bool julia) const { /* Writing initialisation for M_.lead_lag_incidence matrix M_.lead_lag_incidence is a matrix with as many columns as there are @@ -2956,7 +2967,7 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de } //In case of sparse model, writes the block_decomposition structure of the model - if (block_decomposition) + if (block_decomposition || linear_decomposition) { vector state_equ; int count_lead_lag_incidence = 0; @@ -4179,7 +4190,7 @@ DynamicModel::substitutePacExpectation() void DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivatives, int paramsDerivsOrder, const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll, - bool bytecode, const bool nopreprocessoroutput) + bool bytecode, const bool nopreprocessoroutput, bool linear_decomposition) { assert(jacobianExo || !(hessian || thirdDerivatives || paramsDerivsOrder)); @@ -4203,8 +4214,15 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative // Launch computations if (!nopreprocessoroutput) - cout << "Computing dynamic model derivatives:" << endl - << " - order 1" << endl; + { + if (linear_decomposition) + cout << "Computing nonlinear dynamic model derivatives:" << endl + << " - order 1" << endl; + else + cout << "Computing dynamic model derivatives:" << endl + << " - order 1" << endl; + } + computeJacobian(vars); if (hessian) @@ -4231,13 +4249,51 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative computeThirdDerivatives(vars); } + jacob_map_t contemporaneous_jacobian, static_jacobian; + map >, expr_t> first_order_endo_derivatives; + // for each block contains pair + vector > blocks; + vector n_static, n_forward, n_backward, n_mixed; + + if (linear_decomposition) + { + map>, expr_t> first_order_endo_derivatives = collect_first_order_derivatives_endogenous(); + is_equation_linear = equationLinear(first_order_endo_derivatives); + + evaluateAndReduceJacobian(eval_context, contemporaneous_jacobian, static_jacobian, dynamic_jacobian, cutoff, false); + + if (!computeNaturalNormalization()) + computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian, dynamic_jacobian); + + lag_lead_vector_t equation_lag_lead, variable_lag_lead; + + blocks = select_non_linear_equations_and_variables(is_equation_linear, dynamic_jacobian, equation_reordered, variable_reordered, + inv_equation_reordered, inv_variable_reordered, + equation_lag_lead, variable_lag_lead, + n_static, n_forward, n_backward, n_mixed); + + equation_type_and_normalized_equation = equationTypeDetermination(first_order_endo_derivatives, variable_reordered, equation_reordered, 0); + prologue = 0; + epilogue = 0; + + block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(dynamic_jacobian, blocks, equation_type_and_normalized_equation, variable_reordered, equation_reordered, n_static, n_forward, n_backward, n_mixed, block_col_type, linear_decomposition); + + computeChainRuleJacobian(blocks_derivatives); + + blocks_linear = BlockLinear(blocks_derivatives, variable_reordered); + + collect_block_first_order_derivatives(); + + collectBlockVariables(); + + global_temporary_terms = true; + if (!no_tmp_terms) + computeTemporaryTermsOrdered(); + + } + if (block) { - vector n_static, n_forward, n_backward, n_mixed; - jacob_map_t contemporaneous_jacobian, static_jacobian; - - // for each block contains pair - vector> blocks; evaluateAndReduceJacobian(eval_context, contemporaneous_jacobian, static_jacobian, dynamic_jacobian, cutoff, false); @@ -4245,7 +4301,7 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative computePrologueAndEpilogue(static_jacobian, equation_reordered, variable_reordered); - map>, expr_t> first_order_endo_derivatives = collect_first_order_derivatives_endogenous(); + first_order_endo_derivatives = collect_first_order_derivatives_endogenous(); equation_type_and_normalized_equation = equationTypeDetermination(first_order_endo_derivatives, variable_reordered, equation_reordered, mfs); @@ -4256,7 +4312,7 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, dynamic_jacobian, equation_reordered, variable_reordered, blocks, equation_type_and_normalized_equation, false, true, mfs, inv_equation_reordered, inv_variable_reordered, equation_lag_lead, variable_lag_lead, n_static, n_forward, n_backward, n_mixed); - block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(dynamic_jacobian, blocks, equation_type_and_normalized_equation, variable_reordered, equation_reordered, n_static, n_forward, n_backward, n_mixed, block_col_type); + block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(dynamic_jacobian, blocks, equation_type_and_normalized_equation, variable_reordered, equation_reordered, n_static, n_forward, n_backward, n_mixed, block_col_type, linear_decomposition); printBlockDecomposition(blocks); @@ -4685,12 +4741,16 @@ DynamicModel::collectBlockVariables() } void -DynamicModel::writeDynamicFile(const string &basename, bool block, bool bytecode, bool use_dll, int order, bool julia) const +DynamicModel::writeDynamicFile(const string &basename, bool block, bool linear_decomposition, bool bytecode, bool use_dll, int order, bool julia) const { if (block && bytecode) - writeModelEquationsCode_Block(basename, map_idx); + writeModelEquationsCode_Block(basename, map_idx, linear_decomposition); else if (!block && bytecode) - writeModelEquationsCode(basename, map_idx); + { + if (linear_decomposition) + writeModelEquationsCode_Block(basename, map_idx, linear_decomposition); + writeModelEquationsCode(basename, map_idx); + } else if (block && !bytecode) writeSparseDynamicMFile(basename); else if (use_dll) @@ -4891,6 +4951,15 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel &static_model, const boo addEquation(neweq, -1); } +void +DynamicModel::toNonlinearPart(DynamicModel &non_linear_equations_dynamic_model) const +{ + // Convert model local variables (need to be done first) + for (map::const_iterator it = local_variables_table.begin(); + it != local_variables_table.end(); it++) + non_linear_equations_dynamic_model.AddLocalVariable(it->first, it->second); +} + void DynamicModel::toStatic(StaticModel &static_model) const { diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh index 2ec940ab..1e4094a0 100644 --- a/src/DynamicModel.hh +++ b/src/DynamicModel.hh @@ -118,7 +118,7 @@ private: //! Writes the Block reordred structure of the model in M output void writeModelEquationsOrdered_M(const string &basename) const; //! Writes the code of the Block reordred structure of the model in virtual machine bytecode - void writeModelEquationsCode_Block(const string &basename, const map_idx_t &map_idx) const; + void writeModelEquationsCode_Block(const string &basename, const map_idx_t &map_idx, const bool linear_decomposition) const; //! Writes the code of the model in virtual machine bytecode void writeModelEquationsCode(const string &basename, const map_idx_t &map_idx) const; @@ -280,9 +280,9 @@ public: \param no_tmp_terms if true, no temporary terms will be computed in the dynamic files */ void computingPass(bool jacobianExo, bool hessian, bool thirdDerivatives, int paramsDerivsOrder, - const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll, bool bytecode, const bool nopreprocessoroutput); + const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll, bool bytecode, const bool nopreprocessoroutput, bool linear_decomposition); //! Writes model initialization and lead/lag incidence matrix to output - void writeOutput(ostream &output, const string &basename, bool block, bool byte_code, bool use_dll, int order, bool estimation_present, bool compute_xrefs, bool julia) const; + void writeOutput(ostream &output, const string &basename, bool block, bool linear_decomposition, bool byte_code, bool use_dll, int order, bool estimation_present, bool compute_xrefs, bool julia) const; //! Write JSON AST void writeJsonAST(ostream &output) const; @@ -351,9 +351,9 @@ public: //! Adds informations for simulation in a binary file void Write_Inf_To_Bin_File_Block(const string &basename, - const int &num, int &u_count_int, bool &file_open, bool is_two_boundaries) const; + const int &num, int &u_count_int, bool &file_open, bool is_two_boundaries, const bool linear_decomposition) const; //! Writes dynamic model file - void writeDynamicFile(const string &basename, bool block, bool bytecode, bool use_dll, int order, bool julia) const; + void writeDynamicFile(const string &basename, bool block, bool linear_decomposition, bool bytecode, bool use_dll, int order, bool julia) const; //! Writes file containing parameters derivatives void writeParamsDerivativesFile(const string &basename, bool julia) const; @@ -361,6 +361,11 @@ public: /*! It assumes that the static model given in argument has just been allocated */ void toStatic(StaticModel &static_model) const; + + //! Converts to nonlinear model (only the equations) + /*! It assumes that the nonlinear model given in argument has just been allocated */ + void toNonlinearPart(DynamicModel &non_linear_equations_dynamic_model) const; + //! Find endogenous variables not used in model set findUnusedEndogenous(); //! Find exogenous variables not used in model diff --git a/src/DynareBison.yy b/src/DynareBison.yy index c9c77610..5a16ad12 100644 --- a/src/DynareBison.yy +++ b/src/DynareBison.yy @@ -81,7 +81,7 @@ class ParsingDriver; %token BVAR_REPLIC BYTECODE ALL_VALUES_REQUIRED PROPOSAL_DISTRIBUTION REALTIME VINTAGE %token CALIB_SMOOTHER CHANGE_TYPE CHECK CONDITIONAL_FORECAST CONDITIONAL_FORECAST_PATHS CONF_SIG CONSTANT CONTROLLED_VAREXO CORR CUTOFF CYCLE_REDUCTION LOGARITHMIC_REDUCTION %token COMMA CONSIDER_ALL_ENDOGENOUS CONSIDER_ONLY_OBSERVED INITIAL_CONDITION_DECOMPOSITION -%token DATAFILE FILE SERIES DOUBLING DR_CYCLE_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_MAXITER DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION DIFFERENTIATE_FORWARD_VARS +%token DATAFILE FILE SERIES DET_COND_FORECAST DOUBLING DR_CYCLE_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_MAXITER DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION DIFFERENTIATE_FORWARD_VARS %token END ENDVAL EQUAL ESTIMATION ESTIMATED_PARAMS ESTIMATED_PARAMS_BOUNDS ESTIMATED_PARAMS_INIT EXTENDED_PATH ENDOGENOUS_PRIOR %token FILENAME DIRNAME FILTER_STEP_AHEAD FILTERED_VARS FIRST_OBS LAST_OBS SET_TIME OSR_PARAMS_BOUNDS KEEP_KALMAN_ALGO_IF_SINGULARITY_IS_DETECTED %token FLOAT_NUMBER DATES @@ -93,7 +93,7 @@ class ParsingDriver; %token INT_NUMBER %token INV_GAMMA_PDF INV_GAMMA1_PDF INV_GAMMA2_PDF IRF IRF_SHOCKS IRF_PLOT_THRESHOLD IRF_CALIBRATION %token FAST_KALMAN_FILTER KALMAN_ALGO KALMAN_TOL DIFFUSE_KALMAN_TOL SUBSAMPLES OPTIONS TOLF TOLX PLOT_INIT_DATE PLOT_END_DATE -%token LAPLACE LIK_ALGO LIK_INIT LINEAR LOAD_IDENT_FILES LOAD_MH_FILE LOAD_RESULTS_AFTER_LOAD_MH LOAD_PARAMS_AND_STEADY_STATE LOGLINEAR LOGDATA LYAPUNOV LINEAR_APPROXIMATION +%token LAPLACE LIK_ALGO LIK_INIT LINEAR LINEAR_DECOMPOSITION LOAD_IDENT_FILES LOAD_MH_FILE LOAD_RESULTS_AFTER_LOAD_MH LOAD_PARAMS_AND_STEADY_STATE LOGLINEAR LOGDATA LYAPUNOV LINEAR_APPROXIMATION %token LYAPUNOV_FIXED_POINT_TOL LYAPUNOV_DOUBLING_TOL LOG_DEFLATOR LOG_TREND_VAR LOG_GROWTH_FACTOR MARKOWITZ MARGINAL_DENSITY MAX MAXIT %token MFS MH_CONF_SIG MH_DROP MH_INIT_SCALE MH_JSCALE MH_TUNE_JSCALE MH_MODE MH_NBLOCKS MH_REPLIC MH_RECOVER POSTERIOR_MAX_SUBSAMPLE_DRAWS MIN MINIMAL_SOLVING_PERIODS %token MODE_CHECK MODE_CHECK_NEIGHBOURHOOD_SIZE MODE_CHECK_SYMMETRIC_PLOTS MODE_CHECK_NUMBER_OF_POINTS MODE_COMPUTE MODE_FILE MODEL MODEL_COMPARISON MODEL_INFO MSHOCKS ABS SIGN @@ -301,6 +301,7 @@ statement : parameters | gmm_estimation | smm_estimation | shock_groups + | det_cond_forecast | var_expectation_model ; @@ -970,6 +971,7 @@ model_options : BLOCK { driver.block(); } | DIFFERENTIATE_FORWARD_VARS EQUAL '(' symbol_list ')' { driver.differentiate_forward_vars_some(); } | o_linear | PARALLEL_LOCAL_FILES EQUAL '(' parallel_local_filename_list ')' + | LINEAR_DECOMPOSITION { driver.linear_decomposition(); } ; model_options_list : model_options_list COMMA model_options @@ -1477,6 +1479,16 @@ prior_posterior_function_options : o_function | o_sampling_draws ; +det_cond_forecast : DET_COND_FORECAST '(' symbol ')' ';' + { driver.det_cond_forecast_linear_decomposition($3); } + | DET_COND_FORECAST '(' symbol COMMA symbol ')' ';' + { driver.det_cond_forecast_linear_decomposition($3,$5); } + | DET_COND_FORECAST '(' symbol COMMA simul_options_list ')' ';' + { driver.det_cond_forecast_linear_decomposition($3); } + | DET_COND_FORECAST '(' symbol COMMA symbol COMMA simul_options_list ')' ';' + { driver.det_cond_forecast_linear_decomposition($3,$5); } + ; + simul : SIMUL ';' { driver.simul(); } | SIMUL '(' simul_options_list ')' ';' diff --git a/src/DynareFlex.ll b/src/DynareFlex.ll index ae3b0631..2f3be1e0 100644 --- a/src/DynareFlex.ll +++ b/src/DynareFlex.ll @@ -186,6 +186,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 smoother2histval {BEGIN DYNARE_STATEMENT; return token::SMOOTHER2HISTVAL;} perfect_foresight_setup {BEGIN DYNARE_STATEMENT; return token::PERFECT_FORESIGHT_SETUP;} perfect_foresight_solver {BEGIN DYNARE_STATEMENT; return token::PERFECT_FORESIGHT_SOLVER;} +det_cond_forecast {BEGIN DYNARE_STATEMENT; return token::DET_COND_FORECAST;} ; { if (!sigma_e) @@ -776,6 +777,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 use_dll {return token::USE_DLL;} block {return token::BLOCK;} bytecode {return token::BYTECODE;} +linear_decomposition {return token::LINEAR_DECOMPOSITION;} all_values_required {return token::ALL_VALUES_REQUIRED;} no_static {return token::NO_STATIC;} differentiate_forward_vars {return token::DIFFERENTIATE_FORWARD_VARS;} diff --git a/src/ModFile.cc b/src/ModFile.cc index 47246a9a..ad7da6a2 100644 --- a/src/ModFile.cc +++ b/src/ModFile.cc @@ -43,6 +43,8 @@ ModFile::ModFile(WarningConsolidation &warnings_arg) trend_component_model_table, var_model_table), orig_ramsey_dynamic_model(symbol_table, num_constants, external_functions_table, trend_component_model_table, var_model_table), + non_linear_equations_dynamic_model(symbol_table, num_constants, external_functions_table, + trend_component_model_table, var_model_table), epilogue(symbol_table, num_constants, external_functions_table, trend_component_model_table, var_model_table), static_model(symbol_table, num_constants, external_functions_table), @@ -191,9 +193,20 @@ ModFile::checkPass(bool nostrict, bool stochastic) exit(EXIT_FAILURE); } - if (use_dll && (block || byte_code)) + if (use_dll && (block || byte_code || linear_decomposition)) { - cerr << "ERROR: In 'model' block, 'use_dll' option is not compatible with 'block' or 'bytecode'" << endl; + cerr << "ERROR: In 'model' block, 'use_dll' option is not compatible with 'block', 'bytecode' or 'linear_decomposition'" << endl; + exit(EXIT_FAILURE); + } + if (block && linear_decomposition) + { + cerr << "ERROR: In 'model' block, 'block' option is not compatible with 'linear_decomposition'" << endl; + exit(EXIT_FAILURE); + } + + if (!byte_code && linear_decomposition) + { + cerr << "ERROR: For the moment in 'model' block, 'linear_decomposition' option is compatible only with 'bytecode' option" << endl; exit(EXIT_FAILURE); } @@ -683,6 +696,12 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_deri // Compute static model and its derivatives dynamic_model.toStatic(static_model); + if (linear_decomposition) + { + dynamic_model.cloneDynamic(non_linear_equations_dynamic_model); + non_linear_equations_dynamic_model.set_cutoff_to_zero(); + non_linear_equations_dynamic_model.computingPass(true, false, false, 0, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition); + } if (!no_static) { if (mod_file_struct.stoch_simul_present @@ -707,7 +726,7 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_deri || mod_file_struct.calib_smoother_present) { if (mod_file_struct.perfect_foresight_solver_present) - dynamic_model.computingPass(true, false, false, 0, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput); + dynamic_model.computingPass(true, false, false, 0, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition); else { if (mod_file_struct.stoch_simul_present @@ -732,13 +751,13 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_deri int paramsDerivsOrder = 0; if (mod_file_struct.identification_present || mod_file_struct.estimation_analytic_derivation) paramsDerivsOrder = params_derivs_order; - dynamic_model.computingPass(true, hessian, thirdDerivatives, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput); + dynamic_model.computingPass(true, hessian, thirdDerivatives, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition); if (linear && mod_file_struct.ramsey_model_present) - orig_ramsey_dynamic_model.computingPass(true, true, false, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput); + orig_ramsey_dynamic_model.computingPass(true, true, false, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition); } } else // No computing task requested, compute derivatives up to 2nd order by default - dynamic_model.computingPass(true, true, false, 0, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput); + dynamic_model.computingPass(true, true, false, 0, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition); map eqs; if (mod_file_struct.ramsey_model_present) @@ -764,7 +783,7 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_deri for (auto & statement : statements) statement->computingPass(); - epilogue.computingPass(true, true, false, 0, global_eval_context, true, false, false, false, true); + epilogue.computingPass(true, true, false, 0, global_eval_context, true, false, false, false, true, false); } void @@ -903,7 +922,8 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo mOutputFile << "options_.block=" << block << ";" << endl << "options_.bytecode=" << byte_code << ";" << endl - << "options_.use_dll=" << use_dll << ";" << endl; + << "options_.use_dll=" << use_dll << ";" << endl + << "options_.linear_decomposition=" << linear_decomposition << ";" << endl; if (parallel_local_files.size() > 0) { @@ -992,7 +1012,9 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo if (dynamic_model.equation_number() > 0) { - dynamic_model.writeOutput(mOutputFile, basename, block, byte_code, use_dll, mod_file_struct.order_option, mod_file_struct.estimation_present, compute_xrefs, false); + if (linear_decomposition) + non_linear_equations_dynamic_model.writeOutput(mOutputFile, basename, block, true, byte_code, use_dll, mod_file_struct.order_option, mod_file_struct.estimation_present, compute_xrefs, false); + dynamic_model.writeOutput(mOutputFile, basename, block, false, byte_code, use_dll, mod_file_struct.order_option, mod_file_struct.estimation_present, compute_xrefs, false); if (!no_static) static_model.writeOutput(mOutputFile, block); } @@ -1066,7 +1088,14 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo static_model.writeParamsDerivativesFile(basename, false); } - dynamic_model.writeDynamicFile(basename, block, byte_code, use_dll, mod_file_struct.order_option, false); + if (linear_decomposition) + { + non_linear_equations_dynamic_model.writeDynamicFile(basename, block, linear_decomposition, byte_code, use_dll, mod_file_struct.order_option, false); + non_linear_equations_dynamic_model.writeParamsDerivativesFile(basename, false); + } + + dynamic_model.writeDynamicFile(basename, block, false, byte_code, use_dll, mod_file_struct.order_option, false); + dynamic_model.writeParamsDerivativesFile(basename, false); } @@ -1179,7 +1208,7 @@ ModFile::writeExternalFilesJulia(const string &basename, FileOutputType output, if (dynamic_model.equation_number() > 0) { - dynamic_model.writeOutput(jlOutputFile, basename, false, false, false, + dynamic_model.writeOutput(jlOutputFile, basename, false, false, false, false, mod_file_struct.order_option, mod_file_struct.estimation_present, false, true); if (!no_static) @@ -1187,7 +1216,7 @@ ModFile::writeExternalFilesJulia(const string &basename, FileOutputType output, static_model.writeStaticFile(basename, false, false, false, true); static_model.writeParamsDerivativesFile(basename, true); } - dynamic_model.writeDynamicFile(basename, block, byte_code, use_dll, + dynamic_model.writeDynamicFile(basename, block, linear_decomposition, byte_code, use_dll, mod_file_struct.order_option, true); dynamic_model.writeParamsDerivativesFile(basename, true); } diff --git a/src/ModFile.hh b/src/ModFile.hh index 5e47f392..72410555 100644 --- a/src/ModFile.hh +++ b/src/ModFile.hh @@ -67,6 +67,8 @@ public: DynamicModel ramsey_FOC_equations_dynamic_model; //! A copy of the original model, used to test model linearity under ramsey problem DynamicModel orig_ramsey_dynamic_model; + //! A copy of the Dynamic model, containing only the non-linear equations w.r. to endogenous variables + DynamicModel non_linear_equations_dynamic_model; //! Epilogue model, as declared in the "epilogue" block Epilogue epilogue; //! Static model, as derived from the "model" block when leads and lags have been removed @@ -100,6 +102,9 @@ public: with a lead */ vector differentiate_forward_vars_subset; + //! Is the model is block-decomposed according the linear and the non-linear equations + bool linear_decomposition; + //! Are nonstationary variables present ? bool nonstationary_variables; diff --git a/src/ModelTree.cc b/src/ModelTree.cc index 91e126f6..1dde43ac 100644 --- a/src/ModelTree.cc +++ b/src/ModelTree.cc @@ -316,6 +316,105 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, jacob_m } } +vector > +ModelTree::select_non_linear_equations_and_variables(vector is_equation_linear, const dynamic_jacob_map_t &dynamic_jacobian, vector &equation_reordered, vector &variable_reordered, + vector &inv_equation_reordered, vector &inv_variable_reordered, + lag_lead_vector_t &equation_lag_lead, lag_lead_vector_t &variable_lag_lead, + vector &n_static, vector &n_forward, vector &n_backward, vector &n_mixed) +{ + vector eq2endo(equations.size(), 0); + /*equation_reordered.resize(equations.size()); + variable_reordered.resize(equations.size());*/ + unsigned int num = 0; + for (vector::const_iterator it = endo2eq.begin(); it != endo2eq.end(); it++) + if (!is_equation_linear[*it]) + num++; + vector endo2block = vector(endo2eq.size(), 1); + vector, pair, vector > > > components_set(num); + int i = 0, j = 0; + for (vector::const_iterator it = endo2eq.begin(); it != endo2eq.end(); it++, j++) + { + if (!is_equation_linear[*it]) + { + equation_reordered[i] = *it; + variable_reordered[i] = j; + endo2block[j] = 0; + components_set[endo2block[j]].first.insert(i); + /*cout << " -----------------------------------------" << endl; + cout.flush(); + cout << "equation_reordered[" << *it << "] = " << i << endl; + cout.flush(); + cout << "variable_reordered[" << j << "] = " << i << endl; + cout.flush(); + cout << "components_set[" << endo2block[j] << "].first[" << i << "] = " << i << endl; + cout << "endo2block[" << j << "] = " << 0 << endl; + cout.flush(); + */ + i++; + } + } +/* for (unsigned int j = 0; j < is_equation_linear.size() ; j++) + cout << "endo2block[" << j << "] = " << endo2block[j] << endl;*/ + /*cout << "before getVariableLeadLAgByBlock\n"; + cout.flush();*/ + getVariableLeadLagByBlock(dynamic_jacobian, endo2block, endo2block.size(), equation_lag_lead, variable_lag_lead, equation_reordered, variable_reordered); + n_static = vector(endo2eq.size(), 0); + n_forward = vector(endo2eq.size(), 0); + n_backward = vector(endo2eq.size(), 0); + n_mixed = vector(endo2eq.size(), 0); + for (unsigned int i = 0; i < endo2eq.size(); i++) + { + if (variable_lag_lead[variable_reordered[i]].first != 0 && variable_lag_lead[variable_reordered[i]].second != 0) + n_mixed[i]++; + else if (variable_lag_lead[variable_reordered[i]].first == 0 && variable_lag_lead[variable_reordered[i]].second != 0) + n_forward[i]++; + else if (variable_lag_lead[variable_reordered[i]].first != 0 && variable_lag_lead[variable_reordered[i]].second == 0) + n_backward[i]++; + else if (variable_lag_lead[variable_reordered[i]].first == 0 && variable_lag_lead[variable_reordered[i]].second == 0) + n_static[i]++; + } + cout.flush(); + int nb_endo = is_equation_linear.size(); + vector > blocks = vector >(1, make_pair(i, i)); + inv_equation_reordered = vector(nb_endo); + inv_variable_reordered = vector(nb_endo); + for (int i = 0; i < nb_endo; i++) + { + inv_variable_reordered[variable_reordered[i]] = i; + inv_equation_reordered[equation_reordered[i]] = i; + } + return blocks; +} + +bool +ModelTree::computeNaturalNormalization() +{ + bool bool_result = true; + set > result; + endo2eq.resize(equations.size()); + for (int eq = 0; eq < (int) equations.size(); eq++) + if (!is_equation_linear[eq]) + { + BinaryOpNode *eq_node = equations[eq]; + expr_t lhs = eq_node->get_arg1(); + result.clear(); + lhs->collectDynamicVariables(SymbolType::endogenous, result); + if (result.size() == 1 && result.begin()->second == 0) + { + //check if the endogenous variable has not been already used in an other match ! + vector::iterator it = find(endo2eq.begin(), endo2eq.end(), result.begin()->first); + if (it == endo2eq.end()) + endo2eq[result.begin()->first] = eq; + else + { + bool_result = false; + break; + } + } + } + return bool_result; +} + void ModelTree::computePrologueAndEpilogue(const jacob_map_t &static_jacobian_arg, vector &equation_reordered, vector &variable_reordered) { @@ -792,7 +891,7 @@ ModelTree::printBlockDecomposition(const vector> &blocks) const } block_type_firstequation_size_mfs_t -ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_jacobian, vector> &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector &variable_reordered, const vector &equation_reordered, vector &n_static, vector &n_forward, vector &n_backward, vector &n_mixed, vector, pair>> &block_col_type) +ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_jacobian, vector> &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector &variable_reordered, const vector &equation_reordered, vector &n_static, vector &n_forward, vector &n_backward, vector &n_mixed, vector, pair>> &block_col_type, bool linear_decomposition) { int i = 0; int count_equ = 0, blck_count_simult = 0; @@ -836,14 +935,27 @@ ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_j int curr_variable = it.first; int curr_lag = it.second; auto it1 = find(variable_reordered.begin()+first_count_equ, variable_reordered.begin()+(first_count_equ+Blck_Size), curr_variable); - if (it1 != variable_reordered.begin()+(first_count_equ+Blck_Size)) - if (dynamic_jacobian.find({ curr_lag, { equation_reordered[count_equ], curr_variable } }) != dynamic_jacobian.end()) - { - if (curr_lag > Lead) - Lead = curr_lag; - else if (-curr_lag > Lag) - Lag = -curr_lag; - } + if (linear_decomposition) + { + if (dynamic_jacobian.find(make_pair(curr_lag, make_pair(equation_reordered[count_equ], curr_variable))) != dynamic_jacobian.end()) + { + if (curr_lag > Lead) + Lead = curr_lag; + else if (-curr_lag > Lag) + Lag = -curr_lag; + } + } + else + { + if (it1 != variable_reordered.begin()+(first_count_equ+Blck_Size)) + if (dynamic_jacobian.find({ curr_lag, { equation_reordered[count_equ], curr_variable } }) != dynamic_jacobian.end()) + { + if (curr_lag > Lead) + Lead = curr_lag; + else if (-curr_lag > Lag) + Lag = -curr_lag; + } + } } } if ((Lag > 0) && (Lead > 0)) @@ -939,6 +1051,24 @@ ModelTree::reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_j return (block_type_size_mfs); } +vector +ModelTree::equationLinear(map >, expr_t> first_order_endo_derivatives) const +{ + vector is_linear(symbol_table.endo_nbr(), true); + for (map >, expr_t>::const_iterator it = first_order_endo_derivatives.begin(); it != first_order_endo_derivatives.end(); it++) + { + expr_t Id = it->second; + set > endogenous; + Id->collectEndogenous(endogenous); + if (endogenous.size() > 0) + { + int eq = it->first.first; + is_linear[eq] = false; + } + } + return is_linear; +} + vector ModelTree::BlockLinear(const blocks_derivatives_t &blocks_derivatives, const vector &variable_reordered) const { diff --git a/src/ModelTree.hh b/src/ModelTree.hh index 92adc710..44ea9790 100644 --- a/src/ModelTree.hh +++ b/src/ModelTree.hh @@ -171,6 +171,9 @@ protected: //! the file containing the model and the derivatives code ofstream code_file; + + //! Vector indicating if the equation is linear in endogenous variable (true) or not (false) + vector is_equation_linear; //! Computes 1st derivatives /*! \param vars the derivation IDs w.r. to which compute the derivatives */ @@ -250,11 +253,17 @@ protected: If no matching is found with a zero cutoff close to zero an error message is printout. */ void computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian, double cutoff, jacob_map_t &static_jacobian, dynamic_jacob_map_t &dynamic_jacobian); - + //! Try to find a natural normalization if all equations are matched to an endogenous variable on the LHS + bool computeNaturalNormalization(); //! Try to normalized each unnormalized equation (matched endogenous variable only on the LHS) void computeNormalizedEquations(multimap &endo2eqs) const; //! Evaluate the jacobian and suppress all the elements below the cutoff void evaluateAndReduceJacobian(const eval_context_t &eval_context, jacob_map_t &contemporaneous_jacobian, jacob_map_t &static_jacobian, dynamic_jacob_map_t &dynamic_jacobian, double cutoff, bool verbose); + //! Select and reorder the non linear equations of the model + vector > select_non_linear_equations_and_variables(vector is_equation_linear, const dynamic_jacob_map_t &dynamic_jacobian, vector &equation_reordered, vector &variable_reordered, + vector &inv_equation_reordered, vector &inv_variable_reordered, + lag_lead_vector_t &equation_lag_lead, lag_lead_vector_t &variable_lag_lead, + vector &n_static, vector &n_forward, vector &n_backward, vector &n_mixed); //! Search the equations and variables belonging to the prologue and the epilogue of the model void computePrologueAndEpilogue(const jacob_map_t &static_jacobian, vector &equation_reordered, vector &variable_reordered); //! Determine the type of each equation of model and try to normalized the unnormalized equation using computeNormalizedEquations @@ -262,9 +271,11 @@ protected: //! Compute the block decomposition and for a non-recusive block find the minimum feedback set void computeBlockDecompositionAndFeedbackVariablesForEachBlock(const jacob_map_t &static_jacobian, const dynamic_jacob_map_t &dynamic_jacobian, vector &equation_reordered, vector &variable_reordered, vector> &blocks, const equation_type_and_normalized_equation_t &Equation_Type, bool verbose_, bool select_feedback_variable, int mfs, vector &inv_equation_reordered, vector &inv_variable_reordered, lag_lead_vector_t &equation_lag_lead, lag_lead_vector_t &variable_lag_lead_t, vector &n_static, vector &n_forward, vector &n_backward, vector &n_mixed) const; //! Reduce the number of block merging the same type equation in the prologue and the epilogue and determine the type of each block - block_type_firstequation_size_mfs_t reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_jacobian, vector> &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector &variable_reordered, const vector &equation_reordered, vector &n_static, vector &n_forward, vector &n_backward, vector &n_mixed, vector, pair>> &block_col_type); + block_type_firstequation_size_mfs_t reduceBlocksAndTypeDetermination(const dynamic_jacob_map_t &dynamic_jacobian, vector> &blocks, const equation_type_and_normalized_equation_t &Equation_Type, const vector &variable_reordered, const vector &equation_reordered, vector &n_static, vector &n_forward, vector &n_backward, vector &n_mixed, vector, pair>> &block_col_type, bool linear_decomposition); //! Determine the maximum number of lead and lag for the endogenous variable in a bloc void getVariableLeadLagByBlock(const dynamic_jacob_map_t &dynamic_jacobian, const vector &components_set, int nb_blck_sim, lag_lead_vector_t &equation_lead_lag, lag_lead_vector_t &variable_lead_lag, const vector &equation_reordered, const vector &variable_reordered) const; + //! For each equation determine if it is linear or not + vector equationLinear(map >, expr_t> first_order_endo_derivatives) const; //! Print an abstract of the block structure of the model void printBlockDecomposition(const vector> &blocks) const; //! Determine for each block if it is linear or not diff --git a/src/ParsingDriver.cc b/src/ParsingDriver.cc index c6046e4b..9652260d 100644 --- a/src/ParsingDriver.cc +++ b/src/ParsingDriver.cc @@ -771,6 +771,12 @@ ParsingDriver::block() mod_file->block = true; } +void +ParsingDriver::linear_decomposition() +{ + mod_file->linear_decomposition = true; +} + void ParsingDriver::no_static() { @@ -2405,6 +2411,23 @@ ParsingDriver::conditional_forecast_paths() det_shocks.clear(); } +void +ParsingDriver::det_cond_forecast_linear_decomposition(const string & plan) +{ + symbol_list.clear(); + symbol_list.addSymbol(plan); + mod_file->addStatement(make_unique(symbol_list, options_list, mod_file->linear_decomposition)); +} + +void +ParsingDriver::det_cond_forecast_linear_decomposition(const string &plan, const string &dset) +{ + symbol_list.clear(); + symbol_list.addSymbol(plan); + symbol_list.addSymbol(dset); + mod_file->addStatement(make_unique(symbol_list, options_list, mod_file->linear_decomposition)); +} + void ParsingDriver::calib_smoother() { diff --git a/src/ParsingDriver.hh b/src/ParsingDriver.hh index 439ebed7..8ff2dac4 100644 --- a/src/ParsingDriver.hh +++ b/src/ParsingDriver.hh @@ -323,6 +323,9 @@ public: void use_dll(); //! the modelis block decomposed void block(); + //! the model is decomposed according to the linearity of its equations + void linear_decomposition(); + //! the model is stored in a binary file void byte_code(); //! the static model is not computed @@ -673,6 +676,9 @@ public: void conditional_forecast_paths(); //! Plot conditional forecast statement void plot_conditional_forecast(const string &periods = ""); + //! Deterministic conditional forecast statement + void det_cond_forecast_linear_decomposition(const string &plan); + void det_cond_forecast_linear_decomposition(const string &plan, const string &dset); //! Smoother on calibrated models void calib_smoother(); //! Extended path diff --git a/src/StaticModel.cc b/src/StaticModel.cc index 22870c3f..788f53e4 100644 --- a/src/StaticModel.cc +++ b/src/StaticModel.cc @@ -1123,7 +1123,7 @@ StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms computeBlockDecompositionAndFeedbackVariablesForEachBlock(static_jacobian, dynamic_jacobian, equation_reordered, variable_reordered, blocks, equation_type_and_normalized_equation, false, false, mfs, inv_equation_reordered, inv_variable_reordered, equation_lag_lead, variable_lag_lead, n_static, n_forward, n_backward, n_mixed); - block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(dynamic_jacobian, blocks, equation_type_and_normalized_equation, variable_reordered, equation_reordered, n_static, n_forward, n_backward, n_mixed, block_col_type); + block_type_firstequation_size_mfs = reduceBlocksAndTypeDetermination(dynamic_jacobian, blocks, equation_type_and_normalized_equation, variable_reordered, equation_reordered, n_static, n_forward, n_backward, n_mixed, block_col_type, false); printBlockDecomposition(blocks);