DynamicModel.hh 20.3 KB
Newer Older
sebastien's avatar
sebastien committed
1
/*
Sébastien Villemot's avatar
Sébastien Villemot committed
2
 * Copyright (C) 2003-2012 Dynare Team
sebastien's avatar
sebastien committed
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
 *
 * 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 _DYNAMICMODEL_HH
#define _DYNAMICMODEL_HH

using namespace std;
24
#define ZERO_BAND 1e-8
sebastien's avatar
sebastien committed
25

26
27
#include <fstream>

sebastien's avatar
sebastien committed
28
29
30
31
32
33
#include "StaticModel.hh"

//! Stores a dynamic model
class DynamicModel : public ModelTree
{
private:
34
35
36
37
  //! Stores equations declared as [static]
  /*! They will be used in toStatic() to replace equations marked as [dynamic] */
  vector<BinaryOpNode *> static_only_equations;

38
39
40
41
42
43
44
45
46
47
48
  typedef map<pair<int, int>, int> deriv_id_table_t;
  //! Maps a pair (symbol_id, lag) to a deriv ID
  deriv_id_table_t deriv_id_table;
  //! Maps a deriv ID to a pair (symbol_id, lag)
  vector<pair<int, int> > inv_deriv_id_table;

  //! Maps a deriv_id to the column index of the dynamic Jacobian
  /*! Contains only endogenous, exogenous and exogenous deterministic */
  map<int, int> dyn_jacobian_cols_table;

  //! Maximum lag and lead over all types of variables (positive values)
sebastien's avatar
sebastien committed
49
  /*! Set by computeDerivIDs() */
50
51
  int max_lag, max_lead;
  //! Maximum lag and lead over endogenous variables (positive values)
sebastien's avatar
sebastien committed
52
  /*! Set by computeDerivIDs() */
53
54
  int max_endo_lag, max_endo_lead;
  //! Maximum lag and lead over exogenous variables (positive values)
sebastien's avatar
sebastien committed
55
  /*! Set by computeDerivIDs() */
56
57
  int max_exo_lag, max_exo_lead;
  //! Maximum lag and lead over deterministic exogenous variables (positive values)
sebastien's avatar
sebastien committed
58
  /*! Set by computeDerivIDs() */
59
60
  int max_exo_det_lag, max_exo_det_lead;

61
  //! Number of columns of dynamic jacobian
sebastien's avatar
sebastien committed
62
  /*! Set by computeDerivID()s and computeDynJacobianCols() */
63
  int dynJacobianColsNbr;
64
  //! Temporary terms for block decomposed models
65
  vector< vector<temporary_terms_t> > v_temporary_terms;
66

67
  vector<temporary_terms_inuse_t> v_temporary_terms_inuse;
68

69
70
  //! Store the derivatives or the chainrule derivatives:map<pair< equation, pair< variable, lead_lag >, expr_t>
  typedef map< pair< int, pair< int, int> >, expr_t> first_chain_rule_derivatives_t;
71
  first_chain_rule_derivatives_t first_chain_rule_derivatives;
72

sebastien's avatar
sebastien committed
73
74
75
76
  //! Writes dynamic model file (Matlab version)
  void writeDynamicMFile(const string &dynamic_basename) const;
  //! Writes dynamic model file (C version)
  /*! \todo add third derivatives handling */
77
  void writeDynamicCFile(const string &dynamic_basename, const int order) const;
sebastien's avatar
sebastien committed
78
  //! Writes dynamic model file when SparseDLL option is on
79
  void writeSparseDynamicMFile(const string &dynamic_basename, const string &basename) const;
sebastien's avatar
sebastien committed
80
81
  //! Writes the dynamic model equations and its derivatives
  /*! \todo add third derivatives handling in C output */
82
  void writeDynamicModel(ostream &DynamicOutput, bool use_dll) const;
sebastien's avatar
sebastien committed
83
  //! Writes the Block reordred structure of the model in M output
84
  void writeModelEquationsOrdered_M(const string &dynamic_basename) const;
sebastien's avatar
sebastien committed
85
  //! Writes the code of the Block reordred structure of the model in virtual machine bytecode
86
  void writeModelEquationsCode_Block(string &file_name, const string &bin_basename, const map_idx_t &map_idx) const;
87
  //! Writes the code of the model in virtual machine bytecode
88
  void writeModelEquationsCode(string &file_name, const string &bin_basename, const map_idx_t &map_idx) const;
89

