stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Implemented burn-in periods in simulations, irfs, and in real-time simulations.

time-shift
Ondra Kamenik 2011-03-11 23:34:15 +01:00 committed by Sébastien Villemot
parent f48aa28e39
commit b05c7356ee
5 changed files with 64 additions and 45 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

50
dynare++/kord/decision_rule.cweb
View File

@ -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

37
dynare++/kord/decision_rule.hweb
View File

@ -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();@+}

12
dynare++/src/dynare_params.cpp
View File

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

4
dynare++/src/dynare_params.h
View File

@ -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,

6
dynare++/src/main.cpp
View File

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