From 8698b4c540f2f9881f2201a30b6d08344b23068f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Tue, 5 Feb 2019 17:14:59 +0100
Subject: [PATCH] Dynare++: fix computation of IRFs
Since the introduction of the --burn option (in Dynare++ shipped with Dynare
4.3.0), the IRFs reported by Dynare++ were wrong.
The IRFs are computed using a generalized IRF method: the result is
the (average) difference between a simulation with shock and a simulation
without shock. The problem was that the two simulations were not using the same
starting point.
Closes #1634
---
dynare++/kord/decision_rule.cc | 14 ++++++++------
dynare++/kord/decision_rule.hh | 12 ++++++++++--
2 files changed, 18 insertions(+), 8 deletions(-)
diff --git a/dynare++/kord/decision_rule.cc b/dynare++/kord/decision_rule.cc
index fc5dbce2e..7afc2b97a 100644
--- a/dynare++/kord/decision_rule.cc
+++ b/dynare++/kord/decision_rule.cc
@@ -106,7 +106,7 @@ SimResults::simulate(int num_sim, const DecisionRule &dr, const Vector &start,
and shocks are thrown away. */
bool
-SimResults::addDataSet(TwoDMatrix *d, ExplicitShockRealization *sr)
+SimResults::addDataSet(TwoDMatrix *d, ExplicitShockRealization *sr, const ConstVector &st)
{
KORD_RAISE_IF(d->nrows() != num_y,
"Incompatible number of rows for SimResults::addDataSets");
@@ -118,6 +118,10 @@ SimResults::addDataSet(TwoDMatrix *d, ExplicitShockRealization *sr)
data.push_back(new TwoDMatrix((const TwoDMatrix &) (*d), num_burn, num_per));
shocks.push_back(new ExplicitShockRealization(
ConstTwoDMatrix(sr->getShocks(), num_burn, num_per)));
+ if (num_burn == 0)
+ start.emplace_back(st);
+ else
+ start.emplace_back(d->getCol(num_burn-1));
ret = true;
}
@@ -486,7 +490,7 @@ SimulationWorker::operator()(std::mutex &mut)
TwoDMatrix *m = dr.simulate(em, np, st, *esr);
{
std::unique_lock<std::mutex> lk{mut};
- res.addDataSet(m, esr);
+ res.addDataSet(m, esr, st);
}
}
@@ -499,13 +503,11 @@ SimulationIRFWorker::operator()(std::mutex &mut)
auto *esr
= new ExplicitShockRealization(res.control.getShocks(idata));
esr->addToShock(ishock, 0, imp);
- const TwoDMatrix &data = res.control.getData(idata);
- Vector st{data.getCol(res.control.getNumBurn())};
- TwoDMatrix *m = dr.simulate(em, np, st, *esr);
+ TwoDMatrix *m = dr.simulate(em, np, res.control.getStart(idata), *esr);
m->add(-1.0, res.control.getData(idata));
{
std::unique_lock<std::mutex> lk{mut};
- res.addDataSet(m, esr);
+ res.addDataSet(m, esr, res.control.getStart(idata));
}
}
diff --git a/dynare++/kord/decision_rule.hh b/dynare++/kord/decision_rule.hh
index ef7210354..532c5ec22 100644
--- a/dynare++/kord/decision_rule.hh
+++ b/dynare++/kord/decision_rule.hh
@@ -718,7 +718,8 @@ DRFixPoint<t>::calcFixPoint(emethod em, Vector &out)
/* This is a basically a number of matrices of the same dimensions,
which can be obtained as simulation results from a given decision rule
- and shock realizations. We also store the realizations of shocks. */
+ and shock realizations. We also store the realizations of shocks and the
+ starting point of each simulation. */
class ExplicitShockRealization;
class SimResults
@@ -729,6 +730,7 @@ protected:
int num_burn;
vector<TwoDMatrix *> data;
vector<ExplicitShockRealization *> shocks;
+ vector<ConstVector> start;
public:
SimResults(int ny, int nper, int nburn = 0)
: num_y(ny), num_per(nper), num_burn(nburn)
@@ -765,7 +767,13 @@ public:
{
return *(shocks[i]);
}
- bool addDataSet(TwoDMatrix *d, ExplicitShockRealization *sr);
+ const ConstVector &
+ getStart(int i) const
+ {
+ return start[i];
+ }
+
+ bool addDataSet(TwoDMatrix *d, ExplicitShockRealization *sr, const ConstVector &st);
void writeMat(const char *base, const char *lname) const;
void writeMat(mat_t *fd, const char *lname) const;
};
--
GitLab