@q $Id: decision_rule.cweb 1896 2008-06-24 04:01:05Z kamenik $ @>
@q Copyright 2004, Ondra Kamenik @>

@ Start of {\tt decision\_rule.cpp} file.
@c

#include "kord_exception.h"
#include "decision_rule.h"
#include "dynamic_model.h"

#include "SymSchurDecomp.h"

#include <dynlapack.h>

#include <limits>

template <>
int DRFixPoint<KOrder::fold>::max_iter = 10000;
template <>
int DRFixPoint<KOrder::unfold>::max_iter = 10000;
template <>
double DRFixPoint<KOrder::fold>::tol = 1.e-10;
template <>
double DRFixPoint<KOrder::unfold>::tol = 1.e-10;
template <>
int DRFixPoint<KOrder::fold>::max_newton_iter = 50;
template <>
int DRFixPoint<KOrder::unfold>::max_newton_iter = 50;
template <>
int DRFixPoint<KOrder::fold>::newton_pause = 100;
template <>
int DRFixPoint<KOrder::unfold>::newton_pause = 100;
@#
@<|FoldDecisionRule| conversion from |UnfoldDecisionRule|@>;
@<|UnfoldDecisionRule| conversion from |FoldDecisionRule|@>;
@<|SimResults| destructor@>;
@<|SimResults::simulate| code1@>;
@<|SimResults::simulate| code2@>;
@<|SimResults::addDataSet| code@>;
@<|SimResults::writeMat4| code1@>;
@<|SimResults::writeMat4| code2@>;
@<|SimResultsStats::simulate| code@>;
@<|SimResultsStats::writeMat4| code@>;
@<|SimResultsStats::calcMean| code@>;
@<|SimResultsStats::calcVcov| code@>;
@<|SimResultsDynamicStats::simulate| code@>;
@<|SimResultsDynamicStats::writeMat4| code@>;
@<|SimResultsDynamicStats::calcMean| code@>;
@<|SimResultsDynamicStats::calcVariance| code@>;
@<|SimResultsIRF::simulate| code1@>;
@<|SimResultsIRF::simulate| code2@>;
@<|SimResultsIRF::calcMeans| code@>;
@<|SimResultsIRF::calcVariances| code@>;
@<|SimResultsIRF::writeMat4| code@>;
@<|RTSimResultsStats::simulate| code1@>;
@<|RTSimResultsStats::simulate| code2@>;
@<|RTSimResultsStats::writeMat4| code@>;
@<|IRFResults| constructor@>;
@<|IRFResults| destructor@>;
@<|IRFResults::writeMat4| code@>;
@<|SimulationWorker::operator()()| code@>;
@<|SimulationIRFWorker::operator()()| code@>;
@<|RTSimulationWorker::operator()()| code@>;
@<|RandomShockRealization::choleskyFactor| code@>;
@<|RandomShockRealization::schurFactor| code@>;
@<|RandomShockRealization::get| code@>;
@<|ExplicitShockRealization| constructor code@>;
@<|ExplicitShockRealization::get| code@>;
@<|ExplicitShockRealization::addToShock| code@>;
@<|GenShockRealization::get| code@>;

@ 
@<|FoldDecisionRule| conversion from |UnfoldDecisionRule|@>=
FoldDecisionRule::FoldDecisionRule(const UnfoldDecisionRule& udr)
	: DecisionRuleImpl<KOrder::fold>(ctraits<KOrder::fold>::Tpol(udr.nrows(), udr.nvars()),
									 udr.ypart, udr.nu, udr.ysteady)
{
	for (ctraits<KOrder::unfold>::Tpol::const_iterator it = udr.begin();
		 it != udr.end(); ++it) {
		insert(new ctraits<KOrder::fold>::Ttensym(*((*it).second)));
	}
}

