Shocks.cc 9.93 KB
Newer Older
sebastien's avatar
sebastien committed
1
/*
Sébastien Villemot's avatar
Sébastien Villemot committed
2
 * Copyright (C) 2003-2010 Dynare Team
sebastien's avatar
sebastien committed
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
using namespace std;
sebastien's avatar
sebastien committed
21

22
23
#include <cassert>
#include <cstdlib>
sebastien's avatar
sebastien committed
24
25
#include <iostream>

26
27
#include "Shocks.hh"

sebastien's avatar
sebastien committed
28
AbstractShocksStatement::AbstractShocksStatement(bool mshocks_arg,
29
30
                                                 const det_shocks_type &det_shocks_arg,
                                                 const SymbolTable &symbol_table_arg) :
sebastien's avatar
sebastien committed
31
32
33
  mshocks(mshocks_arg),
  det_shocks(det_shocks_arg),
  symbol_table(symbol_table_arg)
34
35
36
{
}

sebastien's avatar
sebastien committed
37
38
void
AbstractShocksStatement::writeDetShocks(ostream &output) const
39
{
sebastien's avatar
sebastien committed
40
  int exo_det_length = 0;
41

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

sebastien's avatar
sebastien committed
49
50
51
52
      for (unsigned int i = 0; i < it->second.size(); i++)
        {
          const int &period1 = it->second[i].period1;
          const int &period2 = it->second[i].period2;
sebastien's avatar
sebastien committed
53
          const NodeID value = it->second[i].value;
sebastien's avatar
sebastien committed
54
55

          if (period1 == period2)
sebastien's avatar
sebastien committed
56
57
58
59
60
61
            {
              output << "set_shocks(" << set_shocks_index << "," << period1
                     << ", " << id << ", ";
              value->writeOutput(output);
              output << ");" << endl;
            }
sebastien's avatar
sebastien committed
62
          else
sebastien's avatar
sebastien committed
63
64
65
66
67
68
            {
              output << "set_shocks(" << set_shocks_index << "," << period1
                     << ":" << period2 << ", " << id << ", ";
              value->writeOutput(output);
              output << ");" << endl;
            }
sebastien's avatar
sebastien committed
69
70
71
72
73
74

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

77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
ShocksStatement::ShocksStatement(const det_shocks_type &det_shocks_arg,
                                 const var_and_std_shocks_type &var_shocks_arg,
                                 const var_and_std_shocks_type &std_shocks_arg,
                                 const covar_and_corr_shocks_type &covar_shocks_arg,
                                 const covar_and_corr_shocks_type &corr_shocks_arg,
                                 const SymbolTable &symbol_table_arg) :
  AbstractShocksStatement(false, det_shocks_arg, symbol_table_arg),
  var_shocks(var_shocks_arg),
  std_shocks(std_shocks_arg),
  covar_shocks(covar_shocks_arg),
  corr_shocks(corr_shocks_arg)
{
}

void
ShocksStatement::writeOutput(ostream &output, const string &basename) const
{
  output << "%" << endl
         << "% SHOCKS instructions" << endl
         << "%" << endl;

  // Write instruction that initializes a shock
  output << "make_ex_;" << endl;

  writeDetShocks(output);
  writeVarAndStdShocks(output);
  writeCovarAndCorrShocks(output);
  if (covar_shocks.size()+corr_shocks.size() > 0)
    output << "M_.sigma_e_is_diagonal = 0;" << endl;
  else
    output << "M_.sigma_e_is_diagonal = 1;" << endl;
}

Sébastien Villemot's avatar
Sébastien Villemot committed
110
111
void
ShocksStatement::writeVarOrStdShock(ostream &output, var_and_std_shocks_type::const_iterator &it,
112
                                    bool stddev) const
Sébastien Villemot's avatar
Sébastien Villemot committed
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
{
  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 << ") = ";
130
  if (stddev)
Sébastien Villemot's avatar
Sébastien Villemot committed
131
132
    output << "(";
  it->second->writeOutput(output);
133
  if (stddev)
Sébastien Villemot's avatar
Sébastien Villemot committed
134
135
136
137
    output << ")^2";
  output << ";" << endl;
}

138
139
void
ShocksStatement::writeVarAndStdShocks(ostream &output) const
140
{
sebastien's avatar
sebastien committed
141
  var_and_std_shocks_type::const_iterator it;
142

143
  for (it = var_shocks.begin(); it != var_shocks.end(); it++)
Sébastien Villemot's avatar
Sébastien Villemot committed
144
    writeVarOrStdShock(output, it, false);
sebastien's avatar
sebastien committed
145

146
  for (it = std_shocks.begin(); it != std_shocks.end(); it++)
Sébastien Villemot's avatar
Sébastien Villemot committed
147
148
149
150
151
152
153
154
155
156
157
158
159
160
    writeVarOrStdShock(output, it, true);
}

void
ShocksStatement::writeCovarOrCorrShock(ostream &output, covar_and_corr_shocks_type::const_iterator &it,
                                       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)));
  string matrix;
  int id1, id2;
  if (type1 == eExogenous)
161
    {
Sébastien Villemot's avatar
Sébastien Villemot committed
162
163
164
165
166
167
168
169
170
      matrix = "M_.Sigma_e";
      id1 = symbol_table.getTypeSpecificID(it->first.first) + 1;
      id2 = symbol_table.getTypeSpecificID(it->first.second) + 1;
    }
  else
    {
      matrix = "M_.H";
      id1 = symbol_table.getObservedVariableIndex(it->first.first) + 1;
      id2 = symbol_table.getObservedVariableIndex(it->first.second) + 1;
171
    }
Sébastien Villemot's avatar
Sébastien Villemot committed
172
173
174
175
176
177
178
179
180

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

sebastien's avatar
sebastien committed
183
void
184
ShocksStatement::writeCovarAndCorrShocks(ostream &output) const
185
{
sebastien's avatar
sebastien committed
186
187
  covar_and_corr_shocks_type::const_iterator it;

188
  for (it = covar_shocks.begin(); it != covar_shocks.end(); it++)
Sébastien Villemot's avatar
Sébastien Villemot committed
189
    writeCovarOrCorrShock(output, it, false);
sebastien's avatar
sebastien committed
190

191
  for (it = corr_shocks.begin(); it != corr_shocks.end(); it++)
Sébastien Villemot's avatar
Sébastien Villemot committed
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
    writeCovarOrCorrShock(output, it, true);
}

void
ShocksStatement::checkPass(ModFileStructure &mod_file_struct)
{
  // Workaround for trac ticket #35
  mod_file_struct.shocks_present = true;

  // Determine if there is a calibrated measurement error
  for (var_and_std_shocks_type::const_iterator it = var_shocks.begin();
       it != var_shocks.end(); it++)
    if (symbol_table.isObservedVariable(it->first))
      mod_file_struct.calibrated_measurement_errors = true;

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

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

  for (covar_and_corr_shocks_type::const_iterator it = corr_shocks.begin();
       it != corr_shocks.end(); it++)
    if (symbol_table.isObservedVariable(it->first.first)
        || symbol_table.isObservedVariable(it->first.second))
      mod_file_struct.calibrated_measurement_errors = true;
223
224
}

sebastien's avatar
sebastien committed
225
226
MShocksStatement::MShocksStatement(const det_shocks_type &det_shocks_arg,
                                   const SymbolTable &symbol_table_arg) :
227
  AbstractShocksStatement(true, det_shocks_arg, symbol_table_arg)
228
229
230
{
}

sebastien's avatar
sebastien committed
231
void
sebastien's avatar
sebastien committed
232
MShocksStatement::writeOutput(ostream &output, const string &basename) const
233
{
sebastien's avatar
sebastien committed
234
  output << "%" << endl
235
         << "% MSHOCKS instructions" << endl
sebastien's avatar
sebastien committed
236
         << "%" << endl;
237

sebastien's avatar
sebastien committed
238
  // Write instruction that initializes a shock
sebastien's avatar
sebastien committed
239
  output << "make_ex_;" << endl;
240

sebastien's avatar
sebastien committed
241
  writeDetShocks(output);
242
}
243

Sébastien Villemot's avatar
Sébastien Villemot committed
244
245
246
247
248
249
250
void
MShocksStatement::checkPass(ModFileStructure &mod_file_struct)
{
  // Workaround for trac ticket #35
  mod_file_struct.shocks_present = true;
}

251
252
253
254
255
256
257
258
259
260
ConditionalForecastPathsStatement::ConditionalForecastPathsStatement(const AbstractShocksStatement::det_shocks_type &paths_arg, const SymbolTable &symbol_table_arg) :
  paths(paths_arg),
  symbol_table(symbol_table_arg),
  path_length(-1)
{
}

void
ConditionalForecastPathsStatement::checkPass(ModFileStructure &mod_file_struct)
{
261
262
  for (AbstractShocksStatement::det_shocks_type::const_iterator it = paths.begin();
       it != paths.end(); it++)
263
264
265
    {
      int this_path_length = 0;
      const vector<AbstractShocksStatement::DetShockElement> &elems = it->second;
266
      for (int i = 0; i < (int) elems.size(); i++)
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
        // 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
ConditionalForecastPathsStatement::writeOutput(ostream &output, const string &basename) const
{
  assert(path_length > 0);
  output << "constrained_vars_ = [];" << endl
         << "constrained_paths_ = zeros(" << paths.size() << ", " << path_length << ");" << endl;

  int k = 1;

288
289
  for (AbstractShocksStatement::det_shocks_type::const_iterator it = paths.begin();
       it != paths.end(); it++)
290
291
292
    {
      output << "constrained_vars_ = strvcat(constrained_vars_, '" << it->first << "');" << endl;
      const vector<AbstractShocksStatement::DetShockElement> &elems = it->second;
293
294
      for (int i = 0; i < (int) elems.size(); i++)
        for (int j = elems[i].period1; j <= elems[i].period2; j++)
295
296
297
298
299
300
301
302
          {
            output << "constrained_paths_(" << k << "," << j << ")=";
            elems[i].value->writeOutput(output);
            output << ";" << endl;
          }
      k++;
    }
}