diff --git a/dynare++/kord/decision_rule.cweb b/dynare++/kord/decision_rule.cweb
index a4a4f59e1f3662a4dd6188543095dfe9587a6088..99536d3de155ea0a27044775c5bb0eecae3cd61f 100644
--- a/dynare++/kord/decision_rule.cweb
+++ b/dynare++/kord/decision_rule.cweb
@@ -112,7 +112,7 @@ void SimResults::simulate(int num_sim, const DecisionRule& dr, const Vector& sta
{
JournalRecordPair paa(journal);
paa << "Performing " << num_sim << " stochastic simulations for "
- << num_per << " periods" << endrec;
+ << num_per << " periods burning " << num_burn << " initial periods" << endrec;
simulate(num_sim, dr, start, vcov);
int thrown = num_sim - data.size();
if (thrown > 0) {
@@ -138,31 +138,34 @@ void SimResults::simulate(int num_sim, const DecisionRule& dr, const Vector& sta
rsrs.push_back(sr);
THREAD* worker = new
SimulationWorker(*this, dr, DecisionRule::horner,
- num_per, start, rsrs.back());
+ num_per+num_burn, start, rsrs.back());
gr.insert(worker);
}
gr.run();
}
-@ This adds the data with the realized shocks. If the data is not
-finite, the both data and shocks are thrown away.
+@ This adds the data with the realized shocks. It takes only periods
+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)
{
KORD_RAISE_IF(d->nrows() != num_y,
"Incompatible number of rows for SimResults::addDataSets");
- KORD_RAISE_IF(d->ncols() != num_per,
+ KORD_RAISE_IF(d->ncols() != num_per+num_burn,
"Incompatible number of cols for SimResults::addDataSets");
+ bool ret = false;
if (d->isFinite()) {
- data.push_back(d);
- shocks.push_back(sr);
- return true;
- } else {
- delete d;
- delete sr;
- return false;
+ data.push_back(new TwoDMatrix((const TwoDMatrix&)(*d),num_burn,num_per));
+ shocks.push_back(new ExplicitShockRealization(
+ ConstTwoDMatrix(sr->getShocks(),num_burn,num_per)));
+ ret = true;
}
+
+ delete d;
+ delete sr;
+ return ret;
}
@
@@ -341,13 +344,12 @@ void SimResultsDynamicStats::calcVariance()
@
@<|SimResultsIRF::simulate| code1@>=
-void SimResultsIRF::simulate(const DecisionRule& dr, const Vector& start,
- Journal& journal)
+void SimResultsIRF::simulate(const DecisionRule& dr, Journal& journal)
{
JournalRecordPair paa(journal);
paa << "Performing " << control.getNumSets() << " IRF simulations for "
<< num_per << " periods; shock=" << ishock << ", impulse=" << imp << endrec;
- simulate(dr, start);
+ simulate(dr);
int thrown = control.getNumSets() - data.size();
if (thrown > 0) {
JournalRecord rec(journal);
@@ -360,13 +362,13 @@ void SimResultsIRF::simulate(const DecisionRule& dr, const Vector& start,
@
@<|SimResultsIRF::simulate| code2@>=
-void SimResultsIRF::simulate(const DecisionRule& dr, const Vector& start)
+void SimResultsIRF::simulate(const DecisionRule& dr)
{
THREAD_GROUP gr;
for (int idata = 0; idata < control.getNumSets(); idata++) {
THREAD* worker = new
SimulationIRFWorker(*this, dr, DecisionRule::horner,
- num_per, start, idata, ishock, imp);
+ num_per, idata, ishock, imp);
gr.insert(worker);
}
gr.run();
@@ -489,8 +491,8 @@ IRFResults::IRFResults(const DynamicModel& mod, const DecisionRule& dr,
}
for (unsigned int ii = 0; ii < irf_list_ind.size(); ii++) {
- irf_res[2*ii]->simulate(dr, model.getSteady(), journal);
- irf_res[2*ii+1]->simulate(dr, model.getSteady(), journal);
+ irf_res[2*ii]->simulate(dr, journal);
+ irf_res[2*ii+1]->simulate(dr, journal);
}
}
@@ -538,6 +540,8 @@ 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);
m->add(-1.0, res.control.getData(idata));
{
@@ -585,16 +589,18 @@ void RTSimulationWorker::operator()()
dy.add(-1.0, ysteady_pred);
sr.get(ip, u);
dr.eval(em, y, dyu);
- nc.update(y);
+ if (ip >= res.num_burn)
+ nc.update(y);
@
@<simulate other real-time periods@>=
-while (y.isFinite() && ip < res.num_per) {
+while (y.isFinite() && ip < res.num_burn + res.num_per) {
ip++;
dy = ypred;
sr.get(ip, u);
dr.eval(em, y, dyu);
- nc.update(y);
+ if (ip >= res.num_burn)
+ nc.update(y);
}
@ This calculates factorization $FF^T=V$ in the Cholesky way. It does
diff --git a/dynare++/kord/decision_rule.hweb b/dynare++/kord/decision_rule.hweb
index caf7315faab80877bd9ac969f39a1c4fdf9a2ae7..f79b6f803e8c949a460e035e57f7e826dde48f7a 100644
--- a/dynare++/kord/decision_rule.hweb
+++ b/dynare++/kord/decision_rule.hweb
@@ -696,11 +696,12 @@ class SimResults {
protected:@;
int num_y;
int num_per;
+ int num_burn;
vector<TwoDMatrix*> data;
vector<ExplicitShockRealization*> shocks;
public:@;
- SimResults(int ny, int nper)
- : num_y(ny), num_per(nper)@+ {}
+ SimResults(int ny, int nper, int nburn = 0)
+ : num_y(ny), num_per(nper), num_burn(nburn)@+ {}
virtual ~SimResults();
void simulate(int num_sim, const DecisionRule& dr, const Vector& start,
const TwoDMatrix& vcov, Journal& journal);
@@ -708,6 +709,8 @@ public:@;
const TwoDMatrix& vcov);
int getNumPer() const
{@+ return num_per;@+}
+ int getNumBurn() const
+ {@+ return num_burn;@+}
int getNumSets() const
{@+ return (int)data.size();@+}
const TwoDMatrix& getData(int i) const
@@ -728,8 +731,8 @@ protected:@;
Vector mean;
TwoDMatrix vcov;
public:@;
- SimResultsStats(int ny, int nper)
- : SimResults(ny, nper), mean(ny), vcov(ny,ny)@+ {}
+ SimResultsStats(int ny, int nper, int nburn = 0)
+ : SimResults(ny, nper, nburn), mean(ny), vcov(ny,ny)@+ {}
void simulate(int num_sim, const DecisionRule& dr, const Vector& start,
const TwoDMatrix& vcov, Journal& journal);
void writeMat4(FILE* fd, const char* lname) const;
@@ -748,8 +751,8 @@ protected:@;
TwoDMatrix mean;
TwoDMatrix variance;
public:@;
- SimResultsDynamicStats(int ny, int nper)
- : SimResults(ny, nper), mean(ny,nper), variance(ny,nper)@+ {}
+ SimResultsDynamicStats(int ny, int nper, int nburn = 0)
+ : SimResults(ny, nper, nburn), mean(ny,nper), variance(ny,nper)@+ {}
void simulate(int num_sim, const DecisionRule& dr, const Vector& start,
const TwoDMatrix& vcov, Journal& journal);
void writeMat4(FILE* fd, const char* lname) const;
@@ -778,12 +781,11 @@ protected:@;
TwoDMatrix variances;
public:@;
SimResultsIRF(const SimResults& cntl, int ny, int nper, int i, double impulse)
- : SimResults(ny, nper), control(cntl),
+ : SimResults(ny, nper, 0), control(cntl),
ishock(i), imp(impulse),
means(ny, nper), variances(ny, nper)@+ {}
- void simulate(const DecisionRule& dr, const Vector& start,
- Journal& journal);
- void simulate(const DecisionRule& dr, const Vector& start);
+ void simulate(const DecisionRule& dr, Journal& journal);
+ void simulate(const DecisionRule& dr);
void writeMat4(FILE* fd, const char* lname) const;
protected:@;
void calcMeans();
@@ -804,13 +806,14 @@ protected:@;
Vector mean;
TwoDMatrix vcov;
int num_per;
+ int num_burn;
NormalConj nc;
int incomplete_simulations;
int thrown_periods;
public:@;
- RTSimResultsStats(int ny, int nper)
+ RTSimResultsStats(int ny, int nper, int nburn = 0)
: mean(ny), vcov(ny, ny),
- num_per(nper), nc(ny),
+ num_per(nper), num_burn(nburn), nc(ny),
incomplete_simulations(0), thrown_periods(0)@+ {}
void simulate(int num_sim, const DecisionRule& dr, const Vector& start,
const TwoDMatrix& vcov, Journal& journal);
@@ -879,7 +882,6 @@ class SimulationIRFWorker : public THREAD {
const DecisionRule& dr;
DecisionRule::emethod em;
int np;
- const Vector& st;
int idata;
int ishock;
double imp;
@@ -887,9 +889,8 @@ public:@;
SimulationIRFWorker(SimResultsIRF& sim_res,
const DecisionRule& dec_rule,
DecisionRule::emethod emet, int num_per,
- const Vector& start, int id,
- int ishck, double impulse)
- : res(sim_res), dr(dec_rule), em(emet), np(num_per), st(start),
+ int id, int ishck, double impulse)
+ : res(sim_res), dr(dec_rule), em(emet), np(num_per),
idata(id), ishock(ishck), imp(impulse)@+ {}
void operator()();
};
@@ -951,12 +952,16 @@ class ExplicitShockRealization : virtual public ShockRealization {
public:@;
ExplicitShockRealization(const TwoDMatrix& sh)
: shocks(sh)@+ {}
+ ExplicitShockRealization(const ConstTwoDMatrix& sh)
+ : shocks(sh)@+ {}
ExplicitShockRealization(const ExplicitShockRealization& sr)
: shocks(sr.shocks)@+ {}
ExplicitShockRealization(ShockRealization& sr, int num_per);
void get(int n, Vector& out);
int numShocks() const
{@+ return shocks.nrows();@+}
+ const TwoDMatrix& getShocks()
+ {@+ return shocks;@+}
void addToShock(int ishock, int iper, double val);
void print() const
{@+ shocks.print();@+}
diff --git a/dynare++/src/dynare_params.cpp b/dynare++/src/dynare_params.cpp
index 3acd42c79b7a4652954ff7b73ac39f800023f791..95a76a674c8c3b4c7e27b0f8dbfe38552fdc0c23 100644
--- a/dynare++/src/dynare_params.cpp
+++ b/dynare++/src/dynare_params.cpp
@@ -13,9 +13,10 @@ const char* help_str =
" --version print version and return\n"
"\n"
"options:\n"
-" --per <num> number of periods simulated [100]\n"
+" --per <num> number of periods simulated after burnt [100]\n"
+" --burn <num> number of periods burnt [0]\n"
" --sim <num> number of simulations [80]\n"
-" --rtper <num> number of RT periods simulated [0]\n"
+" --rtper <num> number of RT periods simulated after burnt [0]\n"
" --rtsim <num> number of RT simulations [0]\n"
" --condper <num> number of periods in cond. simulations [0]\n"
" --condsim <num> number of conditional simulations [0]\n"
@@ -45,7 +46,7 @@ const char* help_str =
const char* dyn_basename(const char* str);
DynareParams::DynareParams(int argc, char** argv)
- : modname(NULL), num_per(100), num_sim(80),
+ : modname(NULL), num_per(100), num_burn(0), num_sim(80),
num_rtper(0), num_rtsim(0),
num_condper(0), num_condsim(0),
num_threads(2), num_steps(0),
@@ -70,6 +71,7 @@ DynareParams::DynareParams(int argc, char** argv)
struct option const opts [] = {
{"periods", required_argument, NULL, opt_per},
{"per", required_argument, NULL, opt_per},
+ {"burn", required_argument, NULL, opt_burn},
{"simulations", required_argument, NULL, opt_sim},
{"sim", required_argument, NULL, opt_sim},
{"rtperiods", required_argument, NULL, opt_rtper},
@@ -108,6 +110,10 @@ DynareParams::DynareParams(int argc, char** argv)
if (1 != sscanf(optarg, "%d", &num_per))
fprintf(stderr, "Couldn't parse integer %s, ignored\n", optarg);
break;
+ case opt_burn:
+ if (1 != sscanf(optarg, "%d", &num_burn))
+ fprintf(stderr, "Couldn't parse integer %s, ignored\n", optarg);
+ break;
case opt_sim:
if (1 != sscanf(optarg, "%d", &num_sim))
fprintf(stderr, "Couldn't parse integer %s, ignored\n", optarg);
diff --git a/dynare++/src/dynare_params.h b/dynare++/src/dynare_params.h
index 3c6a65adc833d47327cd4b5250ff7f2e2623b9cb..594c35ea278a514c0444aea2e08c62a2cd083198 100644
--- a/dynare++/src/dynare_params.h
+++ b/dynare++/src/dynare_params.h
@@ -17,6 +17,7 @@ struct DynareParams {
const char* modname;
std::string basename;
int num_per;
+ int num_burn;
int num_sim;
int num_rtper;
int num_rtsim;
@@ -56,7 +57,8 @@ struct DynareParams {
int getCheckPathPoints() const
{return 10*check_num;}
private:
- enum {opt_per, opt_sim, opt_rtper, opt_rtsim, opt_condper, opt_condsim, opt_prefix, opt_threads,
+ enum {opt_per, opt_burn, opt_sim, opt_rtper, opt_rtsim, opt_condper, opt_condsim,
+ opt_prefix, opt_threads,
opt_steps, opt_seed, opt_order, opt_ss_tol, opt_check,
opt_check_along_path, opt_check_along_shocks, opt_check_on_ellipse,
opt_check_evals, opt_check_scale, opt_check_num, opt_noirfs, opt_irfs,
diff --git a/dynare++/src/main.cpp b/dynare++/src/main.cpp
index fa4eb01dd10e11386b41282a59a1f79077ffc8f1..5ac5d15d1897794a5b2200f2e33e9dc7718fee05 100644
--- a/dynare++/src/main.cpp
+++ b/dynare++/src/main.cpp
@@ -130,7 +130,7 @@ int main(int argc, char** argv)
// simulate conditional
if (params.num_condper > 0 && params.num_condsim > 0) {
- SimResultsDynamicStats rescond(dynare.numeq(), params.num_condper);
+ SimResultsDynamicStats rescond(dynare.numeq(), params.num_condper, 0);
ConstVector det_ss(app.getSS(),0);
rescond.simulate(params.num_condsim, app.getFoldDecisionRule(), det_ss, dynare.getVcov(), journal);
rescond.writeMat4(matfd, params.prefix);
@@ -140,7 +140,7 @@ int main(int argc, char** argv)
//const DecisionRule& dr = app.getUnfoldDecisionRule();
const DecisionRule& dr = app.getFoldDecisionRule();
if (params.num_per > 0 && params.num_sim > 0) {
- SimResultsStats res(dynare.numeq(), params.num_per);
+ SimResultsStats res(dynare.numeq(), params.num_per, params.num_burn);
res.simulate(params.num_sim, dr, dynare.getSteady(), dynare.getVcov(), journal);
res.writeMat4(matfd, params.prefix);
@@ -153,7 +153,7 @@ int main(int argc, char** argv)
// simulate with real-time statistics
if (params.num_rtper > 0 && params.num_rtsim > 0) {
- RTSimResultsStats rtres(dynare.numeq(), params.num_rtper);
+ RTSimResultsStats rtres(dynare.numeq(), params.num_rtper, params.num_burn);
rtres.simulate(params.num_rtsim, dr, dynare.getSteady(), dynare.getVcov(), journal);
rtres.writeMat4(matfd, params.prefix);
}