Commit 8698b4c5 authored by Sébastien Villemot's avatar Sébastien Villemot

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
parent b2fa7dd3
Pipeline #749 passed with stages
in 66 minutes and 52 seconds
......@@ -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));
}
}
......
......@@ -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;
};
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment