NumericalInitialization.cc 14 KB
Newer Older
1
/*
2
 * Copyright (C) 2003-2014 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

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

InitParamStatement::InitParamStatement(int symb_id_arg,
28
                                       const expr_t param_value_arg,
29
                                       const SymbolTable &symbol_table_arg) :
sebastien's avatar
sebastien committed
30
  symb_id(symb_id_arg),
31
32
33
34
35
  param_value(param_value_arg),
  symbol_table(symbol_table_arg)
{
}

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

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

53
54
55
56
57
58
59
60
61
62
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
63
void
64
InitParamStatement::fillEvalContext(eval_context_t &eval_context) const
65
{
sebastien's avatar
sebastien committed
66
67
68
69
  try
    {
      eval_context[symb_id] = param_value->eval(eval_context);
    }
70
  catch (ExprNode::EvalException &e)
sebastien's avatar
sebastien committed
71
72
73
    {
      // Do nothing
    }
74
75
}

76
InitOrEndValStatement::InitOrEndValStatement(const init_values_t &init_values_arg,
77
78
                                             const SymbolTable &symbol_table_arg,
                                             const bool &all_values_required_arg) :
79
  init_values(init_values_arg),
80
81
  symbol_table(symbol_table_arg),
  all_values_required(all_values_required_arg)
82
83
84
{
}

sebastien's avatar
sebastien committed
85
void
86
InitOrEndValStatement::fillEvalContext(eval_context_t &eval_context) const
sebastien's avatar
sebastien committed
87
{
88
  for (init_values_t::const_iterator it = init_values.begin();
89
       it != init_values.end(); it++)
sebastien's avatar
sebastien committed
90
91
92
93
94
    {
      try
        {
          eval_context[it->first] = (it->second)->eval(eval_context);
        }
95
      catch (ExprNode::EvalException &e)
sebastien's avatar
sebastien committed
96
97
98
99
100
101
        {
          // Do nothing
        }
    }
}

102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
set<int>
InitOrEndValStatement::getUninitializedVariables(SymbolType type)
{
  set<int> unused;
  if (!all_values_required)
    return unused;

  if (type == eEndogenous)
    unused = symbol_table.getEndogenous();
  else if (type == eExogenous)
    unused = symbol_table.getExogenous();
  else
    {
      cerr << "ERROR: Shouldn't arrive here." << endl;
      exit(EXIT_FAILURE);
    }

  set<int>::iterator sit;
  for (init_values_t::const_iterator it = init_values.begin();
       it != init_values.end(); it++)
    {
      sit = unused.find(it->first);
      if (sit != unused.end())
        unused.erase(sit);
    }
  return unused;
}

130
131
132
void
InitOrEndValStatement::writeInitValues(ostream &output) const
{
133
  for (init_values_t::const_iterator it = init_values.begin();
134
       it != init_values.end(); it++)
135
    {
sebastien's avatar
sebastien committed
136
      const int symb_id = it->first;
137
      const expr_t expression = it->second;
138

sebastien's avatar
sebastien committed
139
140
      SymbolType type = symbol_table.getType(symb_id);
      int tsid = symbol_table.getTypeSpecificID(symb_id) + 1;
141
142
143
144
145
146
147
148

      if (type == eEndogenous)
        output << "oo_.steady_state";
      else if (type == eExogenous)
        output << "oo_.exo_steady_state";
      else if (type == eExogenousDet)
        output << "oo_.exo_det_steady_state";

sebastien's avatar
sebastien committed
149
      output << "( " << tsid << " ) = ";
150
151
152
153
154
      expression->writeOutput(output);
      output << ";" << endl;
    }
}

155
InitValStatement::InitValStatement(const init_values_t &init_values_arg,
156
157
158
                                   const SymbolTable &symbol_table_arg,
                                   const bool &all_values_required_arg) :
  InitOrEndValStatement(init_values_arg, symbol_table_arg, all_values_required_arg)
159
160
161
{
}

162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
void
InitValStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings)
{
  set<int> exogs = getUninitializedVariables(eExogenous);
  set<int> endogs = getUninitializedVariables(eEndogenous);

  if (endogs.size() > 0)
    {
      cerr << "ERROR: You have not set the following endogenous variables in initval:";
      for (set<int>::const_iterator it = endogs.begin(); it != endogs.end(); it++)
        cerr << " " << symbol_table.getName(*it);
      cerr << endl;
    }

  if (exogs.size() > 0)
    {
      cerr << "ERROR: You have not set the following exogenous variables in initval:";
      for (set<int>::const_iterator it = exogs.begin(); it != exogs.end(); it++)
        cerr << " " << symbol_table.getName(*it);
      cerr << endl;
    }

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

188
189
190
void
InitValStatement::writeOutput(ostream &output, const string &basename) const
{
sebastien's avatar
sebastien committed
191
192
193
  output << "%" << endl
         << "% INITVAL instructions" << endl
         << "%" << endl;
194
  // Writing initval block to set initial values for variables
sebastien's avatar
trunk:    
sebastien committed
195
  output << "options_.initval_file = 0;" << endl;
196
197

  writeInitValues(output);
sebastien's avatar
sebastien committed
198
}
199

sebastien's avatar
sebastien committed
200
201
202
void
InitValStatement::writeOutputPostInit(ostream &output) const
{
203
  output << "if M_.exo_nbr > 0;" << endl
sebastien's avatar
sebastien committed
204
205
206
207
208
         << "\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;
209
210
}

211
EndValStatement::EndValStatement(const init_values_t &init_values_arg,
212
213
214
                                 const SymbolTable &symbol_table_arg,
                                 const bool &all_values_required_arg) :
  InitOrEndValStatement(init_values_arg, symbol_table_arg, all_values_required_arg)
215
216
217
{
}

218
void
219
EndValStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings)
220
{
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
  set<int> exogs = getUninitializedVariables(eExogenous);
  set<int> endogs = getUninitializedVariables(eEndogenous);

  if (endogs.size() > 0)
    {
      cerr << "ERROR: You have not set the following endogenous variables in endval:";
      for (set<int>::const_iterator it = endogs.begin(); it != endogs.end(); it++)
        cerr << " " << symbol_table.getName(*it);
      cerr << endl;
    }

  if (exogs.size() > 0)
    {
      cerr << "ERROR: You have not set the following exogenous variables in endval:";
      for (set<int>::const_iterator it = exogs.begin(); it != exogs.end(); it++)
        cerr << " " << symbol_table.getName(*it);
      cerr << endl;
    }

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

244
245
246
void
EndValStatement::writeOutput(ostream &output, const string &basename) const
{
sebastien's avatar
sebastien committed
247
248
249
  output << "%" << endl
         << "% ENDVAL instructions" << endl
         << "%" << endl;
250
  // Writing endval block to set terminal values for variables
sebastien's avatar
sebastien committed
251
  output << "ys0_= oo_.steady_state;" << endl
252
         << "ex0_ = oo_.exo_steady_state;" << endl;
253
254
255
256

  writeInitValues(output);
}

257
HistValStatement::HistValStatement(const hist_values_t &hist_values_arg,
258
259
260
261
262
263
                                   const SymbolTable &symbol_table_arg) :
  hist_values(hist_values_arg),
  symbol_table(symbol_table_arg)
{
}

264
void
265
HistValStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings)
266
267
268
269
{
  mod_file_struct.histval_present = true;
}

270
271
272
void
HistValStatement::writeOutput(ostream &output, const string &basename) const
{
sebastien's avatar
sebastien committed
273
274
  output << "%" << endl
         << "% HISTVAL instructions" << endl
275
         << "%" << endl
276
         << "M_.endo_histval = zeros(M_.endo_nbr,M_.maximum_endo_lag);" << endl;
277

278
  for (hist_values_t::const_iterator it = hist_values.begin();
279
       it != hist_values.end(); it++)
280
    {
281
282
      int symb_id = it->first.first;
      int lag = it->first.second;
283
      const expr_t expression = it->second;
284

sebastien's avatar
sebastien committed
285
      SymbolType type = symbol_table.getType(symb_id);
Sébastien Villemot's avatar
Sébastien Villemot committed
286

287
288
      // For a lag greater than 1 on endo, or for any exo, lookup for auxiliary variable
      if ((type == eEndogenous && lag < 0) || type == eExogenous)
Sébastien Villemot's avatar
Sébastien Villemot committed
289
290
291
        {
          try
            {
292
              // This function call must remain the 1st statement in this block
Sébastien Villemot's avatar
Sébastien Villemot committed
293
294
              symb_id = symbol_table.searchAuxiliaryVars(symb_id, lag);
              lag = 0;
295
              type = eEndogenous;
Sébastien Villemot's avatar
Sébastien Villemot committed
296
297
298
299
300
301
302
303
304
305
306
307
308
            }
          catch (SymbolTable::SearchFailedException &e)
            {
              if (type == eEndogenous)
                {
                  cerr << "HISTVAL: internal error of Dynare, please contact the developers";
                  exit(EXIT_FAILURE);
                }
              // We don't fail for exogenous, because they are not replaced by
              // auxiliary variables in deterministic mode.
            }
        }

sebastien's avatar
sebastien committed
309
      int tsid = symbol_table.getTypeSpecificID(symb_id) + 1;
310
311

      if (type == eEndogenous)
312
        output << "M_.endo_histval( " << tsid << ", M_.maximum_endo_lag + " << lag << ") = ";
313
      else if (type == eExogenous)
sebastien's avatar
sebastien committed
314
        output << "oo_.exo_simul( M_.maximum_lag + " << lag << ", " << tsid << " ) = ";
315
      else if (type != eExogenousDet)
sebastien's avatar
sebastien committed
316
        output << "oo_.exo_det_simul( M_.maximum_lag + " << lag  << ", " << tsid << " ) = ";
317
318
319
320
321
322

      expression->writeOutput(output);
      output << ";" << endl;
    }
}

sebastien's avatar
sebastien committed
323
324
InitvalFileStatement::InitvalFileStatement(const string &filename_arg) :
  filename(filename_arg)
michel's avatar
michel committed
325
326
327
{
}

sebastien's avatar
sebastien committed
328
329
void
InitvalFileStatement::writeOutput(ostream &output, const string &basename) const
michel's avatar
michel committed
330
{
sebastien's avatar
sebastien committed
331
332
333
334
335
  output << "%" << endl
         << "% INITVAL_FILE statement" << endl
         << "%" << endl
         << "options_.initval_file = 1;" << endl
         << "initvalf('" << filename << "');" << endl;
michel's avatar
michel committed
336
337
}

338
339
340
341
342
343
344
345
346
347
348
HistvalFileStatement::HistvalFileStatement(const string &filename_arg) :
  filename(filename_arg)
{
}

void
HistvalFileStatement::writeOutput(ostream &output, const string &basename) const
{
  output << "histvalf('" << filename << "');" << endl;
}

349
HomotopyStatement::HomotopyStatement(const homotopy_values_t &homotopy_values_arg,
350
                                     const SymbolTable &symbol_table_arg) :
351
352
353
354
355
356
357
358
  homotopy_values(homotopy_values_arg),
  symbol_table(symbol_table_arg)
{
}

void
HomotopyStatement::writeOutput(ostream &output, const string &basename) const
{
sebastien's avatar
sebastien committed
359
360
361
  output << "%" << endl
         << "% HOMOTOPY_SETUP instructions" << endl
         << "%" << endl
362
         << "options_.homotopy_values = [];" << endl;
363

364
  for (homotopy_values_t::const_iterator it = homotopy_values.begin();
365
       it != homotopy_values.end(); it++)
366
    {
sebastien's avatar
sebastien committed
367
      const int &symb_id = it->first;
368
369
      const expr_t expression1 = it->second.first;
      const expr_t expression2 = it->second.second;
370

sebastien's avatar
sebastien committed
371
372
      const SymbolType type = symbol_table.getType(symb_id);
      const int tsid = symbol_table.getTypeSpecificID(symb_id) + 1;
373

sebastien's avatar
sebastien committed
374
      output << "options_.homotopy_values = vertcat(options_.homotopy_values, [ " << type << ", " << tsid << ", ";
375
376
377
378
      if (expression1 != NULL)
        expression1->writeOutput(output);
      else
        output << "NaN";
379
380
      output << ", ";
      expression2->writeOutput(output);
381
      output << "]);" << endl;
382
383
    }
}
sebastien's avatar
sebastien committed
384
385
386
387
388
389
390
391
392
393
394
395

SaveParamsAndSteadyStateStatement::SaveParamsAndSteadyStateStatement(const string &filename_arg) :
  filename(filename_arg)
{
}

void
SaveParamsAndSteadyStateStatement::writeOutput(ostream &output, const string &basename) const
{
  output << "save_params_and_steady_state('" << filename << "');" << endl;
}

396
LoadParamsAndSteadyStateStatement::LoadParamsAndSteadyStateStatement(const string &filename,
397
398
                                                                     const SymbolTable &symbol_table_arg,
                                                                     WarningConsolidation &warnings) :
399
  symbol_table(symbol_table_arg)
sebastien's avatar
sebastien committed
400
{
401
  cout << "Reading " << filename << "." << endl;
sebastien's avatar
sebastien committed
402

403
404
  ifstream f;
  f.open(filename.c_str(), ios::in);
405
  if (f.fail())
406
407
408
409
    {
      cerr << "ERROR: Can't open " << filename << endl;
      exit(EXIT_FAILURE);
    }
sebastien's avatar
sebastien committed
410

411
  while (true)
sebastien's avatar
sebastien committed
412
    {
413
414
415
416
      string symb_name, value;
      f >> symb_name >> value;
      if (f.eof())
        break;
sebastien's avatar
sebastien committed
417
418
419

      try
        {
420
421
          int symb_id = symbol_table.getID(symb_name);
          content[symb_id] = value;
sebastien's avatar
sebastien committed
422
        }
423
      catch (SymbolTable::UnknownSymbolNameException &e)
sebastien's avatar
sebastien committed
424
        {
425
          warnings << "WARNING: Unknown symbol " << symb_name << " in " << filename << endl;
426
427
        }
    }
428
  f.close();
429
430
431
432
433
}

void
LoadParamsAndSteadyStateStatement::writeOutput(ostream &output, const string &basename) const
{
434
435
  for (map<int, string>::const_iterator it = content.begin();
       it != content.end(); it++)
436
    {
437
      switch (symbol_table.getType(it->first))
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
        {
        case eParameter:
          output << "M_.params";
          break;
        case eEndogenous:
          output << "oo_.steady_state";
          break;
        case eExogenous:
          output << "oo_.exo_steady_state";
          break;
        case eExogenousDet:
          output << "oo_.exo_det_steady_state";
          break;
        default:
          cerr << "ERROR: Unsupported variable type for " << symbol_table.getName(it->first) << " in load_params_and_steady_state" << endl;
          exit(EXIT_FAILURE);
sebastien's avatar
sebastien committed
454
        }
455
456
457

      int tsid = symbol_table.getTypeSpecificID(it->first) + 1;
      output << "(" << tsid << ") = " << it->second << ";" << endl;
sebastien's avatar
sebastien committed
458
459
    }
}
460
461

void
462
LoadParamsAndSteadyStateStatement::fillEvalContext(eval_context_t &eval_context) const
463
{
464
465
  for (map<int, string>::const_iterator it = content.begin();
       it != content.end(); it++)
466
467
    eval_context[it->first] = atof(it->second.c_str());
}