diff --git a/DynamicModel.cc b/DynamicModel.cc index cf19468d..896fa47c 100644 --- a/DynamicModel.cc +++ b/DynamicModel.cc @@ -24,6 +24,7 @@ #include #include #include +#include #include "DynamicModel.hh" // For mkdir() and chdir() @@ -579,6 +580,8 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const int lag = it->first.first; unsigned int var = it->first.second.first; unsigned int eq = it->first.second.second; + int eqr = getBlockEquationID(block, eq); + int varr = getBlockVariableID(block, var); if (var != prev_var || lag != prev_lag) { prev_var = var; @@ -590,10 +593,10 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const output << " g1(" << eq+1 << ", " << count_col << ") = "; id->writeOutput(output, local_output_type, local_temporary_terms); - output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, var)) + output << "; % variable=" << symbol_table.getName(symbol_table.getID(eEndogenous, varr)) << "(" << lag - << ") " << var+1 - << ", equation=" << eq+1 << endl; + << ") " << varr+1 << ", " << var+1 + << ", equation=" << eqr+1 << ", " << eq+1 << endl; } prev_var = 999999999; prev_lag = -9999999; @@ -2278,8 +2281,63 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de model at a given period. */ + vector state_var; output << "M_.lead_lag_incidence = ["; // Loop on endogenous variables + int nstatic = 0, + nfwrd = 0, + npred = 0, + nboth = 0; + for (int endoID = 0; endoID < symbol_table.endo_nbr(); endoID++) + { + output << endl; + int sstatic = 1, + sfwrd = 0, + spred = 0, + sboth = 0; + // Loop on periods + for (int lag = -max_endo_lag; lag <= max_endo_lead; lag++) + { + // Print variableID if exists with current period, otherwise print 0 + try + { + int varID = getDerivID(symbol_table.getID(eEndogenous, endoID), lag); + output << " " << getDynJacobianCol(varID) + 1; + if (lag == -1) + { + sstatic = 0; + spred = 1; + } + else if (lag == 1) + { + if (spred == 1) + { + sboth = 1; + spred = 0; + } + else + { + sstatic = 0; + sfwrd = 1; + } + } + } + catch (UnknownDerivIDException &e) + { + output << " 0"; + } + } + nstatic += sstatic; + nfwrd += sfwrd; + npred += spred; + nboth += sboth; + output << ";"; + } + output << "]';" << endl; + output << "M_.nstatic = " << nstatic << ";" << endl + << "M_.nfwrd = " << nfwrd << ";" << endl + << "M_.npred = " << npred << ";" << endl + << "M_.nboth = " << nboth << ";" << endl; for (int endoID = 0; endoID < symbol_table.endo_nbr(); endoID++) { output << endl; @@ -2289,17 +2347,15 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de // Print variableID if exists with current period, otherwise print 0 try { - int varID = getDerivID(symbol_table.getID(eEndogenous, endoID), lag); - output << " " << getDynJacobianCol(varID) + 1; + int varID = getDerivID(variable_reordered[symbol_table.getID(eEndogenous, endoID)], lag); + if (lag < 0 && find(state_var.begin(), state_var.end(), variable_reordered[symbol_table.getID(eEndogenous, endoID)]+1) == state_var.end()) + state_var.push_back(variable_reordered[symbol_table.getID(eEndogenous, endoID)]+1); } catch (UnknownDerivIDException &e) { - output << " 0"; } } - output << ";"; } - output << "]';" << endl; // Write equation tags output << "M_.equations_tags = {" << endl; @@ -2396,6 +2452,47 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de i++; } output << "];\n"; + output << "block_structure.block(" << block+1 << ").other_endogenous_block = ["; + i = 0; + for (set::iterator it_other_endogenous = other_endogenous.begin(); it_other_endogenous != other_endogenous.end(); it_other_endogenous++) + if (*it_other_endogenous >= 0) + { + bool OK = true; + unsigned int j; + for (j = 0; j < block && OK; j++) + for (unsigned int k = 0; k < getBlockSize(j) && OK; k++) + { + //printf("*it_other_endogenous=%d, getBlockVariableID(%d, %d)=%d\n",*it_other_endogenous, j, k, getBlockVariableID(j, k)); + OK = *it_other_endogenous != getBlockVariableID(j, k); + } + if (!OK) + output << " " << j; + i++; + } + output << "];\n"; + + //vector inter_state_var; + output << "block_structure.block(" << block+1 << ").tm1 = zeros(" << i << ", " << state_var.size() << ");\n"; + int count_other_endogenous = 1; + for (set::const_iterator it_other_endogenous = other_endogenous.begin(); it_other_endogenous != other_endogenous.end(); it_other_endogenous++) + { + for (vector::const_iterator it=state_var.begin(); it != state_var.end(); it++) + { + //cout << "block = " << block+1 << " state_var = " << *it << " it_other_endogenous=" << *it_other_endogenous + 1 << "\n"; + if (*it == *it_other_endogenous + 1) + { + output << "block_structure.block(" << block+1 << ").tm1(" + << count_other_endogenous << ", " + << it - state_var.begin()+1 << ") = 1;\n"; + /*output << "block_structure.block(" << block+1 << ").tm1(" + << it - state_var.begin()+1 << ", " + << count_other_endogenous << ") = 1;\n";*/ + //cout << "=>\n"; + } + } + count_other_endogenous++; + } + output << "block_structure.block(" << block+1 << ").other_endo_nbr = " << i << ";\n"; tmp_s.str(""); @@ -2405,13 +2502,16 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de reordered_dynamic_jacobian[make_pair(it->second.first, make_pair(it->first.second, it->first.first))] = it->second.second; output << "block_structure.block(" << block+1 << ").lead_lag_incidence = [];\n"; int last_var = -1; - for (int lag = -max_lag_endo; lag < max_lead_endo+1; lag++) + vector local_state_var; + for (int lag = -1; lag < 1+1; lag++) { last_var = -1; for (dynamic_jacob_map_t::const_iterator it = reordered_dynamic_jacobian.begin(); it != reordered_dynamic_jacobian.end(); it++) { if (lag == it->first.first && last_var != it->first.second.first) { + if (lag == -1) + local_state_var.push_back(getBlockVariableID(block, it->first.second.first)+1); count_lead_lag_incidence++; for (int i = last_var; i < it->first.second.first-1; i++) tmp_s << " 0"; @@ -2426,6 +2526,15 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de output << "block_structure.block(" << block+1 << ").lead_lag_incidence = [ block_structure.block(" << block+1 << ").lead_lag_incidence; " << tmp_s.str() << "]; %lag = " << lag << "\n"; tmp_s.str(""); } + vector inter_state_var; + for (vector::const_iterator it_l=local_state_var.begin(); it_l != local_state_var.end(); it_l++) + for (vector::const_iterator it=state_var.begin(); it != state_var.end(); it++) + if (*it == *it_l) + inter_state_var.push_back(it - state_var.begin()+1); + output << "block_structure.block(" << block+1 << ").sorted_col_dr_ghx = ["; + for (vector::const_iterator it=inter_state_var.begin(); it != inter_state_var.end(); it++) + output << *it << " "; + output << "];\n"; count_lead_lag_incidence = 0; output << "block_structure.block(" << block+1 << ").lead_lag_incidence_other = [];\n"; for (int lag = -1; lag <= 1; lag++) @@ -2497,6 +2606,10 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de output << "];\n"; } // Writing initialization for some other variables + output << "M_.state_var = ["; + for (vector::const_iterator it=state_var.begin(); it != state_var.end(); it++) + output << *it << " "; + output << "];" << endl; output << "M_.exo_names_orig_ord = [1:" << symbol_table.exo_nbr() << "];" << endl << "M_.maximum_lag = " << max_lag << ";" << endl << "M_.maximum_lead = " << max_lead << ";" << endl; @@ -2677,10 +2790,21 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative map >, pair >, int> DynamicModel::get_Derivatives(int block) { + int max_lag, max_lead; map >, pair >, int> Derivatives; Derivatives.clear(); - int max_lag = getBlockMaxLag(block); - int max_lead = getBlockMaxLead(block); + BlockSimulationType simulation_type = getBlockSimulationType(block); + if (simulation_type == EVALUATE_BACKWARD || simulation_type == EVALUATE_FORWARD) + { + max_lag = 1; + max_lead = 1; + setBlockLeadLag(block, max_lag, max_lead); + } + else + { + max_lag = getBlockMaxLag(block); + max_lead = getBlockMaxLead(block); + } int block_size = getBlockSize(block); int block_nb_recursive = block_size - getBlockMfs(block); for (int lag = -max_lag; lag <= max_lead; lag++) @@ -3852,3 +3976,9 @@ DynamicModel::fillEvalContext(eval_context_t &eval_context) const it != trendVars.end(); it++) eval_context[*it] = 2; //not <= 0 bc of log, not 1 bc of powers } + +void +DynamicModel::set_cutoff_to_zero() +{ + cutoff = 0; +} diff --git a/DynamicModel.hh b/DynamicModel.hh index 36606d86..9e738216 100644 --- a/DynamicModel.hh +++ b/DynamicModel.hh @@ -338,6 +338,8 @@ public: //! Fills eval context with values of model local variables and auxiliary variables void fillEvalContext(eval_context_t &eval_context) const; + void set_cutoff_to_zero(); + //! Return the number of blocks virtual unsigned int getNbBlocks() const diff --git a/ModFile.cc b/ModFile.cc index 25fe739d..33c964ca 100644 --- a/ModFile.cc +++ b/ModFile.cc @@ -359,6 +359,10 @@ ModFile::computingPass(bool no_tmp_terms) if (!no_static) { static_model.initializeVariablesAndEquations(); + if (mod_file_struct.stoch_simul_present + || mod_file_struct.estimation_present || mod_file_struct.osr_present + || mod_file_struct.ramsey_policy_present || mod_file_struct.identification_present) + static_model.set_cutoff_to_zero(); static_model.computingPass(global_eval_context, no_tmp_terms, false, block, byte_code); } // Set things to compute for dynamic model @@ -372,6 +376,10 @@ ModFile::computingPass(bool no_tmp_terms) dynamic_model.computingPass(true, false, false, false, global_eval_context, no_tmp_terms, block, use_dll, byte_code); else { + if (mod_file_struct.stoch_simul_present + || mod_file_struct.estimation_present || mod_file_struct.osr_present + || mod_file_struct.ramsey_policy_present || mod_file_struct.identification_present) + dynamic_model.set_cutoff_to_zero(); if (mod_file_struct.order_option < 1 || mod_file_struct.order_option > 3) { cerr << "ERROR: Incorrect order option..." << endl; diff --git a/ModelTree.hh b/ModelTree.hh index 72af80f9..346212df 100644 --- a/ModelTree.hh +++ b/ModelTree.hh @@ -201,6 +201,11 @@ protected: virtual unsigned int getBlockMaxLag(int block_number) const = 0; //! Return the maximum lead in a block virtual unsigned int getBlockMaxLead(int block_number) const = 0; + inline void setBlockLeadLag(int block, int max_lag, int max_lead) + { + block_lag_lead[block] = make_pair(max_lag, max_lead); + }; + //! Return the type of equation (equation_number) belonging to the block block_number virtual EquationType getBlockEquationType(int block_number, int equation_number) const = 0; //! Return true if the equation has been normalized diff --git a/StaticModel.cc b/StaticModel.cc index 93979c94..2488a3cd 100644 --- a/StaticModel.cc +++ b/StaticModel.cc @@ -1609,3 +1609,9 @@ StaticModel::writeAuxVarInitval(ostream &output, ExprNodeOutputType output_type) output << ";" << endl; } } + +void +StaticModel::set_cutoff_to_zero() +{ + cutoff = 0; +} diff --git a/StaticModel.hh b/StaticModel.hh index c8b7fe76..f40b552e 100644 --- a/StaticModel.hh +++ b/StaticModel.hh @@ -208,6 +208,8 @@ public: virtual int getDerivID(int symb_id, int lag) const throw (UnknownDerivIDException); + void set_cutoff_to_zero(); + //! Return the number of blocks virtual unsigned int getNbBlocks() const