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 8698b4c5)

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));
 	}
 }
 

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;
 };