Shocks.cc 7.95 KB
Newer Older
sebastien's avatar
sebastien committed
1
/*
sebastien's avatar
sebastien committed
2
 * Copyright (C) 2003-2009 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
}

sebastien's avatar
sebastien committed
77
void
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
110
111
112
113
114
115
116
117
AbstractShocksStatement::checkPass(ModFileStructure &mod_file_struct)
{
  mod_file_struct.shocks_present = true;
}

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

void
ShocksStatement::writeVarAndStdShocks(ostream &output) const
118
{
sebastien's avatar
sebastien committed
119
  var_and_std_shocks_type::const_iterator it;
120

121
  for (it = var_shocks.begin(); it != var_shocks.end(); it++)
122
    {
sebastien's avatar
sebastien committed
123
      int id = symbol_table.getTypeSpecificID(it->first) + 1;
sebastien's avatar
sebastien committed
124
125
126
127
      const NodeID value = it->second;
      output << "M_.Sigma_e(" << id << ", " << id << ") = ";
      value->writeOutput(output);
      output << ";" << endl;
128
    }
sebastien's avatar
sebastien committed
129

130
  for (it = std_shocks.begin(); it != std_shocks.end(); it++)
131
    {
sebastien's avatar
sebastien committed
132
      int id = symbol_table.getTypeSpecificID(it->first) + 1;
sebastien's avatar
sebastien committed
133
134
135
136
      const NodeID value = it->second;
      output << "M_.Sigma_e(" << id << ", " << id << ") = (";
      value->writeOutput(output);
      output << ")^2;" << endl;
137
138
139
    }
}

sebastien's avatar
sebastien committed
140
void
141
ShocksStatement::writeCovarAndCorrShocks(ostream &output) const
142
{
sebastien's avatar
sebastien committed
143
144
  covar_and_corr_shocks_type::const_iterator it;

145
  for (it = covar_shocks.begin(); it != covar_shocks.end(); it++)
146
    {
sebastien's avatar
sebastien committed
147
148
      int id1 = symbol_table.getTypeSpecificID(it->first.first) + 1;
      int id2 = symbol_table.getTypeSpecificID(it->first.second) + 1;
sebastien's avatar
sebastien committed
149
150
151
152
      const NodeID value = it->second;
      output << "M_.Sigma_e(" << id1 << ", " << id2 << ") = ";
      value->writeOutput(output);
      output << "; M_.Sigma_e(" << id2 << ", " << id1 << ") = M_.Sigma_e("
sebastien's avatar
sebastien committed
153
             << id1 << ", " << id2 << ");\n";
154
    }
sebastien's avatar
sebastien committed
155

156
  for (it = corr_shocks.begin(); it != corr_shocks.end(); it++)
157
    {
sebastien's avatar
sebastien committed
158
159
      int id1 = symbol_table.getTypeSpecificID(it->first.first) + 1;
      int id2 = symbol_table.getTypeSpecificID(it->first.second) + 1;
sebastien's avatar
sebastien committed
160
161
162
163
      const NodeID value = it->second;
      output << "M_.Sigma_e(" << id1 << ", " << id2 << ") = ";
      value->writeOutput(output);
      output << "*sqrt(M_.Sigma_e(" << id1 << ", " << id1 << ")*M_.Sigma_e("
164
             << id2 << ", " << id2 << ")); M_.Sigma_e(" << id2 << ", "
sebastien's avatar
sebastien committed
165
             << id1 << ") = M_.Sigma_e(" << id1 << ", " << id2 << ");\n";
166
167
168
    }
}

sebastien's avatar
sebastien committed
169
170
MShocksStatement::MShocksStatement(const det_shocks_type &det_shocks_arg,
                                   const SymbolTable &symbol_table_arg) :
171
  AbstractShocksStatement(true, det_shocks_arg, symbol_table_arg)
172
173
174
{
}

sebastien's avatar
sebastien committed
175
void
sebastien's avatar
sebastien committed
176
MShocksStatement::writeOutput(ostream &output, const string &basename) const
177
{
sebastien's avatar
sebastien committed
178
  output << "%" << endl
179
         << "% MSHOCKS instructions" << endl
sebastien's avatar
sebastien committed
180
         << "%" << endl;
181

sebastien's avatar
sebastien committed
182
  // Write instruction that initializes a shock
sebastien's avatar
sebastien committed
183
  output << "make_ex_;" << endl;
184

sebastien's avatar
sebastien committed
185
  writeDetShocks(output);
186
}
187
188
189
190
191
192
193
194
195
196
197

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)
{
198
199
  for (AbstractShocksStatement::det_shocks_type::const_iterator it = paths.begin();
       it != paths.end(); it++)
200
201
202
    {
      int this_path_length = 0;
      const vector<AbstractShocksStatement::DetShockElement> &elems = it->second;
203
      for (int i = 0; i < (int) elems.size(); i++)
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
        // 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;

225
226
  for (AbstractShocksStatement::det_shocks_type::const_iterator it = paths.begin();
       it != paths.end(); it++)
227
228
229
    {
      output << "constrained_vars_ = strvcat(constrained_vars_, '" << it->first << "');" << endl;
      const vector<AbstractShocksStatement::DetShockElement> &elems = it->second;
230
231
      for (int i = 0; i < (int) elems.size(); i++)
        for (int j = elems[i].period1; j <= elems[i].period2; j++)
232
233
234
235
236
237
238
239
          {
            output << "constrained_paths_(" << k << "," << j << ")=";
            elems[i].value->writeOutput(output);
            output << ";" << endl;
          }
      k++;
    }
}