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;
    }
}