sebastien's avatar
sebastien committed
90
91
92
93
94
  //! Computes jacobian and prepares for equation normalization
  /*! Using values from initval/endval blocks and parameter initializations:
    - computes the jacobian for the model w.r. to contemporaneous variables
    - removes edges of the incidence matrix when derivative w.r. to the corresponding variable is too close to zero (below the cutoff)
  */
95
  //void evaluateJacobian(const eval_context_t &eval_context, jacob_map *j_m, bool dynamic);
96
97
98
99

  //! return a map on the block jacobian
  map<pair<pair<int, pair<int, int> >, pair<int, int> >, int> get_Derivatives(int block);
  //! Computes chain rule derivatives of the Jacobian w.r. to endogenous variables
100
  void computeChainRuleJacobian(blocks_derivatives_t &blocks_derivatives);
101

sebastien's avatar
sebastien committed
102
  string reform(string name) const;
103
  map_idx_t map_idx;
sebastien's avatar
sebastien committed
104

105
  //! sorts the temporary terms in the blocks order
106
  void computeTemporaryTermsOrdered();
107
108
109

  //! creates a mapping from the index of temporary terms to a natural index
  void computeTemporaryTermsMapping();
sebastien's avatar
sebastien committed
110
  //! Write derivative code of an equation w.r. to a variable
111
  void compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, int lag, const map_idx_t &map_idx) const;
112
  //! Write chain rule derivative code of an equation w.r. to a variable
113
  void compileChainRuleDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int var, int lag, const map_idx_t &map_idx) const;
sebastien's avatar
sebastien committed
114

115
  //! Get the type corresponding to a derivation ID
116
  virtual SymbolType getTypeByDerivID(int deriv_id) const throw (UnknownDerivIDException);
117
  //! Get the lag corresponding to a derivation ID
118
  virtual int getLagByDerivID(int deriv_id) const throw (UnknownDerivIDException);
119
  //! Get the symbol ID corresponding to a derivation ID
120
  virtual int getSymbIDByDerivID(int deriv_id) const throw (UnknownDerivIDException);
121
  //! Compute the column indices of the dynamic Jacobian
122
  void computeDynJacobianCols(bool jacobianExo);
123
124
  //! Computes derivatives of the Jacobian w.r. to trend vars and tests that they are equal to zero
  void testTrendDerivativesEqualToZero(const eval_context_t &eval_context);
125
  //! Collect only the first derivatives
126
  map<pair<int, pair<int, int> >, expr_t> collect_first_order_derivatives_endogenous();
127

sebastien's avatar
sebastien committed
128
129
130
131
  //! Allocates the derivation IDs for all dynamic variables of the model
  /*! Also computes max_{endo,exo}_{lead_lag}, and initializes dynJacobianColsNbr to the number of dynamic endos */
  void computeDerivIDs();

132
  //! Write chain rule derivative of a recursive equation w.r. to a variable
133
  void writeChainRuleDerivative(ostream &output, int eq, int var, int lag, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms) const;
134

135
136
137
  //! Collecte the derivatives w.r. to endogenous of the block, to endogenous of previouys blocks and to exogenous
  void collect_block_first_order_derivatives();

138
139
140
  //! Collecte the informations about exogenous, deterministic exogenous and endogenous from the previous block for each block
  void collectBlockVariables();

141
  //! Factorized code for substitutions of leads/lags
sebastien's avatar
sebastien committed
142
  /*! \param[in] type determines which type of variables is concerned */
143
  void substituteLeadLagInternal(aux_var_t type, bool deterministic_model);
144

145
146
147
148
private:
  //! Indicate if the temporary terms are computed for the overall model (true) or not (false). Default value true
  bool global_temporary_terms;

149
  //! Vector describing equations: BlockSimulationType, if BlockSimulationType == EVALUATE_s then a expr_t on the new normalized equation
150
  equation_type_and_normalized_equation_t equation_type_and_normalized_equation;
151
152

  //! for each block contains pair< Simulation_Type, pair < Block_Size, Recursive_part_Size > >
153
  block_type_firstequation_size_mfs_t block_type_firstequation_size_mfs;
154
155

  //! for all blocks derivatives description
156
  blocks_derivatives_t blocks_derivatives;
157
158

  //! The jacobian without the elements below the cutoff
159
  dynamic_jacob_map_t dynamic_jacobian;
160
161
162
163

  //! Vector indicating if the block is linear in endogenous variable (true) or not (false)
  vector<bool> blocks_linear;

164
165
  //! Map the derivatives for a block pair<lag, make_pair(make_pair(eq, var)), expr_t>
  typedef map<pair< int, pair<int, int> >, expr_t> derivative_t;
