Shocks.cc 17.1 KB
Newer Older
1
/*
Houtan Bastani's avatar
Houtan Bastani committed
2
 * Copyright (C) 2003-2016 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 <cassert>
#include <cstdlib>
22
23
24
25
26
#include <iostream>

#include "Shocks.hh"

AbstractShocksStatement::AbstractShocksStatement(bool mshocks_arg,
27
                                                 bool overwrite_arg,
28
                                                 const det_shocks_t &det_shocks_arg,
29
                                                 const SymbolTable &symbol_table_arg) :
30
  mshocks(mshocks_arg),
31
  overwrite(overwrite_arg),
32
33
34
35
36
37
38
39
40
41
  det_shocks(det_shocks_arg),
  symbol_table(symbol_table_arg)
{
}

void
AbstractShocksStatement::writeDetShocks(ostream &output) const
{
  int exo_det_length = 0;

42
  for (det_shocks_t::const_iterator it = det_shocks.begin();
43
       it != det_shocks.end(); it++)
44
    {
sebastien's avatar
sebastien committed
45
      int id = symbol_table.getTypeSpecificID(it->first) + 1;
46
47
      bool exo_det = (symbol_table.getType(it->first) == eExogenousDet);

48
      for (size_t i = 0; i < it->second.size(); i++)
49
        {
50
51
52
          const int &period1 = it->second[i].period1;
          const int &period2 = it->second[i].period2;
          const expr_t value = it->second[i].value;
53

54
55
56
57
58
59
          output << "M_.det_shocks = [ M_.det_shocks;" << endl
                 << "struct('exo_det'," << (int) exo_det
                 << ",'exo_id'," << id
                 << ",'multiplicative'," << (int) mshocks
                 << ",'periods'," << period1 << ":" << period2
                 << ",'value',";
60
          value->writeOutput(output);
61
          output << ") ];" << endl;
62
63
64
65
66
67
68
69

          if (exo_det && (period2 > exo_det_length))
            exo_det_length = period2;
        }
    }
  output << "M_.exo_det_length = " << exo_det_length << ";\n";
}

70
71
ShocksStatement::ShocksStatement(bool overwrite_arg,
                                 const det_shocks_t &det_shocks_arg,
72
73
74
75
                                 const var_and_std_shocks_t &var_shocks_arg,
                                 const var_and_std_shocks_t &std_shocks_arg,
                                 const covar_and_corr_shocks_t &covar_shocks_arg,
                                 const covar_and_corr_shocks_t &corr_shocks_arg,
76
                                 const SymbolTable &symbol_table_arg) :
77
  AbstractShocksStatement(false, overwrite_arg, det_shocks_arg, symbol_table_arg),
78
79
80
81
82
83
84
85
  var_shocks(var_shocks_arg),
  std_shocks(std_shocks_arg),
  covar_shocks(covar_shocks_arg),
  corr_shocks(corr_shocks_arg)
{
}

void
86
ShocksStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
87
88
89
90
91
{
  output << "%" << endl
         << "% SHOCKS instructions" << endl
         << "%" << endl;

92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
  if (overwrite)
    {
      output << "M_.det_shocks = [];" << endl;

      output << "M_.Sigma_e = zeros(" << symbol_table.exo_nbr() << ", "
              << symbol_table.exo_nbr() << ");" << endl
              << "M_.Correlation_matrix = eye(" << symbol_table.exo_nbr() << ", "
              << symbol_table.exo_nbr() << ");" << endl;

      if (has_calibrated_measurement_errors())
        output << "M_.H = zeros(" << symbol_table.observedVariablesNbr() << ", "
               << symbol_table.observedVariablesNbr() << ");" << endl
               << "M_.Correlation_matrix_ME = eye(" << symbol_table.observedVariablesNbr() << ", "
               << symbol_table.observedVariablesNbr() << ");" << endl;
      else
        output << "M_.H = 0;" << endl
               << "M_.Correlation_matrix_ME = 1;" << endl;

    }

112
113
114
  writeDetShocks(output);
  writeVarAndStdShocks(output);
  writeCovarAndCorrShocks(output);
115
116
117
118
119

  /* M_.sigma_e_is_diagonal is initialized to 1 by ModFile.cc.
     If there are no off-diagonal elements, and we are not in overwrite mode,
     then we don't reset it to 1, since there might be previous shocks blocks
     with off-diagonal elements. */
120
121
  if (covar_shocks.size()+corr_shocks.size() > 0)
    output << "M_.sigma_e_is_diagonal = 0;" << endl;
122
  else if (overwrite)
123
124
125
    output << "M_.sigma_e_is_diagonal = 1;" << endl;
}

Sébastien Villemot's avatar
Sébastien Villemot committed
126
void
127
ShocksStatement::writeVarOrStdShock(ostream &output, var_and_std_shocks_t::const_iterator &it,
128
                                    bool stddev) const
Sébastien Villemot's avatar
Sébastien Villemot committed
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
{
  SymbolType type = symbol_table.getType(it->first);
  assert(type == eExogenous || symbol_table.isObservedVariable(it->first));

  int id;
  if (type == eExogenous)
    {
      output << "M_.Sigma_e(";
      id = symbol_table.getTypeSpecificID(it->first) + 1;
    }
  else
    {
      output << "M_.H(";
      id = symbol_table.getObservedVariableIndex(it->first) + 1;
    }

  output << id << ", " << id << ") = ";
146
  if (stddev)
Sébastien Villemot's avatar
Sébastien Villemot committed
147
148
    output << "(";
  it->second->writeOutput(output);
149
  if (stddev)
Sébastien Villemot's avatar
Sébastien Villemot committed
150
151
152
153
    output << ")^2";
  output << ";" << endl;
}

154
155
void
ShocksStatement::writeVarAndStdShocks(ostream &output) const
156
{
157
  var_and_std_shocks_t::const_iterator it;
158

159
  for (it = var_shocks.begin(); it != var_shocks.end(); it++)
Sébastien Villemot's avatar
Sébastien Villemot committed
160
    writeVarOrStdShock(output, it, false);
161

162
  for (it = std_shocks.begin(); it != std_shocks.end(); it++)
Sébastien Villemot's avatar
Sébastien Villemot committed
163
164
165
166
    writeVarOrStdShock(output, it, true);
}

void
167
ShocksStatement::writeCovarOrCorrShock(ostream &output, covar_and_corr_shocks_t::const_iterator &it,
Sébastien Villemot's avatar
Sébastien Villemot committed
168
169
170
171
172
173
                                       bool corr) const
{
  SymbolType type1 = symbol_table.getType(it->first.first);
  SymbolType type2 = symbol_table.getType(it->first.second);
  assert((type1 == eExogenous && type2 == eExogenous)
         || (symbol_table.isObservedVariable(it->first.first) && symbol_table.isObservedVariable(it->first.second)));
174
  string matrix, corr_matrix;
Sébastien Villemot's avatar
Sébastien Villemot committed
175
176
  int id1, id2;
  if (type1 == eExogenous)
177
    {
Sébastien Villemot's avatar
Sébastien Villemot committed
178
      matrix = "M_.Sigma_e";
179
      corr_matrix = "M_.Correlation_matrix";
Sébastien Villemot's avatar
Sébastien Villemot committed
180
181
182
183
184
185
      id1 = symbol_table.getTypeSpecificID(it->first.first) + 1;
      id2 = symbol_table.getTypeSpecificID(it->first.second) + 1;
    }
  else
    {
      matrix = "M_.H";
186
      corr_matrix = "M_.Correlation_matrix_ME";
Sébastien Villemot's avatar
Sébastien Villemot committed
187
188
      id1 = symbol_table.getObservedVariableIndex(it->first.first) + 1;
      id2 = symbol_table.getObservedVariableIndex(it->first.second) + 1;
189
    }
Sébastien Villemot's avatar
Sébastien Villemot committed
190
191
192
193
194
195
196
197
198

  output << matrix << "(" << id1 << ", " << id2 << ") = ";
  it->second->writeOutput(output);
  if (corr)
    output << "*sqrt(" << matrix << "(" << id1 << ", " << id1 << ")*"
           << matrix << "(" << id2 << ", " << id2 << "))";
  output << ";" << endl
         << matrix << "(" << id2 << ", " << id1 << ") = "
         << matrix << "(" << id1 << ", " << id2 << ");" << endl;
199
200
201
202
203
204
205
206
207

  if (corr)
    {
      output << corr_matrix << "(" << id1 << ", " << id2 << ") = ";
      it->second->writeOutput(output);
      output << ";" << endl
             << corr_matrix << "(" << id2 << ", " << id1 << ") = "
             << corr_matrix << "(" << id1 << ", " << id2 << ");" << endl;
    }
208
209
210
}

