1 /*
2 * Copyright © 2004 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 "kord_exception.hh"
22 #include "decision_rule.hh"
23 #include "dynamic_model.hh"
24 #include "seed_generator.hh"
25
26 #include "SymSchurDecomp.hh"
27
28 #include <dynlapack.h>
29
30 #include <limits>
31 #include <utility>
32 #include <memory>
33
34 // FoldDecisionRule conversion from UnfoldDecisionRule
FoldDecisionRule(const UnfoldDecisionRule & udr)35 FoldDecisionRule::FoldDecisionRule(const UnfoldDecisionRule &udr)
36 : DecisionRuleImpl<Storage::fold>(ctraits<Storage::fold>::Tpol(udr.nrows(), udr.nvars()),
37 udr.ypart, udr.nu, udr.ysteady)
38 {
39 for (const auto &it : udr)
40 insert(std::make_unique<ctraits<Storage::fold>::Ttensym>(*(it.second)));
41 }
42
43 // UnfoldDecisionRule conversion from FoldDecisionRule
UnfoldDecisionRule(const FoldDecisionRule & fdr)44 UnfoldDecisionRule::UnfoldDecisionRule(const FoldDecisionRule &fdr)
45 : DecisionRuleImpl<Storage::unfold>(ctraits<Storage::unfold>::Tpol(fdr.nrows(), fdr.nvars()),
46 fdr.ypart, fdr.nu, fdr.ysteady)
47 {
48 for (const auto &it : fdr)
49 insert(std::make_unique<ctraits<Storage::unfold>::Ttensym>(*(it.second)));
50 }
51
52 /* This runs simulations with an output to journal file. Note that we
53 report how many simulations had to be thrown out due to Nan or Inf. */
54
55 void
simulate(int num_sim,const DecisionRule & dr,const Vector & start,const TwoDMatrix & vcov,Journal & journal)56 SimResults::simulate(int num_sim, const DecisionRule &dr, const Vector &start,
57 const TwoDMatrix &vcov, Journal &journal)
58 {
59 JournalRecordPair paa(journal);
60 paa << "Performing " << num_sim << " stochastic simulations for "
61 << num_per << " periods burning " << num_burn << " initial periods" << endrec;
62 simulate(num_sim, dr, start, vcov);
63 int thrown = num_sim - data.size();
64 if (thrown > 0)
65 {
66 JournalRecord rec(journal);
67 rec << "I had to throw " << thrown << " simulations away due to Nan or Inf" << endrec;
68 }
69 }
70
71 /* This runs a given number of simulations by creating
72 SimulationWorker for each simulation and inserting them to the
73 thread group. */
74
75 void
simulate(int num_sim,const DecisionRule & dr,const Vector & start,const TwoDMatrix & vcov)76 SimResults::simulate(int num_sim, const DecisionRule &dr, const Vector &start,
77 const TwoDMatrix &vcov)
78 {
79 std::vector<RandomShockRealization> rsrs;
80 rsrs.reserve(num_sim);
81
82 sthread::detach_thread_group gr;
83 for (int i = 0; i < num_sim; i++)
84 {
85 RandomShockRealization sr(vcov, seed_generator::get_new_seed());
86 rsrs.push_back(sr);
87 gr.insert(std::make_unique<SimulationWorker>(*this, dr, DecisionRule::emethod::horner,
88 num_per+num_burn, start, rsrs.back()));
89 }
90 gr.run();
91 }
92
93 /* This adds the data with the realized shocks. It takes only periods
94 which are not to be burnt. If the data is not finite, the both data
95 and shocks are thrown away. */
96
97 bool
addDataSet(const TwoDMatrix & d,const ExplicitShockRealization & sr,const ConstVector & st)98 SimResults::addDataSet(const TwoDMatrix &d, const ExplicitShockRealization &sr, const ConstVector &st)
99 {
100 KORD_RAISE_IF(d.nrows() != num_y,
101 "Incompatible number of rows for SimResults::addDataSets");
102 KORD_RAISE_IF(d.ncols() != num_per+num_burn,
103 "Incompatible number of cols for SimResults::addDataSets");
104 bool ret = false;
105 if (d.isFinite())
106 {
107 data.emplace_back(d, num_burn, num_per);
108 shocks.emplace_back(ConstTwoDMatrix(sr.getShocks(), num_burn, num_per));
109 if (num_burn == 0)
110 start.emplace_back(st);
111 else
112 start.emplace_back(d.getCol(num_burn-1));
113 ret = true;
114 }
115
116 return ret;
117 }
118
119 void
writeMat(const std::string & base,const std::string & lname) const120 SimResults::writeMat(const std::string &base, const std::string &lname) const
121 {
122 std::string matfile_name = base + ".mat";
123 mat_t *matfd = Mat_Create(matfile_name.c_str(), nullptr);
124 if (matfd)
125 {
126 writeMat(matfd, lname);
127 Mat_Close(matfd);
128 }
129 }
130
131 /* This save the results as matrices with given prefix and with index
132 appended. If there is only one matrix, the index is not appended. */
133
134 void
writeMat(mat_t * fd,const std::string & lname) const135 SimResults::writeMat(mat_t *fd, const std::string &lname) const
136 {
137 for (int i = 0; i < getNumSets(); i++)
138 {
139 std::string tmp = lname + "_data";
140 if (getNumSets() > 1)
141 tmp += std::to_string(i+1);
142 data[i].writeMat(fd, tmp);
143 }
144 }
145
146 void
simulate(int num_sim,const DecisionRule & dr,const Vector & start,const TwoDMatrix & vcov,Journal & journal)147 SimResultsStats::simulate(int num_sim, const DecisionRule &dr,
148 const Vector &start,
149 const TwoDMatrix &vcov, Journal &journal)
150 {
151 SimResults::simulate(num_sim, dr, start, vcov, journal);
152 {
153 JournalRecordPair paa(journal);
154 paa << "Calculating means from the simulations." << endrec;
155 calcMean();
156 }
157 {
158 JournalRecordPair paa(journal);
159 paa << "Calculating covariances from the simulations." << endrec;
160 calcVcov();
161 }
162 }
163
164 /* Here we do not save the data itself, we save only mean and vcov. */
165 void
writeMat(mat_t * fd,const std::string & lname) const166 SimResultsStats::writeMat(mat_t *fd, const std::string &lname) const
167 {
168 ConstTwoDMatrix(num_y, 1, mean).writeMat(fd, lname + "_mean");;
169 vcov.writeMat(fd, lname + "_vcov");
170 }
171
172 void
calcMean()173 SimResultsStats::calcMean()
174 {
175 mean.zeros();
176 if (data.size()*num_per > 0)
177 {
178 double mult = 1.0/data.size()/num_per;
179 for (const auto &i : data)
180 {
181 for (int j = 0; j < num_per; j++)
182 {
183 ConstVector col{i.getCol(j)};
184 mean.add(mult, col);
185 }
186 }
187 }
188 }
189
190 void
calcVcov()191 SimResultsStats::calcVcov()
192 {
193 if (data.size()*num_per > 1)
194 {
195 vcov.zeros();
196 double mult = 1.0/(data.size()*num_per - 1);
197 for (const auto &d : data)
198 for (int j = 0; j < num_per; j++)
199 for (int m = 0; m < num_y; m++)
200 for (int n = m; n < num_y; n++)
201 {
202 double s = (d.get(m, j)-mean[m])*(d.get(n, j)-mean[n]);
203 vcov.get(m, n) += mult*s;
204 if (m != n)
205 vcov.get(n, m) += mult*s;
206 }
207 }
208 else
209 vcov.infs();
210 }
211
212 void
simulate(int num_sim,const DecisionRule & dr,const Vector & start,const TwoDMatrix & vcov,Journal & journal)213 SimResultsDynamicStats::simulate(int num_sim, const DecisionRule &dr,
214 const Vector &start,
215 const TwoDMatrix &vcov, Journal &journal)
216 {
217 SimResults::simulate(num_sim, dr, start, vcov, journal);
218 {
219 JournalRecordPair paa(journal);
220 paa << "Calculating means of the conditional simulations." << endrec;
221 calcMean();
222 }
223 {
224 JournalRecordPair paa(journal);
225 paa << "Calculating variances of the conditional simulations." << endrec;
226 calcVariance();
227 }
228 }
229
230 void
writeMat(mat_t * fd,const std::string & lname) const231 SimResultsDynamicStats::writeMat(mat_t *fd, const std::string &lname) const
232 {
233 mean.writeMat(fd, lname + "_cond_mean");
234 variance.writeMat(fd, lname + "_cond_variance");
235 }
236
237 void
calcMean()238 SimResultsDynamicStats::calcMean()
239 {
240 mean.zeros();
241 if (data.size() > 0)
242 {
243 double mult = 1.0/data.size();
244 for (int j = 0; j < num_per; j++)
245 {
246 Vector meanj{mean.getCol(j)};
247 for (const auto &i : data)
248 {
249 ConstVector col{i.getCol(j)};
250 meanj.add(mult, col);
251 }
252 }
253 }
254 }
255
256 void
calcVariance()257 SimResultsDynamicStats::calcVariance()
258 {
259 if (data.size() > 1)
260 {
261 variance.zeros();
262 double mult = 1.0/(data.size()-1);
263 for (int j = 0; j < num_per; j++)
264 {
265 ConstVector meanj{mean.getCol(j)};
266 Vector varj{variance.getCol(j)};
267 for (const auto &i : data)
268 {
269 Vector col{i.getCol(j)};
270 col.add(-1.0, meanj);
271 for (int k = 0; k < col.length(); k++)
272 col[k] = col[k]*col[k];
273 varj.add(mult, col);
274 }
275 }
276 }
277 else
278 variance.infs();
279 }
280
281 void
simulate(const DecisionRule & dr,Journal & journal)282 SimResultsIRF::simulate(const DecisionRule &dr, Journal &journal)
283 {
284 JournalRecordPair paa(journal);
285 paa << "Performing " << control.getNumSets() << " IRF simulations for "
286 << num_per << " periods; shock=" << ishock << ", impulse=" << imp << endrec;
287 simulate(dr);
288 int thrown = control.getNumSets() - data.size();
289 if (thrown > 0)
290 {
291 JournalRecord rec(journal);
292 rec << "I had to throw " << thrown
293 << " simulations away due to Nan or Inf" << endrec;
294 }
295 calcMeans();
296 calcVariances();
297 }
298
299 void
simulate(const DecisionRule & dr)300 SimResultsIRF::simulate(const DecisionRule &dr)
301 {
302 sthread::detach_thread_group gr;
303 for (int idata = 0; idata < control.getNumSets(); idata++)
304 gr.insert(std::make_unique<SimulationIRFWorker>(*this, dr, DecisionRule::emethod::horner,
305 num_per, idata, ishock, imp));
306 gr.run();
307 }
308
309 void
calcMeans()310 SimResultsIRF::calcMeans()
311 {
312 means.zeros();
313 if (data.size() > 0)
314 {
315 for (const auto &i : data)
316 means.add(1.0, i);
317 means.mult(1.0/data.size());
318 }
319 }
320
321 void
calcVariances()322 SimResultsIRF::calcVariances()
323 {
324 if (data.size() > 1)
325 {
326 variances.zeros();
327 for (const auto &i : data)
328 {
329 TwoDMatrix d(i);
330 d.add(-1.0, means);
331 for (int j = 0; j < d.nrows(); j++)
332 for (int k = 0; k < d.ncols(); k++)
333 variances.get(j, k) += d.get(j, k)*d.get(j, k);
334 d.mult(1.0/(data.size()-1));
335 }
336 }
337 else
338 variances.infs();
339 }
340
341 void
writeMat(mat_t * fd,const std::string & lname) const342 SimResultsIRF::writeMat(mat_t *fd, const std::string &lname) const
343 {
344 means.writeMat(fd, lname + "_mean");
345 variances.writeMat(fd, lname + "_var");
346 }
347
348 void
simulate(int num_sim,const DecisionRule & dr,const Vector & start,const TwoDMatrix & v,Journal & journal)349 RTSimResultsStats::simulate(int num_sim, const DecisionRule &dr, const Vector &start,
350 const TwoDMatrix &v, Journal &journal)
351 {
352 JournalRecordPair paa(journal);
353 paa << "Performing " << num_sim << " real-time stochastic simulations for "
354 << num_per << " periods" << endrec;
355 simulate(num_sim, dr, start, v);
356 mean = nc.getMean();
357 mean.add(1.0, dr.getSteady());
358 nc.getVariance(vcov);
359 if (thrown_periods > 0)
360 {
361 JournalRecord rec(journal);
362 rec << "I had to throw " << thrown_periods << " periods away due to Nan or Inf" << endrec;
363 JournalRecord rec1(journal);
364 rec1 << "This affected " << incomplete_simulations << " out of "
365 << num_sim << " simulations" << endrec;
366 }
367 }
368
369 void
simulate(int num_sim,const DecisionRule & dr,const Vector & start,const TwoDMatrix & vcov)370 RTSimResultsStats::simulate(int num_sim, const DecisionRule &dr, const Vector &start,
371 const TwoDMatrix &vcov)
372 {
373 std::vector<RandomShockRealization> rsrs;
374 rsrs.reserve(num_sim);
375
376 sthread::detach_thread_group gr;
377 for (int i = 0; i < num_sim; i++)
378 {
379 RandomShockRealization sr(vcov, seed_generator::get_new_seed());
380 rsrs.push_back(sr);
381 gr.insert(std::make_unique<RTSimulationWorker>(*this, dr, DecisionRule::emethod::horner,
382 num_per, start, rsrs.back()));
383 }
384 gr.run();
385 }
386
387 void
writeMat(mat_t * fd,const std::string & lname)388 RTSimResultsStats::writeMat(mat_t *fd, const std::string &lname)
389 {
390 ConstTwoDMatrix(nc.getDim(), 1, mean).writeMat(fd, lname + "_rt_mean");
391 vcov.writeMat(fd, lname + "_rt_vcov");
392 }
393
IRFResults(const DynamicModel & mod,const DecisionRule & dr,const SimResults & control,std::vector<int> ili,Journal & journal)394 IRFResults::IRFResults(const DynamicModel &mod, const DecisionRule &dr,
395 const SimResults &control, std::vector<int> ili,
396 Journal &journal)
397 : model(mod), irf_list_ind(std::move(ili))
398 {
399 int num_per = control.getNumPer();
400 JournalRecordPair pa(journal);
401 pa << "Calculating IRFs against control for " << static_cast<int>(irf_list_ind.size()) << " shocks and for "
402 << num_per << " periods" << endrec;
403 const TwoDMatrix &vcov = mod.getVcov();
404 for (int ishock : irf_list_ind)
405 {
406 double stderror = sqrt(vcov.get(ishock, ishock));
407 irf_res.emplace_back(control, model.numeq(), num_per,
408 ishock, stderror);
409 irf_res.emplace_back(control, model.numeq(), num_per,
410 ishock, -stderror);
411 }
412
413 for (unsigned int ii = 0; ii < irf_list_ind.size(); ii++)
414 {
415 irf_res[2*ii].simulate(dr, journal);
416 irf_res[2*ii+1].simulate(dr, journal);
417 }
418 }
419
420 void
writeMat(mat_t * fd,const std::string & prefix) const421 IRFResults::writeMat(mat_t *fd, const std::string &prefix) const
422 {
423 for (unsigned int i = 0; i < irf_list_ind.size(); i++)
424 {
425 int ishock = irf_list_ind[i];
426 auto shockname = model.getExogNames().getName(ishock);
427 irf_res[2*i].writeMat(fd, prefix + "_irfp_" + shockname);
428 irf_res[2*i+1].writeMat(fd, prefix + "_irfm_" + shockname);
429 }
430 }
431
432 void
operator ()(std::mutex & mut)433 SimulationWorker::operator()(std::mutex &mut)
434 {
435 ExplicitShockRealization esr(sr, np);
436 TwoDMatrix m{dr.simulate(em, np, st, esr)};
437 {
438 std::unique_lock<std::mutex> lk{mut};
439 res.addDataSet(m, esr, st);
440 }
441 }
442
443 /* Here we create a new instance of ExplicitShockRealization of the
444 corresponding control, add the impulse, and simulate. */
445
446 void
operator ()(std::mutex & mut)447 SimulationIRFWorker::operator()(std::mutex &mut)
448 {
449 ExplicitShockRealization esr(res.control.getShocks(idata));
450 esr.addToShock(ishock, 0, imp);
451 TwoDMatrix m{dr.simulate(em, np, res.control.getStart(idata), esr)};
452 m.add(-1.0, res.control.getData(idata));
453 {
454 std::unique_lock<std::mutex> lk{mut};
455 res.addDataSet(m, esr, res.control.getStart(idata));
456 }
457 }
458
459 void
operator ()(std::mutex & mut)460 RTSimulationWorker::operator()(std::mutex &mut)
461 {
462 NormalConj nc(res.nc.getDim());
463 const PartitionY &ypart = dr.getYPart();
464 int nu = dr.nexog();
465 const Vector &ysteady = dr.getSteady();
466
467 // initialize vectors and subvectors for simulation
468 Vector dyu(ypart.nys()+nu);
469 ConstVector ystart_pred(ystart, ypart.nstat, ypart.nys());
470 ConstVector ysteady_pred(ysteady, ypart.nstat, ypart.nys());
471 Vector dy(dyu, 0, ypart.nys());
472 Vector u(dyu, ypart.nys(), nu);
473 Vector y(nc.getDim());
474 ConstVector ypred(y, ypart.nstat, ypart.nys());
475
476 // simulate the first real-time period
477 int ip = 0;
478 dy = ystart_pred;
479 dy.add(-1.0, ysteady_pred);
480 sr.get(ip, u);
481 dr.eval(em, y, dyu);
482 if (ip >= res.num_burn)
483 nc.update(y);
484
485 // simulate other real-time periods
486 while (y.isFinite() && ip < res.num_burn + res.num_per)
487 {
488 ip++;
489 dy = ypred;
490 sr.get(ip, u);
491 dr.eval(em, y, dyu);
492 if (ip >= res.num_burn)
493 nc.update(y);
494 }
495 {
496 std::unique_lock<std::mutex> lk{mut};
497 res.nc.update(nc);
498 if (res.num_per-ip > 0)
499 {
500 res.incomplete_simulations++;
501 res.thrown_periods += res.num_per-ip;
502 }
503 }
504 }
505
506 /* This calculates factorization FFᵀ=V in the Cholesky way. It does
507 not work for semidefinite matrices. */
508
509 void
choleskyFactor(const ConstTwoDMatrix & v)510 RandomShockRealization::choleskyFactor(const ConstTwoDMatrix &v)
511 {
512 factor = v;
513 lapack_int rows = factor.nrows(), lda = factor.getLD();
514 for (int i = 0; i < rows; i++)
515 for (int j = i+1; j < rows; j++)
516 factor.get(i, j) = 0.0;
517 lapack_int info;
518
519 dpotrf("L", &rows, factor.base(), &lda, &info);
520 KORD_RAISE_IF(info != 0,
521 "Info!=0 in RandomShockRealization::choleskyFactor");
522 }
523
524 /* This calculates FFᵀ=V factorization by symmetric Schur
525 decomposition. It works for semidefinite matrices. */
526
527 void
schurFactor(const ConstTwoDMatrix & v)528 RandomShockRealization::schurFactor(const ConstTwoDMatrix &v)
529 {
530 SymSchurDecomp(v).getFactor(factor);
531 }
532
533 void
get(int n,Vector & out)534 RandomShockRealization::get(int n, Vector &out)
535 {
536 KORD_RAISE_IF(out.length() != numShocks(),
537 "Wrong length of out vector in RandomShockRealization::get");
538 Vector d(out.length());
539 for (int i = 0; i < d.length(); i++)
540 d[i] = dis(mtwister);
541 out.zeros();
542 factor.multaVec(out, ConstVector(d));
543 }
544
ExplicitShockRealization(ShockRealization & sr,int num_per)545 ExplicitShockRealization::ExplicitShockRealization(ShockRealization &sr,
546 int num_per)
547 : shocks(sr.numShocks(), num_per)
548 {
549 for (int j = 0; j < num_per; j++)
550 {
551 Vector jcol{shocks.getCol(j)};
552 sr.get(j, jcol);
553 }
554 }
555
556 void
get(int n,Vector & out)557 ExplicitShockRealization::get(int n, Vector &out)
558 {
559 KORD_RAISE_IF(out.length() != numShocks(),
560 "Wrong length of out vector in ExplicitShockRealization::get");
561 int i = n % shocks.ncols();
562 ConstVector icol{shocks.getCol(i)};
563 out = icol;
564 }
565
566 void
addToShock(int ishock,int iper,double val)567 ExplicitShockRealization::addToShock(int ishock, int iper, double val)
568 {
569 KORD_RAISE_IF(ishock < 0 || ishock > numShocks(),
570 "Wrong index of shock in ExplicitShockRealization::addToShock");
571 int j = iper % shocks.ncols();
572 shocks.get(ishock, j) += val;
573 }
574
575 void
get(int n,Vector & out)576 GenShockRealization::get(int n, Vector &out)
577 {
578 KORD_RAISE_IF(out.length() != numShocks(),
579 "Wrong length of out vector in GenShockRealization::get");
580 ExplicitShockRealization::get(n, out);
581 Vector r(numShocks());
582 RandomShockRealization::get(n, r);
583 for (int j = 0; j < numShocks(); j++)
584 if (!std::isfinite(out[j]))
585 out[j] = r[j];
586 }
587