@ 
@<|UnfoldDecisionRule| conversion from |FoldDecisionRule|@>=
UnfoldDecisionRule::UnfoldDecisionRule(const FoldDecisionRule& fdr)
	: DecisionRuleImpl<KOrder::unfold>(ctraits<KOrder::unfold>::Tpol(fdr.nrows(), fdr.nvars()),
									 fdr.ypart, fdr.nu, fdr.ysteady)
{
	for (ctraits<KOrder::fold>::Tpol::const_iterator it = fdr.begin();
		 it != fdr.end(); ++it) {
		insert(new ctraits<KOrder::unfold>::Ttensym(*((*it).second)));
	}
}

@ 
@<|SimResults| destructor@>=
SimResults::~SimResults()
{
	for (int i = 0; i < getNumSets(); i++) {
		delete data[i];
		delete shocks[i];
	}
}

@ This runs simulations with an output to journal file. Note that we
report how many simulations had to be thrown out due to Nan or Inf.

@<|SimResults::simulate| code1@>=
void SimResults::simulate(int num_sim, const DecisionRule& dr, const Vector& start,
						  const TwoDMatrix& vcov, Journal& journal)
{
	JournalRecordPair paa(journal);
	paa << "Performing " << num_sim << " stochastic simulations for "
		<< num_per << " periods" << endrec;
	simulate(num_sim, dr, start, vcov);
	int thrown = num_sim - data.size();
	if (thrown > 0) {
		JournalRecord rec(journal);
		rec << "I had to throw " << thrown << " simulations away due to Nan or Inf" << endrec;
	}
}

@ This runs a given number of simulations by creating
|SimulationWorker| for each simulation and inserting them to the
thread group.