void
211
ShocksStatement::writeCovarAndCorrShocks(ostream &output) const
212
{
213
  covar_and_corr_shocks_t::const_iterator it;
214

215
  for (it = covar_shocks.begin(); it != covar_shocks.end(); it++)
Sébastien Villemot's avatar
Sébastien Villemot committed
216
    writeCovarOrCorrShock(output, it, false);
217

218
  for (it = corr_shocks.begin(); it != corr_shocks.end(); it++)
Sébastien Villemot's avatar
Sébastien Villemot committed
219
220
221
222
    writeCovarOrCorrShock(output, it, true);
}

void
223
ShocksStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings)
Sébastien Villemot's avatar
Sébastien Villemot committed
224
{
225
226
227
  /* Error out if variables are not of the right type. This must be done here
     and not at parsing time (see #448).
     Also Determine if there is a calibrated measurement error */
228
  for (var_and_std_shocks_t::const_iterator it = var_shocks.begin();
Sébastien Villemot's avatar
Sébastien Villemot committed
229
       it != var_shocks.end(); it++)
230
231
232
233
234
235
236
237
238
    {
      if (symbol_table.getType(it->first) != eExogenous
          && !symbol_table.isObservedVariable(it->first))
        {
          cerr << "shocks: setting a variance on '"
               << symbol_table.getName(it->first) << "' is not allowed, because it is neither an exogenous variable nor an observed endogenous variable" << endl;
          exit(EXIT_FAILURE);
        }
    }
Sébastien Villemot's avatar
Sébastien Villemot committed
239

240
  for (var_and_std_shocks_t::const_iterator it = std_shocks.begin();
Sébastien Villemot's avatar
Sébastien Villemot committed
241
       it != std_shocks.end(); it++)
242
243
244
245
246
247
248
249
250
    {
      if (symbol_table.getType(it->first) != eExogenous
          && !symbol_table.isObservedVariable(it->first))
        {
          cerr << "shocks: setting a standard error on '"
               << symbol_table.getName(it->first) << "' is not allowed, because it is neither an exogenous variable nor an observed endogenous variable" << endl;
          exit(EXIT_FAILURE);
        }
    }
Sébastien Villemot's avatar
Sébastien Villemot committed
251

252
  for (covar_and_corr_shocks_t::const_iterator it = covar_shocks.begin();
Sébastien Villemot's avatar
Sébastien Villemot committed
253
       it != covar_shocks.end(); it++)
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
    {
      int symb_id1 = it->first.first;
      int symb_id2 = it->first.second;
      
      if (!((symbol_table.getType(symb_id1) == eExogenous
             && symbol_table.getType(symb_id2) == eExogenous)
            || (symbol_table.isObservedVariable(symb_id1)
                && symbol_table.isObservedVariable(symb_id2))))
        {
          cerr << "shocks: setting a covariance between '"
               << symbol_table.getName(symb_id1) << "' and '"
               << symbol_table.getName(symb_id2) << "'is not allowed; covariances can only be specified for exogenous or observed endogenous variables of same type" << endl;
          exit(EXIT_FAILURE);
        }
    }
Sébastien Villemot's avatar
Sébastien Villemot committed
269

270
  for (covar_and_corr_shocks_t::const_iterator it = corr_shocks.begin();
Sébastien Villemot's avatar
Sébastien Villemot committed
271
       it != corr_shocks.end(); it++)
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
    {
      int symb_id1 = it->first.first;
      int symb_id2 = it->first.second;
      
      if (!((symbol_table.getType(symb_id1) == eExogenous
             && symbol_table.getType(symb_id2) == eExogenous)
            || (symbol_table.isObservedVariable(symb_id1)
                && symbol_table.isObservedVariable(symb_id2))))
        {
          cerr << "shocks: setting a correlation between '"
               << symbol_table.getName(symb_id1) << "' and '"
               << symbol_table.getName(symb_id2) << "'is not allowed; correlations can only be specified for exogenous or observed endogenous variables of same type" << endl;
          exit(EXIT_FAILURE);
        }
    }
287

288
289
290
  // Determine if there is a calibrated measurement error
  mod_file_struct.calibrated_measurement_errors |= has_calibrated_measurement_errors();

291
292
293
  // Fill in mod_file_struct.parameters_with_shocks_values (related to #469)
  for (var_and_std_shocks_t::const_iterator it = var_shocks.begin();
       it != var_shocks.end(); ++it)
294
    it->second->collectVariables(eParameter, mod_file_struct.parameters_within_shocks_values);
295
296
  for (var_and_std_shocks_t::const_iterator it = std_shocks.begin();
       it != std_shocks.end(); ++it)
297
    it->second->collectVariables(eParameter, mod_file_struct.parameters_within_shocks_values);
298
299
  for (covar_and_corr_shocks_t::const_iterator it = covar_shocks.begin();
       it != covar_shocks.end(); ++it)
300
    it->second->collectVariables(eParameter, mod_file_struct.parameters_within_shocks_values);
301
302
  for (covar_and_corr_shocks_t::const_iterator it = corr_shocks.begin();
       it != corr_shocks.end(); ++it)
303
    it->second->collectVariables(eParameter, mod_file_struct.parameters_within_shocks_values);
304

305
306
}

