NumericalInitialization.cc 19 KB
Newer Older
1
/*
Houtan Bastani's avatar
Houtan Bastani committed
2
 * Copyright (C) 2003-2019 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
21
#include <iostream>
#include <fstream>
22
#include <sstream>
23
#include <cstdlib>
24
#include <utility>
25

26
#include "NumericalInitialization.hh"
sebastien's avatar
sebastien committed
27
28

InitParamStatement::InitParamStatement(int symb_id_arg,
29
                                       const expr_t param_value_arg,
30
                                       const SymbolTable &symbol_table_arg) :
31
32
33
  symb_id{symb_id_arg},
  param_value{param_value_arg},
  symbol_table{symbol_table_arg}
34
35
36
{
}

37
void
38
InitParamStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings)
39
40
41
42
43
{
  if (symbol_table.getName(symb_id) == "dsge_prior_weight")
    mod_file_struct.dsge_prior_weight_initialized = true;
}

44
void
45
InitParamStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
46
{
sebastien's avatar
sebastien committed
47
  int id = symbol_table.getTypeSpecificID(symb_id) + 1;
48
49
50
  output << "M_.params( " << id << " ) = ";
  param_value->writeOutput(output);
  output << ";" << endl;
51
52
  if (!minimal_workspace)
    output << symbol_table.getName(symb_id) << " = M_.params( " << id << " );" << endl;
53
54
}

55
56
57
58
59
60
61
62
63
64
65
66
void
InitParamStatement::writeJuliaOutput(ostream &output, const string &basename)
{
  int id = symbol_table.getTypeSpecificID(symb_id) + 1;
  output << "model_.params[ " << id << " ] = ";
  param_value->writeOutput(output);
  output << endl;
  // Do we really need this?
  // if (!minimal_workspace)
  //   output << symbol_table.getName(symb_id) << " = model_.params[ " << id << " ]" << endl;
}

67
68
69
void
InitParamStatement::writeJsonOutput(ostream &output) const
{
70
  output << "{\"statementName\": \"param_init\", \"name\": \"" << symbol_table.getName(symb_id) << "\", " << "\"value\": \"";
71
  param_value->writeJsonOutput(output, {}, {});
72
  output << "\"}";
73
74
}

75
76
77
78
79
80
81
82
83
84
void
InitParamStatement::writeCOutput(ostream &output, const string &basename)
{
  int id = symbol_table.getTypeSpecificID(symb_id);
  output << "params[ " << id << " ] = ";
  param_value->writeOutput(output);
  output << ";" << endl;
  output << "double " << symbol_table.getName(symb_id) << " = params[ " << id << " ];" << endl;
}

sebastien's avatar
sebastien committed
85
void
86
InitParamStatement::fillEvalContext(eval_context_t &eval_context) const
87
{
sebastien's avatar
sebastien committed
88
89
90
91
  try
    {
      eval_context[symb_id] = param_value->eval(eval_context);
    }
92
  catch (ExprNode::EvalException &e)
sebastien's avatar
sebastien committed
93
94
95
    {
      // Do nothing
    }
96
97
}

98
InitOrEndValStatement::InitOrEndValStatement(init_values_t init_values_arg,
99
100
                                             const SymbolTable &symbol_table_arg,
                                             const bool &all_values_required_arg) :
101
102
103
  init_values{move(init_values_arg)},
  symbol_table{symbol_table_arg},
  all_values_required{all_values_required_arg}
104
105
106
{
}

sebastien's avatar
sebastien committed
107
void
108
InitOrEndValStatement::fillEvalContext(eval_context_t &eval_context) const
sebastien's avatar
sebastien committed
109
{
110
  for (const auto & init_value : init_values)
sebastien's avatar
sebastien committed
111
112
113
    {
      try
        {
114
          eval_context[init_value.first] = (init_value.second)->eval(eval_context);
sebastien's avatar
sebastien committed
115
        }
116
      catch (ExprNode::EvalException &e)
sebastien's avatar
sebastien committed
117
118
119
120
121
122
        {
          // Do nothing
        }
    }
}

123
124
125
126
127
128
129
set<int>
InitOrEndValStatement::getUninitializedVariables(SymbolType type)
{
  set<int> unused;
  if (!all_values_required)
    return unused;

130
  if (type == SymbolType::endogenous)
131
    unused = symbol_table.getEndogenous();
132
  else if (type == SymbolType::exogenous)
133
134
135
136
137
138
139
140
    unused = symbol_table.getExogenous();
  else
    {
      cerr << "ERROR: Shouldn't arrive here." << endl;
      exit(EXIT_FAILURE);
    }

  set<int>::iterator sit;
141
  for (const auto & init_value : init_values)
142
    {
143
      sit = unused.find(init_value.first);
144
145
146
147
148
149
      if (sit != unused.end())
        unused.erase(sit);
    }
  return unused;
}

150
151
152
void
InitOrEndValStatement::writeInitValues(ostream &output) const
{
153
  for (const auto & init_value : init_values)
154
    {
155
156
      const int symb_id = init_value.first;
      const expr_t expression = init_value.second;
157

sebastien's avatar
sebastien committed
158
159
      SymbolType type = symbol_table.getType(symb_id);
      int tsid = symbol_table.getTypeSpecificID(symb_id) + 1;
160

161
      if (type == SymbolType::endogenous)
162
        output << "oo_.steady_state";
163
      else if (type == SymbolType::exogenous)
164
        output << "oo_.exo_steady_state";
165
      else if (type == SymbolType::exogenousDet)
166
167
        output << "oo_.exo_det_steady_state";

sebastien's avatar
sebastien committed
168
      output << "( " << tsid << " ) = ";
169
170
171
172
173
      expression->writeOutput(output);
      output << ";" << endl;
    }
}

174
175
176
void
InitOrEndValStatement::writeJsonInitValues(ostream &output) const
{
177
  for (auto it = init_values.begin();
178
       it != init_values.end(); it++)
179
    {
180
181
      if (it != init_values.begin())
        output << ", ";
182
      output << "{\"name\": \"" << symbol_table.getName(it->first) << "\", " << "\"value\": \"";
183
      it->second->writeJsonOutput(output, {}, {});
184
185
186
187
      output << "\"}";
    }
}

188
InitValStatement::InitValStatement(const init_values_t &init_values_arg,
189
190
                                   const SymbolTable &symbol_table_arg,
                                   const bool &all_values_required_arg) :
191
  InitOrEndValStatement{init_values_arg, symbol_table_arg, all_values_required_arg}
192
193
194
{
}

195
196
197
void
InitValStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings)
{
198
199
  set<int> exogs = getUninitializedVariables(SymbolType::exogenous);
  set<int> endogs = getUninitializedVariables(SymbolType::endogenous);
200
201
202
203

  if (endogs.size() > 0)
    {
      cerr << "ERROR: You have not set the following endogenous variables in initval:";
204
205
      for (int endog : endogs)
        cerr << " " << symbol_table.getName(endog);
206
207
208
209
210
211
      cerr << endl;
    }

  if (exogs.size() > 0)
    {
      cerr << "ERROR: You have not set the following exogenous variables in initval:";
212
213
      for (int exog : exogs)
        cerr << " " << symbol_table.getName(exog);
214
215
216
217
218
219
220
      cerr << endl;
    }

  if (endogs.size() > 0 || exogs.size() > 0)
    exit(EXIT_FAILURE);
}

221
void
222
InitValStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
223
{
sebastien's avatar
sebastien committed
224
225
226
  output << "%" << endl
         << "% INITVAL instructions" << endl
         << "%" << endl;
227
  // Writing initval block to set initial values for variables
sebastien's avatar
trunk:    
sebastien committed
228
  output << "options_.initval_file = 0;" << endl;
229
230

  writeInitValues(output);
sebastien's avatar
sebastien committed
231
}
232

233
234
235
void
InitValStatement::writeJsonOutput(ostream &output) const
{
Houtan Bastani's avatar
Houtan Bastani committed
236
  output << "{\"statementName\": \"initval\", \"vals\": [";
237
238
239
240
  writeJsonInitValues(output);
  output << "]}";
}

sebastien's avatar
sebastien committed
241
242
243
void
InitValStatement::writeOutputPostInit(ostream &output) const
{
244
245
246
247
248
249
  output << "if M_.exo_nbr > 0" << endl
         << "\too_.exo_simul = ones(M_.maximum_lag,1)*oo_.exo_steady_state';" << endl
         <<"end" << endl
         << "if M_.exo_det_nbr > 0" << endl
         << "\too_.exo_det_simul = ones(M_.maximum_lag,1)*oo_.exo_det_steady_state';" << endl
         <<"end" << endl;
250
251
}

252
EndValStatement::EndValStatement(const init_values_t &init_values_arg,
253
254
                                 const SymbolTable &symbol_table_arg,
                                 const bool &all_values_required_arg) :
255
  InitOrEndValStatement{init_values_arg, symbol_table_arg, all_values_required_arg}
256
257
258
{
}

259
void
260
EndValStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings)
261
{
262
263
  set<int> exogs = getUninitializedVariables(SymbolType::exogenous);
  set<int> endogs = getUninitializedVariables(SymbolType::endogenous);
264
265
266
267

  if (endogs.size() > 0)
    {
      cerr << "ERROR: You have not set the following endogenous variables in endval:";
268
269
      for (int endog : endogs)
        cerr << " " << symbol_table.getName(endog);
270
271
272
273
274
275
      cerr << endl;
    }

  if (exogs.size() > 0)
    {
      cerr << "ERROR: You have not set the following exogenous variables in endval:";
276
277
      for (int exog : exogs)
        cerr << " " << symbol_table.getName(exog);
278
279
280
281
282
      cerr << endl;
    }

  if (endogs.size() > 0 || exogs.size() > 0)
    exit(EXIT_FAILURE);
283
284
}

285
void
286
EndValStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
287
{
sebastien's avatar
sebastien committed
288
289
290
  output << "%" << endl
         << "% ENDVAL instructions" << endl
         << "%" << endl;
291
  // Writing endval block to set terminal values for variables
sebastien's avatar
sebastien committed
292
  output << "ys0_= oo_.steady_state;" << endl
293
         << "ex0_ = oo_.exo_steady_state;" << endl;
294
295
296
297

  writeInitValues(output);
}

298
299
300
void
EndValStatement::writeJsonOutput(ostream &output) const
{
Houtan Bastani's avatar
Houtan Bastani committed
301
  output << "{\"statementName\": \"endval\", \"vals\": [";
302
303
304
305
  writeJsonInitValues(output);
  output << "]}";
}

306
HistValStatement::HistValStatement(hist_values_t hist_values_arg,
307
308
                                   const SymbolTable &symbol_table_arg,
                                   const bool &all_values_required_arg) :
309
310
311
  hist_values{move(hist_values_arg)},
  symbol_table{symbol_table_arg},
  all_values_required{all_values_required_arg}
312
313
314
{
}

315
void
316
HistValStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings)
317
{
318
319
320
321
322
  if (all_values_required)
    {
      set<int> unused_endo = symbol_table.getEndogenous();
      set<int> unused_exo = symbol_table.getExogenous();

323
      for (const auto & hist_value : hist_values)
324
        {
325
          auto sit = unused_endo.find(hist_value.first.first);
326
327
328
          if (sit != unused_endo.end())
            unused_endo.erase(sit);

329
          sit = unused_exo.find(hist_value.first.first);
330
331
332
333
334
335
336
          if (sit != unused_exo.end())
            unused_exo.erase(sit);
        }

      if (unused_endo.size() > 0)
        {
          cerr << "ERROR: You have not set the following endogenous variables in histval:";
337
338
          for (int it : unused_endo)
            cerr << " " << symbol_table.getName(it);
339
340
341
342
343
344
          cerr << endl;
        }

      if (unused_exo.size() > 0)
        {
          cerr << "ERROR: You have not set the following exogenous variables in endval:";
345
346
          for (int it : unused_exo)
            cerr << " " << symbol_table.getName(it);
347
348
349
350
351
352
          cerr << endl;
        }

      if (unused_endo.size() > 0 || unused_exo.size() > 0)
        exit(EXIT_FAILURE);
    }
353
354
}

355
void
356
HistValStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
357
{
sebastien's avatar
sebastien committed
358
359
  output << "%" << endl
         << "% HISTVAL instructions" << endl
360
         << "%" << endl
361
362
         << "M_.histval_dseries = dseries(zeros(M_.orig_maximum_lag_with_diffs_expanded, M_.orig_endo_nbr"
         << (symbol_table.AuxVarsSize() > 0 ? "+sum([M_.aux_vars.type]==6)" : "")
363
364
         << (symbol_table.exo_nbr() > 0 ? "+M_.exo_nbr" : "")
         << (symbol_table.exo_det_nbr() > 0 ? "+M_.exo_det_nbr" : "")
365
366
         << "), dates(sprintf('%dY', -M_.orig_maximum_lag_with_diffs_expanded+1)), [ M_.endo_names(1:M_.orig_endo_nbr); "
         << (symbol_table.AuxVarsSize() > 0 ? "M_.endo_names([M_.aux_vars(find([M_.aux_vars.type]==6)).endo_index]); " : "")
367
368
369
         << (symbol_table.exo_nbr() > 0 ? "M_.exo_names; " : "")
         << (symbol_table.exo_det_nbr() > 0 ? "M_.exo_det_names; " : "")
         << "]);" << endl;
370

371
  for (const auto & hist_value : hist_values)
372
    {
373
374
375
      int symb_id = hist_value.first.first;
      int lag = hist_value.first.second;
      const expr_t expression = hist_value.second;
376

377
      output << "M_.histval_dseries{'" << symbol_table.getName(symb_id) << "'}(dates('" << lag << "Y'))=";
378
379
380
      expression->writeOutput(output);
      output << ";" << endl;
    }
381

Sébastien Villemot's avatar
Sébastien Villemot committed
382
  output << "if exist(['+' M_.fname '/dynamic_set_auxiliary_series.m'])" << endl
383
       << "  eval(['M_.histval_dseries = ' M_.fname '.dynamic_set_auxiliary_series(M_.histval_dseries, M_.params);']);" << endl
384
       << "end" << endl
385
386
         << "M_.endo_histval = M_.histval_dseries{M_.endo_names{:}}(dates(sprintf('%dY', 1-M_.maximum_lag)):dates('0Y')).data';" << endl
         << "M_.endo_histval(isnan(M_.endo_histval)) = 0;" << endl; // Ensure that lead aux variables do not have a NaN
387
388
389
390
391

  if (symbol_table.exo_nbr() > 0)
    output << "M_.exo_histval = M_.histval_dseries{M_.exo_names{:}}(dates(sprintf('%dY', 1-M_.maximum_lag)):dates('0Y')).data';" << endl;
  if (symbol_table.exo_det_nbr() > 0)
    output << "M_.exo_det_histval = M_.histval_dseries{M_.exo_det_names{:}}(dates(sprintf('%dY', 1-M_.maximum_lag)):dates('0Y')).data';" << endl;
392
393
}

394
395
396
void
HistValStatement::writeJsonOutput(ostream &output) const
{
Houtan Bastani's avatar
Houtan Bastani committed
397
  output << "{\"statementName\": \"histval\", \"vals\": [";
398
  for (auto it = hist_values.begin();
399
400
       it != hist_values.end(); it++)
    {
401
402
      if (it != hist_values.begin())
        output << ", ";
403
404
405
      output << "{ \"name\": \"" << symbol_table.getName(it->first.first) << "\""
             << ", \"lag\": " << it->first.second
             << ", \"value\": \"";
406
      it->second->writeJsonOutput(output, {}, {});
407
408
409
410
411
      output << "\"}";
    }
  output << "]}";
}

412
InitvalFileStatement::InitvalFileStatement(string filename_arg) :
413
  filename{move(filename_arg)}
michel's avatar
michel committed
414
415
416
{
}

sebastien's avatar
sebastien committed
417
void
418
InitvalFileStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
michel's avatar
michel committed
419
{
sebastien's avatar
sebastien committed
420
421
422
423
424
  output << "%" << endl
         << "% INITVAL_FILE statement" << endl
         << "%" << endl
         << "options_.initval_file = 1;" << endl
         << "initvalf('" << filename << "');" << endl;
michel's avatar
michel committed
425
426
}

427
428
429
430
431
432
433
434
void
InitvalFileStatement::writeJsonOutput(ostream &output) const
{
  output << "{\"statementName\": \"init_val_file\""
         << ", \"filename\": \"" << filename << "\""
         << "}";
}

435
HistvalFileStatement::HistvalFileStatement(string filename_arg) :
436
  filename{move(filename_arg)}
437
438
439
440
{
}

void
441
HistvalFileStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
442
443
444
445
{
  output << "histvalf('" << filename << "');" << endl;
}

446
447
448
449
450
451
452
453
void
HistvalFileStatement::writeJsonOutput(ostream &output) const
{
  output << "{\"statementName\": \"hist_val_file\""
         << ", \"filename\": \"" << filename << "\""
         << "}";
}

454
HomotopyStatement::HomotopyStatement(homotopy_values_t homotopy_values_arg,
455
                                     const SymbolTable &symbol_table_arg) :
456
457
  homotopy_values{move(homotopy_values_arg)},
  symbol_table{symbol_table_arg}
458
459
460
461
{
}

void
462
HomotopyStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
463
{
sebastien's avatar
sebastien committed
464
465
466
  output << "%" << endl
         << "% HOMOTOPY_SETUP instructions" << endl
         << "%" << endl
467
         << "options_.homotopy_values = [];" << endl;
468

469
  for (const auto & homotopy_value : homotopy_values)
470
    {
471
472
473
      int symb_id;
      expr_t expression1, expression2;
      tie(symb_id, expression1, expression2) = homotopy_value;
474

sebastien's avatar
sebastien committed
475
476
      const SymbolType type = symbol_table.getType(symb_id);
      const int tsid = symbol_table.getTypeSpecificID(symb_id) + 1;
477

478
      output << "options_.homotopy_values = vertcat(options_.homotopy_values, [ " << static_cast<int>(type) << ", " << tsid << ", ";
479
      if (expression1)
480
481
482
        expression1->writeOutput(output);
      else
        output << "NaN";
483
484
      output << ", ";
      expression2->writeOutput(output);
485
      output << "]);" << endl;
486
487
    }
}
sebastien's avatar
sebastien committed
488

489
490
491
492
493
void
HomotopyStatement::writeJsonOutput(ostream &output) const
{
  output << "{\"statementName\": \"homotopy\", "
         << "\"values\": [";
494
  for (auto it = homotopy_values.begin();
495
       it != homotopy_values.end(); ++it)
496
497
498
    {
      if (it != homotopy_values.begin())
        output << ", ";
499
500
501
502
503
504

      int symb_id;
      expr_t expression1, expression2;
      tie(symb_id, expression1, expression2) = *it;

      output << "{\"name\": \"" << symbol_table.getName(symb_id) << "\""
505
             << ", \"initial_value\": \"";
506
507
      if (expression1)
        expression1->writeJsonOutput(output, {}, {});
508
509
510
      else
        output << "NaN";
      output << "\", \"final_value\": \"";
511
      expression2->writeJsonOutput(output, {}, {});
512
513
514
515
516
517
      output << "\"}";
    }
  output << "]"
         << "}";
}

518
SaveParamsAndSteadyStateStatement::SaveParamsAndSteadyStateStatement(string filename_arg) :
519
  filename{move(filename_arg)}
sebastien's avatar
sebastien committed
520
521
522
523
{
}

void
524
SaveParamsAndSteadyStateStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
sebastien's avatar
sebastien committed
525
526
527
528
{
  output << "save_params_and_steady_state('" << filename << "');" << endl;
}

529
530
531
532
533
534
535
536
void
SaveParamsAndSteadyStateStatement::writeJsonOutput(ostream &output) const
{
  output << "{\"statementName\": \"save_params_and_steady_state\""
         << ", \"filename\": \"" << filename << "\""
         << "}";
}

537
LoadParamsAndSteadyStateStatement::LoadParamsAndSteadyStateStatement(const string &filename,
538
539
                                                                     const SymbolTable &symbol_table_arg,
                                                                     WarningConsolidation &warnings) :
540
  symbol_table{symbol_table_arg}
sebastien's avatar
sebastien committed
541
{
542
  cout << "Reading " << filename << "." << endl;
sebastien's avatar
sebastien committed
543

544
  ifstream f;
545
  f.open(filename, ios::in);
546
  if (f.fail())
547
548
549
550
    {
      cerr << "ERROR: Can't open " << filename << endl;
      exit(EXIT_FAILURE);
    }
sebastien's avatar
sebastien committed
551

552
  while (true)
sebastien's avatar
sebastien committed
553
    {
554
555
556
557
      string symb_name, value;
      f >> symb_name >> value;
      if (f.eof())
        break;
sebastien's avatar
sebastien committed
558
559
560

      try
        {
561
562
          int symb_id = symbol_table.getID(symb_name);
          content[symb_id] = value;
sebastien's avatar
sebastien committed
563
        }
564
      catch (SymbolTable::UnknownSymbolNameException &e)
sebastien's avatar
sebastien committed
565
        {
566
          warnings << "WARNING: Unknown symbol " << symb_name << " in " << filename << endl;
567
568
        }
    }
569
  f.close();
570
571
572
}

void
573
LoadParamsAndSteadyStateStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
574
{
575
  for (const auto & it : content)
576
    {
577
      switch (symbol_table.getType(it.first))
578
        {
579
        case SymbolType::parameter:
580
581
          output << "M_.params";
          break;
582
        case SymbolType::endogenous:
583
584
          output << "oo_.steady_state";
          break;
585
        case SymbolType::exogenous:
586
587
          output << "oo_.exo_steady_state";
          break;
588
        case SymbolType::exogenousDet:
589
590
591
          output << "oo_.exo_det_steady_state";
          break;
        default:
592
          cerr << "ERROR: Unsupported variable type for " << symbol_table.getName(it.first) << " in load_params_and_steady_state" << endl;
593
          exit(EXIT_FAILURE);
sebastien's avatar
sebastien committed
594
        }
595

596
597
      int tsid = symbol_table.getTypeSpecificID(it.first) + 1;
      output << "(" << tsid << ") = " << it.second << ";" << endl;
sebastien's avatar
sebastien committed
598
599
    }
}
600

601
602
603
604
605
void
LoadParamsAndSteadyStateStatement::writeJsonOutput(ostream &output) const
{
  output << "{\"statementName\": \"load_params_and_steady_state\""
         << "\"values\": [";
606
  for (auto it = content.begin();
607
608
609
610
611
612
613
614
615
616
617
       it != content.end(); it++)
    {
      if (it != content.begin())
        output << ", ";
      output << "{\"name\": \"" << symbol_table.getName(it->first) << "\""
             << ", \"value\": \"" << it->second << "\"}";
    }
  output << "]"
         << "}";
}

618
void
619
LoadParamsAndSteadyStateStatement::fillEvalContext(eval_context_t &eval_context) const
620
{
621
  for (const auto & it : content)
622
    eval_context[it.first] = stod(it.second);
623
}