@<|SimResults::simulate| code2@>=
void SimResults::simulate(int num_sim, const DecisionRule& dr, const Vector& start,
						  const TwoDMatrix& vcov)
{
	std::vector<RandomShockRealization> rsrs;
	rsrs.reserve(num_sim);

	THREAD_GROUP gr;
	for (int i = 0; i < num_sim; i++) {
		RandomShockRealization sr(vcov, system_random_generator.int_uniform());
		rsrs.push_back(sr);
		THREAD* worker = new
			SimulationWorker(*this, dr, DecisionRule::horner,
							 num_per, 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.

@<|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,
				  "Incompatible number of cols for SimResults::addDataSets");
	if (d->isFinite()) {
		data.push_back(d);
		shocks.push_back(sr);
		return true;
	} else {
		delete d;
		delete sr;
		return false;
	}
}

@ 
@<|SimResults::writeMat4| code1@>=
void SimResults::writeMat4(const char* base, const char* lname) const
{
	char matfile_name[100];
	sprintf(matfile_name, "%s.mat", base);
	FILE* out;
	if (NULL != (out=fopen(matfile_name, "wb"))) {
		writeMat4(out, lname);
		fclose(out);
	}
}

@ This save the results as matrices with given prefix and with index
appended. If there is only one matrix, the index is not appended.

@<|SimResults::writeMat4| code2@>=
void SimResults::writeMat4(FILE* fd, const char* lname) const
{
	char tmp[100];
	for (int i = 0; i < getNumSets(); i++) {
		if (getNumSets() > 1)
			sprintf(tmp, "%s_data%d", lname, i+1);
		else
			sprintf(tmp, "%s_data", lname);
		ConstTwoDMatrix m(*(data[i]));
		m.writeMat4(fd, tmp);
	}
}

@ 
@<|SimResultsStats::simulate| code@>=
void SimResultsStats::simulate(int num_sim, const DecisionRule& dr,
							   const Vector& start,
							   const TwoDMatrix& vcov, Journal& journal)
{
	SimResults::simulate(num_sim, dr, start, vcov, journal);
	{
		JournalRecordPair paa(journal);
		paa << "Calculating means from the simulations." << endrec;
		calcMean();
	}
	{
		JournalRecordPair paa(journal);
		paa << "Calculating covariances from the simulations." << endrec;
		calcVcov();
	}
}


@ Here we do not save the data itself, we save only mean and vcov.
@<|SimResultsStats::writeMat4| code@>=
void SimResultsStats::writeMat4(FILE* fd, const char* lname) const
{
	char tmp[100];
	sprintf(tmp, "%s_mean", lname);
	ConstTwoDMatrix m(num_y, 1, mean.base());
	m.writeMat4(fd, tmp);
	sprintf(tmp, "%s_vcov", lname);
	ConstTwoDMatrix(vcov).writeMat4(fd, tmp);
}

@ 
@<|SimResultsStats::calcMean| code@>=
void SimResultsStats::calcMean()
{
	mean.zeros();
	if (data.size()*num_per > 0) {
		double mult = 1.0/data.size()/num_per;
		for (unsigned int i = 0; i < data.size(); i++) {
			for (int j = 0; j < num_per; j++) {
				ConstVector col(*data[i], j);
				mean.add(mult, col);
			}
		}
	}
}

@ 
@<|SimResultsStats::calcVcov| code@>=
void SimResultsStats::calcVcov()
{
	if (data.size()*num_per > 1) {
		vcov.zeros();
		double mult = 1.0/(data.size()*num_per - 1);
		for (unsigned int i = 0; i < data.size(); i++) {
			const TwoDMatrix& d = *(data[i]);
			for (int j = 0; j < num_per; j++) {
				for (int m = 0; m < num_y; m++) {
					for (int n = m; n < num_y; n++) {
						double s = (d.get(m,j)-mean[m])*(d.get(n,j)-mean[n]);
						vcov.get(m,n) += mult*s;
						if (m != n)
							vcov.get(n,m) += mult*s;
					}
				}
			}
		}
	} else {
		vcov.infs();
	}
}

@ 
@<|SimResultsDynamicStats::simulate| code@>=
void SimResultsDynamicStats::simulate(int num_sim, const DecisionRule& dr,
									  const Vector& start,
									  const TwoDMatrix& vcov, Journal& journal)
{
	SimResults::simulate(num_sim, dr, start, vcov, journal);
	{
		JournalRecordPair paa(journal);
		paa << "Calculating means of the conditional simulations." << endrec;
		calcMean();
	}
	{
		JournalRecordPair paa(journal);
		paa << "Calculating variances of the conditional simulations." << endrec;
		calcVariance();
	}
}

@ 
@<|SimResultsDynamicStats::writeMat4| code@>=
void SimResultsDynamicStats::writeMat4(FILE* fd, const char* lname) const
{
	char tmp[100];
	sprintf(tmp, "%s_cond_mean", lname);
	ConstTwoDMatrix(mean).writeMat4(fd, tmp);
	sprintf(tmp, "%s_cond_variance", lname);
	ConstTwoDMatrix(variance).writeMat4(fd, tmp);
}

@ 
@<|SimResultsDynamicStats::calcMean| code@>=
void SimResultsDynamicStats::calcMean()
{
	mean.zeros();
	if (data.size() > 0) {
		double mult = 1.0/data.size();
		for (int j = 0; j < num_per; j++) {
			Vector meanj(mean, j);
			for (unsigned int i = 0; i < data.size(); i++) {
				ConstVector col(*data[i], j);
				meanj.add(mult, col);
			}
		}
	}
}

@ 
@<|SimResultsDynamicStats::calcVariance| code@>=
void SimResultsDynamicStats::calcVariance()
{
	if (data.size() > 1) {
		variance.zeros();
		double mult = 1.0/(data.size()-1);
		for (int j = 0; j < num_per; j++) {
			ConstVector meanj(mean, j);
			Vector varj(variance, j);
			for (int i = 0; i < (int)data.size(); i++) {
				Vector col(ConstVector((*data[i]), j));
				col.add(-1.0, meanj);
				for (int k = 0; k < col.length(); k++)
					col[k] = col[k]*col[k];
				varj.add(mult, col);
			}
		}
	} else {
		variance.infs();
	}
}


@ 
@<|SimResultsIRF::simulate| code1@>=
void SimResultsIRF::simulate(const DecisionRule& dr, const Vector& start,
							 Journal& journal)
{
	JournalRecordPair paa(journal);
	paa << "Performing " << control.getNumSets() << " IRF simulations for "
		<< num_per << " periods; shock=" << ishock << ", impulse=" << imp << endrec;
	simulate(dr, start);
	int thrown = control.getNumSets() - data.size();
	if (thrown > 0) {
		JournalRecord rec(journal);
		rec << "I had to throw " << thrown
			<< " simulations away due to Nan or Inf" << endrec;
	}	
	calcMeans();
	calcVariances();
}

@ 
@<|SimResultsIRF::simulate| code2@>=
void SimResultsIRF::simulate(const DecisionRule& dr, const Vector& start)
{
	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);
		gr.insert(worker);
	}
	gr.run();
}

