From fd323a8217b4041076b8c8a53629471e2eead6c2 Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Sat, 15 Mar 2014 21:07:09 +0100 Subject: [PATCH 01/37] extended-preprocessor: starting a C interface (not finished) --- others/C/dynare_driver.h | 32 +++++++++ preprocessor/DynamicModel.cc | 105 ++++++++++++++++++++++++++++-- preprocessor/DynamicModel.hh | 8 ++- preprocessor/ModFile.cc | 108 ++++++++++++++++++++++++++----- preprocessor/ModFile.hh | 2 +- preprocessor/SteadyStateModel.cc | 4 +- preprocessor/SteadyStateModel.hh | 2 +- preprocessor/SymbolTable.cc | 99 +++++++++++++++++++++++++++- preprocessor/SymbolTable.hh | 2 + 9 files changed, 332 insertions(+), 30 deletions(-) create mode 100644 others/C/dynare_driver.h diff --git a/others/C/dynare_driver.h b/others/C/dynare_driver.h new file mode 100644 index 000000000..dfaa7c60a --- /dev/null +++ b/others/C/dynare_driver.h @@ -0,0 +1,32 @@ +/* + * Copyright (C) 2011-2012 Houtan Bastani, Daniel Waggoner, Tao Zha + * + * This is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This code is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * More information is available at . + */ + +#ifndef _DYNARE_C_DRIVER_H +#define _DYNARE_C_DRIVER_H + +#include +#include +#include +#include +#include + +using namespace std; + +struct aux_vars_t { + int endo_index, type, orig_index, orig_lead_lag; +} ; + +#endif // ! _DYNARE_C_DRIVER_H diff --git a/preprocessor/DynamicModel.cc b/preprocessor/DynamicModel.cc index bd2adaede..0c6e801b8 100644 --- a/preprocessor/DynamicModel.cc +++ b/preprocessor/DynamicModel.cc @@ -2958,6 +2958,99 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de void DynamicModel::writeCOutput(ostream &output, const string &basename, bool block_decomposition, bool byte_code, bool use_dll, int order, bool estimation_present) const +{ + int lag_presence[3]; + // Loop on endogenous variables + vector zeta_back, zeta_mixed, zeta_fwrd, zeta_static; + for (int endoID = 0; endoID < symbol_table.endo_nbr(); endoID++) + { + int varID; + // Loop on periods + for (int lag = 0; lag <= 2; lag++) + { + lag_presence[lag] = 1; + try + { + varID = getDerivID(symbol_table.getID(eEndogenous, endoID), lag-1); + } + catch (UnknownDerivIDException &e) + { + lag_presence[lag] = 0; + } + } + if (lag_presence[0] == 1) + if (lag_presence[2] == 1) + zeta_mixed.push_back(endoID); + else + zeta_back.push_back(endoID); + else if (lag_presence[2] == 1) + zeta_fwrd.push_back(endoID); + else + zeta_static.push_back(endoID); + + } + output << "nstatic = " << zeta_static.size() << ";" << endl + << "nfwrd = " << zeta_fwrd.size() << ";" << endl + << "nback = " << zeta_back.size() << ";" << endl + << "nmixed = " << zeta_mixed.size() << ";" << endl; + output << "zeta_static[" << zeta_static.size() << "] = {"; + for (vector::iterator i = zeta_static.begin(); i != zeta_static.end(); ++i) + { + if ( i != zeta_static.begin() ) + output << ","; + output << *i; + } + output << "};" << endl; + + output << "zeta_back[" << zeta_back.size() << "] = {"; + for (vector::iterator i = zeta_back.begin(); i != zeta_back.end(); ++i) + { + if ( i != zeta_back.begin() ) + output << ","; + output << *i; + } + output << "};" << endl; + + output << "zeta_fwrd[" << zeta_fwrd.size() << "] = {"; + for (vector::iterator i = zeta_fwrd.begin(); i != zeta_fwrd.end(); ++i) + { + if ( i != zeta_fwrd.begin() ) + output << ","; + output << *i; + } + output << "};" << endl; + + output << "zeta_mixed[" << zeta_mixed.size() << "] = {"; + for (vector::iterator i = zeta_mixed.begin(); i != zeta_mixed.end(); ++i) + { + if ( i != zeta_mixed.begin() ) + output << ","; + output << *i; + } + output << "};" << endl; + + // Write number of non-zero derivatives + // Use -1 if the derivatives have not been computed + output << "int *NNZDerivatives[3] = {"; + switch (order) + { + case 0: + output << NNZDerivativs[0] << ",-1,-1};" << endl; + break; + case 1: + output << NNZDerivativs[0] << "," << NNZDerivatives[1] << ",-1};" << endl; + break; + case 2: + output << NNZDerivativs[0] << "," << NNZDerivatives[1] << "," << NNZDerivatives[2] << "};" << endl; + break; + default: + cerr << "Order larger than 3 not implemented" << endl; + exit(EXIT_FAILURE); + } +} + +void +DynamicModel::writeCCOutput(ostream &output, const string &basename, bool block_decomposition, bool byte_code, bool use_dll, int order, bool estimation_present) const { int lag_presence[3]; // Loop on endogenous variables @@ -4339,9 +4432,9 @@ DynamicModel::dynamicOnlyEquationsNbr() const } void -DynamicModel::writeFirstDerivativesCC(const string &basename, bool cuda) const +DynamicModel::writeFirstDerivativesC(const string &basename, bool cuda) const { - string filename = basename + "_first_derivatives.cc"; + string filename = basename + "_first_derivatives.c"; ofstream mDynamicModelFile, mDynamicMexFile; mDynamicModelFile.open(filename.c_str(), ios::out | ios::binary); @@ -4396,10 +4489,10 @@ DynamicModel::writeFirstDerivativesCC(const string &basename, bool cuda) const // using compressed sparse row format (CSR) void -DynamicModel::writeSecondDerivativesCC_csr(const string &basename, bool cuda) const +DynamicModel::writeSecondDerivativesC_csr(const string &basename, bool cuda) const { - string filename = basename + "_second_derivatives.cc"; + string filename = basename + "_second_derivatives.c"; ofstream mDynamicModelFile, mDynamicMexFile; mDynamicModelFile.open(filename.c_str(), ios::out | ios::binary); @@ -4490,9 +4583,9 @@ DynamicModel::writeSecondDerivativesCC_csr(const string &basename, bool cuda) co } void -DynamicModel::writeThirdDerivativesCC_csr(const string &basename, bool cuda) const +DynamicModel::writeThirdDerivativesC_csr(const string &basename, bool cuda) const { - string filename = basename + "_third_derivatives.cc"; + string filename = basename + "_third_derivatives.c"; ofstream mDynamicModelFile, mDynamicMexFile; mDynamicModelFile.open(filename.c_str(), ios::out | ios::binary); diff --git a/preprocessor/DynamicModel.hh b/preprocessor/DynamicModel.hh index 8de6bd392..e1ee490bb 100644 --- a/preprocessor/DynamicModel.hh +++ b/preprocessor/DynamicModel.hh @@ -215,6 +215,8 @@ public: void writeOutput(ostream &output, const string &basename, bool block, bool byte_code, bool use_dll, int order, bool estimation_present) const; //! Writes model initialization and lead/lag incidence matrix to C output void writeCOutput(ostream &output, const string &basename, bool block, bool byte_code, bool use_dll, int order, bool estimation_present) const; + //! Writes model initialization and lead/lag incidence matrix to Cpp output + void writeCCOutput(ostream &output, const string &basename, bool block, bool byte_code, bool use_dll, int order, bool estimation_present) const; //! Adds informations for simulation in a binary file void Write_Inf_To_Bin_File_Block(const string &dynamic_basename, const string &bin_basename, @@ -224,11 +226,11 @@ public: //! Writes file containing parameters derivatives void writeParamsDerivativesFile(const string &basename) const; //! Writes CC file containing first order derivatives of model evaluated at steady state - void writeFirstDerivativesCC(const string &basename, bool cuda) const; + void writeFirstDerivativesC(const string &basename, bool cuda) const; //! Writes CC file containing second order derivatives of model evaluated at steady state (compressed sparse column) - void writeSecondDerivativesCC_csr(const string &basename, bool cuda) const; + void writeSecondDerivativesC_csr(const string &basename, bool cuda) const; //! Writes CC file containing third order derivatives of model evaluated at steady state (compressed sparse column) - void writeThirdDerivativesCC_csr(const string &basename, bool cuda) const; + void writeThirdDerivativesC_csr(const string &basename, bool cuda) const; //! Converts to static model (only the equations) /*! It assumes that the static model given in argument has just been allocated */ void toStatic(StaticModel &static_model) const; diff --git a/preprocessor/ModFile.cc b/preprocessor/ModFile.cc index 98ee8d7d4..ece9a709e 100644 --- a/preprocessor/ModFile.cc +++ b/preprocessor/ModFile.cc @@ -805,9 +805,9 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool no_log, b } void -ModFile::writeModelCC(const string &basename, bool cuda) const +ModFile::writeModelC(const string &basename, bool cuda) const { - string filename = basename + ".cc"; + string filename = basename + ".c"; ofstream mDriverCFile; mDriverCFile.open(filename.c_str(), ios::out | ios::binary); @@ -818,15 +818,15 @@ ModFile::writeModelCC(const string &basename, bool cuda) const } mDriverCFile << "/*" << endl - << " * " << filename << " : Driver file for MS-DSGE code" << endl + << " * " << filename << " : Driver file for Dynare C code" << endl << " *" << endl << " * Warning : this file is generated automatically by Dynare" << endl << " * from model file (.mod)" << endl << " */" << endl << endl - << "#include \"dynare_cpp_driver.hh\"" << endl + << "#include \"dynare_driver.h\"" << endl << endl - << "DynareInfo::DynareInfo(void)" << endl + << "struct" << endl << "{" << endl; // Write basic info @@ -846,7 +846,7 @@ ModFile::writeModelCC(const string &basename, bool cuda) const it != statements.end(); it++) (*it)->writeCOutput(mDriverCFile, basename); - mDriverCFile << "}" << endl; + mDriverCFile << "} DynareInfo;" << endl; mDriverCFile.close(); // Write informational m file @@ -876,6 +876,82 @@ ModFile::writeModelCC(const string &basename, bool cuda) const << "% Warning : this file is generated automatically by Dynare" << endl << "% from model file (.mod)" << endl << endl << "disp('The following C file was successfully created:');" << endl + << "ls preprocessorOutput.c" << endl << endl; + mOutputFile.close(); +} + +void +ModFile::writeModelCC(const string &basename, bool cuda) const +{ + string filename = basename + ".cc"; + + ofstream mDriverCFile; + mDriverCFile.open(filename.c_str(), ios::out | ios::binary); + if (!mDriverCFile.is_open()) + { + cerr << "Error: Can't open file " << filename << " for writing" << endl; + exit(EXIT_FAILURE); + } + + mDriverCFile << "/*" << endl + << " * " << filename << " : Driver file for Dynare C++ code" << endl + << " *" << endl + << " * Warning : this file is generated automatically by Dynare" << endl + << " * from model file (.mod)" << endl + << " */" << endl + << endl + << "#include \"dynare_cpp_driver.hh\"" << endl + << endl + << "DynareInfo::DynareInfo(void)" << endl + << "{" << endl; + + // Write basic info + symbol_table.writeCCOutput(mDriverCFile); + + mDriverCFile << endl << "params.resize(param_nbr);" << endl; + + if (dynamic_model.equation_number() > 0) + { + dynamic_model.writeCOutput(mDriverCFile, basename, block, byte_code, use_dll, mod_file_struct.order_option, mod_file_struct.estimation_present); + // if (!no_static) + // static_model.writeCOutput(mOutputFile, block); + } + + // Print statements + for (vector::const_iterator it = statements.begin(); + it != statements.end(); it++) + (*it)->writeCOutput(mDriverCFile, basename); + + mDriverCFile << "};" << endl; + mDriverCFile.close(); + + // Write informational m file + ofstream mOutputFile; + + if (basename.size()) + { + string fname(basename); + fname += ".m"; + mOutputFile.open(fname.c_str(), ios::out | ios::binary); + if (!mOutputFile.is_open()) + { + cerr << "ERROR: Can't open file " << fname + << " for writing" << endl; + exit(EXIT_FAILURE); + } + } + else + { + cerr << "ERROR: Missing file name" << endl; + exit(EXIT_FAILURE); + } + + mOutputFile << "%" << endl + << "% Status : informational m file" << endl + << "%" << endl + << "% Warning : this file is generated automatically by Dynare" << endl + << "% from model file (.mod)" << endl << endl + << "disp('The following C++ file was successfully created:');" << endl << "ls preprocessorOutput.cc" << endl << endl; mOutputFile.close(); } @@ -883,8 +959,8 @@ ModFile::writeModelCC(const string &basename, bool cuda) const void ModFile::writeExternalFiles(const string &basename, FileOutputType output, bool cuda) const { - writeModelCC(basename, cuda); - steady_state_model.writeSteadyStateFileCC(basename, mod_file_struct.ramsey_model_present, cuda); + writeModelC(basename, cuda); + steady_state_model.writeSteadyStateFileC(basename, mod_file_struct.ramsey_model_present, cuda); dynamic_model.writeDynamicFile(basename, block, byte_code, use_dll, mod_file_struct.order_option); @@ -893,18 +969,18 @@ ModFile::writeExternalFiles(const string &basename, FileOutputType output, bool // static_model.writeStaticCFile(basename, block, byte_code, use_dll); - // static_model.writeParamsDerivativesFileCC(basename, cuda); - // static_model.writeAuxVarInitvalCC(mOutputFile, oMatlabOutsideModel, cuda); + // static_model.writeParamsDerivativesFileC(basename, cuda); + // static_model.writeAuxVarInitvalC(mOutputFile, oMatlabOutsideModel, cuda); - // dynamic_model.writeResidualsCC(basename, cuda); - // dynamic_model.writeParamsDerivativesFileCC(basename, cuda); - dynamic_model.writeFirstDerivativesCC(basename, cuda); + // dynamic_model.writeResidualsC(basename, cuda); + // dynamic_model.writeParamsDerivativesFileC(basename, cuda); + dynamic_model.writeFirstDerivativesC(basename, cuda); if (output == second) - dynamic_model.writeSecondDerivativesCC_csr(basename, cuda); + dynamic_model.writeSecondDerivativesC_csr(basename, cuda); else if (output == third) { - dynamic_model.writeSecondDerivativesCC_csr(basename, cuda); - dynamic_model.writeThirdDerivativesCC_csr(basename, cuda); + dynamic_model.writeSecondDerivativesC_csr(basename, cuda); + dynamic_model.writeThirdDerivativesC_csr(basename, cuda); } } diff --git a/preprocessor/ModFile.hh b/preprocessor/ModFile.hh index aca8ef9d1..9885c6e6d 100644 --- a/preprocessor/ModFile.hh +++ b/preprocessor/ModFile.hh @@ -140,7 +140,7 @@ public: ) const; //! Writes C output files only => No further Matlab processing void writeCOutputFiles(const string &basename) const; - void writeModelCC(const string &basename, bool cuda) const; + void writeModelC(const string &basename, bool cuda) const; void writeExternalFiles(const string &basename, FileOutputType output, bool cuda) const; }; diff --git a/preprocessor/SteadyStateModel.cc b/preprocessor/SteadyStateModel.cc index f1069642b..1c590a98a 100644 --- a/preprocessor/SteadyStateModel.cc +++ b/preprocessor/SteadyStateModel.cc @@ -153,9 +153,9 @@ SteadyStateModel::writeSteadyStateFile(const string &basename, bool ramsey_model } void -SteadyStateModel::writeSteadyStateFileCC(const string &basename, bool ramsey_model, bool cuda) const +SteadyStateModel::writeSteadyStateFileC(const string &basename, bool ramsey_model, bool cuda) const { - string filename = basename + "_steadystate.cc"; + string filename = basename + "_steadystate.c"; ofstream output; output.open(filename.c_str(), ios::out | ios::binary); diff --git a/preprocessor/SteadyStateModel.hh b/preprocessor/SteadyStateModel.hh index 817f84f5a..f8bb051f0 100644 --- a/preprocessor/SteadyStateModel.hh +++ b/preprocessor/SteadyStateModel.hh @@ -49,7 +49,7 @@ public: \param[in] ramsey_model Is there a Ramsey model in the MOD file? If yes, then use the "ys" in argument of the steady state file as initial values */ void writeSteadyStateFile(const string &basename, bool ramsey_model) const; - void writeSteadyStateFileCC(const string &basename, bool ramsey_model, bool cuda) const; + void writeSteadyStateFileC(const string &basename, bool ramsey_model, bool cuda) const; }; #endif diff --git a/preprocessor/SymbolTable.cc b/preprocessor/SymbolTable.cc index 001526843..4e207cfe7 100644 --- a/preprocessor/SymbolTable.cc +++ b/preprocessor/SymbolTable.cc @@ -288,6 +288,103 @@ SymbolTable::writeOutput(ostream &output) const throw (NotYetFrozenException) void SymbolTable::writeCOutput(ostream &output) const throw (NotYetFrozenException) +{ + if (!frozen) + throw NotYetFrozenException(); + + output << endl + << "int exo_nbr = " << exo_nbr() << ";" << endl; + if (exo_nbr() > 0) + { + output << "char *exo_names[" << exo_nbr() << "];" << endl; + for (int id = 0; id < exo_nbr(); id++) + output << "exo_names[" << id << "] = \"" << getName(exo_ids[id]) << "\";" << endl; + } + + output << endl + << "int exo_det_nbr = " << exo_det_nbr() << ";" << endl; + if (exo_det_nbr() > 0) + { + output << "char *exo_det_names[" << exo_det_nbr() << "];" << endl; + for (int id = 0; id < exo_det_nbr(); id++) + output << "exo_det_names[" << id << "] = \"" << getName(exo_det_ids[id]) << "\";" << endl; + } + + output << endl + << "int endo_nbr = " << endo_nbr() << ";" << endl; + if (endo_nbr() > 0) + { + output << "char *endo_names[" << endo_nbr() << "];" << endl; + for (int id = 0; id < endo_nbr(); id++) + output << "endo_names[" << id << "] = \"" << getName(endo_ids[id]) << "\";" << endl; + } + + output << endl + << "int param_nbr = " << param_nbr() << ";" << endl; + if (param_nbr() > 0) + { + output << "char *param_names[" << param_nbr() << "];" << endl; + for (int id = 0; id < param_nbr(); id++) + output << "param_names[" << id << "] = \"" << getName(param_ids[id]) << "\";" << endl; + } + + // Write the auxiliary variable table + output << "int aux_var_nbr = " << aux_vars.size() << ";" << endl; + if (aux_vars.size() > 0) + { + output << "struct aux_vars_t *av[" << aux_vars.size() << "];" << endl; + for (int i = 0; i < (int) aux_vars.size(); i++) + { + output << "av[" << i << "].endo_index = " << getTypeSpecificID(aux_vars[i].get_symb_id()) << ";" << endl + << "av[" << i << "].type = " << aux_vars[i].get_type() << ";" << endl; + switch (aux_vars[i].get_type()) + { + case avEndoLead: + case avExoLead: + case avExpectation: + case avMultiplier: + case avDiffForward: + break; + case avEndoLag: + case avExoLag: + output << "av[" << i << "].orig_index = " << getTypeSpecificID(aux_vars[i].get_orig_symb_id()) << ";" << endl + << "av[" << i << "].orig_lead_lag = " << aux_vars[i].get_orig_lead_lag() << ";" << endl; + break; + } + } + } + + output << "int predeterminedNbr = " << predeterminedNbr() << ";" << endl; + if (predeterminedNbr() > 0) + { + output << "int predetermined_variables[" << predeterminedNbr() << "] = {"; + for (set::const_iterator it = predetermined_variables.begin(); + it != predetermined_variables.end(); it++) + { + if ( it != predetermined_variables.begin() ) + output << ","; + output << getTypeSpecificID(*it); + } + output << "};" << endl; + } + + output << "int observedVariablesNbr = " << observedVariablesNbr() << ";" << endl; + if (observedVariablesNbr() > 0) + { + output << "int varobs[" << observedVariablesNbr() << "] = {"; + for (vector::const_iterator it = varobs.begin(); + it != varobs.end(); it++) + { + if ( it != varobs.begin() ) + output << ","; + output << getTypeSpecificID(*it); + } + output << "};" << endl; + } +} + +void +SymbolTable::writeCCOutput(ostream &output) const throw (NotYetFrozenException) { if (!frozen) throw NotYetFrozenException(); @@ -596,7 +693,7 @@ SymbolTable::getEndogenous() const bool SymbolTable::isAuxiliaryVariable(int symb_id) const { - for (int i = 0; i < aux_vars.size(); i++) + for (int i = 0; i < (int) aux_vars.size(); i++) if (aux_vars[i].get_symb_id() == symb_id) return true; return false; diff --git a/preprocessor/SymbolTable.hh b/preprocessor/SymbolTable.hh index 8e0903e5c..3e4c47e12 100644 --- a/preprocessor/SymbolTable.hh +++ b/preprocessor/SymbolTable.hh @@ -278,6 +278,8 @@ public: void writeOutput(ostream &output) const throw (NotYetFrozenException); //! Write C output of this class void writeCOutput(ostream &output) const throw (NotYetFrozenException); + //! Write CC output of this class + void writeCCOutput(ostream &output) const throw (NotYetFrozenException); //! Mark a symbol as predetermined variable void markPredetermined(int symb_id) throw (UnknownSymbolIDException, FrozenException); //! Test if a given symbol is a predetermined variable From 8f06a3134d59f6197c9acb23695b7f04bcad0246 Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Fri, 17 Jan 2014 21:49:04 +0100 Subject: [PATCH 02/37] added unscented nodes generator put nodes generators in setup_integration_nodes.m added stochastic_solvers to homotopic_steps --- matlab/ep/extended_path.m | 18 +++++-- matlab/ep/homotopic_steps.m | 7 ++- matlab/ep/setup_integration_nodes.m | 39 +++++++++++++++ ...tochastic_perfect_foresight_model_solver.m | 33 +----------- ...lve_stochastic_perfect_foresight_model_1.m | 50 ++++++++++--------- 5 files changed, 83 insertions(+), 64 deletions(-) create mode 100644 matlab/ep/setup_integration_nodes.m rename matlab/{ => ep}/setup_stochastic_perfect_foresight_model_solver.m (60%) diff --git a/matlab/ep/extended_path.m b/matlab/ep/extended_path.m index bc98ed25a..6b3ced3ad 100644 --- a/matlab/ep/extended_path.m +++ b/matlab/ep/extended_path.m @@ -35,8 +35,11 @@ global M_ options_ oo_ options_.verbosity = options_.ep.verbosity; verbosity = options_.ep.verbosity+options_.ep.debug; +% Set maximum number of iterations for the deterministic solver. +options_.simul.maxit = options_.ep.maxit; + % Prepare a structure needed by the matlab implementation of the perfect foresight model solver -pfm = setup_stochastic_perfect_foresight_model_solver(M_,options_,oo_,'Tensor-Gaussian-Quadrature'); +pfm = setup_stochastic_perfect_foresight_model_solver(M_,options_,oo_); exo_nbr = M_.exo_nbr; periods = options_.periods; @@ -53,8 +56,6 @@ if isempty(initial_conditions) end end -% Set maximum number of iterations for the deterministic solver. -options_.simul.maxit = options_.ep.maxit; % Set the number of periods for the perfect foresight model periods = options_.ep.periods; @@ -185,7 +186,7 @@ while (t3 if verbose diff --git a/matlab/ep/setup_integration_nodes.m b/matlab/ep/setup_integration_nodes.m new file mode 100644 index 000000000..da8e129b9 --- /dev/null +++ b/matlab/ep/setup_integration_nodes.m @@ -0,0 +1,39 @@ +function [nodes,weights,nnodes] = setup_integration_nodes(EpOptions,pfm) + if EpOptions.stochastic.order + % Compute weights and nodes for the stochastic version of the extended path. + switch EpOptions.IntegrationAlgorithm + case 'Tensor-Gaussian-Quadrature' + % Get the nodes and weights from a univariate Gauss-Hermite quadrature. + [nodes0,weights0] = gauss_hermite_weights_and_nodes(EpOptions.stochastic.quadrature.nodes); + % Replicate the univariate nodes for each innovation and dates, and, if needed, correlate them. + nodes0 = repmat(nodes0,1,pfm.number_of_shocks*pfm.stochastic_order)*kron(eye(pfm.stochastic_order),pfm.Omega); + % Put the nodes and weights in cells + for i=1:pfm.number_of_shocks + rr(i) = {nodes0(:,i)}; + ww(i) = {weights0}; + end + % Build the tensorial grid + nodes = cartesian_product_of_sets(rr{:}); + weights = prod(cartesian_product_of_sets(ww{:}),2); + nnodes = length(pfm.weights); + case 'Stroud-Cubature-3' + [nodes,weights] = cubature_with_gaussian_weight(pfm.number_of_shocks*pfm.stochastic_order,3,'Stroud') + nodes = kron(eye(pfm.stochastic_order),transpose(pfm.Omega))*nodes; + weights = weights; + nnodes = length(pfm.weights); + case 'Stroud-Cubature-5' + [nodes,weights] = cubature_with_gaussian_weight(pfm.number_of_shocks*pfm.stochastic_order,5,'Stroud') + nodes = kron(eye(pfm.stochastic_order),transpose(pfm.Omega))*nodes; + weights = weights; + nnodes = length(weights); + case 'UT_2p+1' + p = pfm.number_of_shocks; + k = EpOptions.ut.k; + C = sqrt(pfm.number_of_shocks + k)*pfm.Omega'; + nodes = [zeros(1,p); -C; C]; + weights = [k/(p+k); (1/(2*(p+k)))*ones(2*p,1)]; + nnodes = 2*p+1; + otherwise + error('Stochastic extended path:: Unknown integration algorithm!') + end +end diff --git a/matlab/setup_stochastic_perfect_foresight_model_solver.m b/matlab/ep/setup_stochastic_perfect_foresight_model_solver.m similarity index 60% rename from matlab/setup_stochastic_perfect_foresight_model_solver.m rename to matlab/ep/setup_stochastic_perfect_foresight_model_solver.m index 42cba7c8a..f6c455d32 100644 --- a/matlab/setup_stochastic_perfect_foresight_model_solver.m +++ b/matlab/ep/setup_stochastic_perfect_foresight_model_solver.m @@ -1,4 +1,4 @@ -function pfm = setup_stochastic_perfect_foresight_model_solver(DynareModel,DynareOptions,DynareOutput,IntegrationMethod) +function pfm = setup_stochastic_perfect_foresight_model_solver(DynareModel,DynareOptions,DynareOutput) % Copyright (C) 2013 Dynare Team % @@ -69,34 +69,3 @@ pfm.verbose = DynareOptions.ep.verbosity; pfm.maxit_ = DynareOptions.simul.maxit; pfm.tolerance = DynareOptions.dynatol.f; -if nargin>3 && DynareOptions.ep.stochastic.order - % Compute weights and nodes for the stochastic version of the extended path. - switch IntegrationMethod - case 'Tensor-Gaussian-Quadrature' - % Get the nodes and weights from a univariate Gauss-Hermite quadrature. - [nodes,weights] = gauss_hermite_weights_and_nodes(DynareOptions.ep.stochastic.quadrature.nodes); - % Replicate the univariate nodes for each innovation and dates, and, if needed, correlate them. - nodes = repmat(nodes,1,pfm.number_of_shocks*pfm.stochastic_order)*kron(eye(pfm.stochastic_order),pfm.Omega); - % Put the nodes and weights in cells - for i=1:pfm.number_of_shocks - rr(i) = {nodes(:,i)}; - ww(i) = {weights}; - end - % Build the tensorial grid - pfm.nodes = cartesian_product_of_sets(rr{:}); - pfm.weights = prod(cartesian_product_of_sets(ww{:}),2); - pfm.nnodes = length(pfm.weights); - case 'Stroud-Cubature-3' - [nodes,weights] = cubature_with_gaussian_weight(pfm.number_of_shocks*pfm.stochastic_order,3,'Stroud') - pfm.nodes = kron(eye(pfm.stochastic_order),transpose(Omega))*nodes; - pfm.weights = weights; - pfm.nnodes = length(pfm.weights); - case 'Stroud-Cubature-5' - [nodes,weights] = cubature_with_gaussian_weight(pfm.number_of_shocks*pfm.stochastic_order,5,'Stroud') - pfm.nodes = kron(eye(pfm.stochastic_order),transpose(Omega))*nodes; - pfm.weights = weights; - pfm.nnodes = length(weights); - otherwise - error('setup_stochastic_perfect_foresight_model_solver:: Unknown integration algorithm!') - end -end diff --git a/matlab/ep/solve_stochastic_perfect_foresight_model_1.m b/matlab/ep/solve_stochastic_perfect_foresight_model_1.m index 01d7bcacc..7c7c3e5aa 100644 --- a/matlab/ep/solve_stochastic_perfect_foresight_model_1.m +++ b/matlab/ep/solve_stochastic_perfect_foresight_model_1.m @@ -1,4 +1,4 @@ -function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo_simul,exo_simul,pfm,nnodes,order) +function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo_simul,exo_simul,EpOptions,pfm,order) % Copyright (C) 2012-2013 Dynare Team % @@ -35,6 +35,7 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo i_cols_T = nonzeros(lead_lag_incidence(1:2,:)'); hybrid_order = pfm.hybrid_order; dr = pfm.dr; + [nodes,weights,nnodes] = setup_integration_nodes(EpOptions,pfm); maxit = pfm.maxit_; tolerance = pfm.tolerance; @@ -42,32 +43,16 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo number_of_shocks = size(exo_simul,2); - [nodes,weights] = gauss_hermite_weights_and_nodes(nnodes); - % make sure that there is a node equal to zero % and permute nodes and weights to have zero first - k = find(abs(nodes) < 1e-12); + k = find(sum(abs(nodes),2) < 1e-12); if ~isempty(k) - nodes = [nodes(k); nodes(1:k-1); nodes(k+1:end)]; + nodes = [nodes(k,:); nodes(1:k-1,:); nodes(k+1:end,:)]; weights = [weights(k); weights(1:k-1); weights(k+1:end)]; else error('there is no nodes equal to zero') end - if number_of_shocks>1 - nodes = repmat(nodes,1,number_of_shocks)*chol(pfm.Sigma); - % to be fixed for Sigma ~= I - for i=1:number_of_shocks - rr(i) = {nodes(:,i)}; - ww(i) = {weights}; - end - nodes = cartesian_product_of_sets(rr{:}); - weights = prod(cartesian_product_of_sets(ww{:}),2); - nnodes = nnodes^number_of_shocks; - else - nodes = nodes*sqrt(pfm.Sigma); - end - if hybrid_order > 0 if hybrid_order == 2 h_correction = 0.5*dr.ghs2(dr.inv_order_var); @@ -90,7 +75,7 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo % The third row block is ny x nnodes^2 % and so on until size ny x nnodes^order world_nbr = 1+(nnodes-1)*order; - Y = repmat(endo_simul(:),1,world_nbr); + Y = repmat(endo_simul(:),1,world_nbr); % The columns of A map the elements of Y such that % each block of Y with ny rows are unfolded column wise @@ -110,8 +95,8 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo for i=2:periods k = n1:n2; for j=1:(1+(nnodes-1)*min(i-1,order)) - i_upd_r(i1:i2) = (n1:n2)+(j-1)*ny*periods; - i_upd_y(i1:i2) = (n1:n2)+ny+(j-1)*ny*(periods+2); + i_upd_r(i1:i2) = k+(j-1)*ny*periods; + i_upd_y(i1:i2) = k+ny+(j-1)*ny*(periods+2); i1 = i2+1; i2 = i2+ny; end @@ -174,7 +159,8 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo Y(i_cols_s,1); Y(i_cols_f,k1)]; end - [d1,jacobian] = dynamic_model(y,innovation,params,steady_state,i+1); + [d1,jacobian] = dynamic_model(y,innovation, ... + params,steady_state,i+1); if i == 1 % in first period we don't keep track of % predetermined variables @@ -248,6 +234,12 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo end end err = max(abs(res(i_upd_r))); + if verbose + [err1, k1] = max(abs(res)); + [err2, k2] = max(abs(err1)); + [err3, k3] = max(abs(err2)); + disp([iter err k1(:,k2(:,:,k3),k3) k2(:,:,k3) k3]) + end if err < tolerance stop = 1; if verbose @@ -264,7 +256,18 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo break end A2 = [nzA{:}]'; + if any(isnan(A2(:,3))) || any(any(any(isnan(res)))) + if verbose + disp(['solve_stochastic_foresight_model_1 encountered ' ... + 'NaN']) + end + flag = 1; + return + end A = [A1; sparse(A2(:,1),A2(:,2),A2(:,3),ny*(periods-order-1)*world_nbr,dimension)]; + if verbose + disp(sprintf('condest %g',condest(A))) + end dy = -A\res(i_upd_r); Y(i_upd_y) = Y(i_upd_y) + dy; end @@ -276,6 +279,7 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo fprintf('\n') ; disp(['WARNING : maximum number of iterations is reached (modify options_.simul.maxit).']) ; fprintf('\n') ; + disp(sprintf('err: %f',err)); end flag = 1;% more iterations are needed. endo_simul = 1; From f5ae90bcc82a7a30a8c70fc6503f781844f4f913 Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Sun, 16 Mar 2014 18:49:07 +0100 Subject: [PATCH 03/37] extended-path: changed homotopy so that integration nodes change too (only for for case 1) --- matlab/ep/extended_path.m | 3 + matlab/ep/homotopic_steps.m | 12 +++- matlab/ep/setup_integration_nodes.m | 70 +++++++++---------- ...lve_stochastic_perfect_foresight_model_1.m | 24 ++++--- 4 files changed, 64 insertions(+), 45 deletions(-) diff --git a/matlab/ep/extended_path.m b/matlab/ep/extended_path.m index 6b3ced3ad..3cf1a08c4 100644 --- a/matlab/ep/extended_path.m +++ b/matlab/ep/extended_path.m @@ -169,6 +169,9 @@ while (t3 @@ -152,7 +159,8 @@ if weight<1 solve_stochastic_perfect_foresight_model(endo_simul0,exo_simul0,pfm,options_.ep.stochastic.quadrature.nodes,options_.ep.stochastic.order); case 1 [flag,tmp] = ... - solve_stochastic_perfect_foresight_model_1(endo_simul0,exo_simul0,options_.ep,pfm,options_.ep.stochastic.order); + solve_stochastic_perfect_foresight_model_1(endo_simul0,exxo_simul,options_.ep,pfm,options_.ep.stochastic.order,weight); + % solve_stochastic_perfect_foresight_model_1(endo_simul0,exo_simul0,options_.ep,pfm,options_.ep.stochastic.order); end end end diff --git a/matlab/ep/setup_integration_nodes.m b/matlab/ep/setup_integration_nodes.m index da8e129b9..8907daac9 100644 --- a/matlab/ep/setup_integration_nodes.m +++ b/matlab/ep/setup_integration_nodes.m @@ -1,39 +1,39 @@ function [nodes,weights,nnodes] = setup_integration_nodes(EpOptions,pfm) if EpOptions.stochastic.order - % Compute weights and nodes for the stochastic version of the extended path. - switch EpOptions.IntegrationAlgorithm - case 'Tensor-Gaussian-Quadrature' - % Get the nodes and weights from a univariate Gauss-Hermite quadrature. - [nodes0,weights0] = gauss_hermite_weights_and_nodes(EpOptions.stochastic.quadrature.nodes); - % Replicate the univariate nodes for each innovation and dates, and, if needed, correlate them. - nodes0 = repmat(nodes0,1,pfm.number_of_shocks*pfm.stochastic_order)*kron(eye(pfm.stochastic_order),pfm.Omega); - % Put the nodes and weights in cells - for i=1:pfm.number_of_shocks - rr(i) = {nodes0(:,i)}; - ww(i) = {weights0}; + % Compute weights and nodes for the stochastic version of the extended path. + switch EpOptions.IntegrationAlgorithm + case 'Tensor-Gaussian-Quadrature' + % Get the nodes and weights from a univariate Gauss-Hermite quadrature. + [nodes0,weights0] = gauss_hermite_weights_and_nodes(EpOptions.stochastic.quadrature.nodes); + % Replicate the univariate nodes for each innovation and dates, and, if needed, correlate them. + nodes0 = repmat(nodes0,1,pfm.number_of_shocks*pfm.stochastic_order)*kron(eye(pfm.stochastic_order),pfm.Omega); + % Put the nodes and weights in cells + for i=1:pfm.number_of_shocks + rr(i) = {nodes0(:,i)}; + ww(i) = {weights0}; + end + % Build the tensorial grid + nodes = cartesian_product_of_sets(rr{:}); + weights = prod(cartesian_product_of_sets(ww{:}),2); + nnodes = length(weights); + case 'Stroud-Cubature-3' + [nodes,weights] = cubature_with_gaussian_weight(pfm.number_of_shocks*pfm.stochastic_order,3,'Stroud') + nodes = kron(eye(pfm.stochastic_order),transpose(pfm.Omega))*nodes; + weights = weights; + nnodes = length(weights); + case 'Stroud-Cubature-5' + [nodes,weights] = cubature_with_gaussian_weight(pfm.number_of_shocks*pfm.stochastic_order,5,'Stroud') + nodes = kron(eye(pfm.stochastic_order),transpose(pfm.Omega))*nodes; + weights = weights; + nnodes = length(weights); + case 'UT_2p+1' + p = pfm.number_of_shocks; + k = EpOptions.ut.k; + C = sqrt(pfm.number_of_shocks + k)*pfm.Omega'; + nodes = [zeros(1,p); -C; C]; + weights = [k/(p+k); (1/(2*(p+k)))*ones(2*p,1)]; + nnodes = 2*p+1; + otherwise + error('Stochastic extended path:: Unknown integration algorithm!') end - % Build the tensorial grid - nodes = cartesian_product_of_sets(rr{:}); - weights = prod(cartesian_product_of_sets(ww{:}),2); - nnodes = length(pfm.weights); - case 'Stroud-Cubature-3' - [nodes,weights] = cubature_with_gaussian_weight(pfm.number_of_shocks*pfm.stochastic_order,3,'Stroud') - nodes = kron(eye(pfm.stochastic_order),transpose(pfm.Omega))*nodes; - weights = weights; - nnodes = length(pfm.weights); - case 'Stroud-Cubature-5' - [nodes,weights] = cubature_with_gaussian_weight(pfm.number_of_shocks*pfm.stochastic_order,5,'Stroud') - nodes = kron(eye(pfm.stochastic_order),transpose(pfm.Omega))*nodes; - weights = weights; - nnodes = length(weights); - case 'UT_2p+1' - p = pfm.number_of_shocks; - k = EpOptions.ut.k; - C = sqrt(pfm.number_of_shocks + k)*pfm.Omega'; - nodes = [zeros(1,p); -C; C]; - weights = [k/(p+k); (1/(2*(p+k)))*ones(2*p,1)]; - nnodes = 2*p+1; - otherwise - error('Stochastic extended path:: Unknown integration algorithm!') end -end diff --git a/matlab/ep/solve_stochastic_perfect_foresight_model_1.m b/matlab/ep/solve_stochastic_perfect_foresight_model_1.m index 7c7c3e5aa..b844f8ce8 100644 --- a/matlab/ep/solve_stochastic_perfect_foresight_model_1.m +++ b/matlab/ep/solve_stochastic_perfect_foresight_model_1.m @@ -1,4 +1,4 @@ -function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo_simul,exo_simul,EpOptions,pfm,order) +function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo_simul,exo_simul,EpOptions,pfm,order,varargin) % Copyright (C) 2012-2013 Dynare Team % @@ -17,6 +17,11 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . + if nargin < 6 + homotopy_parameter = 1; + else + homotopy_parameter = varargin{1}; + end flag = 0; err = 0; stop = 0; @@ -159,7 +164,7 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo Y(i_cols_s,1); Y(i_cols_f,k1)]; end - [d1,jacobian] = dynamic_model(y,innovation, ... + [d1,jacobian] = dynamic_model(y,homotopy_parameter*innovation, ... params,steady_state,i+1); if i == 1 % in first period we don't keep track of @@ -182,7 +187,7 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo y = [Y(i_cols_p,1); Y(i_cols_s,j); Y(i_cols_f,j)]; - [d1,jacobian] = dynamic_model(y,innovation,params,steady_state,i+1); + [d1,jacobian] = dynamic_model(y,homotopy_parameter*innovation,params,steady_state,i+1); i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af]; A1(i_rows,i_cols_A) = jacobian(:,i_cols_j); res(:,i,j) = d1; @@ -194,7 +199,7 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo y = [Y(i_cols_p,j); Y(i_cols_s,j); Y(i_cols_f,j)]; - [d1,jacobian] = dynamic_model(y,innovation,params,steady_state,i+1); + [d1,jacobian] = dynamic_model(y,homotopy_parameter*innovation,params,steady_state,i+1); i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af]; A1(i_rows,i_cols_A) = jacobian(:,i_cols_j); res(:,i,j) = d1; @@ -242,17 +247,16 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo end if err < tolerance stop = 1; + flag = 0;% Convergency obtained. + endo_simul = reshape(Y(:,1),ny,periods+2);%Y(ny+(1:ny),1); if verbose + save ep_test_s1 exo_simul endo_simul Y fprintf('\n') ; disp([' Total time of simulation :' num2str(etime(clock,h1))]) ; fprintf('\n') ; disp([' Convergency obtained.']) ; fprintf('\n') ; end - flag = 0;% Convergency obtained. - endo_simul = reshape(Y(:,1),ny,periods+2);%Y(ny+(1:ny),1); - % figure;plot(Y(16:ny:(periods+2)*ny,:)) - % pause break end A2 = [nzA{:}]'; @@ -260,6 +264,8 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo if verbose disp(['solve_stochastic_foresight_model_1 encountered ' ... 'NaN']) + save ep_test_s2 exo_simul endo_simul + pause end flag = 1; return @@ -280,6 +286,8 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo disp(['WARNING : maximum number of iterations is reached (modify options_.simul.maxit).']) ; fprintf('\n') ; disp(sprintf('err: %f',err)); + save ep_test_s2 exo_simul endo_simul + pause end flag = 1;% more iterations are needed. endo_simul = 1; From 32f7f211d7bc5510ef4d6b2d7c26934f6ea65405 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Tue, 29 Apr 2014 17:56:45 +0200 Subject: [PATCH 04/37] doc: dseries and reporting presentation (draft) --- configure.ac | 1 + doc/Makefile.am | 2 +- doc/dseries-and-reporting/Makefile.am | 16 ++ .../dseriesReporting.tex | 260 ++++++++++++++++++ 4 files changed, 278 insertions(+), 1 deletion(-) create mode 100644 doc/dseries-and-reporting/Makefile.am create mode 100644 doc/dseries-and-reporting/dseriesReporting.tex diff --git a/configure.ac b/configure.ac index 60b946a3c..379b92433 100755 --- a/configure.ac +++ b/configure.ac @@ -181,6 +181,7 @@ AC_CONFIG_FILES([Makefile doc/parallel/Makefile doc/internals/Makefile doc/gsa/Makefile + doc/dseries-and-reporting/Makefile tests/Makefile matlab/dynare_version.m windows/dynare-version.nsi diff --git a/doc/Makefile.am b/doc/Makefile.am index 23b3ae432..6eb8fdcc2 100644 --- a/doc/Makefile.am +++ b/doc/Makefile.am @@ -1,4 +1,4 @@ -SUBDIRS = preprocessor macroprocessor userguide parallel internals gsa +SUBDIRS = preprocessor macroprocessor userguide parallel internals gsa dseries-and-reporting info_TEXINFOS = dynare.texi diff --git a/doc/dseries-and-reporting/Makefile.am b/doc/dseries-and-reporting/Makefile.am new file mode 100644 index 000000000..5d1ff915c --- /dev/null +++ b/doc/dseries-and-reporting/Makefile.am @@ -0,0 +1,16 @@ +if HAVE_PDFLATEX +if HAVE_BEAMER +pdf-local: dseriesReporting.pdf +endif +endif + +SRC = dseriesReporting.tex + +EXTRA_DIST = $(SRC) + +dseriesReporting.pdf: $(SRC) + $(PDFLATEX) dseriesReporting + $(PDFLATEX) dseriesReporting + +clean-local: + rm -f dseriesReporting.pdf *.toc *.aux *.log *.nav *.snm *.vrb *.out *~ diff --git a/doc/dseries-and-reporting/dseriesReporting.tex b/doc/dseries-and-reporting/dseriesReporting.tex new file mode 100644 index 000000000..8871ed949 --- /dev/null +++ b/doc/dseries-and-reporting/dseriesReporting.tex @@ -0,0 +1,260 @@ +\documentclass[10pt]{beamer} +\usepackage[utf8]{inputenc} +\usepackage{color} +\usepackage{amsmath} +\usepackage{epsf} +\usepackage{graphicx} +\usepackage{wasysym} +\usepackage{tikz} +\usetikzlibrary{positioning,shapes,shadows,arrows} +\definecolor{links}{HTML}{0000CC} +\hypersetup{colorlinks,linkcolor=,urlcolor=links} + + +\mode +{ + \usepackage{pgfpages} + \pgfpagesuselayout{4 on 1}[a4paper,border shrink=3mm,landscape] + \usetheme{CambridgeUS} + \usecolortheme{lily} +} + +\mode +{ + \usetheme{CambridgeUS} +} + +\title{Dynare Time Series \& Reporting} +\author{Houtan Bastani} +\institute{CEPREMAP} +\date{13 June 2014} + +\AtBeginSection[] +{ + \begin{frame} + \frametitle{Outline} + \tableofcontents[currentsection] + \end{frame} +} + +\setbeamerfont{frametitle}{family=\rmfamily,series=\bfseries,size={\fontsize{10}{10}}} +\setbeamertemplate{frametitle continuation}[from second] + +\tikzstyle{abstract}=[rectangle, rounded corners, draw=black, anchor=north, fill=blue!10, text centered, minimum height={height("Gp")+2pt}, minimum width=3cm, font=\footnotesize] + + +\begin{document} + +\begin{frame} + \titlepage +\end{frame} + +\begin{frame}[t] + \frametitle{Outline} + \tableofcontents +\end{frame} + + + + +% +% DSERIES +% +\section{Time Series} + +\subsection{Overview} +\begin{frame}[fragile,t] + \frametitle{Overview} + \begin{itemize} + \item Provide support for time series in Dynare + \item Introduced in Dynare 4.4 + \item Currently only used for reporting. + \item Use will increase with time (\textit{e.g.,} to be included in new estimation code) + \end{itemize} +\end{frame} + + +\subsection{Syntax} +\begin{frame}[fragile,t] + \frametitle{Syntax} + \begin{itemize} + \item Two components of a time series: + \begin{itemize} + \item A time component. In Dynare: \texttt{dates} + \item A data component mapped to time. In Dynare: \texttt{dseries} + \end{itemize} + \end{itemize} +\end{frame} + +\subsubsection{\texttt{dates} Syntax} +\begin{frame}[fragile,t] + \frametitle{\texttt{dates} Syntax} + \begin{itemize} + \item The \texttt{dates} command creates an object that represents at least one date at a given frequency + \begin{itemize} + \item Yearly: \texttt{`Y', `y', 1} + \item Quarterly: \texttt{`Q', `q', 4} + \item Monthly: \texttt{`M', `m', 12} + \item Weekly: \texttt{`W', `w', 52} + \end{itemize} + \item It has two slightly different syntaxes + \begin{itemize} + \item One for inclusion in \texttt{.m} files + \item One for inclusion in \texttt{.mod} files (simplified, taking advantage of the preprocessor) + \end{itemize} + \item Minimal restrictions on dates. Can be + \begin{itemize} + \item Negative + \item Empty + \item Noncontiguous + \end{itemize} + \end{itemize} +\end{frame} + + +\begin{frame}[fragile,t] + \frametitle{Creating a new \texttt{dates} object} + \begin{itemize} + \item A single date: + \begin{itemize} + \item In a \texttt{.m} file: \texttt{t = dates(`1999y');} + \item In a \texttt{.mod} file: \texttt{t = 1999y;} + \end{itemize} + \item A date range: + \begin{itemize} + \item In a \texttt{.m} file: \texttt{t = dates(`1999y'):dates(`2020y');} + \item In a \texttt{.mod} file: \texttt{t = 1999y:2020y;} + \end{itemize} + \end{itemize} +\end{frame} + + +\begin{frame}[fragile,t] + \frametitle{Modifying \texttt{dates}} + \begin{itemize} + \item \texttt{append}: appends a date to the date + \begin{itemize} + \item \texttt{t.append(dates(`2021y'));} + \end{itemize} + \item \texttt{pop}: + \item \texttt{sort}: + \end{itemize} +\end{frame} + + +\begin{frame}[fragile,t] + \frametitle{Getting info about \texttt{dates}} + \begin{itemize} + \item \texttt{double}: returns a floating point representation of the date + \begin{itemize} + \item \texttt{t.double;} + \end{itemize} + \item \texttt{freq}: returns the frequency + \begin{itemize} + \item \texttt{t.freq;} + \end{itemize} + \end{itemize} +\end{frame} + + +\begin{frame}[fragile,t] + \frametitle{Comparing \texttt{dates}} + \begin{itemize} + \item + \end{itemize} +\end{frame} + +\subsubsection{\texttt{dseries} Syntax} + + +\subsection{Examples} + + +% +% REPORTING +% +\section{Reporting} +\subsection{Overview} +\begin{frame} + \frametitle{Overview} + \begin{itemize} + \item Introduced in Dynare 4.4 + \item Introduce reporting functionality to Dynare + \begin{itemize} + \item Input: \texttt{dseries} + \item Output: \LaTeX\ report \& compiled \texttt{.pdf} + \end{itemize} + \item Graphs and Tables are modular + \begin{itemize} + \item Can easily be included in another document + \end{itemize} + \item Graphs are produced in Ti$k$Z + \begin{itemize} + \item Scales well + \item Formating follows that of enclosing document + \end{itemize} + \item Works with Matlab \& Octave + \item Works approximately 5 times faster than Iris reporting + \end{itemize} +\end{frame} + +\begin{frame} + \frametitle{Reporting Class Hierarchy} + \centering { + \begin{tikzpicture}[ + node distance = .45cm, + auto, + line/.style={->, >=stealth'}, + ] + \node (Report) [abstract, rectangle split, rectangle split parts=2] + { + \textbf{Report} + \nodepart{second}\texttt{report(...);} + }; + \node (Page) [abstract, rectangle split, rectangle split parts=2, below=of Report] + { + \textbf{Page} + \nodepart{second}\texttt{addPage(...);} + }; + \node (Section) [abstract, rectangle split, rectangle split parts=2, below=of Page] + { + \textbf{Section} + \nodepart{second}\texttt{addSection(...);} + }; + \node (Vspace) [abstract, rectangle split, rectangle split parts=2, below=of Section] + { + \textbf{Vspace} + \nodepart{second}\texttt{addVspace(...);} + }; + \node (Graph) [abstract, rectangle split, rectangle split parts=2, left=of Vspace] + { + \textbf{Graph} + \nodepart{second}\texttt{addGraph(...);} + }; + \node (Table) [abstract, rectangle split, rectangle split parts=2, right=of Vspace, text height=] + { + \textbf{Table} + \nodepart{second}\texttt{addTable(...);} + }; + \node (Series) [abstract, rectangle split, rectangle split parts=2, below=of Vspace] + { + \textbf{Series} + \nodepart{second}\texttt{addSeries(...);} + }; + \draw [line] (Series) to node { } (Table); + \draw [line] (Series) to node { } (Graph); + \draw [line] (Table) to node { } (Section); + \draw [line] (Graph) to node { } (Section); + \draw [line] (Vspace) to node { } (Section); + \draw [line] (Section) to node { } (Page); + \draw [line] (Page) to node { } (Report); + \end{tikzpicture} + } +\end{frame} + +\subsection{Syntax} + + +\subsection{Examples} + +\end{document} From 0cd6c9917e3e377fe90bf3d49d0378f60e455ff9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Scylla=29?= Date: Thu, 27 Mar 2014 10:27:49 +0100 Subject: [PATCH 05/37] Changed name of the clean routine (otherwise if a user types clean instead of clear all the generated files and data would be erased). --- matlab/utilities/general/{clean.m => clean_current_folder.m} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename matlab/utilities/general/{clean.m => clean_current_folder.m} (100%) diff --git a/matlab/utilities/general/clean.m b/matlab/utilities/general/clean_current_folder.m similarity index 100% rename from matlab/utilities/general/clean.m rename to matlab/utilities/general/clean_current_folder.m From f183c047d5d084539092d9ad80ae0dd1c5a03eb6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Scylla=29?= Date: Tue, 22 Apr 2014 15:19:16 +0200 Subject: [PATCH 06/37] Test the existence of the files and folder before deletion. --- matlab/utilities/general/clean_current_folder.m | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/matlab/utilities/general/clean_current_folder.m b/matlab/utilities/general/clean_current_folder.m index 5f1967b30..e3408bec1 100644 --- a/matlab/utilities/general/clean_current_folder.m +++ b/matlab/utilities/general/clean_current_folder.m @@ -19,11 +19,18 @@ function clean() a = dir('*.mod'); + for i = 1:length(a) [junk,basename,extension] = fileparts(a(i).name); - delete([basename '.m']); - delete([basename '.log']); - rmdir(basename,'s'); + if exist([basename '.m']) + delete([basename '.m']); + end + if exist([basename '.log']) + delete([basename '.log']); + end + if exist(basename,'dir') + rmdir(basename,'s'); + end if exist([basename '_steadystate.m']) movefile([basename '_steadystate.m'],['protect_' basename '_steadystate.m']); end From f03dd893a82db4185fdd93d8df55d975b8941514 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Scylla=29?= Date: Wed, 30 Apr 2014 11:11:45 +0200 Subject: [PATCH 07/37] Allow the syntax o.disp() for @dates and @dseries objects. --- matlab/@dates/subsref.m | 5 ++++- matlab/@dseries/subsref.m | 3 +++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/matlab/@dates/subsref.m b/matlab/@dates/subsref.m index 9073f2f46..293ce41c6 100644 --- a/matlab/@dates/subsref.m +++ b/matlab/@dates/subsref.m @@ -37,7 +37,7 @@ switch S(1).type error(['dates::subsref: ' S(1).subs ' is not a method but a member!']) end B = builtin('subsref', A, S(1)); - case {'sort','unique','double','isempty','length'}% Public methods (without arguments) + case {'sort','unique','double','isempty','length'}% Public methods (without input arguments) B = feval(S(1).subs,A); if length(S)>1 && isequal(S(2).type,'()') && isempty(S(2).subs) S = shiftS(S,1); @@ -49,6 +49,9 @@ switch S(1).type else error('dates::subsref: Something is wrong in your syntax!') end + case {'disp'} + feval(S(1).subs,A); + return otherwise error('dates::subsref: Unknown public member or method!') end diff --git a/matlab/@dseries/subsref.m b/matlab/@dseries/subsref.m index 765dc6996..ab3b08dc5 100644 --- a/matlab/@dseries/subsref.m +++ b/matlab/@dseries/subsref.m @@ -157,6 +157,9 @@ switch S(1).type case {'set_names','rename','tex_rename'} B = feval(S(1).subs,A,S(2).subs{:}); S = shiftS(S,1); + case {'disp'} + feval(S(1).subs,A); + return otherwise % Extract a sub-object by selecting one variable. ndx = find(strcmp(S(1).subs,A.name)); if ~isempty(ndx) From 6ca649c4bfdb9c548b64ee113078337545ba984a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Scylla=29?= Date: Fri, 2 May 2014 12:15:58 +0200 Subject: [PATCH 08/37] Added char method for @dates objects. --- matlab/@dates/char.m | 35 +++++++++++++++++++++++++++++++++++ matlab/@dates/subsref.m | 2 +- 2 files changed, 36 insertions(+), 1 deletion(-) create mode 100644 matlab/@dates/char.m diff --git a/matlab/@dates/char.m b/matlab/@dates/char.m new file mode 100644 index 000000000..bdde032dc --- /dev/null +++ b/matlab/@dates/char.m @@ -0,0 +1,35 @@ +function s = char(dd) + +% Given a one element dates object, returns a string with the formatted date. +% +% INPUTS +% o dd dates object with one element +% +% OUTPUTS +% o s a string +% +% SPECIAL REQUIREMENTS +% none + +% Copyright (C) 2014 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +if length(dd)>1 + error('The input argument must be a singleton dates object!') +end + +s = date2string(dd.time, dd.freq); \ No newline at end of file diff --git a/matlab/@dates/subsref.m b/matlab/@dates/subsref.m index 293ce41c6..0cd0cd401 100644 --- a/matlab/@dates/subsref.m +++ b/matlab/@dates/subsref.m @@ -37,7 +37,7 @@ switch S(1).type error(['dates::subsref: ' S(1).subs ' is not a method but a member!']) end B = builtin('subsref', A, S(1)); - case {'sort','unique','double','isempty','length'}% Public methods (without input arguments) + case {'sort','unique','double','isempty','length','char'}% Public methods (without input arguments) B = feval(S(1).subs,A); if length(S)>1 && isequal(S(2).type,'()') && isempty(S(2).subs) S = shiftS(S,1); From 9a6ed32ff3e4c2212e4ebfaafd699708c5e9032a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Tue, 6 May 2014 16:05:47 +0200 Subject: [PATCH 09/37] Use correct texinfo commands for *_plan functions. --- doc/dynare.texi | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/doc/dynare.texi b/doc/dynare.texi index c4f1deb31..06b567827 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -6097,37 +6097,37 @@ can be computed using an extended path method. The forecast scenario describing The first step is the forecast scenario initialization using the function @code{init_plan}: @anchor{init_plan} -@deffn {MATLAB/Octave command} HANDLE = init_plan(DATES) ; +@deftypefn {MATLAB/Octave command} {HANDLE =} init_plan (DATES) ; Creates a new forecast scenario for a forecast period (indicated as a dates class, see @ref{dates class members}). This function return a handle on the new forecast scenario. -@end deffn +@end deftypefn The forecast scenario can contain some simple shocks on the exogenous variables. This shocks are described using the function @code{basic_plan}: @anchor{basic_plan} -@deffn {MATLAB/Octave command} HANDLE = basic_plan(HANDLE, 'VARIABLE_NAME', 'SHOCK_TYPE', DATES, MATLAB VECTOR OF DOUBLE | [DOUBLE | EXPRESSION [DOUBLE | | EXPRESSION] ] ) ; +@deftypefn {MATLAB/Octave command} {HANDLE =} basic_plan (HANDLE, 'VARIABLE_NAME', 'SHOCK_TYPE', DATES, MATLAB VECTOR OF DOUBLE | [DOUBLE | EXPRESSION [DOUBLE | | EXPRESSION] ] ) ; Adds to the forecast scenario a shock on the exogenous variable indicated between quotes in the second argument. The shock type has to be specified in the third argument between quotes: 'surprise' in case of an unexpected shock or 'perfect_foresight' for a perfectly anticipated shock. The fourth argument indicates the period of the shock using a dates class (see @ref{dates class members}). The last argument is the shock path indicated as a Matlab vector of double. This function return the handle of the updated forecast scenario. -@end deffn +@end deftypefn The forecast scenario can also contain a constrained path on an endogenous variable. The values of the related exogenous variable compatible with the constrained path are in this case computed. In other words, a conditional forecast is performed. This kind of shock is described with the function @code{flip_plan}: @anchor{flip_plan} -@deffn {MATLAB/Octave command} HANDLE = flip_plan(HANDLE, 'VARIABLE_NAME, 'VARIABLE_NAME', 'SHOCK_TYPE', DATES, MATLAB VECTOR OF DOUBLE | [DOUBLE | EXPRESSION [DOUBLE | | EXPRESSION] ] ) ; +@deftypefn {MATLAB/Octave command} {HANDLE =} flip_plan (HANDLE, 'VARIABLE_NAME, 'VARIABLE_NAME', 'SHOCK_TYPE', DATES, MATLAB VECTOR OF DOUBLE | [DOUBLE | EXPRESSION [DOUBLE | | EXPRESSION] ] ) ; Adds to the forecast scenario a constrained path on the endogenous variable specified between quotes in the second argument. The associated exogenous variable provided in the third argument between quotes, is considered as an endogenous variable and its values compatible with the constrained path on the endogenous variable will be computed. The nature of the expectation on the constrained path has to be specified in the fourth argument between quotes: 'surprise' in case of an unexpected path or 'perfect_foresight' for a perfectly anticipated path. The fifth argument indicates the period where the path of the endogenous variable is constrained using a dates class (see @ref{dates class members}). The last argument contains the constrained path as a Matlab vector of double. This function return the handle of the updated forecast scenario. -@end deffn +@end deftypefn Once the forecast scenario if fully described, the forecast is computed with the command @code{det_cond_forecast}: @anchor{det_cond_forecast} -@deffn {MATLAB/Octave command} DSERIES = det_cond_forecast(HANDLE[, DSERIES [, DATES]]) ; +@deftypefn {MATLAB/Octave command} {DSERIES =} det_cond_forecast (HANDLE[, DSERIES [, DATES]]) ; Computes the forecast or the conditional forecast using an extended path method for the given forecast scenario (first argument). The past values of the endogenous and exogenous variables provided with a dseries class (see @ref{dseries class members}) can be indicated in the second argument. By default, the past values of the variables are equal to their steady-state values. The initial date of the forecast can be provided in the third argument. By default, the forecast will start at the first date indicated in the @code{init_plan} command. This function returns a dset containing the historical and forecast values for the endogenous and exogenous variables. -@end deffn +@end deftypefn From b3f2822c61afeda00cd636b27bdd2b23170a44b6 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Tue, 6 May 2014 18:05:14 +0200 Subject: [PATCH 10/37] doc: complete dates portion of dseries/reporting slides --- .../dseriesReporting.tex | 135 +++++++++++++----- 1 file changed, 102 insertions(+), 33 deletions(-) diff --git a/doc/dseries-and-reporting/dseriesReporting.tex b/doc/dseries-and-reporting/dseriesReporting.tex index 8871ed949..5b2643749 100644 --- a/doc/dseries-and-reporting/dseriesReporting.tex +++ b/doc/dseries-and-reporting/dseriesReporting.tex @@ -1,16 +1,12 @@ \documentclass[10pt]{beamer} -\usepackage[utf8]{inputenc} -\usepackage{color} -\usepackage{amsmath} -\usepackage{epsf} -\usepackage{graphicx} -\usepackage{wasysym} + \usepackage{tikz} \usetikzlibrary{positioning,shapes,shadows,arrows} +\tikzstyle{abstract}=[rectangle, rounded corners, draw=black, anchor=north, fill=blue!10, text centered, minimum height={height("Gp")+2pt}, minimum width=3cm, font=\footnotesize] + \definecolor{links}{HTML}{0000CC} \hypersetup{colorlinks,linkcolor=,urlcolor=links} - \mode { \usepackage{pgfpages} @@ -24,11 +20,6 @@ \usetheme{CambridgeUS} } -\title{Dynare Time Series \& Reporting} -\author{Houtan Bastani} -\institute{CEPREMAP} -\date{13 June 2014} - \AtBeginSection[] { \begin{frame} @@ -40,8 +31,10 @@ \setbeamerfont{frametitle}{family=\rmfamily,series=\bfseries,size={\fontsize{10}{10}}} \setbeamertemplate{frametitle continuation}[from second] -\tikzstyle{abstract}=[rectangle, rounded corners, draw=black, anchor=north, fill=blue!10, text centered, minimum height={height("Gp")+2pt}, minimum width=3cm, font=\footnotesize] - +\title{Dynare Time Series \& Reporting} +\author{Houtan Bastani} +\institute{CEPREMAP} +\date{13 June 2014} \begin{document} @@ -54,11 +47,8 @@ \tableofcontents \end{frame} - - - % -% DSERIES +% DATES % \section{Time Series} @@ -100,7 +90,10 @@ \item It has two slightly different syntaxes \begin{itemize} \item One for inclusion in \texttt{.m} files - \item One for inclusion in \texttt{.mod} files (simplified, taking advantage of the preprocessor) + \item One for inclusion in \texttt{.mod} files (simplified using the preprocessor) + \begin{itemize} + \item To prevent date translation, escape the date with `\texttt{\$}' (\textit{e.g.,} \texttt{\$2020y}) + \end{itemize} \end{itemize} \item Minimal restrictions on dates. Can be \begin{itemize} @@ -122,8 +115,8 @@ \end{itemize} \item A date range: \begin{itemize} - \item In a \texttt{.m} file: \texttt{t = dates(`1999y'):dates(`2020y');} - \item In a \texttt{.mod} file: \texttt{t = 1999y:2020y;} + \item In a \texttt{.m} file: \texttt{dr = dates(`1999y'):dates(`2020y');} + \item In a \texttt{.mod} file: \texttt{dr = 1999y:2020y;} \end{itemize} \end{itemize} \end{frame} @@ -134,10 +127,42 @@ \begin{itemize} \item \texttt{append}: appends a date to the date \begin{itemize} - \item \texttt{t.append(dates(`2021y'));} + \item \texttt{t=t.append(dates(`1900y')); \% } + \end{itemize} + \item \texttt{horzcat}: horizontal concatenation + \begin{itemize} + \item \texttt{[t t]; \% }; + \end{itemize} + \item \texttt{minus}: either the distance between two \texttt{dates} or lag one \texttt{dates} + \begin{itemize} + \item \texttt{t-t \% [0 0]'} + \item \texttt{t-[3 3]' \% } + \end{itemize} + \item \texttt{plus}: either combine two \texttt{dates} or forward one \texttt{dates} + \begin{itemize} + \item \texttt{t+t \% } + \item \texttt{t+[3 3]' \% } + \end{itemize} + \item \texttt{pop}: remove last element + \begin{itemize} + \item \texttt{t.pop(); \% } + \end{itemize} + \item \texttt{sort}: sort dates in ascending order + \begin{itemize} + \item \texttt{t=t.sort(); \% } + \end{itemize} + \item \texttt{uminus}: shifts dates back one period + \begin{itemize} + \item \texttt{-t; \% } + \end{itemize} + \item \texttt{unique}: removes repetitions + \begin{itemize} + \item \texttt{t.append(dates(`1999y')).unique() \% } + \end{itemize} + \item \texttt{uplus}: shifts dates forward one period + \begin{itemize} + \item \texttt{++t; \% } \end{itemize} - \item \texttt{pop}: - \item \texttt{sort}: \end{itemize} \end{frame} @@ -145,25 +170,69 @@ \begin{frame}[fragile,t] \frametitle{Getting info about \texttt{dates}} \begin{itemize} + \item \texttt{char}: returns a single date as a string + \begin{itemize} + \item \texttt{t(1).char() \% 1999Y} + \end{itemize} \item \texttt{double}: returns a floating point representation of the date \begin{itemize} - \item \texttt{t.double;} + \item \texttt{t.double() \% [1999 1900]'} \end{itemize} \item \texttt{freq}: returns the frequency \begin{itemize} - \item \texttt{t.freq;} + \item \texttt{t.freq; \% 1} + \end{itemize} + \item \texttt{isequal}: returns true if the two arguments are equal + \begin{itemize} + \item \texttt{isequal(t,t) \% 1} + \end{itemize} + \item \texttt{length}: returns the number of dates + \begin{itemize} + \item \texttt{t.length() \% 2} + \end{itemize} + \item \texttt{max}: returns the maximum \texttt{dates} in the arguments + \begin{itemize} + \item \texttt{max(t,dr) \% } + \end{itemize} + \item \texttt{min}: returns the minimum \texttt{dates} in the arguments + \begin{itemize} + \item \texttt{min(t,dr) \% } + \end{itemize} + \item \texttt{eq, ge, gt, le, lt, ne}: returns boolean value of comparison + \begin{itemize} + \item \texttt{t==t \% [1 1]'} + \item \texttt{t>=dates(`1950y') \% [1 0]'} + \item \texttt{t$\thicksim$=dates(`1999y') \% [0 1]'} + \end{itemize} + \end{itemize} +\end{frame} + +\begin{frame}[fragile,t] + \frametitle{Set Operations on \texttt{dates}} + \begin{itemize} + \item \texttt{intersect}: returns the intersection of the arguments + \begin{itemize} + \item \texttt{intersect(t,dr) \% } + \end{itemize} + \item \texttt{isempty}: returns true if the argument is empty + \begin{itemize} + \item \texttt{isempty(t) \% 0} + \end{itemize} + \item \texttt{setdiff}: returns dates present in first arg but not in second + \begin{itemize} + \item \texttt{setdiff(t,dr) \% } + \end{itemize} + \item \texttt{union}: + \begin{itemize} + \item \texttt{union(dr,t) \% } \end{itemize} \end{itemize} \end{frame} -\begin{frame}[fragile,t] - \frametitle{Comparing \texttt{dates}} - \begin{itemize} - \item - \end{itemize} -\end{frame} - +% +% DSERIES +% \subsubsection{\texttt{dseries} Syntax} From c4a4e6eeb3a443d19a2045e31e82796b6f9ba1b1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Wed, 7 May 2014 10:42:28 +0200 Subject: [PATCH 11/37] Octave testsuite: force gnuplot toolkit for text graphics. This is necessary since FLTK is now the default toolkit in octave 3.8. --- tests/run_test_octave.m | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/tests/run_test_octave.m b/tests/run_test_octave.m index fb67fd5f7..b78e25ee5 100644 --- a/tests/run_test_octave.m +++ b/tests/run_test_octave.m @@ -1,4 +1,4 @@ -## Copyright (C) 2009-2012 Dynare Team +## Copyright (C) 2009-2014 Dynare Team ## ## This file is part of Dynare. ## @@ -31,9 +31,8 @@ if !strcmp(dynare_version(), getenv("DYNARE_VERSION")) endif ## Ask gnuplot to create graphics in text mode -## Note that setenv() was introduced in Octave 3.0.2, for compatibility -## with MATLAB -putenv("GNUTERM", "dumb") +graphics_toolkit gnuplot; +setenv("GNUTERM", "dumb"); ## Test MOD files listed in Makefile.am name = getenv("FILESTEM"); From d44825723b7ea820c3ba2c12c9c016858b6cf3b6 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Wed, 7 May 2014 10:49:50 +0200 Subject: [PATCH 12/37] Adds test whether all parameters are set to model_diagnostics.m --- matlab/model_diagnostics.m | 6 ++++++ matlab/test_for_deep_parameters_calibration.m | 9 ++++++--- 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/matlab/model_diagnostics.m b/matlab/model_diagnostics.m index 92122514f..29b56b318 100644 --- a/matlab/model_diagnostics.m +++ b/matlab/model_diagnostics.m @@ -63,6 +63,12 @@ if M.exo_nbr == 0 oo.exo_steady_state = [] ; end + +info=test_for_deep_parameters_calibration(M); +if info + problem_dummy=1; +end; + % check if ys is steady state options.debug=1; %locally set debug option to 1 [dr.ys,params,check1]=evaluate_steady_state(oo.steady_state,M,options,oo,1); diff --git a/matlab/test_for_deep_parameters_calibration.m b/matlab/test_for_deep_parameters_calibration.m index 0f5ab7aec..c5b88b40f 100644 --- a/matlab/test_for_deep_parameters_calibration.m +++ b/matlab/test_for_deep_parameters_calibration.m @@ -1,11 +1,11 @@ -function test_for_deep_parameters_calibration(M_) +function info=test_for_deep_parameters_calibration(M_) % Issues a warning is some of the parameters are NaNs. % % INPUTS % M_ [structure] Description of the (simulated or estimated) model. % % OUTPUTS -% none +% info [scalar] 0 if no problems detected, 1 otherwise % % ALGORITHM % none @@ -13,7 +13,7 @@ function test_for_deep_parameters_calibration(M_) % SPECIAL REQUIREMENTS % none -% Copyright (C) 2010 Dynare Team +% Copyright (C) 2010-2014 Dynare Team % % This file is part of Dynare. % @@ -31,6 +31,7 @@ function test_for_deep_parameters_calibration(M_) % along with Dynare. If not, see . plist = list_of_parameters_calibrated_as_NaN(M_); if ~isempty(plist) + info=1; message = ['Some of the parameters have no value (' ]; for i=1:size(plist,1) if i Date: Wed, 7 May 2014 11:14:54 +0200 Subject: [PATCH 13/37] Ref. manual: fix default number of periods of forecast command. --- doc/dynare.texi | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/dynare.texi b/doc/dynare.texi index 06b567827..5f18125f3 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -5820,7 +5820,7 @@ confidence interval. @table @code @item periods = @var{INTEGER} -Number of periods of the forecast. Default: @code{40} +Number of periods of the forecast. Default: @code{5}. @item conf_sig = @var{DOUBLE} @anchor{conf_sig} Level of significance for confidence From f31d060a4f5e40b4bfd586b1e835a28bbd81a023 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Wed, 7 May 2014 11:16:47 +0200 Subject: [PATCH 14/37] Add all-in-one-file example for smoother2histval. --- tests/Makefile.am | 3 +- .../fs2000_smooth_stoch_simul.mod | 82 +++++++++++++++++++ 2 files changed, 84 insertions(+), 1 deletion(-) create mode 100644 tests/smoother2histval/fs2000_smooth_stoch_simul.mod diff --git a/tests/Makefile.am b/tests/Makefile.am index b53d1f5a6..2600f1295 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -176,7 +176,8 @@ MODFILES = \ filter_step_ahead/fs2000_filter_step_ahead_ML.mod \ loglinear/example4_loglinear.mod \ smoother2histval/fs2000_simul.mod \ - smoother2histval/fs2000_smooth.mod + smoother2histval/fs2000_smooth.mod \ + smoother2histval/fs2000_smooth_stoch_simul.mod XFAIL_MODFILES = ramst_xfail.mod \ estim_param_in_shock_value.mod diff --git a/tests/smoother2histval/fs2000_smooth_stoch_simul.mod b/tests/smoother2histval/fs2000_smooth_stoch_simul.mod new file mode 100644 index 000000000..f599f37b7 --- /dev/null +++ b/tests/smoother2histval/fs2000_smooth_stoch_simul.mod @@ -0,0 +1,82 @@ +// Test that smoother2histval works (everything in a single file) +// Note that an observation equation has been modified in order to have an aux var for lagged endo + +var m P c e W R k d n l gy_obs gp_obs y dA; +varexo e_a e_m; + +parameters alp bet gam mst rho psi del; + +alp = 0.33; +bet = 0.99; +gam = 0.003; +mst = 1.011; +rho = 0.7; +psi = 0.787; +del = 0.02; + +model; +dA = exp(gam+e_a); +log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m; +-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0; +W = l/n; +-(psi/(1-psi))*(c*P/(1-n))+l/n = 0; +R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W; +1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0; +c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1); +P*c = m; +m-1+d = l; +e = exp(e_a); +y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a(-1))); +gy_obs = dA*y/y(-2); +gp_obs = (P/P(-1))*m(-1)/dA; +end; + +initval; +k = 6; +m = mst; +P = 2.25; +c = 0.45; +e = 1; +W = 4; +R = 1.02; +d = 0.85; +n = 0.19; +l = 0.86; +y = 0.6; +gy_obs = exp(gam); +gp_obs = exp(-gam); +dA = exp(gam); +end; + +shocks; +var e_a; stderr 0.014; +var e_m; stderr 0.005; +end; + +steady; + +check; + +estimated_params; +alp, beta_pdf, 0.356, 0.02; +bet, beta_pdf, 0.993, 0.002; +gam, normal_pdf, 0.0085, 0.003; +mst, normal_pdf, 1.0002, 0.007; +rho, beta_pdf, 0.129, 0.223; +psi, beta_pdf, 0.65, 0.05; +del, beta_pdf, 0.01, 0.005; +stderr e_a, inv_gamma_pdf, 0.035449, inf; +stderr e_m, inv_gamma_pdf, 0.008862, inf; +end; + +varobs gp_obs gy_obs; + +options_.solve_tolf = 1e-12; + +estimation(order=1,datafile=fsdat_simul,nobs=192,mh_replic=1500,mh_nblocks=1,mh_jscale=0.8,smoother,consider_all_endogenous); + +smoother2histval(period = 5); + +stoch_simul(nomoments); + +forecast; From 58ba964bb7e286e2fbf679eb39168e1c4e281b4b Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Wed, 7 May 2014 13:34:19 +0200 Subject: [PATCH 15/37] trust_region/dogleg: fixing sign error --- matlab/trust_region.m | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/matlab/trust_region.m b/matlab/trust_region.m index 804ba7ce0..1edb87732 100644 --- a/matlab/trust_region.m +++ b/matlab/trust_region.m @@ -181,6 +181,7 @@ while (niter < maxiter && ~info) % FIXME -- why tolf*n*xn? If abs (e) ~ abs(x) * eps is a vector % of perturbations of x, then norm (fjac*e) <= eps*n*xn, i.e. by % tolf ~ eps we demand as much accuracy as we can expect. + disp([niter fn ratio]) if (fn <= tolf*n*xn) info = 1; % The following tests done only after successful step. @@ -228,7 +229,7 @@ if (xn > delta) bn = norm (b); dxn = delta/xn; snmd = snm/delta; t = (bn/sn) * (bn/xn) * snmd; - t = t - dxn * snmd^2 - sqrt ((t-dxn)^2 + (1-dxn^2)*(1-snmd^2)); + t = t - dxn * snmd^2 + sqrt ((t-dxn)^2 + (1-dxn^2)*(1-snmd^2)); alpha = dxn*(1-snmd^2) / t; else alpha = 0; From fad9aa684636347ac997bd6e5b328a9cf575b86d Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Wed, 7 May 2014 14:23:31 +0200 Subject: [PATCH 16/37] removing debugging printing --- matlab/trust_region.m | 1 - 1 file changed, 1 deletion(-) diff --git a/matlab/trust_region.m b/matlab/trust_region.m index 1edb87732..5fb93c5c6 100644 --- a/matlab/trust_region.m +++ b/matlab/trust_region.m @@ -181,7 +181,6 @@ while (niter < maxiter && ~info) % FIXME -- why tolf*n*xn? If abs (e) ~ abs(x) * eps is a vector % of perturbations of x, then norm (fjac*e) <= eps*n*xn, i.e. by % tolf ~ eps we demand as much accuracy as we can expect. - disp([niter fn ratio]) if (fn <= tolf*n*xn) info = 1; % The following tests done only after successful step. From a6106f07c794ce91172f12613c8e33ced7682cd2 Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Wed, 7 May 2014 21:12:56 +0200 Subject: [PATCH 17/37] making solve1.m robust to sparse Jacobian --- matlab/solve1.m | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/matlab/solve1.m b/matlab/solve1.m index 4c787ddc3..fd57b2618 100644 --- a/matlab/solve1.m +++ b/matlab/solve1.m @@ -90,9 +90,14 @@ for its = 1:maxit g = (fvec'*fjac)'; if debug - disp(['cond(fjac) ' num2str(cond(fjac))]) + disp(['cond(fjac) ' num2str(condest(fjac))]) end - if rcond(fjac) < sqrt(eps) + if issparse(fjac) + rcond_fjac = 1/condest(fjac); + else + rcond_fjac = rcond(fjac); + end + if rcond_fjac < sqrt(eps) fjac2=fjac'*fjac; p=-(fjac2+1e6*sqrt(nn*eps)*max(sum(abs(fjac2)))*eye(nn))\(fjac'*fvec); else From 57ae7d7fd983bea2bbf041c7331cf28e9080c119 Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Wed, 7 May 2014 21:25:10 +0200 Subject: [PATCH 18/37] adding solve_algo=9: using trust_region algorithm on the entire model --- doc/dynare.texi | 6 +++++- matlab/dynare_solve.m | 6 +++++- 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/doc/dynare.texi b/doc/dynare.texi index 5f18125f3..7ddd04dd0 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -2758,7 +2758,8 @@ using the same solver as value @code{1} Use Chris Sims' solver @item 4 -Trust-region solver with autoscaling. +Splits the model into recursive blocks and solves each block in turn +using a trust-region solver with autoscaling. @item 5 Newton algorithm with a sparse Gaussian elimination (SPE) (requires @@ -2777,6 +2778,9 @@ each iteration (requires @code{bytecode} and/or @code{block} option, Newton algorithm with a Stabilized Bi-Conjugate Gradient (BICGSTAB) solver at each iteration (requires @code{bytecode} and/or @code{block} option, @pxref{Model declaration}) + +@item 9 +Trust-region algorithm on the entire model. @end table @noindent diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m index 744d174e2..f98980d97 100644 --- a/matlab/dynare_solve.m +++ b/matlab/dynare_solve.m @@ -109,7 +109,11 @@ if options_.solve_algo == 0 info = 1; end elseif options_.solve_algo == 1 - [x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag,options_.gstep, ... + [x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag,options_.gstep, ... + tolf,options_.solve_tolx, ... + options_.steady.maxit,options_.debug,varargin{:}); +elseif options_.solve_algo == 5 + [x,info]=trust_region(func,x,1:nn,1:nn,jacobian_flag,options_.gstep, ... tolf,options_.solve_tolx, ... options_.steady.maxit,options_.debug,varargin{:}); elseif options_.solve_algo == 2 || options_.solve_algo == 4 From 8bf484a140f8a1d31f474c93e4201740444dcc5c Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Thu, 8 May 2014 12:23:09 +0200 Subject: [PATCH 19/37] fixed expectation model test cases. They had a crazy steady state revealed by the now more accurate nonlinear solver. --- tests/expectations/expectation.mod | 2 +- tests/expectations/expectation_nested.mod | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/expectations/expectation.mod b/tests/expectations/expectation.mod index fde5eada5..c05d7528d 100644 --- a/tests/expectations/expectation.mod +++ b/tests/expectations/expectation.mod @@ -15,7 +15,7 @@ theta = 2.95; phi = 0.1; model; -c*theta*h^(1+psi)=expectation(1)((1-alpha)*y)+expectation(-2)((1-alpha)*y); +c*theta*h^(1+psi)=(expectation(1)((1-alpha)*y)+expectation(-2)((1-alpha)*y))/2; k = beta*(((exp(b)*c)/(exp(b(+1))*c(+1))) *(exp(b(+1))*alpha*y(+1)+(1-delta)*k)); y = exp(a)*(k(-1)^alpha)*(h^(1-alpha)); diff --git a/tests/expectations/expectation_nested.mod b/tests/expectations/expectation_nested.mod index 0bc4069d8..d6bd10aaf 100644 --- a/tests/expectations/expectation_nested.mod +++ b/tests/expectations/expectation_nested.mod @@ -15,7 +15,7 @@ theta = 2.95; phi = 0.1; model; -c*theta*h^(1+psi)=expectation(1)(((1-alpha)*y)+expectation(3)((1-alpha)*y)); +c*theta*h^(1+psi)=expectation(1)(((1-alpha)*y)+expectation(3)((1-alpha)*y))/2; k = beta*(((exp(b)*c)/(exp(b(+1))*c(+1))) *(exp(b(+1))*alpha*y(+1)+(1-delta)*k)); y = exp(a)*(k(-1)^alpha)*(h^(1-alpha)); From e8ac5da9ea73b9af09414f133626d4f5b65c72ec Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Fri, 9 May 2014 15:33:26 +0200 Subject: [PATCH 20/37] Fix trust region on entire model. --- matlab/dynare_solve.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m index f98980d97..3033f1143 100644 --- a/matlab/dynare_solve.m +++ b/matlab/dynare_solve.m @@ -112,7 +112,7 @@ elseif options_.solve_algo == 1 [x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag,options_.gstep, ... tolf,options_.solve_tolx, ... options_.steady.maxit,options_.debug,varargin{:}); -elseif options_.solve_algo == 5 +elseif options_.solve_algo == 9 [x,info]=trust_region(func,x,1:nn,1:nn,jacobian_flag,options_.gstep, ... tolf,options_.solve_tolx, ... options_.steady.maxit,options_.debug,varargin{:}); @@ -165,5 +165,5 @@ elseif options_.solve_algo == 3 [x,info] = csolve(func,x,[],1e-6,500,varargin{:}); end else - error('DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4]') + error('DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4,9]') end From 72950de0a09bbd854ce7996c41f56cfada0d436e Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Sun, 11 May 2014 14:44:06 +0200 Subject: [PATCH 21/37] solve1: initialize fjac only if it needs to be computed numerically. --- matlab/solve1.m | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/matlab/solve1.m b/matlab/solve1.m index fd57b2618..17bb0a3ab 100644 --- a/matlab/solve1.m +++ b/matlab/solve1.m @@ -41,7 +41,6 @@ function [x,check] = solve1(func,x,j1,j2,jacobian_flag,gstep,tolf,tolx,maxit,deb % along with Dynare. If not, see . nn = length(j1); -fjac = zeros(nn,nn) ; g = zeros(nn,1) ; tolmin = tolx ; @@ -71,6 +70,9 @@ end stpmax = stpmx*max([sqrt(x'*x);nn]) ; first_time = 1; +if ~jacobian_flag + fjac = zeros(nn,nn) ; +end for its = 1:maxit if jacobian_flag [fvec,fjac] = feval(func,x,varargin{:}); From b5ca118bfb3a062daf0fdbd086e57c020e4e011a Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Sun, 11 May 2014 14:44:44 +0200 Subject: [PATCH 22/37] trust_region: replace a loop with a matrix expression --- matlab/trust_region.m | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/matlab/trust_region.m b/matlab/trust_region.m index 5fb93c5c6..ec83d5df9 100644 --- a/matlab/trust_region.m +++ b/matlab/trust_region.m @@ -65,7 +65,6 @@ info = 0; fvec = fcn (x, varargin{:}); fvec = fvec(j1); fn = norm (fvec); -jcn = nan(n, 1); % Outer loop. while (niter < maxiter && ~info) @@ -88,9 +87,7 @@ while (niter < maxiter && ~info) end % Get column norms, use them as scaling factors. - for j = 1:n - jcn(j) = norm(fjac(:,j)); - end + jcn = sqrt(sum(fjac.*fjac))'; if (niter == 1) dg = jcn; dg(dg == 0) = 1; From add594ab7eafe1601549f5653c2ceb0d053ad957 Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Mon, 12 May 2014 09:40:12 +0200 Subject: [PATCH 23/37] updating extended_path --- matlab/ep/extended_path.m | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/matlab/ep/extended_path.m b/matlab/ep/extended_path.m index 3cf1a08c4..cacf0530c 100644 --- a/matlab/ep/extended_path.m +++ b/matlab/ep/extended_path.m @@ -59,7 +59,7 @@ end % Set the number of periods for the perfect foresight model periods = options_.ep.periods; -pfm.periods = options_.ep.periods; +pfm.periods = periods; pfm.i_upd = pfm.ny+(1:pfm.periods*pfm.ny); % keep a copy of pfm.i_upd @@ -139,6 +139,17 @@ else pfm.dr = []; end +% number of nonzero derivatives +pfm.nnzA = M_.NNZDerivatives(1); + +% setting up integration nodes if order > 0 +if options_.ep.stochastic.order > 0 + [nodes,weights,nnodes] = setup_integration_nodes(options_.ep,pfm); + pfm.nodes = nodes; + pfm.weights = weights; + pfm.nnodes = nnodes; +end + % Main loop. while (t Date: Mon, 12 May 2014 10:18:43 +0200 Subject: [PATCH 24/37] make extended path algorithm 1 as a self contained problem usable by dynare_solve --- matlab/ep/ep_problem_2.m | 194 +++++++++ ...lve_stochastic_perfect_foresight_model_1.m | 398 ++++++------------ matlab/ep/stroud_judd_7.5.8.m | 9 + 3 files changed, 329 insertions(+), 272 deletions(-) create mode 100644 matlab/ep/ep_problem_2.m create mode 100644 matlab/ep/stroud_judd_7.5.8.m diff --git a/matlab/ep/ep_problem_2.m b/matlab/ep/ep_problem_2.m new file mode 100644 index 000000000..f646545a5 --- /dev/null +++ b/matlab/ep/ep_problem_2.m @@ -0,0 +1,194 @@ +function [res,A,info] = ep_problem_2(y,x,pm) + +info = 0; +res = []; +A = []; + +dynamic_model = pm.dynamic_model; +ny = pm.ny; +params = pm.params; +steady_state = pm.steady_state; +order = pm.order; +nodes = pm.nodes; +nnodes = pm.nnodes; +weights = pm.weights; +h_correction = pm.h_correction; +dimension = pm.dimension; +world_nbr = pm.world_nbr; +nnzA = pm.nnzA; +periods = pm.periods; +i_rows = pm.i_rows'; +i_cols = pm.i_cols; +nyp = pm.nyp; +nyf = pm.nyf; +hybrid_order = pm.hybrid_order; +i_cols_1 = pm.i_cols_1; +i_cols_j = pm.i_cols_j; +icA = pm.icA; +i_cols_T = pm.i_cols_T; + +i_cols_p = i_cols(1:nyp); +i_cols_s = i_cols(nyp+(1:ny)); +i_cols_f = i_cols(nyp+ny+(1:nyf)); +i_cols_A = i_cols; +i_cols_Ap0 = i_cols_p; +i_cols_As = i_cols_s; +i_cols_Af0 = i_cols_f - ny; +i_hc = i_cols_f - 2*ny; + +nzA = cell(periods,world_nbr); +res = zeros(ny,periods,world_nbr); +Y = zeros(ny*(periods+2),world_nbr); +Y(1:ny,1) = pm.y0; +Y(end-ny+1:end,:) = repmat(steady_state,1,world_nbr); +Y(pm.i_upd_y) = y; +offset_r0 = 0; +for i = 1:order+1 + i_w_p = 1; + for j = 1:(1+(nnodes-1)*(i-1)) + innovation = x; + if i <= order && j == 1 + % first world, integrating future shocks + if nargin > 1 + A1 = sparse([],[],[],i*ny,dimension,nnzA*world_nbr); + end + for k=1:nnodes + if nargin > 1 + if i == 2 + i_cols_Ap = i_cols_Ap0; + elseif i > 2 + i_cols_Ap = i_cols_Ap0 + ny*(i-2+(nnodes- ... + 1)*(i-2)*(i-3)/2); + end + if k == 1 + i_cols_Af = i_cols_Af0 + ny*(i-1+(nnodes-1)*i*(i-1)/ ... + 2); + else + i_cols_Af = i_cols_Af0 + ny*(i-1+(nnodes-1)*i*(i-1)/ ... + 2+(i-1)*(nnodes-1) ... + + k - 1); + end + end + if i > 1 + innovation(i+1,:) = nodes(k,:); + end + if k == 1 + k1 = 1; + else + k1 = (nnodes-1)*(i-1)+k; + end + if hybrid_order == 2 && (k > 1 || i == order) + z = [Y(i_cols_p,1); + Y(i_cols_s,1); + Y(i_cols_f,k1)+h_correction(i_hc)]; + else + z = [Y(i_cols_p,1); + Y(i_cols_s,1); + Y(i_cols_f,k1)]; + end + if nargin > 1 + [d1,jacobian] = dynamic_model(z,innovation,params,steady_state,i+1); + if i == 1 + % in first period we don't keep track of + % predetermined variables + i_cols_A = [i_cols_As - ny; i_cols_Af]; + A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(:,i_cols_1); + else + i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af]; + A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(:,i_cols_j); + end + else + d1 = dynamic_model(z,innovation,params,steady_state,i+1); + end + res(:,i,1) = res(:,i,1)+weights(k)*d1; + end + if nargin > 1 + [ir,ic,v] = find(A1); + nzA{i,j} = [ir,ic,v]'; + end + elseif j > 1 + (nnodes-1)*(i-2) + % new world, using previous state of world 1 and hit + % by shocks from nodes + if nargin > 1 + i_cols_Ap = i_cols_Ap0 + ny*(i-2+(nnodes-1)*(i-2)*(i-3)/2); + i_cols_Af = i_cols_Af0 + ny*(i+(nnodes-1)*i*(i-1)/2+j-2); + end + k = j - (nnodes-1)*(i-2); + innovation(i+1,:) = nodes(k,:); + z = [Y(i_cols_p,1); + Y(i_cols_s,j); + Y(i_cols_f,j)]; + if nargin > 1 + [d1,jacobian] = dynamic_model(z,innovation,params,steady_state,i+1); + i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af]; + [ir,ic,v] = find(jacobian(:,i_cols_j)); + nzA{i,j} = [i_rows(ir),i_cols_A(ic), v]'; + else + d1 = dynamic_model(z,innovation,params,steady_state,i+1); + end + res(:,i,j) = d1; + if nargin > 1 + i_cols_Af = i_cols_Af + ny; + end + else + % existing worlds other than 1 + if nargin > 1 + i_cols_Ap = i_cols_Ap0 + ny*(i-2+(nnodes-1)*(i-2)*(i-3)/2+j-1); + i_cols_Af = i_cols_Af0 + ny*(i+(nnodes-1)*i*(i-1)/2+j-2); + end + z = [Y(i_cols_p,j); + Y(i_cols_s,j); + Y(i_cols_f,j)]; + if nargin > 1 + [d1,jacobian] = dynamic_model(z,innovation,params,steady_state,i+1); + i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af]; + [ir,ic,v] = find(jacobian(:,i_cols_j)); + nzA{i,j} = [i_rows(ir),i_cols_A(ic),v]'; + else + d1 = dynamic_model(z,innovation,params,steady_state,i+1); + end + res(:,i,j) = d1; + i_cols_Af = i_cols_Af + ny; + end + i_rows = i_rows + ny; + if mod(j,nnodes) == 0 + i_w_p = i_w_p + 1; + end + if nargin > 1 && i > 1 + i_cols_As = i_cols_As + ny; + end + offset_r0 = offset_r0 + ny; + end + i_cols_p = i_cols_p + ny; + i_cols_s = i_cols_s + ny; + i_cols_f = i_cols_f + ny; +end +for j=1:world_nbr + i_rows_y = i_cols+(order+1)*ny; + offset_c = ny*(order+(nnodes-1)*(order-1)*order/2+j-1); + offset_r = offset_r0+(j-1)*ny; + for i=order+2:periods + if nargin > 1 + [d1,jacobian] = dynamic_model(Y(i_rows_y,j),x,params, ... + steady_state,i+1); + if i < periods + [ir,ic,v] = find(jacobian(:,i_cols_j)); + else + [ir,ic,v] = find(jacobian(:,i_cols_T)); + end + nzA{i,j} = [offset_r+ir,offset_c+icA(ic), v]'; + else + d1 = dynamic_model(Y(i_rows_y,j),x,params, ... + steady_state,i+1); + end + res(:,i,j) = d1; + i_rows_y = i_rows_y + ny; + offset_c = offset_c + world_nbr*ny; + offset_r = offset_r + world_nbr*ny; + end +end +if nargin > 1 + iA = [nzA{:}]'; + A = sparse(iA(:,1),iA(:,2),iA(:,3),dimension,dimension); +end +res = res(pm.i_upd_r); \ No newline at end of file diff --git a/matlab/ep/solve_stochastic_perfect_foresight_model_1.m b/matlab/ep/solve_stochastic_perfect_foresight_model_1.m index b844f8ce8..ddebc0dbe 100644 --- a/matlab/ep/solve_stochastic_perfect_foresight_model_1.m +++ b/matlab/ep/solve_stochastic_perfect_foresight_model_1.m @@ -1,4 +1,4 @@ -function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo_simul,exo_simul,EpOptions,pfm,order,varargin) +function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo_simul,exo_simul,Options,pfm,order,varargin) % Copyright (C) 2012-2013 Dynare Team % @@ -17,281 +17,135 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . - if nargin < 6 - homotopy_parameter = 1; - else - homotopy_parameter = varargin{1}; +global options_ + +if nargin < 6 + homotopy_parameter = 1; +else + homotopy_parameter = varargin{1}; +end +flag = 0; +err = 0; +stop = 0; + +EpOptions = Options.ep; + +params = pfm.params; +steady_state = pfm.steady_state; +ny = pfm.ny; +periods = pfm.periods; +dynamic_model = pfm.dynamic_model; +lead_lag_incidence = pfm.lead_lag_incidence; +nyp = pfm.nyp; +nyf = pfm.nyf; +i_cols_1 = pfm.i_cols_1; +i_cols_A1 = pfm.i_cols_A1; +i_cols_j = pfm.i_cols_j; +i_cols_T = nonzeros(lead_lag_incidence(1:2,:)'); +hybrid_order = pfm.hybrid_order; +dr = pfm.dr; +nodes = pfm.nodes; +weights = pfm.weights; +nnodes = pfm.nnodes; + +maxit = pfm.maxit_; +tolerance = pfm.tolerance; +verbose = pfm.verbose; + +number_of_shocks = size(exo_simul,2); + +% make sure that there is a node equal to zero +% and permute nodes and weights to have zero first +k = find(sum(abs(nodes),2) < 1e-12); +if ~isempty(k) + nodes = [nodes(k,:); nodes(1:k-1,:); nodes(k+1:end,:)]; + weights = [weights(k); weights(1:k-1); weights(k+1:end)]; +else + error('there is no nodes equal to zero') +end + +if hybrid_order > 0 + if hybrid_order == 2 + h_correction = 0.5*dr.ghs2(dr.inv_order_var); end - flag = 0; - err = 0; - stop = 0; +else + h_correction = 0; +end - params = pfm.params; - steady_state = pfm.steady_state; - ny = pfm.ny; - periods = pfm.periods; - dynamic_model = pfm.dynamic_model; - lead_lag_incidence = pfm.lead_lag_incidence; - nyp = pfm.nyp; - nyf = pfm.nyf; - i_cols_1 = pfm.i_cols_1; - i_cols_A1 = pfm.i_cols_A1; - i_cols_j = pfm.i_cols_j; - i_cols_T = nonzeros(lead_lag_incidence(1:2,:)'); - hybrid_order = pfm.hybrid_order; - dr = pfm.dr; - [nodes,weights,nnodes] = setup_integration_nodes(EpOptions,pfm); - - maxit = pfm.maxit_; - tolerance = pfm.tolerance; - verbose = pfm.verbose; +if verbose + disp ([' -----------------------------------------------------']); + disp (['MODEL SIMULATION :']); + fprintf('\n'); +end - number_of_shocks = size(exo_simul,2); +% Each column of Y represents a different world +% The upper right cells are unused +% The first row block is ny x 1 +% The second row block is ny x nnodes +% The third row block is ny x nnodes^2 +% and so on until size ny x nnodes^order +world_nbr = 1+(nnodes-1)*order; +Y = endo_simul(:,2:end-1); +Y = repmat(Y,1,world_nbr); +pfm.y0 = endo_simul(:,1); - % make sure that there is a node equal to zero - % and permute nodes and weights to have zero first - k = find(sum(abs(nodes),2) < 1e-12); - if ~isempty(k) - nodes = [nodes(k,:); nodes(1:k-1,:); nodes(k+1:end,:)]; - weights = [weights(k); weights(1:k-1); weights(k+1:end)]; - else - error('there is no nodes equal to zero') - end - - if hybrid_order > 0 - if hybrid_order == 2 - h_correction = 0.5*dr.ghs2(dr.inv_order_var); - end +% The columns of A map the elements of Y such that +% each block of Y with ny rows are unfolded column wise +% number of blocks +block_nbr = (order+(nnodes-1)*(order-1)*order/2+(periods-order)*world_nbr); +% dimension of the problem +dimension = ny*block_nbr; +pfm.block_nbr = block_nbr; +pfm.dimension = dimension; +if order == 0 + i_upd_r = (1:ny*periods); + i_upd_y = i_upd_r + ny; +else + i_upd_r = zeros(dimension,1); + i_upd_y = i_upd_r; + i_upd_r(1:ny) = (1:ny); + i_upd_y(1:ny) = ny+(1:ny); + i1 = ny+1; + i2 = 2*ny; + n1 = ny+1; + n2 = 2*ny; + for i=2:periods + k = n1:n2; + for j=1:(1+(nnodes-1)*min(i-1,order)) + i_upd_r(i1:i2) = k+(j-1)*ny*periods; + i_upd_y(i1:i2) = k+ny+(j-1)*ny*(periods+2); + i1 = i2+1; + i2 = i2+ny; + end + n1 = n2+1; + n2 = n2+ny; end +end +icA = [find(lead_lag_incidence(1,:)) find(lead_lag_incidence(2,:))+world_nbr*ny ... + find(lead_lag_incidence(3,:))+2*world_nbr*ny]'; +h1 = clock; +pfm.order = order; +pfm.world_nbr = world_nbr; +pfm.nodes = nodes; +pfm.nnodes = nnodes; +pfm.weights = weights; +pfm.h_correction = h_correction; +pfm.i_rows = 1:ny; +i_cols = find(lead_lag_incidence'); +pfm.i_cols = i_cols; +pfm.nyp = nyp; +pfm.nyf = nyf; +pfm.hybrid_order = hybrid_order; +pfm.i_cols_1 = i_cols_1; +pfm.i_cols_h = i_cols_j; +pfm.icA = icA; +pfm.i_cols_T = i_cols_T; +pfm.i_upd_r = i_upd_r; +pfm.i_upd_y = i_upd_y; - if verbose - disp ([' -----------------------------------------------------']); - disp (['MODEL SIMULATION :']); - fprintf('\n'); - end - z = endo_simul(find(lead_lag_incidence')); - [d1,jacobian] = dynamic_model(z,exo_simul,params,steady_state,2); - - % Each column of Y represents a different world - % The upper right cells are unused - % The first row block is ny x 1 - % The second row block is ny x nnodes - % The third row block is ny x nnodes^2 - % and so on until size ny x nnodes^order - world_nbr = 1+(nnodes-1)*order; - Y = repmat(endo_simul(:),1,world_nbr); - - % The columns of A map the elements of Y such that - % each block of Y with ny rows are unfolded column wise - dimension = ny*(order+(nnodes-1)*(order-1)*order/2+(periods-order)*world_nbr); - if order == 0 - i_upd_r = (1:ny*periods); - i_upd_y = i_upd_r + ny; - else - i_upd_r = zeros(dimension,1); - i_upd_y = i_upd_r; - i_upd_r(1:ny) = (1:ny); - i_upd_y(1:ny) = ny+(1:ny); - i1 = ny+1; - i2 = 2*ny; - n1 = ny+1; - n2 = 2*ny; - for i=2:periods - k = n1:n2; - for j=1:(1+(nnodes-1)*min(i-1,order)) - i_upd_r(i1:i2) = k+(j-1)*ny*periods; - i_upd_y(i1:i2) = k+ny+(j-1)*ny*(periods+2); - i1 = i2+1; - i2 = i2+ny; - end - n1 = n2+1; - n2 = n2+ny; - end - end - icA = [find(lead_lag_incidence(1,:)) find(lead_lag_incidence(2,:))+world_nbr*ny ... - find(lead_lag_incidence(3,:))+2*world_nbr*ny]'; - h1 = clock; - for iter = 1:maxit - h2 = clock; - A1 = sparse([],[],[],ny*(order+(nnodes-1)*(order-1)*order/2),dimension,(order+1)*world_nbr*nnz(jacobian)); - res = zeros(ny,periods,world_nbr); - i_rows = 1:ny; - i_cols = find(lead_lag_incidence'); - i_cols_p = i_cols(1:nyp); - i_cols_s = i_cols(nyp+(1:ny)); - i_cols_f = i_cols(nyp+ny+(1:nyf)); - i_cols_A = i_cols; - i_cols_Ap0 = i_cols_p; - i_cols_As = i_cols_s; - i_cols_Af0 = i_cols_f - ny; - i_hc = i_cols_f - 2*ny; - for i = 1:order+1 - i_w_p = 1; - for j = 1:(1+(nnodes-1)*(i-1)) - innovation = exo_simul; - if i <= order && j == 1 - % first world, integrating future shocks - for k=1:nnodes - if i == 2 - i_cols_Ap = i_cols_Ap0; - elseif i > 2 - i_cols_Ap = i_cols_Ap0 + ny*(i-2+(nnodes- ... - 1)*(i-2)*(i-3)/2); - end - if k == 1 - i_cols_Af = i_cols_Af0 + ny*(i-1+(nnodes-1)*i*(i-1)/ ... - 2); - else - i_cols_Af = i_cols_Af0 + ny*(i-1+(nnodes-1)*i*(i-1)/ ... - 2+(i-1)*(nnodes-1) ... - + k - 1); - end - if i > 1 - innovation(i+1,:) = nodes(k,:); - end - if k == 1 - k1 = 1; - else - k1 = (nnodes-1)*(i-1)+k; - end - if hybrid_order == 2 && (k > 1 || i == order) - y = [Y(i_cols_p,1); - Y(i_cols_s,1); - Y(i_cols_f,k1)+h_correction(i_hc)]; - else - y = [Y(i_cols_p,1); - Y(i_cols_s,1); - Y(i_cols_f,k1)]; - end - [d1,jacobian] = dynamic_model(y,homotopy_parameter*innovation, ... - params,steady_state,i+1); - if i == 1 - % in first period we don't keep track of - % predetermined variables - i_cols_A = [i_cols_As - ny; i_cols_Af]; - A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(:,i_cols_1); - else - i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af]; - A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(:,i_cols_j); - end - res(:,i,1) = res(:,i,1)+weights(k)*d1; - end - elseif j > 1 + (nnodes-1)*(i-2) - % new world, using previous state of world 1 and hit - % by shocks from nodes - i_cols_Ap = i_cols_Ap0 + ny*(i-2+(nnodes-1)*(i-2)*(i-3)/2); - i_cols_Af = i_cols_Af0 + ny*(i+(nnodes-1)*i*(i-1)/2+j-2); - k = j - (nnodes-1)*(i-2); - innovation(i+1,:) = nodes(k,:); - y = [Y(i_cols_p,1); - Y(i_cols_s,j); - Y(i_cols_f,j)]; - [d1,jacobian] = dynamic_model(y,homotopy_parameter*innovation,params,steady_state,i+1); - i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af]; - A1(i_rows,i_cols_A) = jacobian(:,i_cols_j); - res(:,i,j) = d1; - i_cols_Af = i_cols_Af + ny; - else - % existing worlds other than 1 - i_cols_Ap = i_cols_Ap0 + ny*(i-2+(nnodes-1)*(i-2)*(i-3)/2+j-1); - i_cols_Af = i_cols_Af0 + ny*(i+(nnodes-1)*i*(i-1)/2+j-2); - y = [Y(i_cols_p,j); - Y(i_cols_s,j); - Y(i_cols_f,j)]; - [d1,jacobian] = dynamic_model(y,homotopy_parameter*innovation,params,steady_state,i+1); - i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af]; - A1(i_rows,i_cols_A) = jacobian(:,i_cols_j); - res(:,i,j) = d1; - i_cols_Af = i_cols_Af + ny; - end - i_rows = i_rows + ny; - if mod(j,nnodes) == 0 - i_w_p = i_w_p + 1; - end - if i > 1 - i_cols_As = i_cols_As + ny; - end - end - i_cols_p = i_cols_p + ny; - i_cols_s = i_cols_s + ny; - i_cols_f = i_cols_f + ny; - end - nzA = cell(periods,world_nbr); - parfor j=1:world_nbr - i_rows_y = find(lead_lag_incidence')+(order+1)*ny; - offset_c = ny*(order+(nnodes-1)*(order-1)*order/2+j-1); - offset_r = (j-1)*ny; - for i=order+2:periods - [d1,jacobian] = dynamic_model(Y(i_rows_y,j), ... - exo_simul,params, ... - steady_state,i+1); - if i == periods - [ir,ic,v] = find(jacobian(:,i_cols_T)); - else - [ir,ic,v] = find(jacobian(:,i_cols_j)); - end - nzA{i,j} = [offset_r+ir,offset_c+icA(ic), v]'; - res(:,i,j) = d1; - i_rows_y = i_rows_y + ny; - offset_c = offset_c + world_nbr*ny; - offset_r = offset_r + world_nbr*ny; - end - end - err = max(abs(res(i_upd_r))); - if verbose - [err1, k1] = max(abs(res)); - [err2, k2] = max(abs(err1)); - [err3, k3] = max(abs(err2)); - disp([iter err k1(:,k2(:,:,k3),k3) k2(:,:,k3) k3]) - end - if err < tolerance - stop = 1; - flag = 0;% Convergency obtained. - endo_simul = reshape(Y(:,1),ny,periods+2);%Y(ny+(1:ny),1); - if verbose - save ep_test_s1 exo_simul endo_simul Y - fprintf('\n') ; - disp([' Total time of simulation :' num2str(etime(clock,h1))]) ; - fprintf('\n') ; - disp([' Convergency obtained.']) ; - fprintf('\n') ; - end - break - end - A2 = [nzA{:}]'; - if any(isnan(A2(:,3))) || any(any(any(isnan(res)))) - if verbose - disp(['solve_stochastic_foresight_model_1 encountered ' ... - 'NaN']) - save ep_test_s2 exo_simul endo_simul - pause - end - flag = 1; - return - end - A = [A1; sparse(A2(:,1),A2(:,2),A2(:,3),ny*(periods-order-1)*world_nbr,dimension)]; - if verbose - disp(sprintf('condest %g',condest(A))) - end - dy = -A\res(i_upd_r); - Y(i_upd_y) = Y(i_upd_y) + dy; - end - - if ~stop - if verbose - fprintf('\n') ; - disp([' Total time of simulation :' num2str(etime(clock,h1))]) ; - fprintf('\n') ; - disp(['WARNING : maximum number of iterations is reached (modify options_.simul.maxit).']) ; - fprintf('\n') ; - disp(sprintf('err: %f',err)); - save ep_test_s2 exo_simul endo_simul - pause - end - flag = 1;% more iterations are needed. - endo_simul = 1; - end - if verbose - disp (['-----------------------------------------------------']) ; - end \ No newline at end of file +options_.solve_algo = 9; +options_.steady.maxit = 100; +y = repmat(steady_state,block_nbr,1); +y = dynare_solve(@ep_problem_2,y,1,exo_simul,pfm); +endo_simul(:,2) = y(1:ny); \ No newline at end of file diff --git a/matlab/ep/stroud_judd_7.5.8.m b/matlab/ep/stroud_judd_7.5.8.m new file mode 100644 index 000000000..40c9bc08d --- /dev/null +++ b/matlab/ep/stroud_judd_7.5.8.m @@ -0,0 +1,9 @@ +function [X,w]=stroud_judd_7.5.8(d) + + E = eye(d); + X = cell(2*d,1); + m = 1; + for i=1:d + X{m} = E(:,i); + m = m+1; + X{m} = -E(:,i); From 979101fa55643eb4a713c87f36752cf91e7793a9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Mon, 12 May 2014 12:14:28 +0200 Subject: [PATCH 25/37] sim1: when there are NaNs/Infs, raise a warning rather than an error. Otherwise the homotopy procedure can fail prematurely. Thanks to Tom Holden for the suggestion. --- matlab/sim1.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/sim1.m b/matlab/sim1.m index e08f57a17..87517461a 100644 --- a/matlab/sim1.m +++ b/matlab/sim1.m @@ -171,7 +171,7 @@ if stop skipline(); fprintf('\nSimulation terminated after %d iterations.\n',iter); fprintf('Total time of simulation: %16.13f\n',etime(clock,h1)); - error('Simulation terminated with NaN or Inf in the residuals or endogenous variables. There is most likely something wrong with your model.'); + fprintf('WARNING: Simulation terminated with NaN or Inf in the residuals or endogenous variables. There is most likely something wrong with your model.\n'); else skipline(); fprintf('\nSimulation concluded successfully after %d iterations.\n',iter); From c83eac26bb5ddb48645e1810be6428dc2919d489 Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Mon, 12 May 2014 12:34:08 +0200 Subject: [PATCH 26/37] extended_preprocessor C API: removing C++ statement from header file (C API implementation not finished) --- others/C/dynare_driver.h | 10 +--------- 1 file changed, 1 insertion(+), 9 deletions(-) diff --git a/others/C/dynare_driver.h b/others/C/dynare_driver.h index dfaa7c60a..56091c7d7 100644 --- a/others/C/dynare_driver.h +++ b/others/C/dynare_driver.h @@ -1,5 +1,5 @@ /* - * Copyright (C) 2011-2012 Houtan Bastani, Daniel Waggoner, Tao Zha + * Copyright (C) 2014 DynareTeam * * This is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -17,14 +17,6 @@ #ifndef _DYNARE_C_DRIVER_H #define _DYNARE_C_DRIVER_H -#include -#include -#include -#include -#include - -using namespace std; - struct aux_vars_t { int endo_index, type, orig_index, orig_lead_lag; } ; From 23134353e353be74c29384d2d34755aa7186d5d0 Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Mon, 12 May 2014 14:17:19 +0200 Subject: [PATCH 27/37] extended_path: update tests models; correct bugs introduced in previous commit --- matlab/ep/ep_problem_2.m | 26 +++++++++---------- matlab/ep/homotopic_steps.m | 6 ++--- ...lve_stochastic_perfect_foresight_model_1.m | 6 ++++- tests/ep/ar.mod | 1 + tests/ep/burnside.mod | 1 + tests/ep/linearmodel.mod | 1 + tests/ep/rbc.mod | 1 + tests/ep/rbc2.mod | 1 + tests/ep/rbcii.mod | 1 + 9 files changed, 26 insertions(+), 18 deletions(-) diff --git a/matlab/ep/ep_problem_2.m b/matlab/ep/ep_problem_2.m index f646545a5..634ae918c 100644 --- a/matlab/ep/ep_problem_2.m +++ b/matlab/ep/ep_problem_2.m @@ -49,11 +49,11 @@ for i = 1:order+1 innovation = x; if i <= order && j == 1 % first world, integrating future shocks - if nargin > 1 + if nargout > 1 A1 = sparse([],[],[],i*ny,dimension,nnzA*world_nbr); end for k=1:nnodes - if nargin > 1 + if nargout > 1 if i == 2 i_cols_Ap = i_cols_Ap0; elseif i > 2 @@ -86,7 +86,7 @@ for i = 1:order+1 Y(i_cols_s,1); Y(i_cols_f,k1)]; end - if nargin > 1 + if nargout > 1 [d1,jacobian] = dynamic_model(z,innovation,params,steady_state,i+1); if i == 1 % in first period we don't keep track of @@ -102,14 +102,14 @@ for i = 1:order+1 end res(:,i,1) = res(:,i,1)+weights(k)*d1; end - if nargin > 1 + if nargout > 1 [ir,ic,v] = find(A1); nzA{i,j} = [ir,ic,v]'; end elseif j > 1 + (nnodes-1)*(i-2) % new world, using previous state of world 1 and hit % by shocks from nodes - if nargin > 1 + if nargout > 1 i_cols_Ap = i_cols_Ap0 + ny*(i-2+(nnodes-1)*(i-2)*(i-3)/2); i_cols_Af = i_cols_Af0 + ny*(i+(nnodes-1)*i*(i-1)/2+j-2); end @@ -118,7 +118,7 @@ for i = 1:order+1 z = [Y(i_cols_p,1); Y(i_cols_s,j); Y(i_cols_f,j)]; - if nargin > 1 + if nargout > 1 [d1,jacobian] = dynamic_model(z,innovation,params,steady_state,i+1); i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af]; [ir,ic,v] = find(jacobian(:,i_cols_j)); @@ -127,34 +127,34 @@ for i = 1:order+1 d1 = dynamic_model(z,innovation,params,steady_state,i+1); end res(:,i,j) = d1; - if nargin > 1 + if nargout > 1 i_cols_Af = i_cols_Af + ny; end else % existing worlds other than 1 - if nargin > 1 + if nargout > 1 i_cols_Ap = i_cols_Ap0 + ny*(i-2+(nnodes-1)*(i-2)*(i-3)/2+j-1); i_cols_Af = i_cols_Af0 + ny*(i+(nnodes-1)*i*(i-1)/2+j-2); end z = [Y(i_cols_p,j); Y(i_cols_s,j); Y(i_cols_f,j)]; - if nargin > 1 + if nargout > 1 [d1,jacobian] = dynamic_model(z,innovation,params,steady_state,i+1); i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af]; [ir,ic,v] = find(jacobian(:,i_cols_j)); nzA{i,j} = [i_rows(ir),i_cols_A(ic),v]'; + i_cols_Af = i_cols_Af + ny; else d1 = dynamic_model(z,innovation,params,steady_state,i+1); end res(:,i,j) = d1; - i_cols_Af = i_cols_Af + ny; end i_rows = i_rows + ny; if mod(j,nnodes) == 0 i_w_p = i_w_p + 1; end - if nargin > 1 && i > 1 + if nargout > 1 && i > 1 i_cols_As = i_cols_As + ny; end offset_r0 = offset_r0 + ny; @@ -168,7 +168,7 @@ for j=1:world_nbr offset_c = ny*(order+(nnodes-1)*(order-1)*order/2+j-1); offset_r = offset_r0+(j-1)*ny; for i=order+2:periods - if nargin > 1 + if nargout > 1 [d1,jacobian] = dynamic_model(Y(i_rows_y,j),x,params, ... steady_state,i+1); if i < periods @@ -187,7 +187,7 @@ for j=1:world_nbr offset_r = offset_r + world_nbr*ny; end end -if nargin > 1 +if nargout > 1 iA = [nzA{:}]'; A = sparse(iA(:,1),iA(:,2),iA(:,3),dimension,dimension); end diff --git a/matlab/ep/homotopic_steps.m b/matlab/ep/homotopic_steps.m index 14bf7ebb9..c4839c46c 100644 --- a/matlab/ep/homotopic_steps.m +++ b/matlab/ep/homotopic_steps.m @@ -73,7 +73,7 @@ while weight<1 solve_stochastic_perfect_foresight_model(endo_simul0,exo_simul0,pfm,options_.ep.stochastic.quadrature.nodes,options_.ep.stochastic.order); case 1 [flag,tmp] = ... - solve_stochastic_perfect_foresight_model_1(endo_simul0,exxo_simul,options_.ep,pfm,options_.ep.stochastic.order,weight); + solve_stochastic_perfect_foresight_model_1(endo_simul0,exxo_simul,options_,pfm,options_.ep.stochastic.order,weight); % solve_stochastic_perfect_foresight_model_1(endo_simul0,exo_simul0,options_.ep,pfm,options_.ep.stochastic.order); end end @@ -114,7 +114,6 @@ while weight<1 endo_simul0 = endo_simul; exo_simul0 = exxo_simul; info.convergence = 0; - info.depth = d; tmp = []; return end @@ -159,7 +158,7 @@ if weight<1 solve_stochastic_perfect_foresight_model(endo_simul0,exo_simul0,pfm,options_.ep.stochastic.quadrature.nodes,options_.ep.stochastic.order); case 1 [flag,tmp] = ... - solve_stochastic_perfect_foresight_model_1(endo_simul0,exxo_simul,options_.ep,pfm,options_.ep.stochastic.order,weight); + solve_stochastic_perfect_foresight_model_1(endo_simul0,exxo_simul,options_,pfm,options_.ep.stochastic.order,weight); % solve_stochastic_perfect_foresight_model_1(endo_simul0,exo_simul0,options_.ep,pfm,options_.ep.stochastic.order); end end @@ -173,7 +172,6 @@ if weight<1 endo_simul0 = endo_simul; exo_simul0 = exxo_simul; info.convergence = 0; - info.depth = d; tmp = []; return else diff --git a/matlab/ep/solve_stochastic_perfect_foresight_model_1.m b/matlab/ep/solve_stochastic_perfect_foresight_model_1.m index ddebc0dbe..50fb79299 100644 --- a/matlab/ep/solve_stochastic_perfect_foresight_model_1.m +++ b/matlab/ep/solve_stochastic_perfect_foresight_model_1.m @@ -147,5 +147,9 @@ pfm.i_upd_y = i_upd_y; options_.solve_algo = 9; options_.steady.maxit = 100; y = repmat(steady_state,block_nbr,1); -y = dynare_solve(@ep_problem_2,y,1,exo_simul,pfm); +[y,info] = dynare_solve(@ep_problem_2,y,1,exo_simul,pfm); +if info + flag = 1; + err = info; +end endo_simul(:,2) = y(1:ny); \ No newline at end of file diff --git a/tests/ep/ar.mod b/tests/ep/ar.mod index a973d1e77..846d74a86 100644 --- a/tests/ep/ar.mod +++ b/tests/ep/ar.mod @@ -40,6 +40,7 @@ ts = extended_path([],1000); options_.ep.verbosity = 0; options_.ep.stochastic.order = 1; +options_.ep.IntegrationAlgorithm='Tensor-Gaussian-Quadrature'; options_.ep.stochastic.nodes = 3; options_.console_mode = 0; diff --git a/tests/ep/burnside.mod b/tests/ep/burnside.mod index 0bad36d87..7253d2f0b 100644 --- a/tests/ep/burnside.mod +++ b/tests/ep/burnside.mod @@ -53,6 +53,7 @@ set_dynare_seed('default'); ts = extended_path([],5000); options_.ep.stochastic.order = 2; +options_.ep.IntegrationAlgorithm='Tensor-Gaussian-Quadrature'; set_dynare_seed('default'); ts1_4 = extended_path([],5000); diff --git a/tests/ep/linearmodel.mod b/tests/ep/linearmodel.mod index f773e651c..fde29d0db 100644 --- a/tests/ep/linearmodel.mod +++ b/tests/ep/linearmodel.mod @@ -36,6 +36,7 @@ options_.console_mode = 0; ts = extended_path([],100); options_.ep.stochastic.status = 1; +options_.ep.IntegrationAlgorithm='Tensor-Gaussian-Quadrature'; options_.ep.order = 1; options_.ep.nnodes = 3; sts = extended_path([],100); diff --git a/tests/ep/rbc.mod b/tests/ep/rbc.mod index cc091e911..9bdbf6a7e 100644 --- a/tests/ep/rbc.mod +++ b/tests/ep/rbc.mod @@ -79,6 +79,7 @@ ts0 = extended_path([],100); options_.ep.stochastic.order = 1; options_.ep.stochastic.nodes = 3; +options_.ep.IntegrationAlgorithm='Tensor-Gaussian-Quadrature'; ts1_3 = extended_path([],100); options_.ep.stochastic.nodes = 5; diff --git a/tests/ep/rbc2.mod b/tests/ep/rbc2.mod index 6384afadc..ff8242843 100644 --- a/tests/ep/rbc2.mod +++ b/tests/ep/rbc2.mod @@ -76,5 +76,6 @@ end; steady; +options_.ep.IntegrationAlgorithm='Tensor-Gaussian-Quadrature'; extended_path(periods=100); diff --git a/tests/ep/rbcii.mod b/tests/ep/rbcii.mod index 186679731..4966f47e6 100644 --- a/tests/ep/rbcii.mod +++ b/tests/ep/rbcii.mod @@ -75,6 +75,7 @@ copyfile('rbcii_steady_state.m','rbcii_steadystate2.m'); ts = extended_path([],100); options_.ep.stochastic.order = 1; + options_.ep.IntegrationAlgorithm='Tensor-Gaussian-Quadrature'; // profile on ts1_4 = extended_path([],100); // profile off From c45f45f58ed144ede58211614101ba763fdbf785 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Tue, 13 May 2014 11:46:19 +0200 Subject: [PATCH 28/37] Forbid model local variables in planner_objective. Otherwise the preprocessor crashes. --- preprocessor/ParsingDriver.cc | 3 +++ 1 file changed, 3 insertions(+) diff --git a/preprocessor/ParsingDriver.cc b/preprocessor/ParsingDriver.cc index 70ebe8712..517b0d882 100644 --- a/preprocessor/ParsingDriver.cc +++ b/preprocessor/ParsingDriver.cc @@ -323,6 +323,9 @@ ParsingDriver::add_model_variable(int symb_id, int lag) if (dynamic_cast(model_tree) != NULL && lag != 0) error("Leads and lags on variables are forbidden in 'planner_objective'."); + if (dynamic_cast(model_tree) != NULL && type == eModelLocalVariable) + error("Model local variable " + mod_file->symbol_table.getName(symb_id) + " cannot be used in 'planner_objective'."); + // It makes sense to allow a lead/lag on parameters: during steady state calibration, endogenous and parameters can be swapped return model_tree->AddVariable(symb_id, lag); } From f21bed4f1f12a3fbd06841aee50d90aa751171f7 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Tue, 13 May 2014 13:00:28 +0200 Subject: [PATCH 29/37] doc: add default values for baxter_king_filter --- doc/dynare.texi | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/dynare.texi b/doc/dynare.texi index 7ddd04dd0..ec7ab1a24 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -9871,7 +9871,7 @@ ts1 is a dseries object: @deftypefn {dseries} {@var{B} = } baxter_king_filter (@var{A}, @var{hf}, @var{lf}, @var{K}) -Implementation of the Baxter and King (1999) band pass filter for @dseries objects. This filter isolates business cycle fluctuations with a period of length ranging between @var{hf} (high frequency) to @var{lf} (low frequency) using a symmetric moving average smoother with @math{2K+1} points, so that K observations at the beginning and at the end of the sample are lost in the computation of the filter. +Implementation of the Baxter and King (1999) band pass filter for @dseries objects. This filter isolates business cycle fluctuations with a period of length ranging between @var{hf} (high frequency) to @var{lf} (low frequency) using a symmetric moving average smoother with @math{2K+1} points, so that K observations at the beginning and at the end of the sample are lost in the computation of the filter. The default value for @var{hf} is @math{6}, for @var{lf} is @math{32}, and for @var{K} is 12. @examplehead @example From 9e86c760b3e5a07c47ca5fbee28453488ccdce16 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Tue, 13 May 2014 18:30:50 +0200 Subject: [PATCH 30/37] doc: finish dseries syntax for slides --- .../dseriesReporting.tex | 210 +++++++++++++++++- 1 file changed, 207 insertions(+), 3 deletions(-) diff --git a/doc/dseries-and-reporting/dseriesReporting.tex b/doc/dseries-and-reporting/dseriesReporting.tex index 5b2643749..694e9fd11 100644 --- a/doc/dseries-and-reporting/dseriesReporting.tex +++ b/doc/dseries-and-reporting/dseriesReporting.tex @@ -58,8 +58,9 @@ \begin{itemize} \item Provide support for time series in Dynare \item Introduced in Dynare 4.4 - \item Currently only used for reporting. + \item Currently only used for reporting \item Use will increase with time (\textit{e.g.,} to be included in new estimation code) + \item NB: More complete information is included in the Dynare manual \end{itemize} \end{frame} @@ -101,6 +102,12 @@ \item Empty \item Noncontiguous \end{itemize} + \item Can index a \texttt{dates} object. If \texttt{t} is a \texttt{dates} object, then + \begin{itemize} + \item \texttt{t(1) \% refers to the first date} + \item \texttt{t(1:3) \% refers to the first three dates} + \item \texttt{t([1,4,5]) \% refers to the first, fourth, and fifth dates} + \end{itemize} \end{itemize} \end{frame} @@ -118,6 +125,7 @@ \item In a \texttt{.m} file: \texttt{dr = dates(`1999y'):dates(`2020y');} \item In a \texttt{.mod} file: \texttt{dr = 1999y:2020y;} \end{itemize} + \item In what follows, \texttt{t} and \texttt{dr} are as above. Operations that assign to \texttt{t} and \texttt{dr} keep the value thereafter. \end{itemize} \end{frame} @@ -149,13 +157,13 @@ \end{itemize} \item \texttt{sort}: sort dates in ascending order \begin{itemize} - \item \texttt{t=t.sort(); \% } + \item \texttt{t.sort(); \% } \end{itemize} \item \texttt{uminus}: shifts dates back one period \begin{itemize} \item \texttt{-t; \% } \end{itemize} - \item \texttt{unique}: removes repetitions + \item \texttt{unique}: removes repetitions (keeping last unique value) \begin{itemize} \item \texttt{t.append(dates(`1999y')).unique() \% } \end{itemize} @@ -217,6 +225,7 @@ \item \texttt{isempty}: returns true if the argument is empty \begin{itemize} \item \texttt{isempty(t) \% 0} + \item \texttt{isempty(dates()) \% 1} \end{itemize} \item \texttt{setdiff}: returns dates present in first arg but not in second \begin{itemize} @@ -234,6 +243,198 @@ % DSERIES % \subsubsection{\texttt{dseries} Syntax} +\begin{frame}[fragile,t] + \frametitle{\texttt{dseries} Syntax} + \begin{itemize} + \item A \texttt{dseries} is composed of one or more individual series + \item All time series in a dseries must have the same frequency + \item A \texttt{dseries} runs from the earliest date to the latest date, with \texttt{NaN}'s inserted to pad the shorter series + \end{itemize} +\end{frame} + + +\begin{frame}[fragile,t] + \frametitle{Creating a new \texttt{dseries} object} + \begin{itemize} + \item Load series from data file (\texttt{.csv, .xls}). Dates in first column, optional variable names in first row + \begin{itemize} + \item \texttt{ts=dseries(`data.csv');} + \end{itemize} + \item Load series from a Matlab/Octave file (\texttt{.m, .mat}). Must define the variables \texttt{INIT\_\_}, \texttt{NAMES\_\_}, and, optionally, \texttt{TEX\_\_}. + \begin{itemize} + \item \texttt{ts=dseries(`data.m');} + \end{itemize} + \item Load data directly. Here: one variable, `MyVar1', with three annual observations starting in 1999. Only first argument is required. + \begin{itemize} + \item \texttt{ts=dseries([1;2;3], `1999y', \{`MyVar1'\}, \{`MyVar\_1'\});} + \end{itemize} + \item Create an empty time series. Usefull for programatically creating a series. + \begin{itemize} + \item \texttt{ts=dseries();} + \item \texttt{ts=dseries(dates(`1999y'));} + \end{itemize} + \end{itemize} +\end{frame} + + +\begin{frame}[fragile,t] + \frametitle{Referencing data from a \texttt{dseries}} + \begin{itemize} + \item For the following, let + \begin{itemize} + \item \texttt{ts1=dseries([1;2;3], `1999y', \{`MyVar1'\}, \{`MyVar\_1'\})} + \item \texttt{ts2=dseries([4;5;6], `2000y', \{`MyVar2'\}, \{`MyVar\_2'\})} + \item \texttt{ts3=[ts1 ts2]} + \end{itemize} + \item To get \texttt{MyVar1} from \texttt{t3}, you have three options + \begin{itemize} + \item \texttt{ts3\{1\}} + \item \texttt{ts3.MyVar1} + \item \texttt{ts3\{`MyVar1'\}} + \end{itemize} + \item To select a subsample of MyVar2 from 2001 to 2002: + \begin{itemize} + \item \texttt{ts3\{2\}(dates(`2001y'):dates(`2002y'))} + \item \texttt{ts3.MyVar2(dates(`2001y'):dates(`2002y'))} + \item \texttt{ts3\{`MyVar2'\}(dates(`2001y'):dates(`2002y'))} + \end{itemize} + \item To see the data for both variables in 2000: + \begin{itemize} + \item \texttt{ts3(`2000y')} + \item \texttt{ts3\{[1 2]\}(`2000y')} + \item \texttt{ts3\{`MyVar1',`MyVar2'\}(`2000y')} + \item \texttt{ts3\{`MyVar@1,2@'\}(`2000y')} + \item \texttt{ts3\{`MyVar[1-2]'\}(`2000y')} + \end{itemize} + \end{itemize} +\end{frame} + + +\begin{frame}[fragile,t] + \frametitle{Getting info about \texttt{dseries}} + \begin{itemize} + \item \texttt{eq, ne}: returns boolean value of element-wise comparison; series must have same start dates and number of observations + \begin{itemize} + \item \texttt{ts1==ts3.MyVar1(dates(`1999y'):dates(`2001y')) \% [1 1 1]'} + \item \texttt{ts1$\thicksim$=ts1 \% [0 0 0]'} + \end{itemize} + \item \texttt{isempty}: returns true if series is empty + \begin{itemize} + \item \texttt{isempty(dseries()) \% 1} + \end{itemize} + \item \texttt{isequal}: returns true if the series are equal + \begin{itemize} + \item \texttt{isequal(ts1,ts1) \% 1} + \end{itemize} + \item \texttt{size}: returns number of observations by number of variables + \begin{itemize} + \item \texttt{ts3.size() \% [4 2]} + \end{itemize} + \end{itemize} +\end{frame} + + +\begin{frame}[fragile,t] + \frametitle{Working with \texttt{dseries}} + \begin{itemize} + \item \texttt{baxter\_king\_filter}: the Baxter and King (1999) band pass filter + \begin{itemize} + \item \texttt{ts1.baxter\_king\_filter(high freq, low freq, K)} + \end{itemize} + \item \texttt{hpcycle}: HP Filters the \texttt{dseries}, returning the business cycle component + \begin{itemize} + \item \texttt{ts1.hptrend(lambda) \% lambda = 1600 by default} + \end{itemize} + \item \texttt{hptrend}: HP Filters the \texttt{dseries}, returning the trend component + \begin{itemize} + \item \texttt{ts1.hptrend(lambda) \% lambda = 1600 by default} + \end{itemize} + \item \texttt{qdiff}: Quarterly difference; works on quarterly, monthly, and weekly series + \begin{itemize} + \item \texttt{ts1.qdiff()} + \end{itemize} + \item \texttt{qgrowth}: Quarterly growth rate: $\frac{ts_t-ts_{t-1}}{ts_{t-1}}$. Works on quarterly, monthly, and weekly series + \begin{itemize} + \item \texttt{ts1.qgrowth()} + \end{itemize} + \item \texttt{yidff}: Yearly difference; works on yearly, quarterly, monthly, and weekly series + \begin{itemize} + \item \texttt{ts1.ydiff()} + \end{itemize} + \item \texttt{ygrowth}: Annual growth rate: $\frac{ts_t-ts_{t-1}}{ts_{t-1}}$. Works on yearly ($t-1$), quarterly ($t-4$), monthly ($t-12$), and weekly ($t-52$) series + \begin{itemize} + \item \texttt{ts1.ygrowth()} + \end{itemize} + \end{itemize} +\end{frame} + + +\begin{frame}[fragile,t] + \frametitle{Operations on \texttt{dseries}} + \begin{itemize} + \item \texttt{abs}: The absolute value of each data point + \item \texttt{cumprod}/\texttt{cumsum}: Cumulative product/sum + \item \texttt{exp}/\texttt{log}: The exponential/log for each data point + \item \texttt{mrdivide}/\texttt{mtimes}: Division/Multiplication + \item \texttt{minus}/\texttt{plus}: Subtraction/Addition + \item \texttt{mpower}: Power + \begin{itemize} + \item \texttt{ts1-ts2} + \item \texttt{ts1-3} + \item \texttt{ts1-[1:3]'} + \end{itemize} + \item \texttt{lag}/\texttt{lead}: lag/lead the series + \begin{itemize} + \item \texttt{ts1(-1)}/\texttt{ts1(2)} + \item \texttt{ts1.lag()}/\texttt{ts1.lead(2)} + \end{itemize} + \item \texttt{uminus}: Equivalent to multiplying by $-1$. + \begin{itemize} + \item \texttt{-ts1} + \end{itemize} + \end{itemize} +\end{frame} + + +\begin{frame}[fragile,t] + \frametitle{Modifying \texttt{dates}} + \begin{itemize} + \item \texttt{align}: Makes both series cover the same time range by padding with \texttt{NaN}'s + \begin{itemize} + \item \texttt{[ts1,ts2]=align(ts1,ts2)} + \end{itemize} + \item \texttt{chain}: + \item \texttt{horzcat}: Join two or more \texttt{dseries} + \begin{itemize} + \item \texttt{ts3=[ts1 ts2]} + \end{itemize} + \item \texttt{insert}: Inserts variables from one \texttt{dseries} into another + \begin{itemize} + \item \texttt{ts1.insert(ts2, 2)} + \end{itemize} + \item \texttt{merge}: Merge two series + \item \texttt{pop}: remove the last variable from \texttt{dseries} + \item \texttt{plot}: plot a \texttt{dseries} + \item \texttt{rename}: Rename a variable in a dseries + \begin{itemize} + \item \texttt{ts1.rename(`MyVar1',`MyFirstVar')} + \end{itemize} + \item \texttt{save}: Save a \texttt{dseries} in either \texttt{.csv} (default), \texttt{.m}, or \texttt{.mat} + \item \texttt{set\_names}: Rename all variables in a \texttt{dseries} + \begin{itemize} + \item \texttt{ts3.set\_names(`NewName1',`NewName2')} + \end{itemize} + \item \texttt{tex\_rename}: Rename the \LaTeX name for a given variable + \begin{itemize} + \item \texttt{ts1.tex\_rename(`MyVar1',`MyVar\textbackslash\_1')} + \end{itemize} + \item \texttt{vertcat}: Add more observations to existing \texttt{dseries} + \begin{itemize} + \item \texttt{[ts1; dseries(4,`2002y',`MyVar1')]} + \end{itemize} + \end{itemize} +\end{frame} + \subsection{Examples} @@ -326,4 +527,7 @@ \subsection{Examples} + +\section{Putting it All Together} + \end{document} From a9e90f5d57b45492e3a275a542ff9151aa509931 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Wed, 14 May 2014 16:54:58 +0200 Subject: [PATCH 31/37] doc: add reporting slides to dseries presentation --- .../dseriesReporting.tex | 108 +++++++++++++++++- 1 file changed, 103 insertions(+), 5 deletions(-) diff --git a/doc/dseries-and-reporting/dseriesReporting.tex b/doc/dseries-and-reporting/dseriesReporting.tex index 694e9fd11..e555dd513 100644 --- a/doc/dseries-and-reporting/dseriesReporting.tex +++ b/doc/dseries-and-reporting/dseriesReporting.tex @@ -424,7 +424,7 @@ \begin{itemize} \item \texttt{ts3.set\_names(`NewName1',`NewName2')} \end{itemize} - \item \texttt{tex\_rename}: Rename the \LaTeX name for a given variable + \item \texttt{tex\_rename}: Rename the \LaTeX\ name for a given variable \begin{itemize} \item \texttt{ts1.tex\_rename(`MyVar1',`MyVar\textbackslash\_1')} \end{itemize} @@ -440,6 +440,8 @@ \subsection{Examples} + + % % REPORTING % @@ -458,19 +460,54 @@ \begin{itemize} \item Can easily be included in another document \end{itemize} - \item Graphs are produced in Ti$k$Z + \item Graphs are produced in Ti$k$Z/PGFPlots (standard in a TeX distribution) \begin{itemize} \item Scales well \item Formating follows that of enclosing document \end{itemize} + \item Dynare provides a subset of the many Ti$k$Z options + \begin{itemize} + \item You can easily modify the Ti$k$Z graph if the option you want is not in Dynare + \end{itemize} \item Works with Matlab \& Octave \item Works approximately 5 times faster than Iris reporting + \item NB: Must install a \LaTeX\ distribution to compile reports + \begin{itemize} + \item Windows: MiKTeX \url{http://miktex.org} + \item Mac OS X: MacTeX \url{http://tug.org/mactex} + \item Linux: TeX Live (from your package manager) + \end{itemize} \end{itemize} \end{frame} + +\begin{frame} + \frametitle{How Reporting Works} + \begin{itemize} + \item Reports are created command by command + \begin{itemize} + \item Hence the order of commands matters + \end{itemize} + \item All reporting commands act on the previously added object until an object of greater or equal hierarchy is added (see next slide) + \begin{itemize} + \item \textit{e.g.,} Once you add a \texttt{Page} to your report with the \texttt{addPage()} command, every \texttt{Section} you add via the \texttt{addSection()} command will be placed on this page. Only when you add another \texttt{Page} will items go on a new page. + \item This will become more clear with an example + \end{itemize} + \item Options to reporting commands are passed in option name/value pairs + \begin{itemize} + \item \textit{e.g.,} \texttt{addPage(`title', \{`Page Title', `Page Subtitle'\})} + \end{itemize} + \end{itemize} +\end{frame} + + \begin{frame} \frametitle{Reporting Class Hierarchy} - \centering { + \begin{itemize} + \item Class names on the top half of the box, constructor names on the bottom + \item Arrows represent what the new object can be added to + \end{itemize} + \begin{center} \begin{tikzpicture}[ node distance = .45cm, auto, @@ -519,15 +556,76 @@ \draw [line] (Section) to node { } (Page); \draw [line] (Page) to node { } (Report); \end{tikzpicture} - } + \end{center} \end{frame} + \subsection{Syntax} +\begin{frame} + \frametitle{Reporting Syntax} + \begin{itemize} + \item \texttt{report(\ldots)}: Create a report + \begin{itemize} + \item Options: \texttt{compiler}, \texttt{showDate}, \texttt{fileName}, \texttt{margin}, \texttt{marginUnit}, \texttt{orientation}, \texttt{paper}, \texttt{title} + \end{itemize} + \item \texttt{addPage(\ldots)}: Add a page to the \texttt{Report} + \begin{itemize} + \item Options: \texttt{footnote}, \texttt{orientation}, \texttt{paper}, \texttt{title}, \texttt{titleFormat} + \end{itemize} + \item \texttt{addSection(\ldots)}: Add a section to the current \texttt{Page} + \begin{itemize} + \item You can think of a section as a matrix. As graphs and/or tables are added section, it fills up from left to right. Once you have added \texttt{cols} objects, a new row is started. + \item Options: \texttt{cols}, \texttt{height} + \end{itemize} + \end{itemize} +\end{frame} + + +\begin{frame} + \frametitle{Reporting Syntax (continued)} + \begin{itemize} + \item \texttt{addGraph(\ldots)}: Add a graph to the current \texttt{Section} + \begin{itemize} + \item Options: \texttt{data}, \texttt{graphDirName}, \texttt{graphName}, \texttt{graphSize}, \texttt{height}, \texttt{showGrid}, \texttt{showLegend}, \texttt{showLegendBox}, \texttt{legendLocation}, \texttt{legendOrientation}, \texttt{legendFontSize}, \texttt{seriesToUse}, \texttt{shade}, \texttt{shadeColor}, \texttt{shadeOpacity}, \texttt{title}, \texttt{titleFormat}, \texttt{width}, \texttt{xlabel}, \texttt{ylabel}, \texttt{xAxisTight}, \texttt{xrange}, \texttt{xTicks}, \texttt{xTickLabels}, \texttt{xTickLabelAnchor}, \texttt{xTickLabelRotation}, \texttt{yAxisTight}, \texttt{yrange}, \texttt{showZeroLine} + \end{itemize} + \item \texttt{addTable(\ldots)}: Add a table to the current \texttt{Section} + \begin{itemize} + \item Options: \texttt{data}, \texttt{showHlines}, \texttt{precision}, \texttt{range}, \texttt{seriesToUse}, \texttt{tableDirName}, \texttt{tableName}, \texttt{title}, \texttt{titleFormat}, \texttt{vlineAfter}, \texttt{vlineAfterEndOfPeriod}, \texttt{showVlines} + \end{itemize} + \item \texttt{addSeries(\ldots)}: Add a series to the current \texttt{Graph} or \texttt{Table} + \begin{itemize} + \item \texttt{data}, \texttt{graphLineStyle}, \texttt{graphLineWidth}, \texttt{graphMarker}, \texttt{graphMarkerEdgeColor}, \texttt{graphMarkerFaceColor}, \texttt{graphMarkerSize}, \texttt{tableDataRhs}, \texttt{tableRowColor}, \texttt{tableRowIndent}, \texttt{tableShowMarkers}, \texttt{tableAlignRight}, \texttt{tableMarkerLimit}, \texttt{tableNegColor}, \texttt{tablePosColor}, \texttt{tableSubSectionHeader}, \texttt{zeroTol} + \end{itemize} + \item Options: \texttt{addVspace(\ldots)}: Add a vertical space to the current \texttt{Section} + \begin{itemize} + \item \texttt{hline}, \texttt{number} + \end{itemize} + \end{itemize} +\end{frame} + + +\begin{frame} + \frametitle{Output} + To create a report: + \begin{itemize} + \item \texttt{write()}: Writes the report to a \LaTeX\ file + \item \texttt{compile(\ldots)}: Compiles the report + \begin{itemize} + \item Options: \texttt{compiler} + \end{itemize} + \end{itemize} + Report Output + \begin{itemize} + \item Unless you pass the \texttt{fileName} option to \texttt{report(\ldots)}, the report will be located in your working directory with the name \texttt{report.tex}. The compiled version will be called \texttt{report.pdf}. + \item Unless you pass the \texttt{graphDirName} or \texttt{graphName} options to \texttt{addGraph(\ldots)}, your graphs will be in a subdirectory of your working directory called \texttt{tmpRepDir}. The default name will take the form \texttt{graph\_pg9\_sec1\_row1\_col5.tex} + \item The same holds for the tables (substituting `table' for `graph' above). + \item Thus you can easily modify these files and include them in another report. + \end{itemize} +\end{frame} \subsection{Examples} - \section{Putting it All Together} \end{document} From d4972e988fbbb19d67eb4bf2396c1821ba2e6dcf Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Fri, 16 May 2014 15:00:43 +0200 Subject: [PATCH 32/37] reporting: fix bug when data option is passed to addGraph --- matlab/reports/@graph/addSeries.m | 22 ---------------------- matlab/reports/@graph/graph.m | 4 ++-- 2 files changed, 2 insertions(+), 24 deletions(-) delete mode 100644 matlab/reports/@graph/addSeries.m diff --git a/matlab/reports/@graph/addSeries.m b/matlab/reports/@graph/addSeries.m deleted file mode 100644 index 6a0f6d23e..000000000 --- a/matlab/reports/@graph/addSeries.m +++ /dev/null @@ -1,22 +0,0 @@ -function o = addSeries(o, varargin) -% function o = addSeries(o, varargin) - -% Copyright (C) 2013-2014 Dynare Team -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see . - -o.series{end+1} = report_series(varargin{:}); -end \ No newline at end of file diff --git a/matlab/reports/@graph/graph.m b/matlab/reports/@graph/graph.m index e79b402f6..5be0ee926 100644 --- a/matlab/reports/@graph/graph.m +++ b/matlab/reports/@graph/graph.m @@ -177,11 +177,11 @@ end if ~isempty(o.data) if isempty(o.seriesToUse) for i=1:o.data.vobs - o = o.addSeries('data', o.data{o.data.name{i}}); + o.series{end+1} = report_series('data', o.data{o.data.name{i}}); end else for i=1:length(o.seriesToUse) - o = o.addSeries('data', o.data{o.seriesToUse{i}}); + o.series{end+1} = report_series('data', o.data{o.seriesToUse{i}}); end end end From 35e56bad7271f31b32366942541fb5516ac53922 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Fri, 16 May 2014 15:29:07 +0200 Subject: [PATCH 33/37] Add back function that was erroneously removed --- matlab/reports/@graph/addSeries.m | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) create mode 100644 matlab/reports/@graph/addSeries.m diff --git a/matlab/reports/@graph/addSeries.m b/matlab/reports/@graph/addSeries.m new file mode 100644 index 000000000..6a0f6d23e --- /dev/null +++ b/matlab/reports/@graph/addSeries.m @@ -0,0 +1,22 @@ +function o = addSeries(o, varargin) +% function o = addSeries(o, varargin) + +% Copyright (C) 2013-2014 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +o.series{end+1} = report_series(varargin{:}); +end \ No newline at end of file From 3516dd59893fd5716015b5f11e3feae0840e8fc7 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Fri, 16 May 2014 16:35:18 +0200 Subject: [PATCH 34/37] reporting: simplify check --- matlab/reports/@report_table/writeTableFile.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/reports/@report_table/writeTableFile.m b/matlab/reports/@report_table/writeTableFile.m index ca6c46aac..8fad2a559 100644 --- a/matlab/reports/@report_table/writeTableFile.m +++ b/matlab/reports/@report_table/writeTableFile.m @@ -33,7 +33,7 @@ function o = writeTableFile(o, pg, sec, row, col) % along with Dynare. If not, see . ne = length(o.series); -if length(o.series) == 0 +if ne == 0 warning('@report_table.write: no series to plot, returning'); return; end From 3af8ab90bb9593191e726249c748822e8437601b Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Fri, 16 May 2014 16:43:45 +0200 Subject: [PATCH 35/37] reporting: bug fix --- matlab/reports/@report_table/writeTableFile.m | 1 + 1 file changed, 1 insertion(+) diff --git a/matlab/reports/@report_table/writeTableFile.m b/matlab/reports/@report_table/writeTableFile.m index 8fad2a559..3aeaaecdc 100644 --- a/matlab/reports/@report_table/writeTableFile.m +++ b/matlab/reports/@report_table/writeTableFile.m @@ -56,6 +56,7 @@ nlhc = 1; if isempty(o.range) dates = getMaxRange(o.series); + o.range = {dates}; else dates = o.range{1}; end From e6b7a5b74a1904fb2002405de019eb912ad9e559 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Fri, 16 May 2014 16:59:56 +0200 Subject: [PATCH 36/37] reporting: fix spacing between sections --- matlab/reports/@section/write.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/reports/@section/write.m b/matlab/reports/@section/write.m index 888e7a9f0..7a9808e8a 100644 --- a/matlab/reports/@section/write.m +++ b/matlab/reports/@section/write.m @@ -71,6 +71,6 @@ for i=1:ne end end end -fprintf(fid, '\\end{tabular}}\n'); +fprintf(fid, '\\end{tabular}}\\\\\n'); fprintf(fid, '%% End Section Object\n\n'); end \ No newline at end of file From e97cf3d1dbacc1adc9b6dc2fe9056b675a3f046e Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Fri, 16 May 2014 17:12:12 +0200 Subject: [PATCH 37/37] doc: add example to dseries/reporting slides --- .../dseriesReporting.tex | 53 +++++++++++++++++++ 1 file changed, 53 insertions(+) diff --git a/doc/dseries-and-reporting/dseriesReporting.tex b/doc/dseries-and-reporting/dseriesReporting.tex index e555dd513..32155c32c 100644 --- a/doc/dseries-and-reporting/dseriesReporting.tex +++ b/doc/dseries-and-reporting/dseriesReporting.tex @@ -627,5 +627,58 @@ \subsection{Examples} \section{Putting it All Together} +\begin{frame}[fragile=singleslide] + \frametitle{Create Report of IRFs from \texttt{example1.mod}} + \begin{itemize} + \item \texttt{example1.mod} is located in the Dynare \texttt{examples} directory + \end{itemize} + \begin{block}{Create \texttt{dseries} from IRFs} +\begin{verbatim} +@#define endovars=["y", "c", "k", "a", "h", "b"] +@#for var in endovars + shocke.@{var} = dseries(@{var}_e, 2014q3, `@{var}'); + shocku.@{var} = dseries(@{var}_u, 2014q3, `@{var}'); +@#endfor +\end{verbatim} + \end{block} +\end{frame} +\begin{frame}[fragile=singleslide] + \frametitle{Create Report of IRFs from \texttt{example1.mod}} + \begin{block}{Populate Report} +\small{ +\begin{verbatim} +@#for shock in ["e", "u"] + report = report.addPage(`title', {`Dseries \& Report Example', ... + `Shock to @{shock}'}, ... + `titleFormat', {`\Large\bfseries', ... + `\large\bfseries'}); + report = report.addSection(`cols', 2); +@# for var in endovars + report = report.addGraph(`data', shock@{shock}.@{var}, ... + `title', `@{var}', ... + `showGrid', false, ... + `showZeroLine', true); +@# endfor + report = report.addSection(`cols', 1); + report = report.addTable(); +@# for var in endovars + report = report.addSeries(`data', shock@{shock}.@{var}); +@# endfor +@#endfor +\end{verbatim} +} + \end{block} +\end{frame} + + +\begin{frame}[fragile=singleslide] + \frametitle{Create Report of IRFs from \texttt{example1.mod}} + \begin{block}{Compile Report} +\begin{verbatim} +report.write(); +report.compile(); +\end{verbatim} + \end{block} +\end{frame} \end{document}