Commit 37abe362 authored by Sébastien Villemot's avatar Sébastien Villemot
Browse files

Preprocessor: added new 'steady_state_model' keyword for automatically creating steady state file

parent ed84e273
......@@ -1075,6 +1075,9 @@ end;
<para>
Note that the model written in the TeX file will differ from the model declared by the user in the some dimensions, see <xref linkend="write_latex_dynamic_model"/>.
</para>
<para>
Also note that this command will not output the contents of the optional <xref linkend="steady_state_model"/> block; it will rather output a static version (<foreignphrase>i.e.</foreignphrase> without leads and lags) of the dynamic model declared in the <xref linkend="model"/> block.
</para>
</refsect1>
</refentry>
......@@ -1660,6 +1663,7 @@ periods 100;
<itemizedlist>
<listitem><para><xref linkend='steady'/></para></listitem>
<listitem><para><xref linkend='homotopy_setup'/></para></listitem>
<listitem><para><xref linkend='steady_state_model'/></para></listitem>
<listitem><para><xref linkend='check'/></para></listitem>
<listitem><para><xref linkend='model_info'/></para></listitem>
<listitem><para><xref linkend='simul'/></para></listitem>
......@@ -1810,6 +1814,85 @@ steady(homotopy_mode = 1, homotopy_steps = 50);
</refentry>
<refentry id="steady_state_model">
<refmeta>
<refentrytitle>steady_state_model</refentrytitle>
</refmeta>
<refnamediv>
<refname>steady_state_model</refname>
<refpurpose>declares the solution of the static model, when analytically known</refpurpose>
</refnamediv>
<refsynopsisdiv>
<cmdsynopsis>
<command>steady_state_model</command>
<arg choice="plain">;</arg>
<sbr/>
<arg choice="plain" rep="repeat">
<replaceable>VARIABLE_NAME</replaceable> = <replaceable>EXPRESSION</replaceable> ;
</arg>
<sbr/>
<command>end</command><arg choice="plain">;</arg>
</cmdsynopsis>
</refsynopsisdiv>
<refsect1><title>Description</title>
<para>
When the analytical solution of the model is known, this command can be used to help Dynare find the steady state in a more efficient and reliable way, especially during estimation where the steady state has to be recomputed for every point in the parameter space.
</para>
<para>
Each line of this block consists of a variable (either an endogenous or a temporary variable) which is assigned an expression (which can contain parameters, exogenous at the steady state, or any endogenous or temporary variable already declared above).
</para>
<para>
Internally, Dynare will create a steady state file called <filename><replaceable>FILENAME</replaceable>_steadystate.m</filename>, using the information provided in this block.
</para>
</refsect1>
<refsect1><title>Example</title>
<informalexample>
<programlisting>
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;
...
// parameter calibration, (dynamic) model declaration, shock calibration...
...
steady_state_model;
dA = exp(gam);
gst = 1/dA; // A temporary variable
m = mst;
// Three other temporary variables
khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1);
nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp );
n = xist/(nust+xist);
P = xist + nust;
k = khst*n;
l = psi*mst*n/( (1-psi)*(1-n) );
c = mst/P;
d = l - mst + 1;
y = k^alp*n^(1-alp)*gst^alp;
R = mst/bet;
W = l/n;
e = 1;
gp_obs = m/dA;
gy_obs = dA;
end;
steady;
</programlisting>
</informalexample>
</refsect1>
</refentry>
<refentry id="check">
<refmeta>
<refentrytitle>check</refentrytitle>
......
......@@ -119,7 +119,7 @@ class ParsingDriver;
%token <string_val> QUOTED_STRING
%token QZ_CRITERIUM FULL
%token RELATIVE_IRF REPLIC RPLOT SAVE_PARAMS_AND_STEADY_STATE
%token SHOCKS SHOCK_DECOMPOSITION SIGMA_E SIMUL SIMUL_ALGO SIMUL_SEED SMOOTHER STACK_SOLVE_ALGO SOLVE_ALGO
%token SHOCKS SHOCK_DECOMPOSITION SIGMA_E SIMUL SIMUL_ALGO SIMUL_SEED SMOOTHER STACK_SOLVE_ALGO STEADY_STATE_MODEL SOLVE_ALGO
%token STDERR STEADY STOCH_SIMUL
%token TEX RAMSEY_POLICY PLANNER_DISCOUNT
%token <string_val> TEX_NAME
......@@ -240,6 +240,7 @@ statement : parameters
| markov_switching
| svar
| external_function
| steady_state_model
;
dsample : DSAMPLE INT_NUMBER ';'
......@@ -1649,6 +1650,18 @@ conditional_forecast_paths_shock_elem : VAR symbol ';' PERIODS period_list ';' V
{ driver.add_det_shock($2, true); }
;
steady_state_model : STEADY_STATE_MODEL ';' { driver.begin_steady_state_model(); }
steady_state_equation_list END { driver.reset_data_tree(); }
;
steady_state_equation_list : steady_state_equation_list steady_state_equation
| steady_state_equation
;
steady_state_equation : symbol EQUAL expression ';'
{ driver.add_steady_state_model_equal($1, $3); }
;
o_dr_algo : DR_ALGO EQUAL INT_NUMBER {
if (*$3 == string("0"))
driver.warning("dr_algo option is now deprecated, and may be removed in a future version of Dynare");
......
......@@ -154,6 +154,7 @@ int sigma_e = 0;
/* Begin of a Dynare block */
<INITIAL>model {BEGIN DYNARE_BLOCK; return token::MODEL;}
<INITIAL>steady_state_model {BEGIN DYNARE_BLOCK; return token::STEADY_STATE_MODEL;}
<INITIAL>initval {BEGIN DYNARE_BLOCK; return token::INITVAL;}
<INITIAL>endval {BEGIN DYNARE_BLOCK; return token::ENDVAL;}
<INITIAL>histval {BEGIN DYNARE_BLOCK; return token::HISTVAL;}
......
......@@ -523,14 +523,13 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
switch (type)
{
case eParameter:
if (output_type == oMatlabOutsideModel)
if (output_type == oMatlabOutsideModel || output_type == oSteadyStateFile)
output << "M_.params" << "(" << tsid + 1 << ")";
else
output << "params" << LEFT_ARRAY_SUBSCRIPT(output_type) << tsid + ARRAY_SUBSCRIPT_OFFSET(output_type) << RIGHT_ARRAY_SUBSCRIPT(output_type);
break;
case eModelLocalVariable:
case eModFileLocalVariable:
if (output_type == oMatlabDynamicModelSparse || output_type == oMatlabStaticModelSparse || output_type == oMatlabDynamicModelSparseLocalTemporaryTerms)
{
output << "(";
......@@ -541,6 +540,10 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
output << datatree.symbol_table.getName(symb_id);
break;
case eModFileLocalVariable:
output << datatree.symbol_table.getName(symb_id);
break;
case eEndogenous:
switch (output_type)
{
......@@ -570,6 +573,9 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
case oMatlabDynamicSteadyStateOperator:
output << "oo_.steady_state(" << tsid + 1 << ")";
break;
case oSteadyStateFile:
output << "ys_(" << tsid + 1 << ")";
break;
default:
assert(false);
}
......@@ -608,6 +614,9 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
case oMatlabDynamicSteadyStateOperator:
output << "oo_.exo_steady_state(" << i << ")";
break;
case oSteadyStateFile:
output << "exo_(" << i << ")";
break;
default:
assert(false);
}
......@@ -646,6 +655,9 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
case oMatlabDynamicSteadyStateOperator:
output << "oo_.exo_det_steady_state(" << tsid + 1 << ")";
break;
case oSteadyStateFile:
output << "exo_(" << i << ")";
break;
default:
assert(false);
}
......
......@@ -71,7 +71,8 @@ enum ExprNodeOutputType
oLatexDynamicSteadyStateOperator, //!< LaTeX code, dynamic model steady state declarations
oMatlabDynamicSteadyStateOperator, //!< Matlab code, dynamic model steady state declarations
oMatlabDynamicModelSparseSteadyStateOperator, //!< Matlab code, dynamic block decomposed model steady state declarations
oMatlabDynamicModelSparseLocalTemporaryTerms //!< Matlab code, dynamic block decomposed model local temporary_terms
oMatlabDynamicModelSparseLocalTemporaryTerms, //!< Matlab code, dynamic block decomposed model local temporary_terms
oSteadyStateFile //!< Matlab code, in the generated steady state file
};
#define IS_MATLAB(output_type) ((output_type) == oMatlabStaticModel \
......@@ -81,7 +82,8 @@ enum ExprNodeOutputType
|| (output_type) == oMatlabDynamicModelSparse \
|| (output_type) == oMatlabDynamicModelSparseLocalTemporaryTerms \
|| (output_type) == oMatlabDynamicSteadyStateOperator \
|| (output_type) == oMatlabDynamicModelSparseSteadyStateOperator)
|| (output_type) == oMatlabDynamicModelSparseSteadyStateOperator \
|| (output_type) == oSteadyStateFile)
#define IS_C(output_type) ((output_type) == oCDynamicModel)
......
......@@ -47,7 +47,9 @@ dynare_m_SOURCES = \
CodeInterpreter.hh \
FlexLexer.h \
ExternalFunctionsTable.cc \
ExternalFunctionsTable.hh
ExternalFunctionsTable.hh \
SteadyStateModel.hh \
SteadyStateModel.cc
# The -I. is for <FlexLexer.h>
dynare_m_CPPFLAGS = $(BOOST_CPPFLAGS) -I.
......
......@@ -24,8 +24,9 @@
#include "ModFile.hh"
ModFile::ModFile() : expressions_tree(symbol_table, num_constants, external_functions_table),
static_model(symbol_table, num_constants, external_functions_table),
steady_state_model(symbol_table, num_constants, external_functions_table),
dynamic_model(symbol_table, num_constants, external_functions_table),
static_model(symbol_table, num_constants, external_functions_table),
linear(false), block(false), byte_code(false),
use_dll(false), no_static(false)
{
......@@ -439,5 +440,8 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all
}
}
// Create steady state file
steady_state_model.writeSteadyStateFile(basename);
cout << "done" << endl;
}
......@@ -30,6 +30,7 @@ using namespace std;
#include "NumericalInitialization.hh"
#include "StaticModel.hh"
#include "DynamicModel.hh"
#include "SteadyStateModel.hh"
#include "Statement.hh"
#include "ExternalFunctionsTable.hh"
......@@ -47,10 +48,12 @@ public:
NumericalConstants num_constants;
//! Expressions outside model block
DataTree expressions_tree;
//! Static Dll model
StaticModel static_model;
//! Dynamic model
//! Static model, as declared in the "steady_state_model" block if present
SteadyStateModel steady_state_model;
//! Dynamic model, as declared in the "model" block
DynamicModel dynamic_model;
//! Static model, as derived from the "model" block when leads and lags have been removed
StaticModel static_model;
//! Option linear
bool linear;
......
......@@ -1391,17 +1391,27 @@ ParsingDriver::add_model_equal_with_zero_rhs(NodeID arg)
void
ParsingDriver::declare_and_init_model_local_variable(string *name, NodeID rhs)
{
int symb_id;
try
{
mod_file->symbol_table.addSymbol(*name, eModelLocalVariable);
symb_id = mod_file->symbol_table.addSymbol(*name, eModelLocalVariable);
}
catch (SymbolTable::AlreadyDeclaredException &e)
{
error("Local model variable " + *name + " declared twice.");
// It can have already been declared in a steady_state_model block, check that it is indeed a ModelLocalVariable
symb_id = mod_file->symbol_table.getID(*name);
if (mod_file->symbol_table.getType(symb_id) != eModelLocalVariable)
error(*name + " has wrong type, you cannot use it within as left-hand side of a pound ('#') expression");
}
int symb_id = mod_file->symbol_table.getID(*name);
model_tree->AddLocalVariable(symb_id, rhs);
try
{
model_tree->AddLocalVariable(symb_id, rhs);
}
catch (DataTree::LocalVariableException &e)
{
error("Local model variable " + *name + " declared twice.");
}
delete name;
}
......@@ -1834,3 +1844,43 @@ ParsingDriver::add_native(const char *s)
{
mod_file->addStatement(new NativeStatement(s));
}
void
ParsingDriver::begin_steady_state_model()
{
set_current_data_tree(&mod_file->steady_state_model);
}
void
ParsingDriver::add_steady_state_model_equal(string *varname, NodeID expr)
{
int id;
try
{
id = mod_file->symbol_table.getID(*varname);
}
catch (SymbolTable::UnknownSymbolNameException &e)
{
// Unknown symbol, declare it as a ModFileLocalVariable
id = mod_file->symbol_table.addSymbol(*varname, eModFileLocalVariable);
}
SymbolType type = mod_file->symbol_table.getType(id);
if (type != eEndogenous && type != eModFileLocalVariable)
error(*varname + " has incorrect type");
try
{
mod_file->steady_state_model.addDefinition(id, expr);
}
catch(SteadyStateModel::AlreadyDefinedException &e)
{
error(*varname + " has already been defined in the steady state block");
}
catch(SteadyStateModel::UndefinedVariableException &e)
{
error(e.varname + " is not yet initialized at this point");
}
delete varname;
}
......@@ -493,6 +493,10 @@ public:
void add_native(const char *s);
//! Resets data_tree and model_tree pointers to default (i.e. mod_file->expressions_tree)
void reset_data_tree();
//! Begin a steady_state_model block
void begin_steady_state_model();
//! Add an assignment equation in steady_state_model block
void add_steady_state_model_equal(string *varname, NodeID expr);
};
#endif // ! PARSING_DRIVER_HH
......@@ -26,7 +26,7 @@ using namespace std;
#include "ModelTree.hh"
//! Stores a static model
//! Stores a static model, as derived from the "model" block when leads and lags have been removed
class StaticModel : public ModelTree
{
private:
......
/*
* Copyright (C) 2010 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 <http://www.gnu.org/licenses/>.
*/
#include <cassert>
#include <algorithm>
#include "SteadyStateModel.hh"
SteadyStateModel::SteadyStateModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants, ExternalFunctionsTable &external_functions_table_arg) :
DataTree(symbol_table_arg, num_constants, external_functions_table)
{
}
void
SteadyStateModel::addDefinition(int symb_id, NodeID expr) throw (UndefinedVariableException, AlreadyDefinedException)
{
assert(symbol_table.getType(symb_id) == eEndogenous
|| symbol_table.getType(symb_id) == eModFileLocalVariable);
// Check that symbol is not already defined
if (find(recursive_order.begin(), recursive_order.end(), symb_id)
!= recursive_order.end())
throw AlreadyDefinedException(symbol_table.getName(symb_id));
// Check that expression has no undefined symbol
set<pair<int, int> > used_symbols;
expr->collectVariables(eEndogenous, used_symbols);
expr->collectVariables(eModFileLocalVariable, used_symbols);
for(set<pair<int, int> >::const_iterator it = used_symbols.begin();
it != used_symbols.end(); it++)
if (find(recursive_order.begin(), recursive_order.end(), it->first)
== recursive_order.end())
throw UndefinedVariableException(symbol_table.getName(it->first));
// Add the variable
recursive_order.push_back(symb_id);
def_table[symb_id] = AddEqual(AddVariable(symb_id), expr);
}
void
SteadyStateModel::writeSteadyStateFile(const string &basename) const
{
if (recursive_order.size() == 0)
return;
string filename = basename + "_steadystate.m";
ofstream output;
output.open(filename.c_str(), ios::out | ios::binary);
if (!output.is_open())
{
cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
exit(EXIT_FAILURE);
}
output << "function [ys_, check_] = " << basename << "_steadystate(ys_orig_, exo_)" << endl
<< "% Steady state generated by Dynare preprocessor" << endl
<< " global M_" << endl
<< " ys_=zeros(" << symbol_table.orig_endo_nbr() << ",1);" << endl;
for(size_t i = 0; i < recursive_order.size(); i++)
{
output << " ";
map<int, NodeID>::const_iterator it = def_table.find(recursive_order[i]);
it->second->writeOutput(output, oSteadyStateFile);
output << ";" << endl;
}
output << " check_=0;" << endl
<< "end" << endl;
}
/*
* Copyright (C) 2010 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 <http://www.gnu.org/licenses/>.
*/
#ifndef _STEADY_STATE_MODEL_HH
#define _STEADY_STATE_MODEL_HH
#include "DataTree.hh"
class SteadyStateModel : public DataTree
{
private:
//! Associates a symbol ID to an expression of the form "var = expr"
map<int, NodeID> def_table;
vector<int> recursive_order;
public:
class AlreadyDefinedException
{
public:
const string &varname;
AlreadyDefinedException(const string &varname_arg) : varname(varname_arg) {}
};
class UndefinedVariableException
{
public:
const string &varname;
UndefinedVariableException(const string &varname_arg) : varname(varname_arg) {}
};
SteadyStateModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants, ExternalFunctionsTable &external_functions_table_arg);
//! Add an expression of the form "var = expr;"
void addDefinition(int symb_id, NodeID expr) throw (UndefinedVariableException, AlreadyDefinedException);
//! Write the steady state file
void writeSteadyStateFile(const string &basename) const;
};
#endif
......@@ -15,6 +15,7 @@ OCTAVE_MODS = \
ramst_steady_state.mod \
example1_varexo_det.mod \
predetermined_variables.mod \
fs2000_ssfile.mod \
block_bytecode/fs2000_simk.mod \
block_bytecode/fs2000_lu.mod \
block_bytecode/fs2000_bicgstab.mod \
......@@ -123,4 +124,6 @@ clean-local:
rm -rf ramsey_objective
rm -f fs2000_ssfile_steadystate.m
rm -f $(shell find -name '*~')
// Test file for "steady_state_model" block.
// Note that this file has leads of 2, so this is also a test for dealing with auxiliary vars.
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));
gy_obs = dA*y/y(-1);
gp_obs = (P/P(-1))*m(-1)/dA;
end;
shocks;
var e_a; stderr 0.014;
var e_m; stderr 0.005;
end;
steady_state_model;
dA = exp(gam);
gst = 1/dA;
m = mst;
khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1);
nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp );
n = xist/(nust+xist);
P = xist + nust;
k = khst*n;
l = psi*mst*n/( (1-psi)*(1-n) );
c = mst/P;
d = l - mst + 1;
y = k^alp*n^(1-alp)*gst^alp;
R = mst/bet;
W = l/n;
e = 1;
gp_obs = m/dA;
gy_obs = dA;
end;
steady;
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment