ModFile.cc 61.3 KB
Newer Older
1
/*
2
 * Copyright (C) 2006-2018 Dynare Team
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
 *
 * 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/>.
 */

20
#include <cstdlib>
21
22
#include <iostream>
#include <fstream>
23
#include <typeinfo>
24
#include <cassert>
25
26

#include <boost/filesystem.hpp>
27

28
#include "ModFile.hh"
29
#include "ConfigFile.hh"
30
#include "ComputingTasks.hh"
31
#include "Shocks.hh"
32

33
ModFile::ModFile(WarningConsolidation &warnings_arg)
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
  : var_model_table{symbol_table},
    trend_component_model_table{symbol_table},
    expressions_tree{symbol_table, num_constants, external_functions_table},
    original_model{symbol_table, num_constants, external_functions_table,
                   trend_component_model_table, var_model_table},
    dynamic_model{symbol_table, num_constants, external_functions_table,
                  trend_component_model_table, var_model_table},
    trend_dynamic_model{symbol_table, num_constants, external_functions_table,
                        trend_component_model_table, var_model_table},
    ramsey_FOC_equations_dynamic_model{symbol_table, num_constants, external_functions_table,
                                       trend_component_model_table, var_model_table},
    orig_ramsey_dynamic_model{symbol_table, num_constants, external_functions_table,
                              trend_component_model_table, var_model_table},
    non_linear_equations_dynamic_model{symbol_table, num_constants, external_functions_table,
                  trend_component_model_table, var_model_table},
    epilogue{symbol_table, num_constants, external_functions_table,
             trend_component_model_table, var_model_table},
    static_model{symbol_table, num_constants, external_functions_table},
    steady_state_model{symbol_table, num_constants, external_functions_table, static_model},
    diff_static_model{symbol_table, num_constants, external_functions_table},
    warnings{warnings_arg}
55
56
57
{
}

58
void
59
ModFile::evalAllExpressions(bool warn_uninit, const bool nopreprocessoroutput)
60
{
61
62
  if (!nopreprocessoroutput)
    cout << "Evaluating expressions...";
63

sebastien's avatar
sebastien committed
64
  // Loop over all statements, and fill global eval context if relevant
65
  for (auto &st : statements)
66
    {
67
      auto ips = dynamic_cast<InitParamStatement *>(st.get());
sebastien's avatar
sebastien committed
68
69
70
      if (ips)
        ips->fillEvalContext(global_eval_context);

71
      auto ies = dynamic_cast<InitOrEndValStatement *>(st.get());
sebastien's avatar
sebastien committed
72
73
74
      if (ies)
        ies->fillEvalContext(global_eval_context);

75
      auto lpass = dynamic_cast<LoadParamsAndSteadyStateStatement *>(st.get());
sebastien's avatar
sebastien committed
76
77
      if (lpass)
        lpass->fillEvalContext(global_eval_context);
78
    }
sebastien's avatar
sebastien committed
79
80

  // Evaluate model local variables
sebastien's avatar
sebastien committed
81
  dynamic_model.fillEvalContext(global_eval_context);
sebastien's avatar
sebastien committed
82

83
84
  if (!nopreprocessoroutput)
    cout << "done" << endl;
sebastien's avatar
sebastien committed
85
86

  // Check if some symbols are not initialized, and give them a zero value then
87
  for (int id = 0; id <= symbol_table.maxID(); id++)
88
    {
sebastien's avatar
sebastien committed
89
      SymbolType type = symbol_table.getType(id);
90
91
      if ((type == SymbolType::endogenous || type == SymbolType::exogenous || type == SymbolType::exogenousDet
           || type == SymbolType::parameter || type == SymbolType::modelLocalVariable)
sebastien's avatar
sebastien committed
92
          && global_eval_context.find(id) == global_eval_context.end())
93
        {
94
          if (warn_uninit)
95
96
            warnings << "WARNING: Can't find a numeric initial value for "
                     << symbol_table.getName(id) << ", using zero" << endl;
sebastien's avatar
sebastien committed
97
          global_eval_context[id] = 0;
98
99
100
101
        }
    }
}

102
void
103
ModFile::addStatement(unique_ptr<Statement> st)
104
{
105
  statements.push_back(move(st));
106
107
}

108
void
109
ModFile::addStatementAtFront(unique_ptr<Statement> st)
110
{
111
  statements.insert(statements.begin(), move(st));
112
113
}

114
void
115
ModFile::checkPass(bool nostrict, bool stochastic)
116
{
117
118
  for (auto & statement : statements)
    statement->checkPass(mod_file_struct, warnings);
119

120
  // Check the steady state block
121
122
  steady_state_model.checkPass(mod_file_struct, warnings);

Houtan Bastani's avatar
Houtan Bastani committed
123
124
125
  // Check epilogue block
  epilogue.checkPass(warnings);

126
127
128
129
130
131
  if (mod_file_struct.write_latex_steady_state_model_present &&
      !mod_file_struct.steady_state_model_present)
    {
      cerr << "ERROR: You cannot have a write_latex_steady_state_model statement without a steady_state_model block." << endl;
      exit(EXIT_FAILURE);
    }
132

133
134
135
136
  // If order option has not been set, default to 2
  if (!mod_file_struct.order_option)
    mod_file_struct.order_option = 2;

137
138
139
140
  param_used_with_lead_lag = dynamic_model.ParamUsedWithLeadLag();
  if (param_used_with_lead_lag)
    warnings << "WARNING: A parameter was used with a lead or a lag in the model block" << endl;

sebastien's avatar
sebastien committed
141
142
143
  bool stochastic_statement_present = mod_file_struct.stoch_simul_present
    || mod_file_struct.estimation_present
    || mod_file_struct.osr_present
144
    || mod_file_struct.ramsey_policy_present
145
    || mod_file_struct.discretionary_policy_present
146
147
    || mod_file_struct.calib_smoother_present
    || stochastic;
sebastien's avatar
sebastien committed
148

149
  // Allow empty model only when doing a standalone BVAR estimation
sebastien's avatar
sebastien committed
150
  if (dynamic_model.equation_number() == 0
151
      && (mod_file_struct.check_present
152
          || mod_file_struct.perfect_foresight_solver_present
sebastien's avatar
sebastien committed
153
154
155
          || stochastic_statement_present))
    {
      cerr << "ERROR: At least one model equation must be declared!" << endl;
156
      exit(EXIT_FAILURE);
sebastien's avatar
sebastien committed
157
158
    }

159
160
161
  if ((mod_file_struct.ramsey_model_present || mod_file_struct.ramsey_policy_present)
      && mod_file_struct.discretionary_policy_present)
    {
Johannes Pfeifer's avatar
Johannes Pfeifer committed
162
      cerr << "ERROR: You cannot use the discretionary_policy command when you use either ramsey_model or ramsey_policy and vice versa" << endl;
163
164
165
      exit(EXIT_FAILURE);
    }

166
  if (((mod_file_struct.ramsey_model_present || mod_file_struct.discretionary_policy_present)
167
       && !mod_file_struct.planner_objective_present)
168
      || (!(mod_file_struct.ramsey_model_present || mod_file_struct.discretionary_policy_present)
169
          && mod_file_struct.planner_objective_present))
170
    {
171
      cerr << "ERROR: A planner_objective statement must be used with a ramsey_model, a ramsey_policy or a discretionary_policy statement and vice versa." << endl;
172
173
174
      exit(EXIT_FAILURE);
    }

175
176
177
178
179
180
181
182
  if ((mod_file_struct.osr_present && (!mod_file_struct.osr_params_present || !mod_file_struct.optim_weights_present))
      || ((!mod_file_struct.osr_present || !mod_file_struct.osr_params_present) && mod_file_struct.optim_weights_present)
      || ((!mod_file_struct.osr_present || !mod_file_struct.optim_weights_present) && mod_file_struct.osr_params_present))
    {
      cerr << "ERROR: The osr statement must be used with osr_params and optim_weights." << endl;
      exit(EXIT_FAILURE);
    }

183
  if (mod_file_struct.perfect_foresight_solver_present && stochastic_statement_present)
184
    {
185
      cerr << "ERROR: A .mod file cannot contain both one of {perfect_foresight_solver,simul} and one of {stoch_simul, estimation, osr, ramsey_policy, discretionary_policy}. This is not possible: one cannot mix perfect foresight context with stochastic context in the same file." << endl;
186
      exit(EXIT_FAILURE);
187
188
    }

189
  if (mod_file_struct.k_order_solver && byte_code)
sebastien's avatar
sebastien committed
190
    {
191
      cerr << "ERROR: 'k_order_solver' (which is implicit if order >= 3), is not yet compatible with 'bytecode'." << endl;
sebastien's avatar
sebastien committed
192
193
194
      exit(EXIT_FAILURE);
    }

195
196
197
198
199
200
201
202
203
204
205
206
  if (use_dll && (block || byte_code || linear_decomposition))
    {
      cerr << "ERROR: In 'model' block, 'use_dll' option is not compatible with 'block', 'bytecode' or 'linear_decomposition'" << endl;
      exit(EXIT_FAILURE);
    }
  if (block && linear_decomposition)
    {
      cerr << "ERROR: In 'model' block, 'block' option is not compatible with 'linear_decomposition'" << endl;
      exit(EXIT_FAILURE);
    }

  if (!byte_code && linear_decomposition)
207
    {
208
      cerr << "ERROR: For the moment in 'model' block, 'linear_decomposition' option is compatible only with 'bytecode' option" << endl;
209
210
211
      exit(EXIT_FAILURE);
    }

212
213
214
215
216
217
218
  if (block || byte_code)
    if (dynamic_model.isModelLocalVariableUsed())
      {
        cerr << "ERROR: In 'model' block, 'block' or 'bytecode' options are not yet compatible with pound expressions" << endl;
        exit(EXIT_FAILURE);
      }

219
  if ((stochastic_statement_present || mod_file_struct.check_present || mod_file_struct.steady_present) && no_static)
220
    {
221
      cerr << "ERROR: no_static option is incompatible with stoch_simul, estimation, osr, ramsey_policy, discretionary_policy, steady and check commands" << endl;
222
223
      exit(EXIT_FAILURE);
    }
224

225
226
227
228
229
230
231
232
  if (mod_file_struct.dsge_var_estimated)
    if (!mod_file_struct.dsge_prior_weight_in_estimated_params)
      {
        cerr << "ERROR: When estimating a DSGE-VAR model and estimating the weight of the prior, dsge_prior_weight must "
             << "be referenced in the estimated_params block." << endl;
        exit(EXIT_FAILURE);
      }

233
234
  if (symbol_table.exists("dsge_prior_weight"))
    {
235
      if (symbol_table.getType("dsge_prior_weight") != SymbolType::parameter)
236
237
238
239
240
        {
          cerr << "ERROR: dsge_prior_weight may only be used as a parameter." << endl;
          exit(EXIT_FAILURE);
        }
      else
241
242
243
        warnings << "WARNING: When estimating a DSGE-Var, declaring dsge_prior_weight as a "
                 << "parameter is deprecated. The preferred method is to do this via "
                 << "the dsge_var option in the estimation statement." << endl;
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260

      if (mod_file_struct.dsge_var_estimated || !mod_file_struct.dsge_var_calibrated.empty())
        {
          cerr << "ERROR: dsge_prior_weight can either be declared as a parameter (deprecated) or via the dsge_var option "
               << "to the estimation statement (preferred), but not both." << endl;
          exit(EXIT_FAILURE);
        }

      if (!mod_file_struct.dsge_prior_weight_initialized && !mod_file_struct.dsge_prior_weight_in_estimated_params)
        {
          cerr << "ERROR: If dsge_prior_weight is declared as a parameter, it must either be initialized or placed in the "
               << "estimated_params block." << endl;
          exit(EXIT_FAILURE);
        }

      if (mod_file_struct.dsge_prior_weight_initialized && mod_file_struct.dsge_prior_weight_in_estimated_params)
        {
261
          cerr << "ERROR: dsge_prior_weight cannot be both initialized and estimated." << endl;
262
263
264
265
          exit(EXIT_FAILURE);
        }
    }

266
  if (mod_file_struct.dsge_prior_weight_in_estimated_params)
267
    if (!mod_file_struct.dsge_var_estimated && !mod_file_struct.dsge_var_calibrated.empty())
268
      {
269
        cerr << "ERROR: If dsge_prior_weight is in the estimated_params block, the prior weight cannot be calibrated "
270
271
272
             << "via the dsge_var option in the estimation statement." << endl;
        exit(EXIT_FAILURE);
      }
273
274
275
276
277
278
    else if (!mod_file_struct.dsge_var_estimated && !symbol_table.exists("dsge_prior_weight"))
      {
        cerr << "ERROR: If dsge_prior_weight is in the estimated_params block, it must either be declared as a parameter "
             << "(deprecated) or the dsge_var option must be passed to the estimation statement (preferred)." << endl;
        exit(EXIT_FAILURE);
      }
279
280
281
282
283
284

  if (dynamic_model.staticOnlyEquationsNbr() != dynamic_model.dynamicOnlyEquationsNbr())
    {
      cerr << "ERROR: the number of equations marked [static] must be equal to the number of equations marked [dynamic]" << endl;
      exit(EXIT_FAILURE);
    }
285
286
287

  if (dynamic_model.staticOnlyEquationsNbr() > 0
      && (mod_file_struct.ramsey_model_present || mod_file_struct.discretionary_policy_present))
288
    {
289
      cerr << "ERROR: marking equations as [static] or [dynamic] is not possible with ramsey_model, ramsey_policy or discretionary_policy" << endl;
290
291
      exit(EXIT_FAILURE);
    }
292

293
  if (stochastic_statement_present
294
295
296
297
298
299
300
301
302
303
      && (dynamic_model.isUnaryOpUsed(UnaryOpcode::sign)
          || dynamic_model.isUnaryOpUsed(UnaryOpcode::abs)
          || dynamic_model.isBinaryOpUsed(BinaryOpcode::max)
          || dynamic_model.isBinaryOpUsed(BinaryOpcode::min)
          || dynamic_model.isBinaryOpUsed(BinaryOpcode::greater)
          || dynamic_model.isBinaryOpUsed(BinaryOpcode::less)
          || dynamic_model.isBinaryOpUsed(BinaryOpcode::greaterEqual)
          || dynamic_model.isBinaryOpUsed(BinaryOpcode::lessEqual)
          || dynamic_model.isBinaryOpUsed(BinaryOpcode::equalEqual)
          || dynamic_model.isBinaryOpUsed(BinaryOpcode::different)))
304
    warnings << "WARNING: you are using a function (max, min, abs, sign) or an operator (<, >, <=, >=, ==, !=) which is unsuitable for a stochastic context; see the reference manual, section about \"Expressions\", for more details." << endl;
305

306
  if (linear
307
308
309
310
311
312
313
314
315
316
      && (dynamic_model.isUnaryOpUsed(UnaryOpcode::sign)
          || dynamic_model.isUnaryOpUsed(UnaryOpcode::abs)
          || dynamic_model.isBinaryOpUsed(BinaryOpcode::max)
          || dynamic_model.isBinaryOpUsed(BinaryOpcode::min)
          || dynamic_model.isBinaryOpUsed(BinaryOpcode::greater)
          || dynamic_model.isBinaryOpUsed(BinaryOpcode::less)
          || dynamic_model.isBinaryOpUsed(BinaryOpcode::greaterEqual)
          || dynamic_model.isBinaryOpUsed(BinaryOpcode::lessEqual)
          || dynamic_model.isBinaryOpUsed(BinaryOpcode::equalEqual)
          || dynamic_model.isBinaryOpUsed(BinaryOpcode::different)))
317
318
    warnings << "WARNING: you have declared your model 'linear' but you are using a function (max, min, abs, sign) or an operator (<, >, <=, >=, ==, !=) which potentially makes it non-linear." << endl;

319
320
321
322
323
324
325
326
327
328
329
  // Test if some estimated parameters are used within the values of shocks
  // statements (see issue #469)
  set<int> parameters_intersect;
  set_intersection(mod_file_struct.parameters_within_shocks_values.begin(),
                   mod_file_struct.parameters_within_shocks_values.end(),
                   mod_file_struct.estimated_parameters.begin(),
                   mod_file_struct.estimated_parameters.end(),
                   inserter(parameters_intersect, parameters_intersect.begin()));
  if (parameters_intersect.size() > 0)
    {
      cerr << "ERROR: some estimated parameters (";
330
      for (auto it = parameters_intersect.begin();
331
           it != parameters_intersect.end();)
332
333
334
335
336
337
338
339
        {
          cerr << symbol_table.getName(*it);
          if (++it != parameters_intersect.end())
            cerr << ", ";
        }
      cerr << ") also appear in the expressions defining the variance/covariance matrix of shocks; this is not allowed." << endl;
      exit(EXIT_FAILURE);
    }
340

341
  // Check if some exogenous is not used in the model block, Issue #841
Houtan Bastani's avatar
Houtan Bastani committed
342
343
344
345
346
  set<int> unusedExo0 = dynamic_model.findUnusedExogenous();
  set<int> unusedExo;
  set_difference(unusedExo0.begin(), unusedExo0.end(),
                 mod_file_struct.pac_params.begin(), mod_file_struct.pac_params.end(),
                 inserter(unusedExo, unusedExo.begin()));
347
  if (unusedExo.size() > 0)
348
    {
349
      ostringstream unused_exos;
350
351
      for (int it : unusedExo)
        unused_exos << symbol_table.getName(it) << " ";
352
353
354
355
356

      if (nostrict)
        warnings << "WARNING: " << unused_exos.str()
                 << "not used in model block, removed by nostrict command-line option" << endl;
      else
357
        {
358
359
          cerr << "ERROR: " << unused_exos.str() << "not used in model block. To bypass this error, use the `nostrict` option. This may lead to crashes or unexpected behavior." << endl;
          exit(EXIT_FAILURE);
360
361
        }
    }
sebastien's avatar
sebastien committed
362
363
364
}

void
365
ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const bool nopreprocessoroutput, const bool transform_unary_ops)
sebastien's avatar
sebastien committed
366
{
367
  // Save the original model (must be done before any model transformations by preprocessor)
368
  // - except adl and diff which we always want expanded
369
  dynamic_model.substituteAdl();
Stéphane Adjemian's avatar
Stéphane Adjemian committed
370
  dynamic_model.setLeadsLagsOrig();
371
  original_model = dynamic_model;
372

373
374
  if (nostrict)
    {
375
      set<int> unusedEndogs = dynamic_model.findUnusedEndogenous();
376
      for (int unusedEndog : unusedEndogs)
377
        {
378
          symbol_table.changeType(unusedEndog, SymbolType::unusedEndogenous);
379
          warnings << "WARNING: '" << symbol_table.getName(unusedEndog)
380
381
382
383
                   << "' not used in model block, removed by nostrict command-line option" << endl;
        }
    }

384
  // Get all equation tags associated with VARs and Pac Models
385
  set<string> eqtags;
Houtan Bastani's avatar
Houtan Bastani committed
386
387
  for (auto const & it : trend_component_model_table.getEqTags())
    for (auto & it1 : it.second)
388
      eqtags.insert(it1);
389

Houtan Bastani's avatar
Houtan Bastani committed
390
391
  for (auto const & it : var_model_table.getEqTags())
    for (auto & it1 : it.second)
Houtan Bastani's avatar
Houtan Bastani committed
392
393
      eqtags.insert(it1);

394
  if (transform_unary_ops)
395
396
    dynamic_model.substituteUnaryOps(diff_static_model);
  else
Houtan Bastani's avatar
Houtan Bastani committed
397
    // substitute only those unary ops that appear in auxiliary model equations
398
    dynamic_model.substituteUnaryOps(diff_static_model, eqtags);
399

400
  // Create auxiliary variable and equations for Diff operators that appear in VAR equations
Houtan Bastani's avatar
Houtan Bastani committed
401
  ExprNode::subst_table_t diff_subst_table;
402
  dynamic_model.substituteDiff(diff_static_model, diff_subst_table);
403

404
405
406
  // Fill Trend Component Model Table
  dynamic_model.fillTrendComponentModelTable();
  original_model.fillTrendComponentModelTableFromOrigModel(diff_static_model);
407
  dynamic_model.fillTrendComponentmodelTableAREC(diff_subst_table);
Houtan Bastani's avatar
Houtan Bastani committed
408
409
  dynamic_model.fillVarModelTable();
  original_model.fillVarModelTableFromOrigModel(diff_static_model);
410

Houtan Bastani's avatar
Houtan Bastani committed
411
  // Pac Model
412
  for (auto & statement : statements)
413
    {
414
      auto pms = dynamic_cast<PacModelStatement *>(statement.get());
415
      if (pms != nullptr)
Houtan Bastani's avatar
Houtan Bastani committed
416
         {
417
418
419
420
421
422
423
424
425
426
           vector<int> lhs;
           vector<bool> nonstationary;
           int growth_symb_id, max_lag;
           string aux_model_name, pac_model_name;
           tie ( pac_model_name, aux_model_name, growth_symb_id ) =
             pms->getPacModelInfoForPacExpectation();
           if (trend_component_model_table.isExistingTrendComponentModelName(aux_model_name))
             {
               max_lag = trend_component_model_table.getMaxLag(aux_model_name) + 1;
               lhs = dynamic_model.getUndiffLHSForPac(aux_model_name, diff_subst_table);
427
428
               // All lhs variables in a trend component model are nonstationary
               nonstationary.insert(nonstationary.end(), trend_component_model_table.getDiff(aux_model_name).size(), true);
429
             }
Houtan Bastani's avatar
Houtan Bastani committed
430
431
           else if (var_model_table.isExistingVarModelName(aux_model_name))
             {
432
               max_lag = var_model_table.getMaxLag(aux_model_name);
Houtan Bastani's avatar
Houtan Bastani committed
433
               lhs = var_model_table.getLhs(aux_model_name);
434
435
               // nonstationary variables in a VAR are those that are in diff
               nonstationary = var_model_table.getDiff(aux_model_name);
Houtan Bastani's avatar
Houtan Bastani committed
436
             }
437
           else
Houtan Bastani's avatar
Houtan Bastani committed
438
             {
Houtan Bastani's avatar
Houtan Bastani committed
439
440
               cerr << "Error: aux_model_name not recognized as VAR model or Trend Component model"
                    << endl;
441
               exit(EXIT_FAILURE);
Houtan Bastani's avatar
Houtan Bastani committed
442
             }
443
444
445
446
447
448
           pms->fillUndiffedLHS(lhs);
           dynamic_model.walkPacParameters();
           int pac_max_lag = original_model.getPacMaxLag(pac_model_name);
           dynamic_model.fillPacExpectationVarInfo(pac_model_name, lhs, max_lag,
                                                   pac_max_lag, nonstationary, growth_symb_id);
           dynamic_model.substitutePacExpectation();
Houtan Bastani's avatar
Houtan Bastani committed
449
450
         }
     }
Houtan Bastani's avatar
Houtan Bastani committed
451
452
453

  dynamic_model.addEquationsForVar();

454
455
456
  if (symbol_table.predeterminedNbr() > 0)
    dynamic_model.transformPredeterminedVariables();

457
458
459
  // Create auxiliary vars for Expectation operator
  dynamic_model.substituteExpectation(mod_file_struct.partial_information);

460
461
462
  if (nonstationary_variables)
    {
      dynamic_model.detrendEquations();
463
      trend_dynamic_model = dynamic_model;
464
465
466
      dynamic_model.removeTrendVariableFromEquations();
    }

467
  mod_file_struct.orig_eq_nbr = dynamic_model.equation_number();
468
  if (mod_file_struct.ramsey_model_present)
469
    {
470
      PlannerObjectiveStatement *pos = nullptr;
471
      for (auto & statement : statements)
472
        {
473
          auto pos2 = dynamic_cast<PlannerObjectiveStatement *>(statement.get());
474
475
476
477
478
479
480
481
          if (pos2 != nullptr)
            if (pos != nullptr)
              {
                cerr << "ERROR: there can only be one planner_objective statement" << endl;
                exit(EXIT_FAILURE);
              }
            else
              pos = pos2;
482
        }
483
484
      assert(pos != nullptr);
      const StaticModel &planner_objective = pos->getPlannerObjective();
485
486
487
488

      /*
        clone the model then clone the new equations back to the original because
        we have to call computeDerivIDs (in computeRamseyPolicyFOCs and computingPass)
489
      */
490
      if (linear)
491
492
        orig_ramsey_dynamic_model = dynamic_model;
      ramsey_FOC_equations_dynamic_model = dynamic_model;
493
      ramsey_FOC_equations_dynamic_model.computeRamseyPolicyFOCs(planner_objective, nopreprocessoroutput);
494
      ramsey_FOC_equations_dynamic_model.replaceMyEquations(dynamic_model);
495
      mod_file_struct.ramsey_eq_nbr = dynamic_model.equation_number() - mod_file_struct.orig_eq_nbr;
496
497
    }

498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
  // Workaround for #1193
  if (!mod_file_struct.hist_vals_wrong_lag.empty())
    {
      bool err = false;
      for (map<int, int>::const_iterator it = mod_file_struct.hist_vals_wrong_lag.begin();
           it != mod_file_struct.hist_vals_wrong_lag.end(); it++)
          if (dynamic_model.minLagForSymbol(it->first) > it->second - 1)
            {
              cerr << "ERROR: histval: variable " << symbol_table.getName(it->first)
                   << " does not appear in the model with the lag " << it->second - 1
                   << " (see the reference manual for the timing convention in 'histval')" << endl;
              err = true;
            }
      if (err)
        exit(EXIT_FAILURE);
    }

515
516
517
518
519
520
521
522
  /* Handle var_expectation_model statements: collect information about them,
     create the new corresponding parameters, and the expressions to replace
     the var_expectation statements.
     TODO: move information collection to checkPass(), within a new
     VarModelTable class */
  map<string, expr_t> var_expectation_subst_table;
  for (auto & statement : statements)
    {
523
      auto vems = dynamic_cast<VarExpectationModelStatement *>(statement.get());
524
525
526
      if (!vems)
        continue;

527
528
      int max_lag;
      vector<int> lhs;
529
      auto &model_name = vems->model_name;
530
      if (var_model_table.isExistingVarModelName(vems->aux_model_name))
531
        {
532
533
          max_lag = var_model_table.getMaxLag(vems->aux_model_name);
          lhs = var_model_table.getLhs(vems->aux_model_name);
534
        }
535
      else if (trend_component_model_table.isExistingTrendComponentModelName(vems->aux_model_name))
536
        {
537
          max_lag = trend_component_model_table.getMaxLag(vems->aux_model_name) + 1;
538
          lhs = dynamic_model.getUndiffLHSForPac(vems->aux_model_name, diff_subst_table);
539
540
        }
      else
541
        {
Houtan Bastani's avatar
Houtan Bastani committed
542
          cerr << "ERROR: var_expectation_model " << model_name
543
               << " refers to nonexistent auxiliary model " << vems->aux_model_name << endl;
544
545
546
          exit(EXIT_FAILURE);
        }

547
      /* Create auxiliary parameters and the expression to be substituted into
548
549
         the var_expectations statement */
      auto subst_expr = dynamic_model.Zero;
550
      for (int lag = 0; lag < max_lag; lag++)
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
        for (auto variable : lhs)
          {
            string param_name = "var_expectation_model_" + model_name + '_' + symbol_table.getName(variable) + '_' + to_string(lag);
            int new_param_id = symbol_table.addSymbol(param_name, SymbolType::parameter);
            vems->aux_params_ids.push_back(new_param_id);

            subst_expr = dynamic_model.AddPlus(subst_expr,
                                               dynamic_model.AddTimes(dynamic_model.AddVariable(new_param_id),
                                                                      dynamic_model.AddVariable(variable, -lag)));
          }

      if (var_expectation_subst_table.find(model_name) != var_expectation_subst_table.end())
        {
          cerr << "ERROR: model name '" << model_name << "' is used by several var_expectation_model statements" << endl;
          exit(EXIT_FAILURE);
        }
      var_expectation_subst_table[model_name] = subst_expr;
    }
  // And finally perform the substitutions
  dynamic_model.substituteVarExpectation(var_expectation_subst_table);

sebastien's avatar
sebastien committed
572
573
574
  if (mod_file_struct.stoch_simul_present
      || mod_file_struct.estimation_present
      || mod_file_struct.osr_present
575
      || mod_file_struct.ramsey_policy_present
576
      || mod_file_struct.discretionary_policy_present
577
578
      || mod_file_struct.calib_smoother_present
      || stochastic )
sebastien's avatar
sebastien committed
579
    {
580
581
582
583
584
585
586
587
588
589
590
      // In stochastic models, create auxiliary vars for leads and lags greater than 2, on both endos and exos
      dynamic_model.substituteEndoLeadGreaterThanTwo(false);
      dynamic_model.substituteExoLead(false);
      dynamic_model.substituteEndoLagGreaterThanTwo(false);
      dynamic_model.substituteExoLag(false);
    }
  else
    {
      // In deterministic models, create auxiliary vars for leads and lags endogenous greater than 2, only on endos (useless on exos)
      dynamic_model.substituteEndoLeadGreaterThanTwo(true);
      dynamic_model.substituteEndoLagGreaterThanTwo(true);
sebastien's avatar
sebastien committed
591
    }
592

593
  dynamic_model.updateVarAndTrendModel();
Houtan Bastani's avatar
Houtan Bastani committed
594

595
  if (differentiate_forward_vars)
596
    dynamic_model.differentiateForwardVars(differentiate_forward_vars_subset);
597

598
  if (mod_file_struct.dsge_var_estimated || !mod_file_struct.dsge_var_calibrated.empty())
599
600
    try
      {
601
        int sid = symbol_table.addSymbol("dsge_prior_weight", SymbolType::parameter);
602
        if (!mod_file_struct.dsge_var_calibrated.empty())
603
604
605
          addStatementAtFront(make_unique<InitParamStatement>(sid,
                                                              expressions_tree.AddNonNegativeConstant(mod_file_struct.dsge_var_calibrated),
                                                              symbol_table));
606
607
608
609
610
611
612
613
      }
    catch (SymbolTable::AlreadyDeclaredException &e)
      {
        cerr << "ERROR: dsge_prior_weight should not be declared as a model variable / parameter "
             << "when the dsge_var option is passed to the estimation statement." << endl;
        exit(EXIT_FAILURE);
      }

sebastien's avatar
sebastien committed
614
615
616
  // Freeze the symbol table
  symbol_table.freeze();

617
618
619
  if (compute_xrefs)
    dynamic_model.computeXrefs();

620
  /*
621
    Enforce the same number of equations and endogenous, except in three cases:
622
    - ramsey_model, ramsey_policy or discretionary_policy is used
623
    - a BVAR command is used and there is no equation (standalone BVAR estimation)
624
    - nostrict option is passed and there are more endogs than equations (dealt with before freeze)
625
  */
626
  if (!(mod_file_struct.ramsey_model_present || mod_file_struct.discretionary_policy_present)
627
      && !(mod_file_struct.bvar_present && dynamic_model.equation_number() == 0)
628
      && !(mod_file_struct.occbin_option)
sebastien's avatar
sebastien committed
629
      && (dynamic_model.equation_number() != symbol_table.endo_nbr()))
630
    {
sebastien's avatar
sebastien committed
631
      cerr << "ERROR: There are " << dynamic_model.equation_number() << " equations but " << symbol_table.endo_nbr() << " endogenous variables!" << endl;
632
      exit(EXIT_FAILURE);
633
    }
634

635
  if (symbol_table.exo_det_nbr() > 0 && mod_file_struct.perfect_foresight_solver_present)
636
    {
Houtan Bastani's avatar
Houtan Bastani committed
637
      cerr << "ERROR: A .mod file cannot contain both one of {perfect_foresight_solver, simul} and varexo_det declaration (all exogenous variables are deterministic in this case)" << endl;
638
639
      exit(EXIT_FAILURE);
    }
640

641
642
643
644
645
646
  if (mod_file_struct.ramsey_policy_present && symbol_table.exo_det_nbr() > 0)
    {
      cerr << "ERROR: ramsey_policy is incompatible with deterministic exogenous variables" << endl;
      exit(EXIT_FAILURE);
    }

647
  if (mod_file_struct.ramsey_policy_present)
648
    for (auto & statement : statements)
649
      {
650
        auto *rps = dynamic_cast<RamseyPolicyStatement *>(statement.get());
651
        if (rps != nullptr)
652
653
654
          rps->checkRamseyPolicyList();
      }

655
656
657
658
659
660
  if (mod_file_struct.identification_present && symbol_table.exo_det_nbr() > 0)
    {
      cerr << "ERROR: identification is incompatible with deterministic exogenous variables" << endl;
      exit(EXIT_FAILURE);
    }

661
662
663
664
665
666
667
668
  if (!nopreprocessoroutput)
    if (!mod_file_struct.ramsey_model_present)
      cout << "Found " << dynamic_model.equation_number() << " equation(s)." << endl;
    else
      {
        cout << "Found " << mod_file_struct.orig_eq_nbr  << " equation(s)." << endl;
        cout << "Found " << dynamic_model.equation_number() << " FOC equation(s) for Ramsey Problem." << endl;
      }
669

670
  if (symbol_table.exists("dsge_prior_weight"))
671
672
673
674
    if (mod_file_struct.bayesian_irf_present)
      {
        if (symbol_table.exo_nbr() != symbol_table.observedVariablesNbr())
          {
675
676
            cerr << "ERROR: When estimating a DSGE-Var and the bayesian_irf option is passed to the estimation "
                 << "statement, the number of shocks must equal the number of observed variables." << endl;
677
678
679
680
681
682
            exit(EXIT_FAILURE);
          }
      }
    else
      if (symbol_table.exo_nbr() < symbol_table.observedVariablesNbr())
        {
683
684
          cerr << "ERROR: When estimating a DSGE-Var, the number of shocks must be "
               << "greater than or equal to the number of observed variables." << endl;
685
686
          exit(EXIT_FAILURE);
        }
687
688
689
}

void
690
ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_derivs_order, const bool nopreprocessoroutput)
691
692
{
  // Mod file may have no equation (for example in a standalone BVAR estimation)
sebastien's avatar
sebastien committed
693
  if (dynamic_model.equation_number() > 0)
694
    {
695
696
      if (nonstationary_variables)
        trend_dynamic_model.runTrendTest(global_eval_context);
697

sebastien's avatar
sebastien committed
698
      // Compute static model and its derivatives
699
      static_model = static_cast<StaticModel>(dynamic_model);
700
701
	  if (linear_decomposition)
        {
702
          non_linear_equations_dynamic_model = dynamic_model;
703
          non_linear_equations_dynamic_model.set_cutoff_to_zero();
704
          non_linear_equations_dynamic_model.computingPass(true, 1, 0, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition);
705
        }
706
      if (!no_static)
707
708
709
710
711
712
713
        {
          if (mod_file_struct.stoch_simul_present
              || mod_file_struct.estimation_present || mod_file_struct.osr_present
              || mod_file_struct.ramsey_model_present || mod_file_struct.identification_present
              || mod_file_struct.calib_smoother_present)
            static_model.set_cutoff_to_zero();

714
          int derivsOrder = 1;
715
716
          int paramsDerivsOrder = 0;
          if (mod_file_struct.identification_present || mod_file_struct.estimation_analytic_derivation)
717
718
719
720
721
            {
              derivsOrder = 2;
              paramsDerivsOrder = params_derivs_order;
            }
          static_model.computingPass(derivsOrder, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, byte_code, nopreprocessoroutput);
722
        }
sebastien's avatar
sebastien committed
723
      // Set things to compute for dynamic model
724
      if (mod_file_struct.perfect_foresight_solver_present || mod_file_struct.check_present
725
726
727
728
729
730
          || mod_file_struct.stoch_simul_present
          || mod_file_struct.estimation_present || mod_file_struct.osr_present
          || mod_file_struct.ramsey_model_present || mod_file_struct.identification_present
          || mod_file_struct.calib_smoother_present)
        {
          if (mod_file_struct.perfect_foresight_solver_present)
731
            dynamic_model.computingPass(true, 1, 0, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition);
732
733
734
735
736
737
738
739
740
741
742
743
          else
            {
              if (mod_file_struct.stoch_simul_present
                  || mod_file_struct.estimation_present || mod_file_struct.osr_present
                  || mod_file_struct.ramsey_model_present || mod_file_struct.identification_present
                  || mod_file_struct.calib_smoother_present)
                dynamic_model.set_cutoff_to_zero();
              if (mod_file_struct.order_option < 1 || mod_file_struct.order_option > 3)
                {
                  cerr << "ERROR: Incorrect order option..." << endl;
                  exit(EXIT_FAILURE);
                }
744
745
746
747
748
              int derivsOrder = mod_file_struct.order_option;
              if (mod_file_struct.identification_present || linear || output == FileOutputType::second)
                derivsOrder = max(derivsOrder, 2);
              if (mod_file_struct.estimation_analytic_derivation || output == FileOutputType::third)
                derivsOrder = max(derivsOrder, 3);
749
750
751
              int paramsDerivsOrder = 0;
              if (mod_file_struct.identification_present || mod_file_struct.estimation_analytic_derivation)
                paramsDerivsOrder = params_derivs_order;
752
              dynamic_model.computingPass(true, derivsOrder, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition);
753
              if (linear && mod_file_struct.ramsey_model_present)
754
                orig_ramsey_dynamic_model.computingPass(true, 2, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition);
755
756
757
            }
        }
      else // No computing task requested, compute derivatives up to 2nd order by default
758
        dynamic_model.computingPass(true, 2, 0, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition);
759

760
761
762
763
764
      map<int, string> eqs;
      if (mod_file_struct.ramsey_model_present)
        orig_ramsey_dynamic_model.setNonZeroHessianEquations(eqs);
      else
        dynamic_model.setNonZeroHessianEquations(eqs);
765

766
767
      if (linear && !eqs.empty())
        {
768
769
770
771
772
773
774
775
776
          cerr << "ERROR: If the model is declared linear the second derivatives must be equal to zero." << endl
               << "       The following equations had non-zero second derivatives:" << endl;
          for (map<int, string >::const_iterator it = eqs.begin(); it != eqs.end(); it++)
            {
              cerr << "       * Eq # " << it->first+1;
              if (!it->second.empty())
                cerr << " [" << it->second << "]";
              cerr << endl;
            }
777
778
          exit(EXIT_FAILURE);
        }
779
    }
780

781
782
  for (auto & statement : statements)
    statement->computingPass();
Houtan Bastani's avatar
Houtan Bastani committed
783

784
  epilogue.computingPass(true, 2, 0, global_eval_context, true, false, false, false, true, false);
785
786
787
}

