From 583e2dfc0af1c43611444be964fdfba0060c9833 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Tue, 5 Feb 2019 18:00:41 +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
(manually cherry-picked from 8698b4c540f2f9881f2201a30b6d08344b23068f)
---
dynare++/kord/decision_rule.cweb | 15 +++++++++------
dynare++/kord/decision_rule.hweb | 8 ++++++--
2 files changed, 15 insertions(+), 8 deletions(-)
diff --git a/dynare++/kord/decision_rule.cweb b/dynare++/kord/decision_rule.cweb
index c5ab3e5800..827f778d9c 100644
--- a/dynare++/kord/decision_rule.cweb
+++ b/dynare++/kord/decision_rule.cweb
@@ -100,6 +100,7 @@ SimResults::~SimResults()
for (int i = 0; i < getNumSets(); i++) {
delete data[i];
delete shocks[i];
+ delete start[i];
}
}
@@ -149,7 +150,7 @@ which are not to be burnt. If the data is not finite, the both data
and shocks are thrown away.
@<|SimResults::addDataSet| code@>=
-bool SimResults::addDataSet(TwoDMatrix* d, ExplicitShockRealization* sr)
+bool SimResults::addDataSet(TwoDMatrix* d, ExplicitShockRealization* sr, const ConstVector& st)
{
KORD_RAISE_IF(d->nrows() != num_y,
"Incompatible number of rows for SimResults::addDataSets");
@@ -160,6 +161,10 @@ bool 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.push_back(new Vector(st));
+ else
+ start.push_back(new Vector(ConstVector(*d, num_burn-1)));
ret = true;
}
@@ -527,7 +532,7 @@ void SimulationWorker::operator()()
TwoDMatrix* m = dr.simulate(em, np, st, *esr);
{
SYNCHRO syn(&res, "simulation");
- res.addDataSet(m, esr);
+ res.addDataSet(m, esr, st);
}
}
@@ -540,13 +545,11 @@ void SimulationIRFWorker::operator()()
ExplicitShockRealization* esr =
new ExplicitShockRealization(res.control.getShocks(idata));
esr->addToShock(ishock, 0, imp);
- const TwoDMatrix& data = res.control.getData(idata);
- ConstVector st(data, 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));
{
SYNCHRO syn(&res, "simulation");
- res.addDataSet(m, esr);
+ res.addDataSet(m, esr, res.control.getStart(idata));
}
}
diff --git a/dynare++/kord/decision_rule.hweb b/dynare++/kord/decision_rule.hweb
index 25cf370313..3d3796a8c2 100644
--- a/dynare++/kord/decision_rule.hweb
+++ b/dynare++/kord/decision_rule.hweb
@@ -694,7 +694,8 @@ bool 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.
@<|SimResults| class declaration@>=
class ExplicitShockRealization;
@@ -705,6 +706,7 @@ protected:@;
int num_burn;
vector<TwoDMatrix*> data;
vector<ExplicitShockRealization*> shocks;
+ vector<Vector *> start;
public:@;
SimResults(int ny, int nper, int nburn = 0)
: num_y(ny), num_per(nper), num_burn(nburn)@+ {}
@@ -723,7 +725,9 @@ public:@;
{@+ return *(data[i]);@+}
const ExplicitShockRealization& getShocks(int i) const
{ @+ return *(shocks[i]);@+}
- bool addDataSet(TwoDMatrix* d, ExplicitShockRealization* sr);
+ const Vector& 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