307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
bool
ShocksStatement::has_calibrated_measurement_errors() const
{
  for (var_and_std_shocks_t::const_iterator it = var_shocks.begin();
       it != var_shocks.end(); it++)
    if (symbol_table.isObservedVariable(it->first))
      return true;

  for (var_and_std_shocks_t::const_iterator it = std_shocks.begin();
       it != std_shocks.end(); it++)
    if (symbol_table.isObservedVariable(it->first))
      return true;

  for (covar_and_corr_shocks_t::const_iterator it = covar_shocks.begin();
       it != covar_shocks.end(); it++)
    if (symbol_table.isObservedVariable(it->first.first)
        || symbol_table.isObservedVariable(it->first.second))
      return true;

  for (covar_and_corr_shocks_t::const_iterator it = corr_shocks.begin();
       it != corr_shocks.end(); it++)
    if (symbol_table.isObservedVariable(it->first.first)
        || symbol_table.isObservedVariable(it->first.second))
      return true;

  return false;
}

MShocksStatement::MShocksStatement(bool overwrite_arg,
                                   const det_shocks_t &det_shocks_arg,
337
                                   const SymbolTable &symbol_table_arg) :
338
  AbstractShocksStatement(true, overwrite_arg, det_shocks_arg, symbol_table_arg)
339
340
341
342
{
}

void
343
MShocksStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
344
{
sebastien's avatar
sebastien committed
345
  output << "%" << endl
346
         << "% MSHOCKS instructions" << endl
sebastien's avatar
sebastien committed
347
         << "%" << endl;
348

349
350
351
  if (overwrite)
    output << "M_.det_shocks = [];" << endl;

352
353
  writeDetShocks(output);
}
354

355
ConditionalForecastPathsStatement::ConditionalForecastPathsStatement(const AbstractShocksStatement::det_shocks_t &paths_arg) :
356
357
358
359
360
361
  paths(paths_arg),
  path_length(-1)
{
}

void
362
ConditionalForecastPathsStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings)
363
{
364
  for (AbstractShocksStatement::det_shocks_t::const_iterator it = paths.begin();
365
       it != paths.end(); it++)
366
367
    {
      int this_path_length = 0;
368
      const vector<AbstractShocksStatement::DetShockElement> &elems = it->second;
369
      for (int i = 0; i < (int) elems.size(); i++)
370
371
372
373
374
375
376
377
378
379
380
381
382
        // Period1 < Period2, as enforced in ParsingDriver::add_period()
        this_path_length = max(this_path_length, elems[i].period2);
      if (path_length == -1)
        path_length = this_path_length;
      else if (path_length != this_path_length)
        {
          cerr << "conditional_forecast_paths: all constrained paths must have the same length!" << endl;
          exit(EXIT_FAILURE);
        }
    }
}

void
383
ConditionalForecastPathsStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
384
385
386
387
388
389
390
{
  assert(path_length > 0);
  output << "constrained_vars_ = [];" << endl
         << "constrained_paths_ = zeros(" << paths.size() << ", " << path_length << ");" << endl;

  int k = 1;

391
  for (AbstractShocksStatement::det_shocks_t::const_iterator it = paths.begin();
392
       it != paths.end(); it++)
393
    {
394
      if (it == paths.begin())
395
396
397
        {
          output << "constrained_vars_ = " << it->first +1 << ";" << endl;
        }
398
      else
399
400
401
402
        {
          output << "constrained_vars_ = [constrained_vars_; " << it->first +1 << "];" << endl;
        }

403

404
      const vector<AbstractShocksStatement::DetShockElement> &elems = it->second;
405
406
      for (int i = 0; i < (int) elems.size(); i++)
        for (int j = elems[i].period1; j <= elems[i].period2; j++)
407
408
409
410
411
412
413
414
          {
            output << "constrained_paths_(" << k << "," << j << ")=";
            elems[i].value->writeOutput(output);
            output << ";" << endl;
          }
      k++;
    }
}
415
416
417
418
419
420
421
422

MomentCalibration::MomentCalibration(const constraints_t &constraints_arg,
                                     const SymbolTable &symbol_table_arg)
  : constraints(constraints_arg), symbol_table(symbol_table_arg)
{
}

void
423
MomentCalibration::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
424
{
425
  output << "options_.endogenous_prior_restrictions.moment = {" << endl;
426
427
428
429
430
  for (size_t i = 0; i < constraints.size(); i++)
    {
      const Constraint &c = constraints[i];
      output << "'" << symbol_table.getName(c.endo1) << "', "
             << "'" << symbol_table.getName(c.endo2) << "', "
431
             << c.lags << ", "
432
433
434
435
436
437
438
             << "[ " << c.lower_bound << ", " << c.upper_bound << " ];"
             << endl;
    }
  output << "};" << endl;
}

IrfCalibration::IrfCalibration(const constraints_t &constraints_arg,
439
440
441
                               const SymbolTable &symbol_table_arg,
                               const OptionsList &options_list_arg)
  : constraints(constraints_arg), symbol_table(symbol_table_arg), options_list(options_list_arg)
442
443
444
445
{
}

void
446
IrfCalibration::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
447
{
448
449
  options_list.writeOutput(output);

450
  output << "options_.endogenous_prior_restrictions.irf = {" << endl;
451
452
453
454
455
  for (size_t i = 0; i < constraints.size(); i++)
    {
      const Constraint &c = constraints[i];
      output << "'" << symbol_table.getName(c.endo) << "', "
             << "'" << symbol_table.getName(c.exo) << "', "
456
             << c.periods << ", "
457
458
459
460
461
             << "[ " << c.lower_bound << ", " << c.upper_bound << " ];"
             << endl;
    }
  output << "};" << endl;
}
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480

ShockGroupsStatement::ShockGroupsStatement(const group_t &shock_groups_arg, const string &name_arg)
  : shock_groups(shock_groups_arg), name(name_arg)
{
}

void
ShockGroupsStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
{
  for (vector<Group>::const_iterator it = shock_groups.begin(); it != shock_groups.end(); it++)
    {
      output << "M_.shock_groups." << name
             << "." << it->name << " = {";
      for ( vector<string>::const_iterator it1 = it->list.begin(); it1 != it->list.end(); it1++)
        output << " '" << *it1 << "'";
      output << "};" << endl;
    }
}