166
  //! Vector of derivative for each blocks
167
  vector<derivative_t> derivative_endo, derivative_other_endo, derivative_exo, derivative_exo_det;
168

169
  //!List for each block and for each lag-lead all the other endogenous variables and exogenous variables
170
171
172
  typedef set<int> var_t;
  typedef map<int, var_t> lag_var_t;
  vector<lag_var_t> other_endo_block, exo_block, exo_det_block;
173

174
175
176
  //!List for each block the exogenous variables
  vector<pair<var_t, int> > block_var_exo;

177
  map< int, map<int, int> > block_exo_index, block_det_exo_index, block_other_endo_index;
178
179
180

  //! for each block described the number of static, forward, backward and mixed variables in the block
  /*! pair< pair<static, forward>, pair<backward,mixed> > */
181
  vector<pair< pair<int, int>, pair<int, int> > > block_col_type;
182
183
184
185
186
187

  //! List for each variable its block number and its maximum lag and lead inside the block
  vector<pair<int, pair<int, int> > > variable_block_lead_lag;
  //! List for each equation its block number
  vector<int> equation_block;

188
189
190
  //!Maximum lead and lag for each block on endogenous of the block, endogenous of the previous blocks, exogenous and deterministic exogenous
  vector<pair<int, int> > endo_max_leadlag_block, other_endo_max_leadlag_block, exo_max_leadlag_block, exo_det_max_leadlag_block, max_leadlag_block;

sebastien's avatar
sebastien committed
191
public:
192
  DynamicModel(SymbolTable &symbol_table_arg, NumericalConstants &num_constants_arg, ExternalFunctionsTable &external_functions_table_argx);
sebastien's avatar
sebastien committed
193
194
  //! Adds a variable node
  /*! This implementation allows for non-zero lag */
sebastien's avatar
sebastien committed
195
  virtual VariableNode *AddVariable(int symb_id, int lag = 0);
196
  
sebastien's avatar
sebastien committed
197
  //! Execute computations (variable sorting + derivation)
198
199
200
201
  /*!
    \param jacobianExo whether derivatives w.r. to exo and exo_det should be in the Jacobian (derivatives w.r. to endo are always computed)
    \param hessian whether 2nd derivatives w.r. to exo, exo_det and endo should be computed (implies jacobianExo = true)
    \param thirdDerivatives whether 3rd derivatives w.r. to endo/exo/exo_det should be computed (implies jacobianExo = true)
sebastien's avatar
sebastien committed
202
    \param paramsDerivatives whether 2nd derivatives w.r. to a pair (endo/exo/exo_det, parameter) should be computed (implies jacobianExo = true)
203
204
205
    \param eval_context evaluation context for normalization
    \param no_tmp_terms if true, no temporary terms will be computed in the dynamic files
  */
sebastien's avatar
sebastien committed
206
  void computingPass(bool jacobianExo, bool hessian, bool thirdDerivatives, bool paramsDerivatives,
207
                     const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll, bool bytecode);
sebastien's avatar
sebastien committed
208
  //! Writes model initialization and lead/lag incidence matrix to output
209
  void writeOutput(ostream &output, const string &basename, bool block, bool byte_code, bool use_dll, int order, bool estimation_present) const;
210

sebastien's avatar
sebastien committed
211
  //! Adds informations for simulation in a binary file
212
  void Write_Inf_To_Bin_File_Block(const string &dynamic_basename, const string &bin_basename,
213
                                   const int &num, int &u_count_int, bool &file_open, bool is_two_boundaries) const;
sebastien's avatar
sebastien committed
214
  //! Writes dynamic model file
215
  void writeDynamicFile(const string &basename, bool block, bool bytecode, bool use_dll, int order) const;
sebastien's avatar
sebastien committed
216
217
  //! Writes file containing parameters derivatives
  void writeParamsDerivativesFile(const string &basename) const;
sebastien's avatar
sebastien committed
218
219
220
  //! 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;
221

222
223
224
225
  //! Copies a dynamic model (only the equations)
  /*! It assumes that the dynamic model given in argument has just been allocated */
  void cloneDynamic(DynamicModel &dynamic_model) const;

226
  //! Replaces model equations with derivatives of Lagrangian w.r.t. endogenous
227
  void computeRamseyPolicyFOCs(const StaticModel &static_model);
Houtan Bastani's avatar
Houtan Bastani committed
228
  //! Replaces the model equations in dynamic_model with those in this model