@ 
@<|SimResultsIRF::calcMeans| code@>=
void SimResultsIRF::calcMeans()
{
	means.zeros();
	if (data.size() > 0) {
		for (unsigned int i = 0; i < data.size(); i++)
			means.add(1.0, *(data[i]));
		means.mult(1.0/data.size());
	}
}

@ 
@<|SimResultsIRF::calcVariances| code@>=
void SimResultsIRF::calcVariances()
{
	if (data.size() > 1) {
		variances.zeros();
		for (unsigned int i = 0; i < data.size(); i++) {
			TwoDMatrix d((const TwoDMatrix&)(*(data[i])));
			d.add(-1.0, means);
			for (int j = 0; j < d.nrows(); j++)
				for (int k = 0;	k < d.ncols(); k++)
					variances.get(j,k) += d.get(j,k)*d.get(j,k);
			d.mult(1.0/(data.size()-1));
		}
	} else {
		variances.infs();
	}
}

@ 
@<|SimResultsIRF::writeMat4| code@>=
void SimResultsIRF::writeMat4(FILE* fd, const char* lname) const
{
	char tmp[100];
	sprintf(tmp, "%s_mean", lname);
	means.writeMat4(fd, tmp);
	sprintf(tmp, "%s_var", lname);
	variances.writeMat4(fd, tmp);
}

@ 
@<|RTSimResultsStats::simulate| code1@>=
void RTSimResultsStats::simulate(int num_sim, const DecisionRule& dr, const Vector& start,
								 const TwoDMatrix& v, Journal& journal)
{
	JournalRecordPair paa(journal);
	paa << "Performing " << num_sim << " real-time stochastic simulations for "
		<< num_per << " periods" << endrec;
	simulate(num_sim, dr, start, v);
	mean = nc.getMean();
	mean.add(1.0, dr.getSteady());
	nc.getVariance(vcov);
	if (thrown_periods > 0) {
		JournalRecord rec(journal);
		rec << "I had to throw " << thrown_periods << " periods away due to Nan or Inf" << endrec;
		JournalRecord rec1(journal);
		rec1 << "This affected " << incomplete_simulations << " out of "
			 << num_sim << " simulations" << endrec;
	}
}

@ 
@<|RTSimResultsStats::simulate| code2@>=
void RTSimResultsStats::simulate(int num_sim, const DecisionRule& dr, const Vector& start,
								 const TwoDMatrix& vcov)
{
	std::vector<RandomShockRealization> rsrs;
	rsrs.reserve(num_sim);

	THREAD_GROUP gr;
	for (int i = 0; i < num_sim; i++) {
		RandomShockRealization sr(vcov, system_random_generator.int_uniform());
		rsrs.push_back(sr);
		THREAD* worker = new
			RTSimulationWorker(*this, dr, DecisionRule::horner,
							   num_per, start, rsrs.back());
		gr.insert(worker);
	}
	gr.run();
}