void
788
789
ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_global, bool no_log, bool no_warn,
                          bool console, bool nograph, bool nointeractive, const ConfigFile &config_file,
790
791
792
                          bool check_model_changes, bool minimal_workspace, bool compute_xrefs,
                          const bool nopreprocessoroutput, const string &mexext,
                          const boost::filesystem::path &matlabroot,
793
                          const boost::filesystem::path &dynareroot, bool onlymodel) const
794
{
795
  bool hasModelChanged = !dynamic_model.isChecksumMatching(basename, block);
796
797
798
799
800
801
  if (!check_model_changes)
    hasModelChanged = true;

  if (hasModelChanged)
    {
      // Erase possible remnants of previous runs
802
803
804
805
806
807
808
809
810
811
812
      /* Under MATLAB+Windows (but not under Octave nor under GNU/Linux or
         macOS), if we directly remove the "+" subdirectory, then the
         preprocessor is not able to recreate it afterwards (presumably because
         MATLAB maintains some sort of lock on it). The workaround is to rename
         it before deleting it. */
      if (boost::filesystem::exists("+" + basename))
	{
	  auto tmp = boost::filesystem::unique_path();
	  boost::filesystem::rename("+" + basename, tmp);
	  boost::filesystem::remove_all(tmp);
	}
813
814
      boost::filesystem::remove_all(basename + "/model/src");
      boost::filesystem::remove_all(basename + "/model/bytecode");
815
816
    }

817
  ofstream mOutputFile;
818

819
820
  if (basename.size())
    {
821
822
      boost::filesystem::create_directory("+" + basename);
      string fname = "+" + basename + "/driver.m";
823
      mOutputFile.open(fname, ios::out | ios::binary);
824
825
      if (!mOutputFile.is_open())
        {
Houtan Bastani's avatar
Houtan Bastani committed
826
          cerr << "ERROR: Can't open file " << fname << " for writing" << endl;
827
          exit(EXIT_FAILURE);
828
829
830
831
        }
    }
  else
    {
sebastien's avatar
sebastien committed
832
      cerr << "ERROR: Missing file name" << endl;
833
      exit(EXIT_FAILURE);
834
835
    }

sebastien's avatar
sebastien committed
836
  mOutputFile << "%" << endl
Houtan Bastani's avatar
Houtan Bastani committed
837
              << "% Status : main Dynare file" << endl
sebastien's avatar
sebastien committed
838
839
840
              << "%" << endl
              << "% Warning : this file is generated automatically by Dynare" << endl
              << "%           from model file (.mod)" << endl << endl;
841

842
843
844
  if (no_warn)
    mOutputFile << "warning off" << endl; // This will be executed *after* function warning_config()

845
  if (clear_all)
846
847
    mOutputFile << "if isoctave || matlab_ver_less_than('8.6')" << endl
                << "    clear all" << endl
848
849
850
851
                << "else" << endl
                << "    clearvars -global" << endl
                << "    clear_persistent_variables(fileparts(which('dynare')), false)" << endl
                << "end" << endl;
852
  else if (clear_global)
853
    mOutputFile << "clear M_ options_ oo_ estim_params_ bayestopt_ dataset_ dataset_info estimation_info ys0_ ex0_;" << endl;
854

855
  mOutputFile << "tic0 = tic;" << endl
856
              << "% Define global variables." << endl
857
              << "global M_ options_ oo_ estim_params_ bayestopt_ dataset_ dataset_info estimation_info ys0_ ex0_" << endl
858
859
              << "options_ = [];" << endl
              << "M_.fname = '" << basename << "';" << endl
860
861
862
              << "M_.dynare_version = '" << PACKAGE_VERSION << "';" << endl
              << "oo_.dynare_version = '" << PACKAGE_VERSION << "';" << endl
              << "options_.dynare_version = '" << PACKAGE_VERSION << "';" << endl
863
864
              << "%" << endl
              << "% Some global variables initialization" << endl
865
              << "%" << endl;
866
867
  if (!onlymodel)
    config_file.writeHooks(mOutputFile);
868
  mOutputFile << "global_initialization;" << endl
869
870
              << "diary off;" << endl;
  if (!no_log)
871
    mOutputFile << "diary('" << basename << ".log');" << endl;
872

873
874
875
  if (minimal_workspace)
    mOutputFile << "options_.minimal_workspace = 1;" << endl;

876
  if (console)
877
878
    mOutputFile << "options_.console_mode = 1;" << endl
                << "options_.nodisplay = 1;" << endl;
879
880
  if (nograph)
    mOutputFile << "options_.nograph = 1;" << endl;
881
882
883

  if (nointeractive)
    mOutputFile << "options_.nointeractive = 1;" << endl;
884
885
886
887

  if (param_used_with_lead_lag)
    mOutputFile << "M_.parameter_used_with_lead_lag = true;" << endl;

888
889
  if (!nopreprocessoroutput)
    cout << "Processing outputs ..." << endl;
890
891
892

  symbol_table.writeOutput(mOutputFile);

Houtan Bastani's avatar
Houtan Bastani committed
893
  var_model_table.writeOutput(basename, mOutputFile);
894
  trend_component_model_table.writeOutput(basename, mOutputFile);
895

896
  // Initialize M_.Sigma_e, M_.Correlation_matrix, M_.H, and M_.Correlation_matrix_ME
Sébastien Villemot's avatar
Sébastien Villemot committed
897
  mOutputFile << "M_.Sigma_e = zeros(" << symbol_table.exo_nbr() << ", "
898
899
              << symbol_table.exo_nbr() << ");" << endl
              << "M_.Correlation_matrix = eye(" << symbol_table.exo_nbr() << ", "
Sébastien Villemot's avatar
Sébastien Villemot committed
900
901
902
903
              << symbol_table.exo_nbr() << ");" << endl;

  if (mod_file_struct.calibrated_measurement_errors)
    mOutputFile << "M_.H = zeros(" << symbol_table.observedVariablesNbr() << ", "
904
905
                << symbol_table.observedVariablesNbr() << ");" << endl
                << "M_.Correlation_matrix_ME = eye(" << symbol_table.observedVariablesNbr() << ", "
Sébastien Villemot's avatar
Sébastien Villemot committed
906
907
                << symbol_table.observedVariablesNbr() << ");" << endl;
  else
908
909
    mOutputFile << "M_.H = 0;" << endl
                << "M_.Correlation_matrix_ME = 1;" << endl;
Sébastien Villemot's avatar
Sébastien Villemot committed
910

911
912
913
  // May be later modified by a shocks block
  mOutputFile << "M_.sigma_e_is_diagonal = 1;" << endl;

914
915
916
  // Initialize M_.det_shocks
  mOutputFile << "M_.det_shocks = [];" << endl;

917
918
919
  if (linear == 1)
    mOutputFile << "options_.linear = 1;" << endl;

920
  mOutputFile << "options_.block=" << block << ";" << endl