229
230
  void replaceMyEquations(DynamicModel &dynamic_model) const;

231
232
233
234
235
236
237
238
239
  //! Adds an equation marked as [static]
  void addStaticOnlyEquation(expr_t eq);

  //! Returns number of static only equations
  size_t staticOnlyEquationsNbr() const;
  
  //! Returns number of dynamic only equations
  size_t dynamicOnlyEquationsNbr() const;

240
241
242
  //! Writes LaTeX file with the equations of the dynamic model
  void writeLatexFile(const string &basename) const;

243
244
  virtual int getDerivID(int symb_id, int lag) const throw (UnknownDerivIDException);
  virtual int getDynJacobianCol(int deriv_id) const throw (UnknownDerivIDException);
245
  virtual void addAllParamDerivId(set<int> &deriv_id_set);
246

247
  //! Returns true indicating that this is a dynamic model
248
249
250
251
252
  virtual bool
  isDynamic() const
  {
    return true;
  };
sebastien's avatar
sebastien committed
253

254
255
256
  //! Drive test of detrended equations
  void runTrendTest(const eval_context_t &eval_context);

257
  //! Transforms the model by removing all leads greater or equal than 2 on endos
sebastien's avatar
sebastien committed
258
  /*! Note that this can create new lags on endos and exos */
259
  void substituteEndoLeadGreaterThanTwo(bool deterministic_model);
sebastien's avatar
sebastien committed
260

261
  //! Transforms the model by removing all lags greater or equal than 2 on endos
262
  void substituteEndoLagGreaterThanTwo(bool deterministic_model);
263

sebastien's avatar
sebastien committed
264
265
  //! Transforms the model by removing all leads on exos
  /*! Note that this can create new lags on endos and exos */
266
  void substituteExoLead(bool deterministic_model);
sebastien's avatar
sebastien committed
267
268

  //! Transforms the model by removing all lags on exos
269
  void substituteExoLag(bool deterministic_model);
sebastien's avatar
sebastien committed
270

271
272
273
  //! Transforms the model by removing all oExpectation
  void substituteExpectation(bool partial_information_model);

274
275
276
  //! Transforms the model by decreasing the lead/lag of predetermined variables in model equations by one
  void transformPredeterminedVariables();

277
278
279
280
281
282
  //! Transforms the model by removing trends specified by the user
  void detrendEquations();

  //! Transforms the model by replacing trend variables with a 1
  void removeTrendVariableFromEquations();

sebastien's avatar
sebastien committed
283
  //! Fills eval context with values of model local variables and auxiliary variables
284
  void fillEvalContext(eval_context_t &eval_context) const;
285
286

  //! Return the number of blocks
287
288
289
290
291
  virtual unsigned int
  getNbBlocks() const
  {
    return (block_type_firstequation_size_mfs.size());
  };
292
  //! Determine the simulation type of each block
293
294
295
296
297
  virtual BlockSimulationType
  getBlockSimulationType(int block_number) const
  {
    return (block_type_firstequation_size_mfs[block_number].first.first);
  };
298
  //! Return the first equation number of a block
299
300
301
302
303
  virtual unsigned int
  getBlockFirstEquation(int block_number) const
  {
    return (block_type_firstequation_size_mfs[block_number].first.second);
  };
304
  //! Return the size of the block block_number
305
306
307
308
309
  virtual unsigned int
  getBlockSize(int block_number) const
  {
    return (block_type_firstequation_size_mfs[block_number].second.first);
  };
310
311
312
313
314
315
316
317
318
319
320
321
  //! Return the number of exogenous variable in the block block_number
  virtual unsigned int
  getBlockExoSize(int block_number) const
  {
    return (block_var_exo[block_number].first.size());
  };
  //! Return the number of colums in the jacobian matrix for exogenous variable in the block block_number
  virtual unsigned int
  getBlockExoColSize(int block_number) const
  {
    return (block_var_exo[block_number].second);
  };
322
  //! Return the number of feedback variable of the block block_number
323
324
325
326
327
  virtual unsigned int
  getBlockMfs(int block_number) const
  {
    return (block_type_firstequation_size_mfs[block_number].second.second);
  };
328
  //! Return the maximum lag in a block
329
330
331
332
333
  virtual unsigned int
  getBlockMaxLag(int block_number) const
  {
    return (block_lag_lead[block_number].first);
  };
334
  //! Return the maximum lead in a block
335
336
337
338
339
  virtual unsigned int
  getBlockMaxLead(int block_number) const
  {
    return (block_lag_lead[block_number].second);
  };
340
  //! Return the type of equation (equation_number) belonging to the block block_number
341
342
343
344
345
  virtual EquationType
  getBlockEquationType(int block_number, int equation_number) const
  {
    return (equation_type_and_normalized_equation[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]].first);
  };
346
  //! Return true if the equation has been normalized
347
348
349
350
351
  virtual bool
  isBlockEquationRenormalized(int block_number, int equation_number) const
  {
    return (equation_type_and_normalized_equation[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]].first == E_EVALUATE_S);
  };
352
353
354
  //! Return the expr_t of the equation equation_number belonging to the block block_number
  virtual expr_t
  getBlockEquationExpr(int block_number, int equation_number) const
355
356
357
  {
    return (equations[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]]);
  };
358
359
360
  //! Return the expr_t of the renormalized equation equation_number belonging to the block block_number
  virtual expr_t
  getBlockEquationRenormalizedExpr(int block_number, int equation_number) const
361
362
363
  {
    return (equation_type_and_normalized_equation[equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]].second);
  };
364
  //! Return the original number of equation equation_number belonging to the block block_number
365
366
367
368
369
  virtual int
  getBlockEquationID(int block_number, int equation_number) const
  {
    return (equation_reordered[block_type_firstequation_size_mfs[block_number].first.second+equation_number]);
  };
370
  //! Return the original number of variable variable_number belonging to the block block_number
371
372
373
374
375
  virtual int
  getBlockVariableID(int block_number, int variable_number) const
  {
    return (variable_reordered[block_type_firstequation_size_mfs[block_number].first.second+variable_number]);
  };
376
377
378
379
380
381
382
  //! Return the original number of the exogenous variable varexo_number belonging to the block block_number
  virtual int
  getBlockVariableExoID(int block_number, int variable_number) const
  {
    map<int, var_t>::const_iterator it = exo_block[block_number].find(variable_number);
    return (it->first);
  };
383
  //! Return the position of equation_number in the block number belonging to the block block_number
384
385
386
387
388
  virtual int
  getBlockInitialEquationID(int block_number, int equation_number) const
  {
    return ((int) inv_equation_reordered[equation_number] - (int) block_type_firstequation_size_mfs[block_number].first.second);
  };
389
  //! Return the position of variable_number in the block number belonging to the block block_number
390
391
392
393
394
  virtual int
  getBlockInitialVariableID(int block_number, int variable_number) const
  {
    return ((int) inv_variable_reordered[variable_number] - (int) block_type_firstequation_size_mfs[block_number].first.second);
  };
395
396
397
398
  //! Return the block number containing the endogenous variable variable_number
  int
  getBlockVariableID(int variable_number) const
  {
399
    return (variable_block_lead_lag[variable_number].first);
400
401
402
403
404
405
406
407
408
  };
  //! Return the position of the exogenous variable_number in the block number belonging to the block block_number
  virtual int
  getBlockInitialExogenousID(int block_number, int variable_number) const
  {
    map< int, map<int, int> >::const_iterator it = block_exo_index.find(block_number);
    if (it != block_exo_index.end())
      {
        map<int, int>::const_iterator it1 = it->second.find(variable_number);
409
        if (it1 != it->second.end())
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
          return it1->second;
        else
          return -1;
      }
    else
      return (-1);
  };
  //! Return the position of the deterministic exogenous variable_number in the block number belonging to the block block_number
  virtual int
  getBlockInitialDetExogenousID(int block_number, int variable_number) const
  {
    map< int, map<int, int> >::const_iterator it = block_det_exo_index.find(block_number);
    if (it != block_det_exo_index.end())
      {
        map<int, int>::const_iterator it1 = it->second.find(variable_number);
425
        if (it1 != it->second.end())
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
          return it1->second;
        else
          return -1;
      }
    else
      return (-1);
  };
  //! Return the position of the other endogenous variable_number in the block number belonging to the block block_number
  virtual int
  getBlockInitialOtherEndogenousID(int block_number, int variable_number) const
  {
    map< int, map<int, int> >::const_iterator it = block_other_endo_index.find(block_number);
    if (it != block_other_endo_index.end())
      {
        map<int, int>::const_iterator it1 = it->second.find(variable_number);
441
        if (it1 != it->second.end())
442
443
444
445
446
447
448
          return it1->second;
        else
          return -1;
      }
    else
      return (-1);
  };
449
  bool isModelLocalVariableUsed() const;
sebastien's avatar
sebastien committed
450
451
452
};

#endif