1 /*
2 * Copyright © 2004-2011 Ondra Kamenik
3 * Copyright © 2019 Dynare Team
4 *
5 * This file is part of Dynare.
6 *
7 * Dynare is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * Dynare is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with Dynare. If not, see <http://www.gnu.org/licenses/>.
19 */
20
21 #include <sstream>
22 #include <fstream>
23
24 #include "dynare3.hh"
25 #include "dynare_exception.hh"
26 #include "planner_builder.hh"
27 #include "forw_subst_builder.hh"
28
29 #include "utils/cc/exception.hh"
30 #include "parser/cc/parser_exception.hh"
31 #include "parser/cc/atom_substitutions.hh"
32 #include "../tl/cc/tl_exception.hh"
33 #include "../kord/kord_exception.hh"
34
35 /**************************************************************************************/
36 /* DynareNameList class */
37 /**************************************************************************************/
38 std::vector<int>
selectIndices(const std::vector<std::string> & ns) const39 DynareNameList::selectIndices(const std::vector<std::string> &ns) const
40 {
41 std::vector<int> res;
42 for (const auto &n : ns)
43 {
44 int j = 0;
45 while (j < getNum() && n != getName(j))
46 j++;
47 if (j == getNum())
48 throw DynareException(__FILE__, __LINE__,
49 std::string("Couldn't find name for ") + n
50 +" in DynareNameList::selectIndices");
51 res.push_back(j);
52 }
53 return res;
54 }
55
56 /**************************************************************************************/
57 /* Dynare class */
58 /**************************************************************************************/
59
Dynare(const std::string & modname,int ord,double sstol,Journal & jr)60 Dynare::Dynare(const std::string &modname, int ord, double sstol, Journal &jr)
61 : journal(jr), md(1), ss_tol(sstol)
62 {
63 std::ifstream f{modname};
64 if (f.fail())
65 throw DynareException(__FILE__, __LINE__, std::string{"Could not open model file "}+modname);
66
67 std::ostringstream buffer;
68 buffer << f.rdbuf();
69 std::string contents{buffer.str()};
70
71 try
72 {
73 model = std::make_unique<ogdyn::DynareParser>(contents, ord);
74 }
75 catch (const ogp::ParserException &pe)
76 {
77 // Compute line and column, given the offset in the file
78 int line = 1;
79 int col = 0;
80 size_t i = 0;
81 while (i < contents.length() && i < static_cast<size_t>(pe.offset()))
82 {
83 if (contents[i] == '\n')
84 {
85 line++;
86 col = 0;
87 }
88 i++;
89 col++;
90 }
91 throw DynareException(pe.message(), modname, line, col);
92 }
93 ysteady = std::make_unique<Vector>(model->getAtoms().ny());
94 dnl = std::make_unique<DynareNameList>(*this);
95 denl = std::make_unique<DynareExogNameList>(*this);
96 dsnl = std::make_unique<DynareStateNameList>(*this, *dnl, *denl);
97 fe = std::make_unique<ogp::FormulaEvaluator>(model->getParser());
98 fde = std::make_unique<ogp::FormulaDerEvaluator>(model->getParser());
99 writeModelInfo(journal);
100 }
101
Dynare(const std::vector<std::string> & endo,const std::vector<std::string> & exo,const std::vector<std::string> & par,const std::string & equations,int ord,double sstol,Journal & jr)102 Dynare::Dynare(const std::vector<std::string> &endo,
103 const std::vector<std::string> &exo,
104 const std::vector<std::string> &par,
105 const std::string &equations, int ord,
106 double sstol, Journal &jr)
107 : journal(jr), md(1), ss_tol(sstol)
108 {
109 try
110 {
111 model = std::make_unique<ogdyn::DynareSPModel>(endo, exo, par, equations, ord);
112 }
113 catch (const ogp::ParserException &pe)
114 {
115 throw DynareException(pe.message(), pe.offset());
116 }
117 ysteady = std::make_unique<Vector>(model->getAtoms().ny());
118 dnl = std::make_unique<DynareNameList>(*this);
119 denl = std::make_unique<DynareExogNameList>(*this);
120 dsnl = std::make_unique<DynareStateNameList>(*this, *dnl, *denl);
121 fe = std::make_unique<ogp::FormulaEvaluator>(model->getParser());
122 fde = std::make_unique<ogp::FormulaDerEvaluator>(model->getParser());
123 writeModelInfo(journal);
124 }
125
Dynare(const Dynare & dynare)126 Dynare::Dynare(const Dynare &dynare)
127 : journal(dynare.journal), md(dynare.md), ss_tol(dynare.ss_tol)
128 {
129 model = dynare.model->clone();
130 ysteady = std::make_unique<Vector>(*(dynare.ysteady));
131 dnl = std::make_unique<DynareNameList>(*this);
132 denl = std::make_unique<DynareExogNameList>(*this);
133 dsnl = std::make_unique<DynareStateNameList>(*this, *dnl, *denl);
134 fe = std::make_unique<ogp::FormulaEvaluator>(model->getParser());
135 fde = std::make_unique<ogp::FormulaDerEvaluator>(model->getParser());
136 }
137
138 void
writeMat(mat_t * fd,const std::string & prefix) const139 Dynare::writeMat(mat_t *fd, const std::string &prefix) const
140 {
141 getAllEndoNames().writeMat(fd, prefix + "_vars");
142 getAllEndoNames().writeMatIndices(fd, prefix);
143 getStateNames().writeMat(fd, prefix + "_state_vars");
144 getExogNames().writeMat(fd, prefix + "_shocks");
145 getExogNames().writeMatIndices(fd, prefix);
146 model->getVcov().writeMat(fd, prefix + "_vcov_exo");
147 TwoDMatrix aux(1, 1);
148 aux.get(0, 0) = nstat();
149 aux.writeMat(fd, prefix + "_nstat");
150 aux.get(0, 0) = npred();
151 aux.writeMat(fd, prefix + "_npred");
152 aux.get(0, 0) = nboth();
153 aux.writeMat(fd, prefix + "_nboth");
154 aux.get(0, 0) = nforw();
155 aux.writeMat(fd, prefix + "_nforw");
156 }
157
158 void
writeDump(const std::string & basename) const159 Dynare::writeDump(const std::string &basename) const
160 {
161 std::string fname(basename + ".dump");
162 std::ofstream out(fname);
163 model->dump_model(out);
164 out.close();
165 }
166
167 void
solveDeterministicSteady(Vector & steady)168 Dynare::solveDeterministicSteady(Vector &steady)
169 {
170 JournalRecordPair pa(journal);
171 pa << "Non-linear solver for deterministic steady state" << endrec;
172 steady = const_cast<const Vector &>(model->getInit());
173 DynareVectorFunction dvf(*this);
174 DynareJacobian dj(*this);
175 ogu::NLSolver nls(dvf, dj, 500, ss_tol, journal);
176 int iter;
177 if (!nls.solve(steady, iter))
178 throw DynareException(__FILE__, __LINE__,
179 "Could not obtain convergence in non-linear solver");
180 }
181
182 // Evaluate system at given yₜ=yₜ₊₁=yₜ₋₁, and given shocks xₜ
183 void
evaluateSystem(Vector & out,const ConstVector & yy,const Vector & xx)184 Dynare::evaluateSystem(Vector &out, const ConstVector &yy, const Vector &xx)
185 {
186 ConstVector yym(yy, nstat(), nys());
187 ConstVector yyp(yy, nstat()+npred(), nyss());
188 evaluateSystem(out, yym, yy, yyp, xx);
189 }
190
191 /* Evaluate system at given y*ₜ₋₁, yₜ, y**ₜ₊₁ and at exogenous xₜ, all three
192 vectors yym, yy, and yyp have the respective lengths of y*ₜ₋₁, yₜ, y**ₜ₊₁ */
193 void
evaluateSystem(Vector & out,const ConstVector & yym,const ConstVector & yy,const ConstVector & yyp,const Vector & xx)194 Dynare::evaluateSystem(Vector &out, const ConstVector &yym, const ConstVector &yy,
195 const ConstVector &yyp, const Vector &xx)
196 {
197 ogdyn::DynareAtomValues dav(model->getAtoms(), model->getParams(), yym, yy, yyp, xx);
198 DynareEvalLoader del(model->getAtoms(), out);
199 fe->eval(dav, del);
200 }
201
202 void
calcDerivatives(const Vector & yy,const Vector & xx)203 Dynare::calcDerivatives(const Vector &yy, const Vector &xx)
204 {
205 ConstVector yym(yy, nstat(), nys());
206 ConstVector yyp(yy, nstat()+npred(), nyss());
207 ogdyn::DynareAtomValues dav(model->getAtoms(), model->getParams(), yym, yy, yyp, xx);
208 DynareDerEvalLoader ddel(model->getAtoms(), md, model->getOrder());
209 for (int iord = 1; iord <= model->getOrder(); iord++)
210 fde->eval(dav, ddel, iord);
211 }
212
213 void
calcDerivativesAtSteady()214 Dynare::calcDerivativesAtSteady()
215 {
216 Vector xx(nexog());
217 xx.zeros();
218 calcDerivatives(*ysteady, xx);
219 }
220
221 void
writeModelInfo(Journal & jr) const222 Dynare::writeModelInfo(Journal &jr) const
223 {
224 // write info on variables
225 {
226 JournalRecordPair rp(journal);
227 rp << "Information on variables" << endrec;
228 JournalRecord rec1(journal);
229 rec1 << "Number of endogenous: " << ny() << endrec;
230 JournalRecord rec2(journal);
231 rec2 << "Number of exogenous: " << nexog() << endrec;
232 JournalRecord rec3(journal);
233 rec3 << "Number of static: " << nstat() << endrec;
234 JournalRecord rec4(journal);
235 rec4 << "Number of predetermined: " << npred()+nboth() << endrec;
236 JournalRecord rec5(journal);
237 rec5 << "Number of forward looking: " << nforw()+nboth() << endrec;
238 JournalRecord rec6(journal);
239 rec6 << "Number of both: " << nboth() << endrec;
240 }
241
242 // write info on planner variables
243 const ogdyn::PlannerInfo *pinfo = model->get_planner_info();
244 if (pinfo)
245 {
246 JournalRecordPair rp(journal);
247 rp << "Information on planner variables" << endrec;
248 JournalRecord rec1(journal);
249 rec1 << "Number of Lagrange multipliers: " << pinfo->num_lagrange_mults << endrec;
250 JournalRecord rec2(journal);
251 rec2 << "Number of auxiliary variables: " << pinfo->num_aux_variables << endrec;
252 JournalRecord rec3(journal);
253 rec3 << "Number of new terms in the tree: " << pinfo->num_new_terms << endrec;
254 }
255
256 // write info on forward substitutions
257 const ogdyn::ForwSubstInfo *finfo = model->get_forw_subst_info();
258 if (finfo)
259 {
260 JournalRecordPair rp(journal);
261 rp << "Information on forward substitutions" << endrec;
262 JournalRecord rec1(journal);
263 rec1 << "Number of affected equations: " << finfo->num_affected_equations << endrec;
264 JournalRecord rec2(journal);
265 rec2 << "Number of substituted terms: " << finfo->num_subst_terms << endrec;
266 JournalRecord rec3(journal);
267 rec3 << "Number of auxiliary variables: " << finfo->num_aux_variables << endrec;
268 JournalRecord rec4(journal);
269 rec4 << "Number of new terms in the tree: " << finfo->num_new_terms << endrec;
270 }
271
272 // write info on substitutions
273 const ogp::SubstInfo *sinfo = model->get_subst_info();
274 if (sinfo)
275 {
276 JournalRecordPair rp(journal);
277 rp << "Information on substitutions" << endrec;
278 JournalRecord rec1(journal);
279 rec1 << "Number of substitutions: " << sinfo->num_substs << endrec;
280 }
281 }
282
DynareNameList(const Dynare & dynare)283 DynareNameList::DynareNameList(const Dynare &dynare)
284 {
285 for (int i = 0; i < dynare.ny(); i++)
286 {
287 int j = dynare.model->getAtoms().y2outer_endo()[i];
288 const std::string &name = dynare.model->getAtoms().get_endovars()[j];
289 names.push_back(name);
290 }
291 }
292
DynareStateNameList(const Dynare & dynare,const DynareNameList & dnl,const DynareExogNameList & denl)293 DynareStateNameList::DynareStateNameList(const Dynare &dynare, const DynareNameList &dnl,
294 const DynareExogNameList &denl)
295 {
296 for (int i = 0; i < dynare.nys(); i++)
297 names.push_back(dnl.getName(i+dynare.nstat()));
298 for (int i = 0; i < dynare.nexog(); i++)
299 names.push_back(denl.getName(i));
300 }
301
DynareExogNameList(const Dynare & dynare)302 DynareExogNameList::DynareExogNameList(const Dynare &dynare)
303 {
304 for (int i = 0; i < dynare.nexog(); i++)
305 {
306 int j = dynare.model->getAtoms().y2outer_exo()[i];
307 const std::string &name = dynare.model->getAtoms().get_exovars()[j];
308 names.push_back(name);
309 }
310 }
311
DynareEvalLoader(const ogp::FineAtoms & a,Vector & out)312 DynareEvalLoader::DynareEvalLoader(const ogp::FineAtoms &a, Vector &out)
313 : Vector(out)
314 {
315 if (a.ny() != out.length())
316 throw DynareException(__FILE__, __LINE__, "Wrong length of out vector in DynareEvalLoader constructor");
317 }
318
319 /* This clears the container of model derivatives and initializes it inserting
320 empty sparse tensors up to the given order. */
DynareDerEvalLoader(const ogp::FineAtoms & a,TensorContainer<FSSparseTensor> & mod_ders,int order)321 DynareDerEvalLoader::DynareDerEvalLoader(const ogp::FineAtoms &a,
322 TensorContainer<FSSparseTensor> &mod_ders,
323 int order)
324 : atoms(a), md(mod_ders)
325 {
326 md.clear();
327 for (int iord = 1; iord <= order; iord++)
328 {
329 auto t = std::make_unique<FSSparseTensor>(iord, atoms.ny()+atoms.nys()+atoms.nyss()+atoms.nexo(), atoms.ny());
330 md.insert(std::move(t));
331 }
332 }
333
334 void
load(int i,int iord,const int * vars,double res)335 DynareDerEvalLoader::load(int i, int iord, const int *vars, double res)
336 {
337 FSSparseTensor &t = md.get(Symmetry{iord});
338 IntSequence s(iord, 0);
339 for (int j = 0; j < iord; j++)
340 s[j] = atoms.get_pos_of_all(vars[j]);
341 t.insert(s, i, res);
342 }
343
DynareJacobian(Dynare & dyn)344 DynareJacobian::DynareJacobian(Dynare &dyn)
345 : Jacobian(dyn.ny()), d(dyn)
346 {
347 zeros();
348 }
349
350 void
eval(const Vector & yy)351 DynareJacobian::eval(const Vector &yy)
352 {
353 ogdyn::DynareSteadyAtomValues
354 dav(d.getModel().getAtoms(), d.getModel().getParams(), yy);
355 zeros();
356 d.fde->eval(dav, *this, 1);
357 }
358
359 void
load(int i,int iord,const int * vars,double res)360 DynareJacobian::load(int i, int iord, const int *vars, double res)
361 {
362 if (iord != 1)
363 throw DynareException(__FILE__, __LINE__,
364 "Derivative order different from order=1 in DynareJacobian::load");
365
366 int t = vars[0];
367 int j = d.getModel().getAtoms().get_pos_of_all(t);
368 if (j < d.nyss())
369 get(i, j+d.nstat()+d.npred()) += res;
370 else if (j < d.nyss()+d.ny())
371 get(i, j-d.nyss()) += res;
372 else if (j < d.nyss()+d.ny()+d.nys())
373 get(i, j-d.nyss()-d.ny()+d.nstat()) += res;
374 }
375
376 void
eval(const ConstVector & in,Vector & out)377 DynareVectorFunction::eval(const ConstVector &in, Vector &out)
378 {
379 check_for_eval(in, out);
380 Vector xx(d.nexog());
381 xx.zeros();
382 d.evaluateSystem(out, in, xx);
383 }
384