@ 
@<|RTSimResultsStats::writeMat4| code@>=
void RTSimResultsStats::writeMat4(FILE* fd, const char* lname)
{
	char tmp[100];
	sprintf(tmp, "%s_rt_mean", lname);
	ConstTwoDMatrix m(nc.getDim(), 1, mean.base());
	m.writeMat4(fd, tmp);
	sprintf(tmp, "%s_rt_vcov", lname);
	ConstTwoDMatrix(vcov).writeMat4(fd, tmp);
}

@ 
@<|IRFResults| constructor@>=
IRFResults::IRFResults(const DynamicModel& mod, const DecisionRule& dr,
					   const SimResults& control, const vector<int>& ili,
					   Journal& journal)
	: model(mod), irf_list_ind(ili)
{
	int num_per = control.getNumPer();
	JournalRecordPair pa(journal);
	pa << "Calculating IRFs against control for " << (int)irf_list_ind.size() << " shocks and for "
	   << num_per << " periods" << endrec;
	const TwoDMatrix& vcov = mod.getVcov();
	for (unsigned int ii = 0; ii < irf_list_ind.size(); ii++) {
		int ishock = irf_list_ind[ii];
		double stderror = sqrt(vcov.get(ishock,ishock));
		irf_res.push_back(new SimResultsIRF(control, model.numeq(), num_per,
											ishock, stderror));
		irf_res.push_back(new SimResultsIRF(control, model.numeq(), num_per,
											ishock, -stderror));
	}

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

@ 
@<|IRFResults| destructor@>=
IRFResults::~IRFResults()
{
	for (unsigned int i = 0; i < irf_res.size(); i++)
		delete irf_res[i];
}

@ 
@<|IRFResults::writeMat4| code@>=
void IRFResults::writeMat4(FILE* fd, const char* prefix) const
{
	for (unsigned int i = 0; i < irf_list_ind.size(); i++) {
		char tmp[100];
		int ishock = irf_list_ind[i];
		const char* shockname = model.getExogNames().getName(ishock);
		sprintf(tmp, "%s_irfp_%s", prefix, shockname);
		irf_res[2*i]->writeMat4(fd, tmp);
		sprintf(tmp, "%s_irfm_%s", prefix, shockname);
		irf_res[2*i+1]->writeMat4(fd, tmp);
	}
}

@ 
@<|SimulationWorker::operator()()| code@>=
void SimulationWorker::operator()()
{
	ExplicitShockRealization* esr = new ExplicitShockRealization(sr, np);
	TwoDMatrix* m = dr.simulate(em, np, st, *esr);
	{
		SYNCHRO syn(&res, "simulation");
		res.addDataSet(m, esr);
	}
}

@ Here we create a new instance of |ExplicitShockRealization| of the
corresponding control, add the impulse, and simulate.

@<|SimulationIRFWorker::operator()()| code@>=
void SimulationIRFWorker::operator()()
{
	ExplicitShockRealization* esr =
	    new ExplicitShockRealization(res.control.getShocks(idata));
	esr->addToShock(ishock, 0, imp);
	TwoDMatrix* m = dr.simulate(em, np, st, *esr);
	m->add(-1.0, res.control.getData(idata));
	{
		SYNCHRO syn(&res, "simulation");
		res.addDataSet(m, esr);
	}
}

@ 
@<|RTSimulationWorker::operator()()| code@>=
void RTSimulationWorker::operator()()
{
	NormalConj nc(res.nc.getDim());
	const PartitionY& ypart = dr.getYPart();
	int nu = dr.nexog();
	const Vector& ysteady = dr.getSteady();

	@<initialize vectors and subvectors for simulation@>;
	@<simulate the first real-time period@>;
    @<simulate other real-time periods@>;
	{
		SYNCHRO syn(&res, "rtsimulation");
		res.nc.update(nc);
		if (res.num_per-ip > 0) {
			res.incomplete_simulations++;
			res.thrown_periods += res.num_per-ip;
		}		
	}
}

@ 
@<initialize vectors and subvectors for simulation@>=
	Vector dyu(ypart.nys()+nu);
	ConstVector ystart_pred(ystart, ypart.nstat, ypart.nys());
	ConstVector ysteady_pred(ysteady, ypart.nstat, ypart.nys());
	Vector dy(dyu, 0, ypart.nys());
	Vector u(dyu, ypart.nys(), nu);
	Vector y(nc.getDim());
	ConstVector ypred(y, ypart.nstat, ypart.nys());

@ 
@<simulate the first real-time period@>=
	int ip = 0;
	dy = ystart_pred;
	dy.add(-1.0, ysteady_pred);
	sr.get(ip, u);
	dr.eval(em, y, dyu);
	nc.update(y);

@
@<simulate other real-time periods@>=
while (y.isFinite() && ip < res.num_per) {
	ip++;
	dy = ypred;
	sr.get(ip, u);
	dr.eval(em, y, dyu);
	nc.update(y);
}

@ This calculates factorization $FF^T=V$ in the Cholesky way. It does
not work for semidefinite matrices.
 
@<|RandomShockRealization::choleskyFactor| code@>=
void RandomShockRealization::choleskyFactor(const TwoDMatrix& v)
{
	factor = v;
	lapack_int rows = factor.nrows();
	for (int i = 0; i < rows; i++)
		for (int j = i+1; j < rows; j++)
			factor.get(i,j) = 0.0;
	lapack_int info;

	dpotrf("L", &rows, factor.base(), &rows, &info);
	KORD_RAISE_IF(info != 0,
				  "Info!=0 in RandomShockRealization::choleskyFactor");
}

@ This calculates $FF^T=V$ factorization by symmetric Schur
decomposition. It works for semidifinite matrices.
 
@<|RandomShockRealization::schurFactor| code@>=
void RandomShockRealization::schurFactor(const TwoDMatrix& v)
{
	SymSchurDecomp ssd(v);
	ssd.getFactor(factor);
}

@ 
@<|RandomShockRealization::get| code@>=
void RandomShockRealization::get(int n, Vector& out)
{
 	KORD_RAISE_IF(out.length() != numShocks(),
				  "Wrong length of out vector in RandomShockRealization::get");
	Vector d(out.length());
	for (int i = 0; i < d.length(); i++) {
		d[i] = mtwister.normal();
	}
	out.zeros();
	factor.multaVec(out, ConstVector(d));
}

@ 
@<|ExplicitShockRealization| constructor code@>=
ExplicitShockRealization::ExplicitShockRealization(ShockRealization& sr,
												   int num_per)
	: shocks(sr.numShocks(), num_per)
{
	for (int j = 0; j < num_per; j++) {
		Vector jcol(shocks, j);
		sr.get(j, jcol);
	}
}

@ 
@<|ExplicitShockRealization::get| code@>=
void ExplicitShockRealization::get(int n, Vector& out)
{
 	KORD_RAISE_IF(out.length() != numShocks(),
				  "Wrong length of out vector in ExplicitShockRealization::get");
	int i = n % shocks.ncols();
	ConstVector icol(shocks, i);
	out = icol;
}

@ 
@<|ExplicitShockRealization::addToShock| code@>=
void ExplicitShockRealization::addToShock(int ishock, int iper, double val)
{
	KORD_RAISE_IF(ishock < 0 || ishock > numShocks(),
				  "Wrong index of shock in ExplicitShockRealization::addToShock");
	int j = iper % shocks.ncols();
	shocks.get(ishock, j) += val;
}


@ 
@<|GenShockRealization::get| code@>=
void GenShockRealization::get(int n, Vector& out)
{
	KORD_RAISE_IF(out.length() != numShocks(),
				  "Wrong length of out vector in GenShockRealization::get");
	ExplicitShockRealization::get(n, out);
	Vector r(numShocks());
	RandomShockRealization::get(n, r);
	for (int j = 0; j < numShocks(); j++)
#ifndef _MSC_VER
		if (! isfinite(out[j]))
#else
		if (!_finite(out[j]))
#endif
			out[j] = r[j];
}


@ End of {\tt decision\_rule.cpp} file.