dynare++/integ/cc/vector_function.cweb

-@q $Id: vector_function.cweb 431 2005-08-16 15:41:01Z kamenik $ @>
-@q Copyright 2005, Ondra Kamenik @>
-
-@ This is {\tt vector\_function.cpp} file
-
-@c
-
-#include "vector_function.h"
-
-#include <dynlapack.h>
-
-#include <cmath>
-
-#include <cstring>
-#include <algorithm>
-
-@<|ParameterSignal| constructor code@>;
-@<|ParameterSignal| copy constructor code@>;
-@<|ParameterSignal::signalAfter| code@>;
-@<|VectorFunctionSet| constructor 1 code@>;
-@<|VectorFunctionSet| constructor 2 code@>;
-@<|VectorFunctionSet| destructor code@>;
-@<|GaussConverterFunction| constructor code 1@>;
-@<|GaussConverterFunction| constructor code 2@>;
-@<|GaussConverterFunction| copy constructor code@>;
-@<|GaussConverterFunction::eval| code@>;
-@<|GaussConverterFunction::multiplier| code@>;
-@<|GaussConverterFunction::calcCholeskyFactor| code@>;
-
-@ Just an easy constructor of sequence of booleans defaulting to
-change everywhere.
-
-@<|ParameterSignal| constructor code@>=
-ParameterSignal::ParameterSignal(int n)
-	: data(new bool[n]), num(n)
-{
-	for (int i = 0; i < num; i++)
-		data[i] = true;
-}
-
-@ 
-@<|ParameterSignal| copy constructor code@>=
-ParameterSignal::ParameterSignal(const ParameterSignal& sig)
-	: data(new bool[sig.num]), num(sig.num)
-{
-	memcpy(data, sig.data, num);
-}
-
-@ This sets |false| (no change) before a given parameter, and |true|
-(change) after the given parameter (including).
-
-@<|ParameterSignal::signalAfter| code@>=
-void ParameterSignal::signalAfter(int l)
-{
-	for (int i = 0; i < std::min(l,num); i++)
-		data[i] = false;
-	for (int i = l; i < num; i++)
-		data[i] = true;
-}
-
-@ This constructs a function set hardcopying also the first.
-@<|VectorFunctionSet| constructor 1 code@>=
-VectorFunctionSet::VectorFunctionSet(const VectorFunction& f, int n)
-	: funcs(n), first_shallow(false)
-{
-	for (int i = 0; i < n; i++)
-		funcs[i] = f.clone();
-}
-
-@ This constructs a function set with shallow copy in the first and
-hard copies in others.
-
-@<|VectorFunctionSet| constructor 2 code@>=
-VectorFunctionSet::VectorFunctionSet(VectorFunction& f, int n)
-	: funcs(n), first_shallow(true)
-{
-	if (n > 0)
-		funcs[0] = &f;
-	for (int i = 1; i < n; i++)
-		funcs[i] = f.clone();
-}
-
-@ This deletes the functions. The first is deleted only if it was not
-a shallow copy.
-
-@<|VectorFunctionSet| destructor code@>=
-VectorFunctionSet::~VectorFunctionSet()
-{
-	unsigned int start = first_shallow ? 1 : 0;
-	for (unsigned int i = start; i < funcs.size(); i++)
-		delete funcs[i];
-}
-
-@ Here we construct the object from the given function $f$ and given
-variance-covariance matrix $\Sigma=$|vcov|. The matrix $A$ is
-calculated as lower triangular and yields $\Sigma=AA^T$.
-
-@<|GaussConverterFunction| constructor code 1@>=
-GaussConverterFunction::GaussConverterFunction(VectorFunction& f, const GeneralMatrix& vcov)
-	: VectorFunction(f), func(&f), delete_flag(false), A(vcov.numRows(), vcov.numRows()),
-	  multiplier(calcMultiplier()) 
-{
-	// todo: raise if |A.numRows() != indim()|
-	calcCholeskyFactor(vcov);
-}
-
-@ Here we construct the object in the same way, however we mark the
-function as to be deleted.
-
-@<|GaussConverterFunction| constructor code 2@>=
-GaussConverterFunction::GaussConverterFunction(VectorFunction* f, const GeneralMatrix& vcov)
-	: VectorFunction(*f), func(f), delete_flag(true), A(vcov.numRows(), vcov.numRows()),
-	  multiplier(calcMultiplier()) 
-{
-	// todo: raise if |A.numRows() != indim()|
-	calcCholeskyFactor(vcov);
-}
-
-
-@ 
-@<|GaussConverterFunction| copy constructor code@>=
-GaussConverterFunction::GaussConverterFunction(const GaussConverterFunction& f)
-	: VectorFunction(f), func(f.func->clone()), delete_flag(true), A(f.A),
-	  multiplier(f.multiplier)
-{
-} 
-
-@ Here we evaluate the function
-$g(y)={1\over\sqrt{\pi^n}}f\left(\sqrt{2}Ay\right)$. Since the matrix $A$ is lower
-triangular, the change signal for the function $f$ will look like
-$(0,\ldots,0,1,\ldots,1)$ where the first $1$ is in the same position
-as the first change in the given signal |sig| of the input
-$y=$|point|.
-
-@<|GaussConverterFunction::eval| code@>=
-void GaussConverterFunction::eval(const Vector& point, const ParameterSignal& sig, Vector& out)
-{
-	ParameterSignal s(sig);
-	int i = 0;
-	while (i < indim() && !sig[i])
-		i++;
-	s.signalAfter(i);
-
-	Vector x(indim());
-	x.zeros();
-	A.multaVec(x, point);
-	x.mult(sqrt(2.0));
-
-	func->eval(x, s, out);
-
-	out.mult(multiplier);
-}
-
-@ This returns $1\over\sqrt{\pi^n}$.
-@<|GaussConverterFunction::multiplier| code@>=
-double GaussConverterFunction::calcMultiplier() const
-{
-	return sqrt(pow(M_PI, -1*indim()));
-}
-
-@ 
-@<|GaussConverterFunction::calcCholeskyFactor| code@>=
-void GaussConverterFunction::calcCholeskyFactor(const GeneralMatrix& vcov)
-{
-	A = vcov;
-
-	lapack_int rows = A.numRows();
-	for (int i = 0; i < rows; i++)
-		for (int j = i+1; j < rows; j++)
-			A.get(i,j) = 0.0;
-
-	lapack_int info;
-	dpotrf("L", &rows, A.base(), &rows, &info);
-	// todo: raise if |info!=1|
-}
-
-
-@ End of {\tt vector\_function.cpp} file

dynare++/integ/cc/vector_function.hweb

-@q $Id: vector_function.hweb 431 2005-08-16 15:41:01Z kamenik $ @>
-@q Copyright 2005, Ondra Kamenik @>
-
-@*2 Vector function. This is {\tt vector\_function.h} file
-
-This file defines interface for functions taking a vector as an input
-and returning a vector (with a different size) as an output. We are
-also introducing a parameter signalling; it is a boolean vector which
-tracks parameters which were changed from the previous call. The
-|VectorFunction| implementation can exploit this information and
-evaluate the function more efficiently. The information can be
-completely ignored.
-
-From the signalling reason, and from other reasons, the function
-evaluation is not |const|.
-
-@s ParameterSignal int
-@s VectorFunction int
-@s VectorFunctionSet int
-@s GaussConverterFunction int
-
-@c
-#ifndef VECTOR_FUNCTION_H
-#define VECTOR_FUNCTION_H
-
-#include "Vector.h"
-#include "GeneralMatrix.h"
-
-#include <vector>
-
-@<|ParameterSignal| class declaration@>;
-@<|VectorFunction| class declaration@>;
-@<|VectorFunctionSet| class declaration@>;
-@<|GaussConverterFunction| class declaration@>;
-
-#endif
-
-@ This is a simple class representing a vector of booleans. The items
-night be retrieved or changed, or can be set |true| after some
-point. This is useful when we multiply the vector with lower
-triangular matrix.
-
-|true| means that a parameter was changed.
-
-@<|ParameterSignal| class declaration@>=
-class ParameterSignal {
-protected:@;
-	bool* data;
-	int num;
-public:@;
-	ParameterSignal(int n);
-	ParameterSignal(const ParameterSignal& sig);
-	~ParameterSignal()
-		{@+ delete [] data;@+}
-	void signalAfter(int l);
-	const bool& operator[](int i) const
-		{@+ return data[i];@+} 
-	bool& operator[](int i)
-		{@+ return data[i];@+} 
-};
-
-@ This is the abstract class for vector function. At this level of
-abstraction we only need to know size of input vector and a size of
-output vector.
-
-The important thing here is a clone method, we will need to make hard
-copies of vector functions since the evaluations are not |const|. The
-hardcopies apply for parallelization.
-
-@<|VectorFunction| class declaration@>=
-class VectorFunction {
-protected:@;
-	int in_dim;
-	int out_dim;
-public:@;
-	VectorFunction(int idim, int odim)
-		: in_dim(idim), out_dim(odim)@+ {}
-	VectorFunction(const VectorFunction& func)
-		: in_dim(func.in_dim), out_dim(func.out_dim)@+ {}
-	virtual ~VectorFunction()@+ {}
-	virtual VectorFunction* clone() const =0;
-	virtual void eval(const Vector& point, const ParameterSignal& sig, Vector& out) =0;
-	int indim() const
-		{@+ return in_dim;@+}
-	int outdim() const
-		{@+ return out_dim;@+}
-};
-
-@ This makes |n| copies of |VectorFunction|. The first constructor
-make exactly |n| new copies, the second constructor copies only the
-pointer to the first and others are hard (real) copies.
-
-The class is useful for making a given number of copies at once, and
-this set can be reused many times if we need mupliple copis of the
-function (for example for paralelizing the code).
-
-@<|VectorFunctionSet| class declaration@>=
-class VectorFunctionSet {
-protected:@;
-	std::vector<VectorFunction*> funcs;
-	bool first_shallow;
-public:@;
-	VectorFunctionSet(const VectorFunction& f, int n);
-	VectorFunctionSet(VectorFunction& f, int n);
-	~VectorFunctionSet();
-	VectorFunction& getFunc(int i)
-		{@+ return *(funcs[i]);@+}
-	int getNum() const
-		{@+ return funcs.size(); @+}
-};
-
-@ This class wraps another |VectorFunction| to allow integration of a
-function through normally distributed inputs. Namely, if one wants to
-integrate
-$${1\over\sqrt{(2\pi)^n\vert\Sigma\vert}}\int f(x)e^{-{1\over2}x^T\Sigma^{-1}x}{\rm d}x$$
-then if we write $\Sigma=AA^T$ and $x=\sqrt{2}Ay$, we get integral
-$${1\over\sqrt{(2\pi)^n\vert\Sigma\vert}}
-\int f\left(\sqrt{2}Ay\right)e^{-y^Ty}\sqrt{2^n}\vert A\vert{\rm d}y=
-{1\over\sqrt{\pi^n}}\int f\left(\sqrt{2}Ay\right)e^{-y^Ty}{\rm d}y,$$
-which means that a given function $f$ we have to wrap to yield a function
-$$g(y)={1\over\sqrt{\pi^n}}f\left(\sqrt{2}Ay\right).$$
-This is exactly what this class is doing. This transformation is
-useful since the Gauss--Hermite points and weights are defined for
-weighting function $e^{-y^2}$, so this transformation allows using
-Gauss--Hermite quadratures seemlessly in a context of integration through
-normally distributed inputs.
-
-The class maintains a pointer to the function $f$. When the object is
-constructed by the first constructor, the $f$ is not copied. If the
-object of this class is copied, then $f$ is copied and we need to
-remember to destroy it in the desctructor; hence |delete_flag|. The
-second constructor takes a pointer to the function and differs from
-the first only by setting |delete_flag| to |true|.
-
-@<|GaussConverterFunction| class declaration@>=
-class GaussConverterFunction : public VectorFunction {
-protected:@;
-	VectorFunction* func;
-	bool delete_flag;
-	GeneralMatrix A;
-	double multiplier;
-public:@;
-	GaussConverterFunction(VectorFunction& f, const GeneralMatrix& vcov);
-	GaussConverterFunction(VectorFunction* f, const GeneralMatrix& vcov);
-	GaussConverterFunction(const GaussConverterFunction& f);
-	virtual ~GaussConverterFunction()
-		{@+ if (delete_flag) delete func; @+}
-	virtual VectorFunction* clone() const
-		{@+ return new GaussConverterFunction(*this);@+}
-	virtual void eval(const Vector& point, const ParameterSignal& sig, Vector& out);	
-private:@;
-	double calcMultiplier() const;
-	void calcCholeskyFactor(const GeneralMatrix& vcov);
-};
-
-@ End of {\tt vector\_function.h} file

dynare++/integ/src/Makefile.am

-noinst_PROGRAMS = quadrature-points
-
-quadrature_points_SOURCES = quadrature-points.cpp
-quadrature_points_CPPFLAGS = -I../.. -I../../sylv/cc -I../../integ/cc -I../../tl/cc
-quadrature_points_CXXFLAGS = $(PTHREAD_CFLAGS)
-quadrature_points_LDADD = ../cc/libinteg.a ../../tl/cc/libtl.a ../../parser/cc/libparser.a ../../sylv/cc/libsylv.a ../../utils/cc/libutils.a $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS) $(PTHREAD_LIBS)

dynare++/integ/src/quadrature-points.cpp

-// Copyright (C) 2008-2011, Ondra Kamenik
-
-#include "parser/cc/matrix_parser.h"
-#include "utils/cc/memory_file.h"
-#include "utils/cc/exception.h"
-#include "sylv/cc/GeneralMatrix.h"
-#include "sylv/cc/Vector.h"
-#include "sylv/cc/SymSchurDecomp.h"
-#include "sylv/cc/SylvException.h"
-#include "integ/cc/quadrature.h"
-#include "integ/cc/smolyak.h"
-#include "integ/cc/product.h"
-
-#include <getopt.h>
-#include <cstdio>
-
-#include <cmath>
-
-struct QuadParams {
-	const char* outname;
-	const char* vcovname;
-	int max_level;
-	double discard_weight;
-	QuadParams(int argc, char** argv);
-	void check_consistency() const;
-private:
-	enum {opt_max_level, opt_discard_weight, opt_vcov};
-};
-
-QuadParams::QuadParams(int argc, char** argv)
-	: outname(NULL), vcovname(NULL), max_level(3), discard_weight(0.0)
-{
-	if (argc == 1) {
-		// print the help and exit
-		exit(1);
-	}
-
-	outname = argv[argc-1];
-	argc--;
-
-	struct option const opts [] = {
-		{"max-level", required_argument, NULL, opt_max_level},
-		{"discard-weight", required_argument, NULL, opt_discard_weight},
-		{"vcov", required_argument, NULL, opt_vcov},
-		{NULL, 0, NULL, 0}
-	};
-
-	int ret;
-	int index;
-	while (-1 != (ret = getopt_long(argc, argv, "", opts, &index))) {
-		switch (ret) {
-		case opt_max_level:
-			if (1 != sscanf(optarg, "%d", &max_level))
-				fprintf(stderr, "Couldn't parse integer %s, ignored\n", optarg);
-			break;
-		case opt_discard_weight:
-			if (1 != sscanf(optarg, "%lf", &discard_weight))
-				fprintf(stderr, "Couldn't parse float %s, ignored\n", optarg);
-			break;
-		case opt_vcov:
-			vcovname = optarg;
-			break;
-		}
-	}
-
-	check_consistency();
-}
-
-void QuadParams::check_consistency() const
-{
-	if (outname == NULL) {
-		fprintf(stderr, "Error: output name not set\n");
-		exit(1);
-	}
-
-	if (vcovname == NULL) {
-		fprintf(stderr, "Error: vcov file name not set\n");
-		exit(1);
-	}
-}
-
-/** Utility class for ordering pointers to vectors according their
- * ordering. */
-struct OrderVec {
-	bool operator()(const Vector* a, const Vector* b) const
-		{return *a < *b;}
-};
-
-int main(int argc, char** argv)
-{
-	QuadParams params(argc, argv);
-
-	// open output file for writing
-	FILE* fout;
-	if (NULL == (fout=fopen(params.outname, "w"))) {
-		fprintf(stderr, "Could not open %s for writing\n", params.outname);
-		exit(1);
-	}
-
-	try {
-
-		// open memory file for vcov
-		ogu::MemoryFile vcov_mf(params.vcovname);
-	
-		// parse the vcov matrix
-		ogp::MatrixParser mp;
-		mp.parse(vcov_mf.length(), vcov_mf.base());
-		if (mp.nrows() != mp.ncols())
-			throw ogu::Exception(__FILE__,__LINE__,
-								 "VCOV matrix not square");
-		// and put to the GeneralMatrix
-		GeneralMatrix vcov(mp.nrows(), mp.ncols());
-		vcov.zeros();
-		for (ogp::MPIterator it = mp.begin(); it != mp.end(); ++it)
-			vcov.get(it.row(), it.col()) = *it;
-	
-		// calculate the factor A of vcov, so that A*A^T=VCOV
-		GeneralMatrix A(vcov.numRows(), vcov.numRows());
-		SymSchurDecomp ssd(vcov);
-		ssd.getFactor(A);
-
-		// construct Gauss-Hermite quadrature
-		GaussHermite ghq;
-		// construct Smolyak quadrature
-		int level = params.max_level;
-		SmolyakQuadrature sq(vcov.numRows(), level, ghq);
-
-		printf("Dimension:                %d\n", vcov.numRows());
-		printf("Maximum level:            %d\n", level);
-		printf("Total number of nodes:    %d\n", sq.numEvals(level));
-
-		// put the points to the vector
-		std::vector<Vector*> points;
-		for (smolpit qit = sq.start(level); qit != sq.end(level); ++qit)
-			points.push_back(new Vector((const Vector&)qit.point()));
-		// sort and uniq
-		OrderVec ordvec;
-		std::sort(points.begin(), points.end(), ordvec);
-		std::vector<Vector*>::iterator new_end = std::unique(points.begin(), points.end());
-		for (std::vector<Vector*>::iterator it = new_end; it != points.end(); ++it)
-			delete *it;
-		points.erase(new_end, points.end());
-
-		printf("Duplicit nodes removed:   %lu\n", (unsigned long) (sq.numEvals(level)-points.size()));
-
-		// calculate weights and mass
-		double mass = 0.0;
-		std::vector<double> weights;
-		for (int i = 0; i < (int)points.size(); i++) {
-			weights.push_back(std::exp(-points[i]->dot(*(points[i]))));
-			mass += weights.back();
-		}
-
-		// calculate discarded mass
-		double discard_mass = 0.0;
-		for (int i = 0; i < (int)weights.size(); i++)
-			if (weights[i]/mass < params.discard_weight)
-				discard_mass += weights[i];
-
-		printf("Total mass discarded:     %f\n", discard_mass/mass);
-	
-		// dump the results
-		int npoints = 0;
-		double upscale_weight = 1/(mass-discard_mass);
-		Vector x(vcov.numRows());
-		for (int i = 0; i < (int)weights.size(); i++)
-			if (weights[i]/mass >= params.discard_weight) {
-				// print the upscaled weight
-				fprintf(fout, "%20.16g", upscale_weight*weights[i]);
-				// multiply point with the factor A and sqrt(2)
-				A.multVec(0.0, x, std::sqrt(2.), *(points[i]));
-				// print the coordinates
-				for (int j = 0; j < x.length(); j++)
-					fprintf(fout, " %20.16g", x[j]);
-				fprintf(fout, "\n");
-				npoints++;
-			}
-
-		printf("Final number of points:   %d\n", npoints);
-
-		fclose(fout);
-
-	} catch (const SylvException& e) {
-		e.printMessage();
-		return 1;
-	} catch (const ogu::Exception& e) {
-		e.print();
-		return 1;
-	}
-
-	return 0;
-}

dynare++/integ/testing/Makefile.am

-check_PROGRAMS = tests
-
-tests_SOURCES = tests.cpp
-tests_CPPFLAGS = -I../cc -I../../tl/cc -I../../sylv/cc -I$(top_srcdir)/mex/sources
-tests_CXXFLAGS = $(PTHREAD_CFLAGS)
-tests_LDFLAGS = $(LDFLAGS_MATIO)
-tests_LDADD = ../../tl/cc/libtl.a ../../sylv/cc/libsylv.a ../cc/libinteg.a $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS) $(PTHREAD_LIBS) $(LIBADD_MATIO)
-
-check-local:
-	./tests

dynare++/integ/testing/tests.cpp

-/* $Id: tests.cpp 431 2005-08-16 15:41:01Z kamenik $ */
-/* Copyright 2005, Ondra Kamenik */
-
-#include "GeneralMatrix.h"
-#include <dynlapack.h>
-#include "SylvException.h"
-
-#include "rfs_tensor.h"
-#include "normal_moments.h"
-
-#include "vector_function.h"
-#include "quadrature.h"
-#include "smolyak.h"
-#include "product.h"
-#include "quasi_mcarlo.h"
-
-#include <cstdio>
-#include <cstring>
-#include <sys/time.h>
-#include <cmath>
-
-const int num_threads = 2; // does nothing if DEBUG defined
-
-// evaluates unfolded (Dx)^k power, where x is a vector, D is a
-// Cholesky factor (lower triangular)
-class MomentFunction : public VectorFunction {
-	GeneralMatrix D;
-	int k;
-public:
-	MomentFunction(const GeneralMatrix& inD, int kk)
-		: VectorFunction(inD.numRows(), UFSTensor::calcMaxOffset(inD.numRows(), kk)),
-		  D(inD), k(kk) {}
-	MomentFunction(const MomentFunction& func)
-		: VectorFunction(func), D(func.D), k(func.k) {}
-	VectorFunction* clone() const
-		{return new MomentFunction(*this);}
-	void eval(const Vector& point, const ParameterSignal& sig, Vector& out);
-};
-
-void MomentFunction::eval(const Vector& point, const ParameterSignal& sig, Vector& out)
-{
-	if (point.length() != indim() || out.length() != outdim()) {
-		printf("Wrong length of vectors in MomentFunction::eval\n");
-		exit(1);
-	}
-	Vector y(point);
-	y.zeros();
-	D.multaVec(y, point);
-	URSingleTensor ypow(y, k);
-	out.zeros();
-	out.add(1.0, ypow.getData());
-}
-
-class TensorPower : public VectorFunction {
-	int k;
-public:
-	TensorPower(int nvar, int kk)
-		: VectorFunction(nvar, UFSTensor::calcMaxOffset(nvar, kk)), k(kk) {}
-	TensorPower(const TensorPower& func)
-		: VectorFunction(func), k(func.k) {}
-	VectorFunction* clone() const
-		{return new TensorPower(*this);}
-	void eval(const Vector& point, const ParameterSignal& sig, Vector& out);
-};
-
-void TensorPower::eval(const Vector& point, const ParameterSignal& sig, Vector& out)
-{
-	if (point.length() != indim() || out.length() != outdim()) {
-		printf("Wrong length of vectors in TensorPower::eval\n");
-		exit(1);
-	}
-	URSingleTensor ypow(point, k);
-	out.zeros();
-	out.add(1.0, ypow.getData());
-}
-
-
-// evaluates (1+1/d)^d*(x_1*...*x_d)^(1/d), its integral over <0,1>^d
-// is 1.0, and its variation grows exponetially
-// d = dim
-class Function1 : public VectorFunction {
-	int dim;
-public:
-	Function1(int d)
-		: VectorFunction(d, 1), dim(d) {}
-	Function1(const Function1& f)
-		: VectorFunction(f.indim(), f.outdim()), dim(f.dim) {}
-	VectorFunction* clone() const
-		{return new Function1(*this);}
-	virtual void eval(const Vector& point, const ParameterSignal& sig, Vector& out);
-};
-
-void Function1::eval(const Vector& point, const ParameterSignal& sig, Vector& out)
-{
-	if (point.length() != dim || out.length() != 1) {
-		printf("Wrong length of vectors in Function1::eval\n");
-		exit(1);
-	}
-	double r = 1;
-	for (int i = 0; i < dim; i++)
-		r *= point[i];
-	r = pow(r, 1.0/dim);
-	r *= pow(1.0 + 1.0/dim, (double)dim);
-	out[0] = r;
-}
-
-// evaluates Function1 but with transformation x_i=0.5(y_i+1)
-// this makes the new function integrate over <-1,1>^d to 1.0
-class Function1Trans : public Function1 {
-public:
-	Function1Trans(int d)
-		: Function1(d) {}
-	Function1Trans(const Function1Trans& func)
-		: Function1(func) {}
-	VectorFunction* clone() const
-		{return new Function1Trans(*this);}
-	virtual void eval(const Vector& point, const ParameterSignal& sig, Vector& out);
-};
-
-void Function1Trans::eval(const Vector& point, const ParameterSignal& sig, Vector& out)
-{
-	Vector p(point.length());
-	for (int i = 0; i < p.length(); i++)
-		p[i] = 0.5*(point[i]+1);
-	Function1::eval(p, sig, out);
-	out.mult(pow(0.5,indim()));
-}
-
-
-// WallTimer class. Constructor saves the wall time, destructor
-// cancels the current time from the saved, and prints the message
-// with time information
-class WallTimer {
-	char mes[100];
-	struct timeval start;
-	bool new_line;
-public:
-	WallTimer(const char* m, bool nl = true)
-		{strcpy(mes, m);new_line = nl; gettimeofday(&start, NULL);}
-	~WallTimer()
-		{
-			struct timeval end;
-			gettimeofday(&end, NULL);
-			printf("%s%8.4g", mes,
-				   end.tv_sec-start.tv_sec + (end.tv_usec-start.tv_usec)*1.0e-6);
-			if (new_line)
-				printf("\n");
-		}
-};
-
-/****************************************************/
-/*     declaration of TestRunnable class            */
-/****************************************************/
-class TestRunnable {
-	char name[100];
-public:
-	int dim; // dimension of the solved problem
-	int nvar; // number of variable of the solved problem
-	TestRunnable(const char* n, int d, int nv)
-		: dim(d), nvar(nv)
-		{strncpy(name, n, 100);}
-	bool test() const;
-	virtual bool run() const =0;
-	const char* getName() const
-		{return name;}
-protected:
-	static bool smolyak_normal_moments(const GeneralMatrix& m, int imom, int level);
-	static bool product_normal_moments(const GeneralMatrix& m, int imom, int level);
-	static bool qmc_normal_moments(const GeneralMatrix& m, int imom, int level);
-	static bool smolyak_product_cube(const VectorFunction& func, const Vector& res,
-									 double tol, int level);
-	static bool qmc_cube(const VectorFunction& func, double res, double tol, int level);
-};
-
-bool TestRunnable::test() const
-{
-	printf("Running test <%s>\n",name);
-	bool passed;
-	{
-		WallTimer tim("Wall clock time ", false);
-		passed = run();
-	}
-	if (passed) {
-		printf("............................ passed\n\n");
-		return passed;
-	} else {
-		printf("............................ FAILED\n\n");
-		return passed;
-	}
-}
-
-
-/****************************************************/
-/*     definition of TestRunnable static methods    */
-/****************************************************/
-bool TestRunnable::smolyak_normal_moments(const GeneralMatrix& m, int imom, int level)
-{
-	// first make m*m' and then Cholesky factor
-	GeneralMatrix mtr(m, "transpose");
-	GeneralMatrix msq(m, mtr);
-
-	// make vector function
-	int dim = m.numRows();
-	TensorPower tp(dim, imom);
-	GaussConverterFunction func(tp, msq);
-
-	// smolyak quadrature
-	Vector smol_out(UFSTensor::calcMaxOffset(dim, imom));
-	{
-		WallTimer tim("\tSmolyak quadrature time:         ");
-		GaussHermite gs;
-		SmolyakQuadrature quad(dim, level, gs);
-		quad.integrate(func, level, num_threads, smol_out);
-		printf("\tNumber of Smolyak evaluations:    %d\n", quad.numEvals(level));
-	}
-
-	// check against theoretical moments
-	UNormalMoments moments(imom, msq);
-	smol_out.add(-1.0, (moments.get(Symmetry(imom)))->getData());
-	printf("\tError:                         %16.12g\n", smol_out.getMax());
-	return smol_out.getMax() < 1.e-7;
-}
-
-bool TestRunnable::product_normal_moments(const GeneralMatrix& m, int imom, int level)
-{
-	// first make m*m' and then Cholesky factor
-	GeneralMatrix mtr(m, "transpose");
-	GeneralMatrix msq(m, mtr);
-
-	// make vector function
-	int dim = m.numRows();
-	TensorPower tp(dim, imom);
-	GaussConverterFunction func(tp, msq);
-
-	// product quadrature
-	Vector prod_out(UFSTensor::calcMaxOffset(dim, imom));
-	{
-		WallTimer tim("\tProduct quadrature time:         ");
-		GaussHermite gs;
-		ProductQuadrature quad(dim, gs);
-		quad.integrate(func, level, num_threads, prod_out);
-		printf("\tNumber of product evaluations:    %d\n", quad.numEvals(level));
-	}
-
-	// check against theoretical moments
-	UNormalMoments moments(imom, msq);
-	prod_out.add(-1.0, (moments.get(Symmetry(imom)))->getData());
-	printf("\tError:                         %16.12g\n", prod_out.getMax());
-	return prod_out.getMax() < 1.e-7;
-}
-
-bool TestRunnable::qmc_normal_moments(const GeneralMatrix& m, int imom, int level)
-{
-	// first make m*m' and then Cholesky factor
-	GeneralMatrix mtr(m, "transpose");
-	GeneralMatrix msq(m, mtr);
-	GeneralMatrix mchol(msq);
-	int rows = mchol.numRows();
-	for (int i = 0; i < rows; i++)
-		for (int j = i+1; j < rows; j++)
-			mchol.get(i,j) = 0.0;
-	int info;
-	dpotrf("L", &rows, mchol.base(), &rows, &info);
-
-	// make vector function
-	MomentFunction func(mchol, imom);
-
-	// permutation schemes
-	WarnockPerScheme wps;
-	ReversePerScheme rps;
-	IdentityPerScheme ips;
-	PermutationScheme* scheme[] = {&wps, &rps, &ips};
-	const char* labs[] = {"Warnock", "Reverse", "Identity"};
-
-	// theoretical result
-	int dim = mchol.numRows();
-	UNormalMoments moments(imom, msq);
-	Vector res((const Vector&)((moments.get(Symmetry(imom)))->getData()));
-
-	// quasi monte carlo normal quadrature
-	double max_error = 0.0;
-	Vector qmc_out(UFSTensor::calcMaxOffset(dim, imom));
-	for (int i = 0; i < 3; i++) {
-		{
-			char mes[100];
-			sprintf(mes, "\tQMC normal quadrature time %8s:         ", labs[i]);
-			WallTimer tim(mes);
-			QMCarloNormalQuadrature quad(dim, level, *(scheme[i]));
-			quad.integrate(func, level, num_threads, qmc_out);
-		}
-		qmc_out.add(-1.0, res);
-		printf("\tError %8s:                         %16.12g\n", labs[i], qmc_out.getMax());
-		if (qmc_out.getMax() > max_error) {
-			max_error = qmc_out.getMax();
-		}
-	}
-
-	return max_error < 1.e-7;
-}
-
-
-bool TestRunnable::smolyak_product_cube(const VectorFunction& func, const Vector& res,
-										double tol, int level)
-{
-	if (res.length() != func.outdim()) {
-		fprintf(stderr, "Incompatible dimensions of check value and function.\n");
-		exit(1);
-	}
-
-	GaussLegendre glq;
-	Vector out(func.outdim());
-	double smol_error;
-	double prod_error;
-	{
-		WallTimer tim("\tSmolyak quadrature time:         ");
-		SmolyakQuadrature quad(func.indim(), level, glq);
-		quad.integrate(func, level, num_threads, out);
-		out.add(-1.0, res);
-		smol_error = out.getMax();
-		printf("\tNumber of Smolyak evaluations:    %d\n", quad.numEvals(level));
-		printf("\tError:                            %16.12g\n", smol_error);
-	}
-	{
-		WallTimer tim("\tProduct quadrature time:         ");
-		ProductQuadrature quad(func.indim(), glq);
-		quad.integrate(func, level, num_threads, out);
-		out.add(-1.0, res);
-		prod_error = out.getMax();
-		printf("\tNumber of product evaluations:    %d\n", quad.numEvals(level));
-		printf("\tError:                            %16.12g\n", prod_error);
-	}
-
-	return smol_error < tol && prod_error < tol;
-}
-
-bool TestRunnable::qmc_cube(const VectorFunction& func, double res, double tol, int level)
-{
-	Vector r(1);
-	double error1;
-	{
-		WallTimer tim("\tQuasi-Monte Carlo (Warnock scrambling) time:  ");
-		WarnockPerScheme wps;
-		QMCarloCubeQuadrature qmc(func.indim(), level, wps);
-//		qmc.savePoints("warnock.txt", level);
-		qmc.integrate(func, level, num_threads, r);
-		error1 = std::max(res - r[0], r[0] - res);
-		printf("\tQuasi-Monte Carlo (Warnock scrambling) error: %16.12g\n",
-			   error1);
-	}
-	double error2;
-	{
-		WallTimer tim("\tQuasi-Monte Carlo (reverse scrambling) time:  ");
-		ReversePerScheme rps;
-		QMCarloCubeQuadrature qmc(func.indim(), level, rps);
-//		qmc.savePoints("reverse.txt", level);
-		qmc.integrate(func, level, num_threads, r);
-		error2 = std::max(res - r[0], r[0] - res);
-		printf("\tQuasi-Monte Carlo (reverse scrambling) error: %16.12g\n",
-			   error2);
-	}
-	double error3;
-	{
-		WallTimer tim("\tQuasi-Monte Carlo (no scrambling) time:       ");
-		IdentityPerScheme ips;
-		QMCarloCubeQuadrature qmc(func.indim(), level, ips);
-//		qmc.savePoints("identity.txt", level);
-		qmc.integrate(func, level, num_threads, r);
-		error3 = std::max(res - r[0], r[0] - res);
-		printf("\tQuasi-Monte Carlo (no scrambling) error:      %16.12g\n",
-			   error3);
-	}
-
-
-	return error1 < tol && error2 < tol && error3 < tol;
-}
-
-/****************************************************/
-/*     definition of TestRunnable subclasses        */
-/****************************************************/
-class SmolyakNormalMom1 : public TestRunnable {
-public:
-	SmolyakNormalMom1()
-		: TestRunnable("Smolyak normal moments (dim=2, level=4, order=4)", 4, 2) {}
-
-	bool run() const
-		{
-			GeneralMatrix m(2,2);
-			m.zeros(); m.get(0,0)=1; m.get(1,1)=1;
-			return smolyak_normal_moments(m, 4, 4);
-		}
-};
-
-class SmolyakNormalMom2 : public TestRunnable {
-public:
-	SmolyakNormalMom2()
-		: TestRunnable("Smolyak normal moments (dim=3, level=8, order=8)", 8, 3) {}
-
-	bool run() const
-		{
-			GeneralMatrix m(3,3);
-			m.zeros();
-			m.get(0,0)=1; m.get(0,2)=0.5; m.get(1,1)=1;
-			m.get(1,0)=0.5;m.get(2,2)=2;m.get(2,1)=4;
-			return smolyak_normal_moments(m, 8, 8);
-		}
-};
-
-class ProductNormalMom1 : public TestRunnable {
-public:
-	ProductNormalMom1()
-		: TestRunnable("Product normal moments (dim=2, level=4, order=4)", 4, 2) {}
-
-	bool run() const
-		{
-			GeneralMatrix m(2,2);
-			m.zeros(); m.get(0,0)=1; m.get(1,1)=1;
-			return product_normal_moments(m, 4, 4);
-		}
-};
-
-class ProductNormalMom2 : public TestRunnable {
-public:
-	ProductNormalMom2()
-		: TestRunnable("Product normal moments (dim=3, level=8, order=8)", 8, 3) {}
-
-	bool run() const
-		{
-			GeneralMatrix m(3,3);
-			m.zeros();
-			m.get(0,0)=1; m.get(0,2)=0.5; m.get(1,1)=1;
-			m.get(1,0)=0.5;m.get(2,2)=2;m.get(2,1)=4;
-			return product_normal_moments(m, 8, 8);
-		}
-};
-
-class QMCNormalMom1 : public TestRunnable {
-public:
-	QMCNormalMom1()
-		: TestRunnable("QMC normal moments (dim=2, level=1000, order=4)", 4, 2) {}
-
-	bool run() const
-		{
-			GeneralMatrix m(2,2);
-			m.zeros(); m.get(0,0)=1; m.get(1,1)=1;
-			return qmc_normal_moments(m, 4, 1000);
-		}
-};
-
-class QMCNormalMom2 : public TestRunnable {
-public:
-	QMCNormalMom2()
-		: TestRunnable("QMC normal moments (dim=3, level=10000, order=8)", 8, 3) {}
-
-	bool run() const
-		{
-			GeneralMatrix m(3,3);
-			m.zeros();
-			m.get(0,0)=1; m.get(0,2)=0.5; m.get(1,1)=1;
-			m.get(1,0)=0.5;m.get(2,2)=2;m.get(2,1)=4;
-			return qmc_normal_moments(m, 8, 10000);
-		}
-};
-
-
-
-// note that here we pass 1,1 to tls since smolyak has its own PascalTriangle
-class F1GaussLegendre : public TestRunnable {
-public:
-	F1GaussLegendre()
-		: TestRunnable("Function1 Gauss-Legendre (dim=6, level=13", 1, 1) {}
-
-	bool run() const
-		{
-			Function1Trans f1(6);
-			Vector res(1); res[0] = 1.0;
-			return smolyak_product_cube(f1, res, 1e-2, 13);
-		}
-};
-
-
-class F1QuasiMCarlo : public TestRunnable {
-public:
-	F1QuasiMCarlo()
-		: TestRunnable("Function1 Quasi-Monte Carlo (dim=6, level=1000000)", 1, 1) {}
-
-	bool run() const
-		{
-			Function1 f1(6);
-			return qmc_cube(f1, 1.0, 1.e-4, 1000000);
-		}
-};
-
-int main()
-{
-	TestRunnable* all_tests[50];
-	// fill in vector of all tests
-	int num_tests = 0;
-	all_tests[num_tests++] = new SmolyakNormalMom1();
-	all_tests[num_tests++] = new SmolyakNormalMom2();
-	all_tests[num_tests++] = new ProductNormalMom1();
-	all_tests[num_tests++] = new ProductNormalMom2();
-	all_tests[num_tests++] = new QMCNormalMom1();
-	all_tests[num_tests++] = new QMCNormalMom2();
-/*
-	all_tests[num_tests++] = new F1GaussLegendre();
-	all_tests[num_tests++] = new F1QuasiMCarlo();
-*/
-	// find maximum dimension and maximum nvar
-	int dmax=0;
-	int nvmax = 0;
-	for (int i = 0; i < num_tests; i++) {
-		if (dmax < all_tests[i]->dim)
-			dmax = all_tests[i]->dim;
-		if (nvmax < all_tests[i]->nvar)
-			nvmax = all_tests[i]->nvar;
-	}
-	tls.init(dmax, nvmax); // initialize library
-	THREAD_GROUP::max_parallel_threads = num_threads;
-
-	// launch the tests
-	int success = 0;
-	for (int i = 0; i < num_tests; i++) {
-		try {
-			if (all_tests[i]->test())
-				success++;
-		} catch (const TLException& e) {
-			printf("Caugth TL exception in <%s>:\n", all_tests[i]->getName());
-			e.print();
-		} catch (SylvException& e) {
-			printf("Caught Sylv exception in <%s>:\n", all_tests[i]->getName());
-			e.printMessage();
-		}
-	}
-
-	printf("There were %d tests that failed out of %d tests run.\n",
-		   num_tests - success, num_tests);
-
-	// destroy
-	for (int i = 0; i < num_tests; i++) {
-		delete all_tests[i];
-	}
-
-	return 0;
-}

dynare++/kord/Makefile.am

-CWEBSRC = \
-	faa_di_bruno.cweb \
-	korder_stoch.cweb \
-	journal.cweb \
-	decision_rule.cweb \
-	dynamic_model.cweb \
-	random.cweb \
-	first_order.cweb \
-	normal_conjugate.cweb \
-	approximation.cweb \
-	global_check.cweb \
-	korder.cweb \
-	kord_exception.hweb \
-	random.hweb \
-	journal.hweb \
-	approximation.hweb \
-	korder_stoch.hweb \
-	dynamic_model.hweb \
-	decision_rule.hweb \
-	korder.hweb \
-	normal_conjugate.hweb \
-	first_order.hweb \
-	mersenne_twister.hweb \
-	global_check.hweb \
-	faa_di_bruno.hweb
-
-GENERATED_FILES = \
-	faa_di_bruno.cpp \
-	korder_stoch.cpp \
-	journal.cpp \
-	decision_rule.cpp \
-	dynamic_model.cpp \
-	random.cpp \
-	first_order.cpp \
-	normal_conjugate.cpp \
-	approximation.cpp \
-	global_check.cpp \
-	korder.cpp \
-	kord_exception.h \
-	random.h \
-	journal.h \
-	approximation.h \
-	korder_stoch.h \
-	dynamic_model.h \
-	decision_rule.h \
-	korder.h \
-	normal_conjugate.h \
-	first_order.h \
-	mersenne_twister.h \
-	global_check.h \
-	faa_di_bruno.h
-
-noinst_LIBRARIES = libkord.a
-
-libkord_a_SOURCES = $(CWEBSRC) $(GENERATED_FILES)
-libkord_a_CPPFLAGS = -I../sylv/cc -I../tl/cc -I../integ/cc -I$(top_srcdir)/mex/sources $(CPPFLAGS_MATIO)
-libkord_a_CXXFLAGS = $(PTHREAD_CFLAGS)
-
-BUILT_SOURCES = $(GENERATED_FILES)
-
-EXTRA_DIST = main.web dummy.ch
-
-check_PROGRAMS = tests
-
-tests_SOURCES = tests.cpp
-tests_CPPFLAGS = -I../sylv/cc -I../tl/cc -I../integ/cc -I$(top_srcdir)/mex/sources
-tests_CXXFLAGS = $(PTHREAD_CFLAGS)
-tests_LDFLAGS = $(LDFLAGS_MATIO)
-tests_LDADD = libkord.a ../tl/cc/libtl.a ../sylv/cc/libsylv.a $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS) $(PTHREAD_LIBS) $(LIBADD_MATIO)
-
-check-local:
-	./tests
-
-%.cpp: %.cweb dummy.ch
-	$(CTANGLE) -bhp $< dummy.ch $@
-
-%.h: %.hweb dummy.ch
-	$(CTANGLE) -bhp $< dummy.ch $@
-
-if HAVE_CWEAVE
-if HAVE_PDFTEX
-if HAVE_EPLAIN
-pdf-local: kord.pdf
-
-kord.pdf: main.web $(CWEBSRC)
-	$(CWEAVE) -bhp main.web
-	$(PDFTEX) main
-	mv main.pdf kord.pdf
-endif
-endif
-endif
-
-CLEANFILES = kord.pdf main.idx main.log main.scn main.tex main.toc out.txt

dynare++/kord/approximation.cweb

-@q $Id: approximation.cweb 2344 2009-02-09 20:36:08Z michel $ @>
-@q Copyright 2005, Ondra Kamenik @>
-
-@ Start of {\tt approximation.cpp} file.
-
-@c
-#include "kord_exception.h"
-#include "approximation.h"
-#include "first_order.h"
-#include "korder_stoch.h"
-
-@<|ZAuxContainer| constructor code@>;
-@<|ZAuxContainer::getType| code@>;
-@<|Approximation| constructor code@>;
-@<|Approximation| destructor code@>;
-@<|Approximation::getFoldDecisionRule| code@>;
-@<|Approximation::getUnfoldDecisionRule| code@>;
-@<|Approximation::approxAtSteady| code@>;
-@<|Approximation::walkStochSteady| code@>;
-@<|Approximation::saveRuleDerivs| code@>;
-@<|Approximation::calcStochShift| code@>;
-@<|Approximation::check| code@>;
-@<|Approximation::calcYCov| code@>;
-
-@ 
-@<|ZAuxContainer| constructor code@>=
-ZAuxContainer::ZAuxContainer(const _Ctype* gss, int ngss, int ng, int ny, int nu)
-	: StackContainer<FGSTensor>(4,1)
-{
-	stack_sizes[0] = ngss; stack_sizes[1] = ng;
-	stack_sizes[2] = ny; stack_sizes[3] = nu;
-	conts[0] = gss;
-	calculateOffsets();
-}
-
-
-@ The |getType| method corresponds to
-$f(g^{**}(y^*,u',\sigma),0,0,0)$. For the first argument we return
-|matrix|, for other three we return |zero|.
-
-@<|ZAuxContainer::getType| code@>=
-ZAuxContainer::itype ZAuxContainer::getType(int i, const Symmetry& s) const
-{
-	if (i == 0)
-		if (s[2] > 0)
-			return zero;
-		else
-			return matrix;
-	return zero;
-}
-
-
-@ 
-@<|Approximation| constructor code@>=
-Approximation::Approximation(DynamicModel& m, Journal& j, int ns, bool dr_centr, double qz_crit)
-	: model(m), journal(j), rule_ders(NULL), rule_ders_ss(NULL), fdr(NULL), udr(NULL),
-	  ypart(model.nstat(), model.npred(), model.nboth(), model.nforw()),
-	  mom(UNormalMoments(model.order(), model.getVcov())), nvs(4), steps(ns),
-	  dr_centralize(dr_centr), qz_criterium(qz_crit), ss(ypart.ny(), steps+1)
-{
-	nvs[0] = ypart.nys(); nvs[1] = model.nexog();
-	nvs[2] = model.nexog(); nvs[3] = 1;
-
-	ss.nans();
-}
-
-@ 
-@<|Approximation| destructor code@>=
-Approximation::~Approximation()
-{
-	if (rule_ders_ss) delete rule_ders_ss;
-	if (rule_ders) delete rule_ders;
-	if (fdr) delete fdr;
-	if (udr) delete udr;
-}
-
-@ This just returns |fdr| with a check that it is created.
-@<|Approximation::getFoldDecisionRule| code@>=
-const FoldDecisionRule& Approximation::getFoldDecisionRule() const
-{
-	KORD_RAISE_IF(fdr == NULL,
-				  "Folded decision rule has not been created in Approximation::getFoldDecisionRule");
-	return *fdr;
-}
-
-
-@ This just returns |udr| with a check that it is created.
-@<|Approximation::getUnfoldDecisionRule| code@>=
-const UnfoldDecisionRule& Approximation::getUnfoldDecisionRule() const
-{
-	KORD_RAISE_IF(udr == NULL,
-				  "Unfolded decision rule has not been created in Approximation::getUnfoldDecisionRule");
-	return *udr;
-}
-
-
-@ This methods assumes that the deterministic steady state is
-|model.getSteady()|. It makes an approximation about it and stores the
-derivatives to |rule_ders| and |rule_ders_ss|. Also it runs a |check|
-for $\sigma=0$.
-
-@<|Approximation::approxAtSteady| code@>=
-void Approximation::approxAtSteady()
-{
-	model.calcDerivativesAtSteady();
-	FirstOrder fo(model.nstat(), model.npred(), model.nboth(), model.nforw(),
-				  model.nexog(), *(model.getModelDerivatives().get(Symmetry(1))),
-				  journal, qz_criterium);
-	KORD_RAISE_IF_X(! fo.isStable(),
-					"The model is not Blanchard-Kahn stable",
-					KORD_MD_NOT_STABLE);
-
-	if (model.order() >= 2) {
-		KOrder korder(model.nstat(), model.npred(), model.nboth(), model.nforw(),
-					  model.getModelDerivatives(), fo.getGy(), fo.getGu(),
-					  model.getVcov(), journal);
-		korder.switchToFolded();
-		for (int k = 2; k <= model.order(); k++)
-			korder.performStep<KOrder::fold>(k);
-		
-		saveRuleDerivs(korder.getFoldDers());
-	} else {
-		FirstOrderDerivs<KOrder::fold> fo_ders(fo);
-		saveRuleDerivs(fo_ders);
-	}
-	check(0.0);
-}
-
-@ This is the core routine of |Approximation| class.
-
-First we solve for the approximation about the deterministic steady
-state. Then we perform |steps| cycles toward the stochastic steady
-state. Each cycle moves the size of shocks by |dsigma=1.0/steps|. At
-the end of a cycle, we have |rule_ders| being the derivatives at
-stochastic steady state for $\sigma=sigma\_so\_far+dsigma$ and
-|model.getSteady()| being the steady state.
-
-If the number of |steps| is zero, the decision rule |dr| at the bottom
-is created from derivatives about deterministic steady state, with
-size of $\sigma=1$. Otherwise, the |dr| is created from the
-approximation about stochastic steady state with $\sigma=0$.
-
-Within each cycle, we first make a backup of the last steady (from
-initialization or from a previous cycle), then we calculate the fix
-point of the last rule with $\sigma=dsigma$. This becomes a new steady
-state at the $\sigma=sigma\_so\_far+dsigma$. We calculate expectations
-of $g^{**}(y,\sigma\eta_{t+1},\sigma$ expressed as a Taylor expansion
-around the new $\sigma$ and the new steady state. Then we solve for
-the decision rule with explicit $g^{**}$ at $t+1$ and save the rule.
-
-After we reached $\sigma=1$, the decision rule is formed.
-
-The biproduct of this method is the matrix |ss|, whose columns are
-steady states for subsequent $\sigma$s. The first column is the
-deterministic steady state, the last column is the stochastic steady
-state for a full size of shocks ($\sigma=1$). There are |steps+1|
-columns.
-
-@<|Approximation::walkStochSteady| code@>=
-void Approximation::walkStochSteady()
-{
-	@<initial approximation at deterministic steady@>;
-	double sigma_so_far = 0.0;
-	double dsigma = (steps == 0)? 0.0 : 1.0/steps;
-	for (int i = 1; i <= steps; i++) {
-		JournalRecordPair pa(journal);
-		pa << "Approximation about stochastic steady for sigma=" << sigma_so_far+dsigma << endrec;
-
-		Vector last_steady((const Vector&)model.getSteady());
-
-		@<calculate fix-point of the last rule for |dsigma|@>;
-		@<calculate |hh| as expectations of the last $g^{**}$@>;
-		@<form |KOrderStoch|, solve and save@>;
-
-		check(sigma_so_far+dsigma);
-		sigma_so_far += dsigma;
-	}
-
-	@<construct the resulting decision rules@>;
-}
-
-@ Here we solve for the deterministic steady state, calculate
-approximation at the deterministic steady and save the steady state
-to |ss|.
-
-@<initial approximation at deterministic steady@>=
-	model.solveDeterministicSteady();
-	approxAtSteady();
-	Vector steady0(ss, 0);
-	steady0 = (const Vector&)model.getSteady();
-
-@ We form the |DRFixPoint| object from the last rule with
-$\sigma=dsigma$. Then we save the steady state to |ss|. The new steady
-is also put to |model.getSteady()|.
-
-@<calculate fix-point of the last rule for |dsigma|@>=
-	DRFixPoint<KOrder::fold> fp(*rule_ders, ypart, model.getSteady(), dsigma);
-	bool converged = fp.calcFixPoint(DecisionRule::horner, model.getSteady());
-	JournalRecord rec(journal);
-	rec << "Fix point calcs: iter=" << fp.getNumIter() << ", newton_iter="
-		<< fp.getNewtonTotalIter() << ", last_newton_iter=" << fp.getNewtonLastIter() << ".";
-	if (converged)
-		rec << " Converged." << endrec;
-	else {
-		rec << " Not converged!!" << endrec;
-		KORD_RAISE_X("Fix point calculation not converged", KORD_FP_NOT_CONV);
-	}
-	Vector steadyi(ss, i);
-	steadyi = (const Vector&)model.getSteady();
-
-@ We form the steady state shift |dy|, which is the new steady state
-minus the old steady state. Then we create |StochForwardDerivs|
-object, which calculates the derivatives of $g^{**}$ expectations at
-new sigma and new steady.
-
-@<calculate |hh| as expectations of the last $g^{**}$@>=
-	Vector dy((const Vector&)model.getSteady());
-	dy.add(-1.0, last_steady);
-
-	StochForwardDerivs<KOrder::fold> hh(ypart, model.nexog(), *rule_ders_ss, mom, dy,
-										dsigma, sigma_so_far);
-	JournalRecord rec1(journal);
-	rec1 << "Calculation of g** expectations done" << endrec;
-
-
-@ We calculate derivatives of the model at the new steady, form
-|KOrderStoch| object and solve, and save the rule.
-
-@<form |KOrderStoch|, solve and save@>=
-	model.calcDerivativesAtSteady();
-	KOrderStoch korder_stoch(ypart, model.nexog(), model.getModelDerivatives(),
-							 hh, journal);
-	for (int d = 1; d <= model.order(); d++) {
-		korder_stoch.performStep<KOrder::fold>(d);
-	}
-	saveRuleDerivs(korder_stoch.getFoldDers());
-
-
-@ 
-@<construct the resulting decision rules@>=
-	if (fdr) {
-		delete fdr;
-		fdr = NULL;
-	}
-	if (udr) {
-		delete udr;
-		udr = NULL;
-	}
-
-	fdr = new FoldDecisionRule(*rule_ders, ypart, model.nexog(),
-							   model.getSteady(), 1.0-sigma_so_far);
-	if (steps == 0 && dr_centralize) {
-		@<centralize decision rule for zero steps@>;
-	}
-
-
-@ 
-@<centralize decision rule for zero steps@>=
-	DRFixPoint<KOrder::fold> fp(*rule_ders, ypart, model.getSteady(), 1.0);
-	bool converged = fp.calcFixPoint(DecisionRule::horner, model.getSteady());
-	JournalRecord rec(journal);
-	rec << "Fix point calcs: iter=" << fp.getNumIter() << ", newton_iter="
-		<< fp.getNewtonTotalIter() << ", last_newton_iter=" << fp.getNewtonLastIter() << ".";
-	if (converged)
-		rec << " Converged." << endrec;
-	else {
-		rec << " Not converged!!" << endrec;
-		KORD_RAISE_X("Fix point calculation not converged", KORD_FP_NOT_CONV);
-	}
-
-	{
-		JournalRecordPair recp(journal);
-		recp << "Centralizing about fix-point." << endrec;
-		FoldDecisionRule* dr_backup = fdr;
-		fdr = new FoldDecisionRule(*dr_backup, model.getSteady());
-		delete dr_backup;
-	}
-
-
-@ Here we simply make a new hardcopy of the given rule |rule_ders|,
-and make a new container of in-place subtensors of the derivatives
-corresponding to forward looking variables. The given container comes
-from a temporary object and will be destroyed.
- 
-@<|Approximation::saveRuleDerivs| code@>=
-void Approximation::saveRuleDerivs(const FGSContainer& g)
-{
-	if (rule_ders) {
-		delete rule_ders;
-		delete rule_ders_ss;
-	}
-	rule_ders = new FGSContainer(g);
-	rule_ders_ss = new FGSContainer(4);
-	for (FGSContainer::iterator run = (*rule_ders).begin(); run != (*rule_ders).end(); ++run) {
-		FGSTensor* ten = new FGSTensor(ypart.nstat+ypart.npred, ypart.nyss(), *((*run).second));
-		rule_ders_ss->insert(ten);
-	}
-}
-
-@ This method calculates a shift of the system equations due to
-integrating shocks at a given $\sigma$ and current steady state. More precisely, if
-$$F(y,u,u',\sigma)=f(g^{**}(g^*(y,u,\sigma),u',\sigma),g(y,u,\sigma),y,u)$$
-then the method returns a vector
-$$\sum_{d=1}{1\over d!}\sigma^d\left[F_{u'^d}\right]_{\alpha_1\ldots\alpha_d}
-\Sigma^{\alpha_1\ldots\alpha_d}$$
-
-For a calculation of $\left[F_{u'^d}\right]$ we use |@<|ZAuxContainer|
-class declaration@>|, so we create its object. In each cycle we
-calculate $\left[F_{u'^d}\right]$@q'@>, and then multiply with the shocks,
-and add the ${\sigma^d\over d!}$ multiple to the result.
-
-@<|Approximation::calcStochShift| code@>=
-void Approximation::calcStochShift(Vector& out, double at_sigma) const
-{
-	KORD_RAISE_IF(out.length() != ypart.ny(),
-				  "Wrong length of output vector for Approximation::calcStochShift");
-	out.zeros();
-
-	ZAuxContainer zaux(rule_ders_ss, ypart.nyss(), ypart.ny(),
-					   ypart.nys(), model.nexog());
-
-	int dfac = 1;
-	for (int d = 1; d <= rule_ders->getMaxDim(); d++, dfac*=d) {
-		if ( KOrder::is_even(d)) {
-			Symmetry sym(0,d,0,0);
-			@<calculate $F_{u'^d}$ via |ZAuxContainer|@>;@q'@>
-			@<multiply with shocks and add to result@>;
-		}
-	}
-}
-
-@ 
-@<calculate $F_{u'^d}$ via |ZAuxContainer|@>=
-	FGSTensor* ten = new FGSTensor(ypart.ny(), TensorDimens(sym, nvs));
-	ten->zeros();
-	for (int l = 1; l <= d; l++) {
-		const FSSparseTensor* f = model.getModelDerivatives().get(Symmetry(l));
-		zaux.multAndAdd(*f, *ten);
-	}
-
-@
-@<multiply with shocks and add to result@>=
-	FGSTensor* tmp = new FGSTensor(ypart.ny(), TensorDimens(Symmetry(0,0,0,0), nvs));
-	tmp->zeros();
-	ten->contractAndAdd(1, *tmp, *(mom.get(Symmetry(d))));
-
-	out.add(pow(at_sigma,d)/dfac, tmp->getData());
-	delete ten;
-	delete tmp;
-
-
-@ This method calculates and reports
-$$f(\bar y)+\sum_{d=1}{1\over d!}\sigma^d\left[F_{u'^d}\right]_{\alpha_1\ldots\alpha_d}
-\Sigma^{\alpha_1\ldots\alpha_d}$$
-at $\bar y$, zero shocks and $\sigma$. This number should be zero.
-
-We evaluate the error both at a given $\sigma$ and $\sigma=1.0$.
-
-@<|Approximation::check| code@>=
-void Approximation::check(double at_sigma) const
-{
-	Vector stoch_shift(ypart.ny());
-	Vector system_resid(ypart.ny());
-	Vector xx(model.nexog());
-	xx.zeros();
-	model.evaluateSystem(system_resid, model.getSteady(), xx);
-	calcStochShift(stoch_shift, at_sigma);
-	stoch_shift.add(1.0, system_resid);
-	JournalRecord rec1(journal);
-	rec1 << "Error of current approximation for shocks at sigma " << at_sigma
-		 << " is " << stoch_shift.getMax() << endrec;
-	calcStochShift(stoch_shift, 1.0);
-	stoch_shift.add(1.0, system_resid);
-	JournalRecord rec2(journal);
-	rec2 << "Error of current approximation for full shocks is " << stoch_shift.getMax() << endrec;
-}
-
-@ The method returns unconditional variance of endogenous variables
-based on the first order. The first order approximation looks like
-$$\hat y_t=g_{y^*}\hat y^*_{t-1}+g_uu_t$$
-where $\hat y$ denotes a deviation from the steady state. It can be written as
-$$\hat y_t=\left[0\, g_{y^*}\, 0\right]\hat y_{t-1}+g_uu_t$$
-which yields unconditional covariance $V$ for which
-$$V=GVG^T + g_u\Sigma g_u^T,$$
-where $G=[0\, g_{y^*}\, 0]$ and $\Sigma$ is the covariance of the shocks. 
-
-For solving this Lyapunov equation we use the Sylvester module, which
-solves equation of the type
-$$AX+BX(C\otimes\cdots\otimes C)=D$$
-So we invoke the Sylvester solver for the first dimension with $A=I$,
-$B=-G$, $C=G^T$ and $D=g_u\Sigma g_u^T$.
-
-
-@<|Approximation::calcYCov| code@>=
-TwoDMatrix* Approximation::calcYCov() const
-{
-	const TwoDMatrix& gy = *(rule_ders->get(Symmetry(1,0,0,0)));
-	const TwoDMatrix& gu = *(rule_ders->get(Symmetry(0,1,0,0)));
-	TwoDMatrix G(model.numeq(), model.numeq());
-	G.zeros();
-	G.place(gy, 0, model.nstat());
-	TwoDMatrix B((const TwoDMatrix&)G);
-	B.mult(-1.0);
-	TwoDMatrix C(G, "transpose");
-	TwoDMatrix A(model.numeq(), model.numeq());
-	A.zeros();
-	for (int i = 0; i < model.numeq(); i++)
-		A.get( i,i)	= 1.0;
-
-	TwoDMatrix guSigma(gu, model.getVcov());
-	TwoDMatrix guTrans(gu, "transpose");
-	TwoDMatrix* X = new TwoDMatrix(guSigma, guTrans);
-	
-	GeneralSylvester gs(1, model.numeq(), model.numeq(), 0,
-						A.base(), B.base(), C.base(), X->base());
-	gs.solve();
-
-	return X;
-}
-
-@ End of {\tt approximation.cpp} file.

dynare++/kord/approximation.hweb

-@q $Id: approximation.hweb 2352 2009-09-03 19:18:15Z michel $ @>
-@q Copyright 2005, Ondra Kamenik @>
-
-@*2 Approximating model solution. Start of {\tt approximation.h} file.
-
-The class |Approximation| in this file is a main interface to the
-algorithms calculating approximations to the decision rule about
-deterministic and stochastic steady states.
-
-The approximation about a deterministic steady state is solved by
-classes |@<|FirstOrder| class declaration@>| and |@<|KOrder| class
-declaration@>|. The approximation about the stochastic steady state is
-solved by class |@<|KOrderStoch| class declaration@>| together with a
-method of |Approximation| class |@<|Approximation::walkStochSteady|
-code@>|.
-
-The approximation about the stochastic steady state is done with
-explicit expression of forward derivatives of $g^{**}$. More formally,
-we have to solve the decision rule $g$ from the implicit system:
-$$E_t(f(g^{**}(g^*(y^*,u_t,\sigma),u_{t+1},\sigma),g(y^*,u_t,\sigma),y_t,u_t))=0$$
-The term within the expectations can be Taylor expanded, and the
-expectation can be driven into the formula. However, when doing this
-at $\sigma\not=0$, the term $g^{**}$ at $\sigma\not=0$ is dependent on
-$u_{t+1}$ and thus the integral of its approximation includes all
-derivatives wrt. $u$ of $g^{**}$. Note that for $\sigma=0$, the
-derivatives of $g^{**}$ in this context are constant. This is the main
-difference between the approximation at deterministic steady
-($\sigma=0$), and stochastic steady ($\sigma\not=0$). This means that
-$k$-order derivative of the above equation at $\sigma\not=0$ depends of
-all derivatives of $g^**$ (including those with order greater than
-$k$).
-
-The explicit expression of the forward $g^{**}$ means that the
-derivatives of $g$ are not solved simultaneously, but that the forward
-derivatives of $g^{**}$ are calculated as an extrapolation based on
-the approximation at lower $\sigma$. This is exactly what does the
-|@<|Approximation::walkStochSteady| code@>|. It starts at the
-deterministic steady state, and in a few steps it adds to $\sigma$
-explicitly expressing forward $g^{**}$ from a previous step.
-
-Further details on the both solution methods are given in (todo: put
-references here when they exist).
-
-Very important note: all classes here used for calculation of decision
-rule approximation are folded. For the time being, it seems that faa
-Di Bruno formula is quicker for folded tensors, and that is why we
-stick to folded tensors here. However, when the calcs are done, we
-calculate also its unfolded versions, to be available for simulations
-and so on.
-
-@s ZAuxContainer int
-@s Approximation int
-@c
-#ifndef APPROXIMATION_H
-#define APPROXIMATION_H
-
-#include "dynamic_model.h"
-#include "decision_rule.h"
-#include "korder.h"
-#include "journal.h"
-
-@<|ZAuxContainer| class declaration@>;
-@<|Approximation| class declaration@>;
-
-#endif
-
-
-@ This class is used to calculate derivatives by faa Di Bruno of the
-$$f(g^{**}(g^*(y^*,u,\sigma),u',\sigma),g(y^*,u,\sigma),y^*,u)$$ with
-respect $u'$. In order to keep it as simple as possible, the class
-represents an equivalent (with respect to $u'$) container for
-$f(g^{**}(y^*,u',\sigma),0,0,0)$. The class is used only for
-evaluation of approximation error in |Approximation| class, which is
-calculated in |Approximation::calcStochShift| method.
-
-Since it is a folded version, we inherit from
-|StackContainer<FGSTensor>| and |FoldedStackContainer|. To construct
-it, we need only the $g^{**}$ container and size of stacks.
-
-@<|ZAuxContainer| class declaration@>=
-class ZAuxContainer : public StackContainer<FGSTensor>, public FoldedStackContainer {
-public:@;
-	typedef StackContainer<FGSTensor>::_Ctype _Ctype;
-	typedef StackContainer<FGSTensor>::itype itype;
-	ZAuxContainer(const _Ctype* gss, int ngss, int ng, int ny, int nu);
-	itype getType(int i, const Symmetry& s) const;
-};
-
-
-
-@ This class provides an interface to approximation algorithms. The
-core method is |walkStochSteady| which calculates the approximation
-about stochastic steady state in a given number of steps. The number
-is given as a parameter |ns| of the constructor. If the number is
-equal to zero, the resulted approximation is about the deterministic
-steady state.
-
-An object is constructed from the |DynamicModel|, and the number of
-steps |ns|. Also, we pass a reference to journal. That's all. The
-result of the core method |walkStochSteady| is a decision rule |dr|
-and a matrix |ss| whose columns are steady states for increasing
-$\sigma$ during the walk. Both can be retrived by public methods. The
-first column of the matrix is the deterministic steady state, the last
-is the stochastic steady state for the full size shocks.
-
-The method |walkStochSteady| calls the following methods:
-|approxAtSteady| calculates an initial approximation about the
-deterministic steady, |saveRuleDerivs| saves derivatives of a rule for
-the following step in |rule_ders| and |rule_ders_ss| (see
-|@<|Approximation::saveRuleDerivs| code@>| for their description),
-|check| reports an error of the current approximation and
-|calcStochShift| (called from |check|) calculates a shift of the
-system equations due to uncertainity.
-
-dr\_centralize is a new option. dynare++ was automatically expressing 
-results around the fixed point instead of the deterministic steady 
-state. dr\_centralize controls this behavior. 
-
-
-@<|Approximation| class declaration@>=
-class Approximation {
-	DynamicModel& model;
-	Journal& journal;
-	FGSContainer* rule_ders;
-	FGSContainer* rule_ders_ss;
-	FoldDecisionRule* fdr;
-	UnfoldDecisionRule* udr;
-	const PartitionY ypart;
-	const FNormalMoments mom;
-	IntSequence nvs;
-	int steps;
-	bool dr_centralize;
-	double qz_criterium;
-	TwoDMatrix ss;
-public:@;
-	Approximation(DynamicModel& m, Journal& j, int ns, bool dr_centr, double qz_crit);
-	virtual ~Approximation();
-
-	const FoldDecisionRule& getFoldDecisionRule() const;
-	const UnfoldDecisionRule& getUnfoldDecisionRule() const;
-	const TwoDMatrix& getSS() const
-		{@+ return ss;@+}
-	const DynamicModel& getModel() const
-		{@+ return model;@+}
-
-	void walkStochSteady();
-	TwoDMatrix* calcYCov() const;
-	const FGSContainer* get_rule_ders() const
-	      	{@+ return rule_ders;@+}	   
-	const FGSContainer* get_rule_ders_ss() const
-	      	{@+ return rule_ders;@+}	   
-protected:@;
-	void approxAtSteady();
-	void calcStochShift(Vector& out, double at_sigma) const;
-	void saveRuleDerivs(const FGSContainer& g);
-	void check(double at_sigma) const;
-};
-
-
-@ End of {\tt approximation.h} file.
-

dynare++/kord/decision_rule.cweb

-@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::writeMat| code1@>;
-@<|SimResults::writeMat| code2@>;
-@<|SimResultsStats::simulate| code@>;
-@<|SimResultsStats::writeMat| code@>;
-@<|SimResultsStats::calcMean| code@>;
-@<|SimResultsStats::calcVcov| code@>;
-@<|SimResultsDynamicStats::simulate| code@>;
-@<|SimResultsDynamicStats::writeMat| code@>;
-@<|SimResultsDynamicStats::calcMean| code@>;
-@<|SimResultsDynamicStats::calcVariance| code@>;
-@<|SimResultsIRF::simulate| code1@>;
-@<|SimResultsIRF::simulate| code2@>;
-@<|SimResultsIRF::calcMeans| code@>;
-@<|SimResultsIRF::calcVariances| code@>;
-@<|SimResultsIRF::writeMat| code@>;
-@<|RTSimResultsStats::simulate| code1@>;
-@<|RTSimResultsStats::simulate| code2@>;
-@<|RTSimResultsStats::writeMat| code@>;
-@<|IRFResults| constructor@>;
-@<|IRFResults| destructor@>;
-@<|IRFResults::writeMat| 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 burning " << num_burn << " initial 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+num_burn, start, rsrs.back());
-		gr.insert(worker);
-	}
-	gr.run();
-}
-
-@ 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+num_burn,
-				  "Incompatible number of cols for SimResults::addDataSets");
-	bool ret = false;
-	if (d->isFinite()) {
-		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;
-}
-
-@ 
-@<|SimResults::writeMat| code1@>=
-void SimResults::writeMat(const char* base, const char* lname) const
-{
-	char matfile_name[100];
-	sprintf(matfile_name, "%s.mat", base);
-  mat_t* matfd = Mat_Create(matfile_name, NULL);
-  if (matfd != NULL) {
-		writeMat(matfd, lname);
-    Mat_Close(matfd);
-	}
-}
-
-@ 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::writeMat| code2@>=
-void SimResults::writeMat(mat_t* 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.writeMat(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::writeMat| code@>=
-void SimResultsStats::writeMat(mat_t* fd, const char* lname) const
-{
-	char tmp[100];
-	sprintf(tmp, "%s_mean", lname);
-	ConstTwoDMatrix m(num_y, 1, mean.base());
-	m.writeMat(fd, tmp);
-	sprintf(tmp, "%s_vcov", lname);
-	ConstTwoDMatrix(vcov).writeMat(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::writeMat| code@>=
-void SimResultsDynamicStats::writeMat(mat_t* fd, const char* lname) const
-{
-	char tmp[100];
-	sprintf(tmp, "%s_cond_mean", lname);
-	ConstTwoDMatrix(mean).writeMat(fd, tmp);
-	sprintf(tmp, "%s_cond_variance", lname);
-	ConstTwoDMatrix(variance).writeMat(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, Journal& journal)
-{
-	JournalRecordPair paa(journal);
-	paa << "Performing " << control.getNumSets() << " IRF simulations for "
-		<< num_per << " periods; shock=" << ishock << ", impulse=" << imp << endrec;
-	simulate(dr);
-	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)
-{
-	THREAD_GROUP gr;
-	for (int idata = 0; idata < control.getNumSets(); idata++) {
-		THREAD* worker = new
-			SimulationIRFWorker(*this, dr, DecisionRule::horner,
-								num_per, 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::writeMat| code@>=
-void SimResultsIRF::writeMat(mat_t* fd, const char* lname) const
-{
-	char tmp[100];
-	sprintf(tmp, "%s_mean", lname);
-	means.writeMat(fd, tmp);
-	sprintf(tmp, "%s_var", lname);
-	variances.writeMat(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::writeMat| code@>=
-void RTSimResultsStats::writeMat(mat_t* fd, const char* lname)
-{
-	char tmp[100];
-	sprintf(tmp, "%s_rt_mean", lname);
-	ConstTwoDMatrix m(nc.getDim(), 1, mean.base());
-	m.writeMat(fd, tmp);
-	sprintf(tmp, "%s_rt_vcov", lname);
-	ConstTwoDMatrix(vcov).writeMat(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, journal);
-		irf_res[2*ii+1]->simulate(dr, journal);
-	}
-}
-
-@ 
-@<|IRFResults| destructor@>=
-IRFResults::~IRFResults()
-{
-	for (unsigned int i = 0; i < irf_res.size(); i++)
-		delete irf_res[i];
-}
-
-@ 
-@<|IRFResults::writeMat| code@>=
-void IRFResults::writeMat(mat_t* 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]->writeMat(fd, tmp);
-		sprintf(tmp, "%s_irfm_%s", prefix, shockname);
-		irf_res[2*i+1]->writeMat(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);
-	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));
-	{
-		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);
-    if (ip >= res.num_burn)
-		nc.update(y);
-
-@
-@<simulate other real-time periods@>=
-while (y.isFinite() && ip < res.num_burn + res.num_per) {
-	ip++;
-	dy = ypred;
-	sr.get(ip, u);
-	dr.eval(em, y, dyu);
-    if (ip >= res.num_burn)
-	    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++)
-		if (! isfinite(out[j]))
-			out[j] = r[j];
-}
-
-
-@ End of {\tt decision\_rule.cpp} file.

dynare++/kord/decision_rule.hweb

-@q $Id: decision_rule.hweb 2336 2009-01-14 10:37:02Z kamenik $ @>
-@q Copyright 2004, Ondra Kamenik @>
-
-@*2 Decision rule and simulation. Start of {\tt decision\_rule.h} file.
-
-The main purpose of this file is a decision rule representation which
-can run a simulation. So we define an interface for classes providing
-realizations of random shocks, and define the class
-|DecisionRule|. The latter basically takes tensor container of
-derivatives of policy rules, and adds them up with respect to
-$\sigma$. The class allows to specify the $\sigma$ different from $1$.
-
-In addition, we provide classes for running simulations and storing
-the results, calculating some statistics and generating IRF. The class
-|DRFixPoint| allows for calculation of the fix point of a given
-decision rule.
-
-@s DecisionRule int
-@s DecisionRuleImpl int
-@s FoldDecisionRule int
-@s UnfoldDecisionRule int
-@s ShockRealization int
-@s DRFixPoint int
-@s SimResults int
-@s SimResultsStats int
-@s SimResultsDynamicStats int
-@s RTSimResultsStats int
-@s SimResultsIRF int
-@s IRFResults int
-@s SimulationWorker int
-@s RTSimulationWorker int
-@s SimulationIRFWorker int
-@s RandomShockRealization int
-@s ExplicitShockRealization int
-@s GenShockRealization int
-@s IRFShockRealization int
-
-@c
-#ifndef DECISION_RULE_H
-#define DECISION_RULE_H
-
-#include <matio.h>
-  
-#include "kord_exception.h"
-#include "korder.h"
-#include "normal_conjugate.h"
-#include "mersenne_twister.h"
-
-@<|ShockRealization| class declaration@>;
-@<|DecisionRule| class declaration@>;
-@<|DecisionRuleImpl| class declaration@>;
-@<|FoldDecisionRule| class declaration@>;
-@<|UnfoldDecisionRule| class declaration@>;
-@<|DRFixPoint| class declaration@>;
-@<|SimResults| class declaration@>;
-@<|SimResultsStats| class declaration@>;
-@<|SimResultsDynamicStats| class declaration@>;
-@<|SimResultsIRF| class declaration@>;
-@<|RTSimResultsStats| class declaration@>;
-@<|IRFResults| class declaration@>;
-@<|SimulationWorker| class declaration@>;
-@<|SimulationIRFWorker| class declaration@>;
-@<|RTSimulationWorker| class declaration@>;
-@<|RandomShockRealization| class declaration@>;
-@<|ExplicitShockRealization| class declaration@>;
-@<|GenShockRealization| class declaration@>;
-
-#endif
-
-@ This is a general interface to a shock realizations. The interface
-has only one method returning the shock realizations at the given
-time. This method is not constant, since it may change a state of the
-object.
-
-@<|ShockRealization| class declaration@>=
-class ShockRealization {
-public:@;
-	virtual ~ShockRealization()@+ {}
-	virtual void get(int n, Vector& out) =0;
-	virtual int numShocks() const =0;
-};
-
-@ This class is an abstract interface to decision rule. Its main
-purpose is to define a common interface for simulation of a decision
-rule. We need only a simulate, evaluate, cetralized clone and output
-method. The |simulate| method simulates the rule for a given
-realization of the shocks. |eval| is a primitive evaluation (it takes
-a vector of state variables (predetermined, both and shocks) and
-returns the next period variables. Both input and output are in
-deviations from the rule's steady. |evaluate| method makes only one
-step of simulation (in terms of absolute values, not
-deviations). |centralizedClone| returns a new copy of the decision
-rule, which is centralized about provided fix-point. And finally
-|writeMat| writes the decision rule to the MAT file.
-
-@<|DecisionRule| class declaration@>=
-class DecisionRule {
-public:@;
-	enum emethod {@+ horner, trad @+};
-	virtual ~DecisionRule()@+ {}
-	virtual TwoDMatrix* simulate(emethod em, int np, const Vector& ystart,
-								 ShockRealization& sr) const =0;
-	virtual void eval(emethod em, Vector& out, const ConstVector& v) const =0;
-	virtual void evaluate(emethod em, Vector& out, const ConstVector& ys,
-						  const ConstVector& u) const =0;
-	virtual void writeMat(mat_t* fd, const char* prefix) const =0;
-	virtual DecisionRule* centralizedClone(const Vector& fixpoint) const =0;
-	virtual const Vector& getSteady() const =0;
-	virtual int nexog() const =0;
-	virtual const PartitionY& getYPart() const =0;
-};
-
-@ The main purpose of this class is to implement |DecisionRule|
-interface, which is a simulation. To be able to do this we have to
-know the partitioning of state vector $y$ since we will need to pick
-only predetermined part $y^*$. Also, we need to know the steady state.
-
-The decision rule will take the form: $$y_t-\bar
-y=\sum_{i=0}^n\left[g_{(yu)^i}\right]_{\alpha_1\ldots\alpha_i}\prod_{m=1}^i
-\left[\matrix{y^*_{t-1}-\bar y^*\cr u_t}\right]^{\alpha_m},$$ where
-the tensors $\left[g_{(yu)^i}\right]$ are tensors of the constructed
-container, and $\bar y$ is the steady state.
-
-If we know the fix point of the rule (conditional zero shocks)
-$\tilde y$, the rule can be transformed to so called ``centralized''
-form. This is very similar to the form above but the zero dimensional
-tensor is zero:
-$$y_t-\tilde y=\sum_{i=1}^n
-\left[\tilde g_{(yu)^i}\right]_{\alpha_1\ldots\alpha_i}\prod_{m=1}^i
-\left[\matrix{y^*_{t-1}-\tilde y^*\cr u_t}\right]^{\alpha_m}.$$
-We provide a method and a constructor to transform a rule to the centralized form.
-
-The class is templated, the template argument is either |KOrder::fold|
-or |KOrder::unfold|. So, there are two implementations of |DecisionRule| interface.
-
-@<|DecisionRuleImpl| class declaration@>=
-template <int t>
-class DecisionRuleImpl : public ctraits<t>::Tpol, public DecisionRule {
-protected:@;
-	typedef typename ctraits<t>::Tpol _Tparent;
-	const Vector ysteady;
-	const PartitionY ypart;
-	const int nu;
-public:@;
-	DecisionRuleImpl(const _Tparent& pol, const PartitionY& yp, int nuu,
-					 const Vector& ys)
-		: ctraits<t>::Tpol(pol), ysteady(ys), ypart(yp), nu(nuu)@+ {}
-	DecisionRuleImpl(_Tparent& pol, const PartitionY& yp, int nuu,
-					 const Vector& ys)
-		: ctraits<t>::Tpol(0, yp.ny(), pol), ysteady(ys), ypart(yp),
-		nu(nuu)@+ {}
-	DecisionRuleImpl(const _Tg& g, const PartitionY& yp, int nuu,
-					 const Vector& ys, double sigma)
-		: ctraits<t>::Tpol(yp.ny(), yp.nys()+nuu), ysteady(ys), ypart(yp), nu(nuu)
-		{@+ fillTensors(g, sigma);@+}
-	DecisionRuleImpl(const DecisionRuleImpl<t>& dr, const ConstVector& fixpoint)
-		: ctraits<t>::Tpol(dr.ypart.ny(), dr.ypart.nys()+dr.nu),
-	  	ysteady(fixpoint), ypart(dr.ypart), nu(dr.nu)
-		{@+ centralize(dr);@+}
-	const Vector& getSteady() const
-		{@+ return ysteady;@+}
-	@<|DecisionRuleImpl::simulate| code@>;
-	@<|DecisionRuleImpl::evaluate| code@>;
-	@<|DecisionRuleImpl::centralizedClone| code@>;
-	@<|DecisionRuleImpl::writeMat| code@>;
-	int nexog() const
-		{@+ return nu;@+}
-	const PartitionY& getYPart() const
-		{@+ return ypart;}
-protected:@;
-	@<|DecisionRuleImpl::fillTensors| code@>;
-	@<|DecisionRuleImpl::centralize| code@>;
-	@<|DecisionRuleImpl::eval| code@>;
-};
-
-@ Here we have to fill the tensor polynomial. This involves two
-separated actions. First is to evaluate the approximation at a given
-$\sigma$, the second is to compile the tensors $[g_{{(yu)}^{i+j}}]$ from
-$[g_{y^iu^j}]$. The first action is done here, the second is done by
-method |addSubTensor| of a full symmetry tensor.
-
-The way how the evaluation is done is described here:
-
-The $q-$order approximation to the solution can be written as:
-
-$$
-\eqalign{
-y_t-\bar y &= \sum_{l=1}^q{1\over l!}\left[\sum_{i+j+k=l}
-\left(\matrix{l\cr i,j,k}\right)\left[g_{y^iu^j\sigma^k}\right]
-_{\alpha_1\ldots\alpha_j\beta_1\ldots\beta_j}
-\prod_{m=1}^i[y^*_{t-1}-\bar y^*]^{\alpha_m}
-\prod_{n=1}^j[u_t]^{\beta_m}\sigma^k\right]\cr
- &= \sum_{l=1}^q\left[\sum_{i+j\leq l}\left(\matrix{i+j\cr i}\right)
-\left[\sum_{k=0}^{l-i-j}{1\over l!}
-\left(\matrix{l\cr k}\right)\left[g_{y^iu^j\sigma^k}\right]\sigma^k\right]
-\prod_{m=1}^i[y^*_{t-1}-\bar y^*]^{\alpha_m}
-\prod_{n=1}^j[u_t]^{\beta_m}\sigma^k\right]
-}
-$$
-
-This means that for each $i+j+k=l$ we have to add
-$${1\over l!}\left(\matrix{l\cr
-k}\right)\left[g_{y^iu^j\sigma^k}\right]\cdot\sigma^k=
-{1\over (i+j)!k!}\left[g_{y^iu^j\sigma^k}\right]\cdot\sigma^k$$ to
-$g_{(yu)^{i+j}}$. In addition, note that the multiplier
-$\left(\matrix{i+j\cr i}\right)$ is applied when the fully symmetric
-tensor $[g_{(yu)^{i+j}}]$ is evaluated.
-
-So we go through $i+j=d=0\ldots q$ and in each loop we form the fully
-symmetric tensor $[g_{(yu)^l}]$ and insert it to the container.
-
-@<|DecisionRuleImpl::fillTensors| code@>=
-void fillTensors(const _Tg& g, double sigma)
-{
-	IntSequence tns(2);
-	tns[0] = ypart.nys(); tns[1] = nu;
-	int dfact = 1;
-	for (int d = 0; d <= g.getMaxDim(); d++, dfact*=d) {
-		_Ttensym* g_yud = new _Ttensym(ypart.ny(), ypart.nys()+nu, d);
-		g_yud->zeros();
-		@<fill tensor of |g_yud| of dimension |d|@>;
-		this->insert(g_yud);
-	}
-}
-
-@ Here we have to fill the tensor $\left[g_{(yu)^d}\right]$. So we go
-through all pairs $(i,j)$ giving $i+j=d$, and through all $k$ from
-zero up to maximal dimension minus $d$. In this way we go through all
-symmetries of $g_{y^iu^j\sigma^k}$ which will be added to $g_{(yu)^d}$.
-
-Note that at the beginning, |dfact| is a factorial of |d|. We
-calculate |kfact| is equal to $k!$. As indicated in
-|@<|DecisionRuleImpl::fillTensors| code@>|, the added tensor is thus
-multiplied with ${1\over d!k!}\sigma^k$.
-
-@<fill tensor of |g_yud| of dimension |d|@>=
-	for (int i = 0; i <= d; i++) {
-		int j = d-i;
-		int kfact = 1;
-		_Ttensor tmp(ypart.ny(),
-					 TensorDimens(Symmetry(i,j), tns));
-		tmp.zeros();
-		for (int k = 0; k+d <= g.getMaxDim(); k++, kfact*=k) {
-			Symmetry sym(i,j,0,k);
-			if (g.check(sym)) {
-				double mult = pow(sigma,k)/dfact/kfact;
-				tmp.add(mult,*(g.get(sym)));
-			}
-		}
-		g_yud->addSubTensor(tmp);
-	}
-
-@ The centralization is straightforward. We suppose here that the
-object's steady state is the fix point $\tilde y$. It is clear that
-the new derivatives $\left[\tilde g_{(yu)^i}\right]$ will be equal to
-the derivatives of the original decision rule |dr| at the new steady
-state $\tilde y$. So, the new derivatives are obtained by derivating the
-given decision rule $dr$ and evaluating its polynomial at
-$$dstate=\left[\matrix{\tilde y^*-\bar y^*\cr 0}\right],$$
-where $\bar y$ is the steady state of the original rule |dr|.
-
-@<|DecisionRuleImpl::centralize| code@>=
-void centralize(const DecisionRuleImpl& dr)
-{
-	Vector dstate(ypart.nys() + nu);
-	dstate.zeros();
-	Vector dstate_star(dstate, 0, ypart.nys());
-	ConstVector newsteady_star(ysteady, ypart.nstat, ypart.nys());
-	ConstVector oldsteady_star(dr.ysteady, ypart.nstat, ypart.nys());
-	dstate_star.add(1.0, newsteady_star);
-	dstate_star.add(-1.0, oldsteady_star);
-
-	_Tpol pol(dr);
-	int dfac = 1;
-	for (int d = 1; d <= dr.getMaxDim(); d++, dfac *= d) {
-		pol.derivative(d-1);
-		_Ttensym* der = pol.evalPartially(d, dstate);
-		der->mult(1.0/dfac);
-		this->insert(der);
-	}
-}
-
-@ Here we evaluate repeatedly the polynomial storing results in the
-created matrix. For exogenous shocks, we use |ShockRealization|
-class, for predetermined variables, we use |ystart| as the first
-state. The |ystart| vector is required to be all state variables
-|ypart.ny()|, although only the predetermined part of |ystart| is
-used.
-
-We simulate in terms of $\Delta y$, this is, at the beginning the
-|ysteady| is canceled from |ystart|, we simulate, and at the end
-|ysteady| is added to all columns of the result.
-
-@<|DecisionRuleImpl::simulate| code@>=
-TwoDMatrix* simulate(emethod em, int np, const Vector& ystart,
-					 ShockRealization& sr) const
-{
-	KORD_RAISE_IF(ysteady.length() != ystart.length(),
-				  "Start and steady lengths differ in DecisionRuleImpl::simulate");
-	TwoDMatrix* res = new TwoDMatrix(ypart.ny(), np);
-
-	@<initialize vectors and subvectors for simulation@>;
-	@<perform the first step of simulation@>;
-	@<perform all other steps of simulations@>;
-	@<add the steady state to columns of |res|@>;
-	return res;
-}
-
-@ Here allocate the stack vector $(\Delta y^*, u)$, define the
-subvectors |dy|, and |u|, then we pickup predetermined parts of
-|ystart| and |ysteady|.
-
-@<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);
-
-
-@ We cancel |ysteady| from |ystart|, get realization to |u|, and
-evaluate the polynomial.
-
-@<perform the first step of simulation@>=
-	dy = ystart_pred;
-	dy.add(-1.0, ysteady_pred);
-	sr.get(0, u);
-	Vector out(*res, 0);
-	eval(em, out, dyu);
-
-@ Also clear. If the result at some period is not finite, we pad the
-rest of the matrix with zeros.
-
-@<perform all other steps of simulations@>=
-        int i=1;
-	while (i < np) {
-		ConstVector ym(*res, i-1);
-		ConstVector dym(ym, ypart.nstat, ypart.nys());
-		dy = dym;
-		sr.get(i, u);
-		Vector out(*res, i);
-		eval(em, out, dyu);
-		if (! out.isFinite()) {
-			if (i+1 < np) {
-				TwoDMatrix rest(*res, i+1, np-i-1);
-				rest.zeros();
-			}
-			break;
-		}
-		i++;
-	}
-
-@ Even clearer. We add the steady state to the numbers computed above
-and leave the padded columns to zero.
-
-@<add the steady state to columns of |res|@>=
-	for (int j = 0; j < i; j++) {
-		Vector col(*res, j);
-		col.add(1.0, ysteady);
-	}
-
-
-@ This is one period evaluation of the decision rule. The simulation
-is a sequence of repeated one period evaluations with a difference,
-that the steady state (fix point) is cancelled and added once. Hence
-we have two special methods.
-
-@<|DecisionRuleImpl::evaluate| code@>=
-void evaluate(emethod em, Vector& out, const ConstVector& ys,
-			  const ConstVector& u) const
-{
-	KORD_RAISE_IF(ys.length() != ypart.nys() || u.length() != nu,
-				  "Wrong dimensions of input vectors in DecisionRuleImpl::evaluate");
-	KORD_RAISE_IF(out.length() != ypart.ny(),
-				  "Wrong dimension of output vector in DecisionRuleImpl::evaluate");
-	ConstVector ysteady_pred(ysteady, ypart.nstat, ypart.nys());
-	Vector ys_u(ypart.nys()+nu);
-	Vector ys_u1(ys_u, 0, ypart.nys());
-	ys_u1 = ys;
-	ys_u1.add(-1.0, ysteady_pred);
-	Vector ys_u2(ys_u, ypart.nys(), nu);
-	ys_u2 = u;
-	eval(em, out, ys_u);
-	out.add(1.0, ysteady);
-}
-
-@ This is easy. We just return the newly created copy using the
-centralized constructor.
-
-@<|DecisionRuleImpl::centralizedClone| code@>=
-DecisionRule* centralizedClone(const Vector& fixpoint) const
-{
-	return new DecisionRuleImpl<t>(*this, fixpoint);
-}
-
-@ Here we only encapsulate two implementations to one, deciding
-according to the parameter.
-
-@<|DecisionRuleImpl::eval| code@>=
-void eval(emethod em, Vector& out, const ConstVector& v) const
-{
-	if (em == DecisionRule::horner)
-		_Tparent::evalHorner(out, v);
-	else
-		_Tparent::evalTrad(out, v);
-}
-
-@ Write the decision rule and steady state to the MAT file.
-@<|DecisionRuleImpl::writeMat| code@>=
-void writeMat(mat_t* fd, const char* prefix) const
-{
-	ctraits<t>::Tpol::writeMat(fd, prefix);
-	TwoDMatrix dum(ysteady.length(), 1);
-	dum.getData() = ysteady;
-	char tmp[100];
-	sprintf(tmp, "%s_ss", prefix);
-	ConstTwoDMatrix(dum).writeMat(fd, tmp);
-}
-
-@ This is exactly the same as |DecisionRuleImpl<KOrder::fold>|. The
-only difference is that we have a conversion from
-|UnfoldDecisionRule|, which is exactly
-|DecisionRuleImpl<KOrder::unfold>|.
-
-@<|FoldDecisionRule| class declaration@>=
-class UnfoldDecisionRule;
-class FoldDecisionRule : public DecisionRuleImpl<KOrder::fold> {
-	friend class UnfoldDecisionRule;
-public:@;
-	FoldDecisionRule(const ctraits<KOrder::fold>::Tpol& pol, const PartitionY& yp, int nuu,
-					 const Vector& ys)
-		: DecisionRuleImpl<KOrder::fold>(pol, yp, nuu, ys) {}
-	FoldDecisionRule(ctraits<KOrder::fold>::Tpol& pol, const PartitionY& yp, int nuu,
-					 const Vector& ys)
-		: DecisionRuleImpl<KOrder::fold>(pol, yp, nuu, ys) {}
-	FoldDecisionRule(const ctraits<KOrder::fold>::Tg& g, const PartitionY& yp, int nuu,
-					 const Vector& ys, double sigma)
-		: DecisionRuleImpl<KOrder::fold>(g, yp, nuu, ys, sigma) {}
-	FoldDecisionRule(const DecisionRuleImpl<KOrder::fold>& dr, const ConstVector& fixpoint)
-		: DecisionRuleImpl<KOrder::fold>(dr, fixpoint) {}
-	FoldDecisionRule(const UnfoldDecisionRule& udr);
-};
-
-@ This is exactly the same as |DecisionRuleImpl<KOrder::unfold>|, but
-with a conversion from |FoldDecisionRule|, which is exactly
-|DecisionRuleImpl<KOrder::fold>|.
-
-@<|UnfoldDecisionRule| class declaration@>=
-class UnfoldDecisionRule : public DecisionRuleImpl<KOrder::unfold> {
-	friend class FoldDecisionRule; 
-public:@;
-	UnfoldDecisionRule(const ctraits<KOrder::unfold>::Tpol& pol, const PartitionY& yp, int nuu,
-					 const Vector& ys)
-		: DecisionRuleImpl<KOrder::unfold>(pol, yp, nuu, ys) {}
-	UnfoldDecisionRule(ctraits<KOrder::unfold>::Tpol& pol, const PartitionY& yp, int nuu,
-					 const Vector& ys)
-		: DecisionRuleImpl<KOrder::unfold>(pol, yp, nuu, ys) {}
-	UnfoldDecisionRule(const ctraits<KOrder::unfold>::Tg& g, const PartitionY& yp, int nuu,
-					 const Vector& ys, double sigma)
-		: DecisionRuleImpl<KOrder::unfold>(g, yp, nuu, ys, sigma) {}
-	UnfoldDecisionRule(const DecisionRuleImpl<KOrder::unfold>& dr, const ConstVector& fixpoint)
-		: DecisionRuleImpl<KOrder::unfold>(dr, fixpoint) {}
-	UnfoldDecisionRule(const FoldDecisionRule& udr);
-};
-
-
-@ This class serves for calculation of the fix point of the decision
-rule given that the shocks are zero. The class is very similar to the
-|DecisionRuleImpl|. Besides the calculation of the fix point, the only
-difference between |DRFixPoint| and |DecisionRuleImpl| is that the
-derivatives wrt. shocks are ignored (since shocks are zero during the
-calculations). That is why have a different |fillTensor| method.
-
-The solution algorithm is Newton and is described in
-|@<|DRFixPoint::solveNewton| code@>|. It solves $F(y)=0$, where
-$F=g(y,0)-y$. The function $F$ is given by its derivatives |bigf|. The
-Jacobian of the solved system is given by derivatives stored in
-|bigfder|.
-
-@<|DRFixPoint| class declaration@>=
-template <int t>
-class DRFixPoint : public ctraits<t>::Tpol {
-	typedef typename ctraits<t>::Tpol _Tparent;
-	static int max_iter;
-	static int max_newton_iter;
-	static int newton_pause;
-	static double tol;
-	const Vector ysteady;
-	const PartitionY ypart;
-	_Tparent* bigf;
-	_Tparent* bigfder;
-public:@;
-	typedef typename DecisionRule::emethod emethod;
-	@<|DRFixPoint| constructor code@>;
-	@<|DRFixPoint| destructor code@>;
-	@<|DRFixPoint::calcFixPoint| code@>;
-	int getNumIter() const
-		{@+ return iter;@+}
-	int getNewtonLastIter() const
-		{@+ return newton_iter_last;@+}
-	int getNewtonTotalIter() const
-		{@+ return newton_iter_total;@+}
-protected:@;
-	@<|DRFixPoint::fillTensors| code@>;
-	@<|DRFixPoint::solveNewton| code@>;
-private:@;
-	int iter;
-	int newton_iter_last;
-	int newton_iter_total;	
-};
-
-
-@ Here we have to setup the function $F=g(y,0)-y$ and ${\partial
-F\over\partial y}$. The former is taken from the given derivatives of
-$g$ where a unit matrix is subtracted from the first derivative
-(|Symmetry(1)|). Then the derivative of the $F$ polynomial is
-calculated.
-
-@<|DRFixPoint| constructor code@>=
-DRFixPoint(const _Tg& g, const PartitionY& yp,
-					   const Vector& ys, double sigma)
-	: ctraits<t>::Tpol(yp.ny(), yp.nys()),
-	  ysteady(ys), ypart(yp), bigf(NULL), bigfder(NULL)
-{
-	fillTensors(g, sigma);
-	_Tparent yspol(ypart.nstat, ypart.nys(), *this);
-	bigf = new _Tparent((const _Tparent&) yspol);
-	_Ttensym* frst = bigf->get(Symmetry(1));
-	for (int i = 0; i < ypart.nys(); i++)
-		frst->get(i,i) = frst->get(i,i) - 1;
-	bigfder = new _Tparent(*bigf, 0);
-}
-
-@ 
-@<|DRFixPoint| destructor code@>=
-virtual ~DRFixPoint()
-{
-	if (bigf)
-		delete bigf;
-	if (bigfder)
-		delete bigfder;
-}
-
-@ Here we fill the tensors for the |DRFixPoint| class. We ignore the
-derivatives $g_{y^iu^j\sigma^k}$ for which $j>0$. So we go through all
-dimensions |d|, and all |k| such that |d+k| is between the maximum
-dimension and |d|, and add ${\sigma^k\over d!k!}g_{y^d\sigma^k}$ to
-the tensor $g_{(y)^d}$.
-
-@<|DRFixPoint::fillTensors| code@>=
-void fillTensors(const _Tg& g, double sigma)
-{
-	int dfact = 1;
-	for (int d = 0; d <= g.getMaxDim(); d++, dfact*=d) {
-		_Ttensym* g_yd = new _Ttensym(ypart.ny(), ypart.nys(), d);
-		g_yd->zeros();
-		int kfact = 1;
-		for (int k = 0; d+k <= g.getMaxDim(); k++, kfact*=k) {
-			if (g.check(Symmetry(d,0,0,k))) {
-				const _Ttensor* ten = g.get(Symmetry(d,0,0,k));
-				double mult = pow(sigma,k)/dfact/kfact;
-				g_yd->add(mult, *ten);
-			}
-		}
-		this->insert(g_yd);
-	}
-}
-
-@ This tries to solve polynomial equation $F(y)=0$, where $F$
-polynomial is |bigf| and its derivative is in |bigfder|. It returns
-true if the Newton converged. The method takes the given vector as
-initial guess, and rewrites it with a solution. The method guarantees
-to return the vector, which has smaller norm of the residual. That is
-why the input/output vector |y| is always changed.
-
-The method proceeds with a Newton step, if the Newton step improves
-the residual error. So we track residual errors in |flastnorm| and
-|fnorm| (former and current). In addition, at each step we search for
-an underrelaxation parameter |urelax|, which improves the residual. If
-|urelax| is less that |urelax_threshold|, we stop searching and stop
-the Newton.
-
-@<|DRFixPoint::solveNewton| code@>=
-bool solveNewton(Vector& y)
-{
-	const double urelax_threshold = 1.e-5;
-	Vector sol((const Vector&) y);
-	Vector delta(y.length());
-	newton_iter_last = 0;
-	bool delta_finite = true;
-	double flastnorm = 0.0;
-	double fnorm = 0.0;
-	bool converged = false;
-	double urelax = 1.0;
-
-	do {
-		_Ttensym* jacob = bigfder->evalPartially(1, sol);
-		bigf->evalHorner(delta, sol);
-		if (newton_iter_last == 0)
-			flastnorm = delta.getNorm();
-		delta_finite = delta.isFinite();
-		if (delta_finite) {
-			ConstTwoDMatrix(*jacob).multInvLeft(delta);
-			@<find |urelax| improving residual@>;
-			sol.add(-urelax, delta);
-			delta_finite = delta.isFinite();
-		}
-		delete jacob;
-		newton_iter_last++;
-		converged = delta_finite && fnorm < tol;
-		flastnorm = fnorm;
-	} while (!converged && newton_iter_last < max_newton_iter
-			 && urelax > urelax_threshold);
-
-	newton_iter_total += newton_iter_last;
-	if (! converged)
-		newton_iter_last = 0;
-	y = (const Vector&)sol; 
-	return converged;
-}
-
-@ Here we find the |urelax|. We cycle as long as the new residual size
-|fnorm| is greater than last residual size |flastnorm|. If the urelax
-is less than |urelax_threshold| we give up. The |urelax| is damped by
-the ratio of |flastnorm| and |fnorm|. It the ratio is close to one, we
-damp by one half.
-
-@<find |urelax| improving residual@>=
-	bool urelax_found = false;
-	urelax = 1.0;
-	while (!urelax_found && urelax > urelax_threshold) {
-		Vector soltmp((const Vector&)sol);
-		soltmp.add(-urelax, delta);
-		Vector f(sol.length());
-		bigf->evalHorner(f, soltmp);
-		fnorm = f.getNorm();
-		if (fnorm <= flastnorm)
-			urelax_found = true;
-		else
-			urelax *= std::min(0.5, flastnorm/fnorm);
-	}
-
-
-@ This method solves the fix point of the no-shocks rule
-$y_{t+1}=f(y_t)$. It combines dull steps with Newton attempts. The
-dull steps correspond to evaluations setting $y_{t+1}=f(y_t)$. For
-reasonable models the dull steps converge to the fix-point but very
-slowly. That is why we make Newton attempt from time to time. The
-frequency of the Newton attempts is given by |newton_pause|. We
-perform the calculations in deviations from the steady state. So, at
-the end, we have to add the steady state.
-
-The method also sets the members |iter|, |newton_iter_last| and
-|newton_iter_total|. These numbers can be examined later.
-
-The |out| vector is not touched if the algorithm has not convered.
-
-@<|DRFixPoint::calcFixPoint| code@>=
-bool calcFixPoint(emethod em, Vector& out)
-{
-	KORD_RAISE_IF(out.length() != ypart.ny(),
-				  "Wrong length of out in DRFixPoint::calcFixPoint");
-
-	Vector delta(ypart.nys());
-	Vector ystar(ypart.nys());
-	ystar.zeros();
-
-	iter = 0;
-	newton_iter_last = 0;
-	newton_iter_total = 0;
-	bool converged = false;
-	do {
-		if ((iter/newton_pause)*newton_pause == iter)
-			converged = solveNewton(ystar);
-		if (! converged) {
-			bigf->evalHorner(delta, ystar);
-			KORD_RAISE_IF_X(! delta.isFinite(),
-							"NaN or Inf asserted in DRFixPoint::calcFixPoint",
-							KORD_FP_NOT_FINITE);
-			ystar.add(1.0, delta);
-			converged = delta.getNorm() < tol;
-		}
-		iter++;
-	} while (iter < max_iter && ! converged);
-
-	if (converged) {
-		_Tparent::evalHorner(out, ystar);
-		out.add(1.0, ysteady);
-	}
-
-	return converged;
-}
-
-
-@ 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.
-
-@<|SimResults| class declaration@>=
-class ExplicitShockRealization;
-class SimResults {
-protected:@;
-	int num_y;
-	int num_per;
-	int num_burn;
-	vector<TwoDMatrix*> data;
-	vector<ExplicitShockRealization*> shocks;
-public:@;
-	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);
-	void simulate(int num_sim, const DecisionRule& dr, const Vector& start,
-				  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
-		{@+ return *(data[i]);@+}
-	const ExplicitShockRealization& getShocks(int i) const
-		{ @+ return *(shocks[i]);@+}
-	bool addDataSet(TwoDMatrix* d, ExplicitShockRealization* sr);
-	void writeMat(const char* base, const char* lname) const;
-	void writeMat(mat_t* fd, const char* lname) const;
-};
-
-@ This does the same as |SimResults| plus it calculates means and
-covariances of the simulated data.
-
-@<|SimResultsStats| class declaration@>=
-class SimResultsStats : public SimResults {
-protected:@;
-	Vector mean;
-	TwoDMatrix vcov;
-public:@;
-	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 writeMat(mat_t* fd, const char* lname) const;
-protected:@;
-	void calcMean();
-	void calcVcov();
-};
-
-@ This does the similar thing as |SimResultsStats| but the statistics are
-not calculated over all periods but only within each period. Then we
-do not calculate covariances with periods but only variances.
-
-@<|SimResultsDynamicStats| class declaration@>=
-class SimResultsDynamicStats : public SimResults {
-protected:@;
-	TwoDMatrix mean;
-	TwoDMatrix variance;
-public:@;
-	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 writeMat(mat_t* fd, const char* lname) const; 
-protected:@;
-	void calcMean();
-	void calcVariance();
-};
-
-
-@ This goes through control simulation results, and for each control
-it adds a given impulse to a given shock and runs a simulation. The
-control simulation is then cancelled and the result is stored. After
-that these results are averaged with variances calculated.
-
-The means and the variances are then written to the MAT-4 file.
-
-@<|SimResultsIRF| class declaration@>=
-class SimulationIRFWorker;
-class SimResultsIRF : public SimResults {
-	friend class SimulationIRFWorker;
-protected:@;
-	const SimResults& control;
-	int ishock;
-	double imp;
-	TwoDMatrix means;
-	TwoDMatrix variances;
-public:@;
-	SimResultsIRF(const SimResults& cntl, int ny, int nper, int i, double impulse)
-		: SimResults(ny, nper, 0), control(cntl),
-		  ishock(i), imp(impulse),
-		  means(ny, nper), variances(ny, nper)@+ {}
-	void simulate(const DecisionRule& dr, Journal& journal);
-	void simulate(const DecisionRule& dr);
-	void writeMat(mat_t* fd, const char* lname) const;
-protected:@;
-	void calcMeans();
-	void calcVariances();
-};
-
-@ This simulates and gathers all statistics from the real time
-simulations. In the |simulate| method, it runs |RTSimulationWorker|s
-which accummulate information from their own estimates. The estimation
-is done by means of |NormalConj| class, which is a conjugate family of
-densities for normal distibutions.
-
-@<|RTSimResultsStats| class declaration@>=
-class RTSimulationWorker;
-class RTSimResultsStats {
-	friend class RTSimulationWorker;
-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, int nburn = 0)
-		: mean(ny), vcov(ny, 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);
-	void simulate(int num_sim, const DecisionRule& dr, const Vector& start,
-				  const TwoDMatrix& vcov);
-	void writeMat(mat_t* fd, const char* lname);
-};
-
-@ For each shock, this simulates plus and minus impulse. The class
-maintains a vector of simulation results, each gets a particular shock
-and sign (positive/negative). The results of type |SimResultsIRF| are
-stored in a vector so that even ones are positive, odd ones are
-negative.
-
-The constructor takes a reference to the control simulations, which
-must be finished before the constructor is called. The control
-simulations are passed to all |SimResultsIRF|s.
-
-The constructor also takes the vector of indices of exogenous
-variables (|ili|) for which the IRFs are generated. The list is kept
-(as |irf_list_ind|) for other methods.
- 
-@<|IRFResults| class declaration@>=
-class DynamicModel;
-class IRFResults {
-	vector<SimResultsIRF*> irf_res;
-	const DynamicModel& model;
-	vector<int> irf_list_ind;
-public:@;
-	IRFResults(const DynamicModel& mod, const DecisionRule& dr,
-			   const SimResults& control, const vector<int>& ili,
-			   Journal& journal);
-	~IRFResults();
-	void writeMat(mat_t* fd, const char* prefix) const;
-};
-
-@ This worker simulates the given decision rule and inserts the result
-to |SimResults|.
-
-@<|SimulationWorker| class declaration@>=
-class SimulationWorker : public THREAD {
-protected:@;
-	SimResults& res;
-	const DecisionRule& dr;
-	DecisionRule::emethod em;
-	int np;
-	const Vector& st;
-	ShockRealization& sr;
-public:@;
-	SimulationWorker(SimResults& sim_res,
-					 const DecisionRule& dec_rule,
-					 DecisionRule::emethod emet, int num_per,
-					 const Vector& start, ShockRealization& shock_r)
-		: res(sim_res), dr(dec_rule), em(emet), np(num_per), st(start), sr(shock_r) {}
-	void operator()();
-};
-
-@ This worker simulates a given impulse |imp| to a given shock
-|ishock| based on a given control simulation with index |idata|. The
-control simulations are contained in |SimResultsIRF| which is passed
-to the constructor.
-
-@<|SimulationIRFWorker| class declaration@>=
-class SimulationIRFWorker : public THREAD {
-	SimResultsIRF& res;
-	const DecisionRule& dr;
-	DecisionRule::emethod em;
-	int np;
-	int idata;
-	int ishock;
-	double imp;	
-public:@;
-	SimulationIRFWorker(SimResultsIRF& sim_res,
-						const DecisionRule& dec_rule,
-						DecisionRule::emethod emet, int num_per,
-						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()();
-};
-
-@ This class does the real time simulation job for
-|RTSimResultsStats|. It simulates the model period by period. It
-accummulates the information in the |RTSimResultsStats::nc|. If NaN or
-Inf is observed, it ends the simulation and adds to the
-|thrown_periods| of |RTSimResultsStats|.
-
-@<|RTSimulationWorker| class declaration@>=
-class RTSimulationWorker : public THREAD {
-protected:@;
-	RTSimResultsStats& res;
-	const DecisionRule& dr;
-	DecisionRule::emethod em;
-	int np;
-	const Vector& ystart;
-	ShockRealization& sr;
-public:@;
-	RTSimulationWorker(RTSimResultsStats& sim_res,
-					   const DecisionRule& dec_rule,
-					   DecisionRule::emethod emet, int num_per,
-					   const Vector& start, ShockRealization& shock_r)
-		: res(sim_res), dr(dec_rule), em(emet), np(num_per), ystart(start), sr(shock_r) {}
-	void operator()();
-};
-
-@ This class generates draws from Gaussian distribution with zero mean
-and the given variance-covariance matrix. It stores the factor of vcov
-$V$ matrix, yielding $FF^T = V$.
- 
-@<|RandomShockRealization| class declaration@>=
-class RandomShockRealization : virtual public ShockRealization {
-protected:@;
-	MersenneTwister mtwister;
-	TwoDMatrix factor;
-public:@;
-	RandomShockRealization(const TwoDMatrix& v, unsigned int iseed)
-		: mtwister(iseed), factor(v.nrows(),v.nrows())
-		{@+schurFactor(v);@+}
-	RandomShockRealization(const RandomShockRealization& sr)
-		: mtwister(sr.mtwister), factor(sr.factor)@+ {}
-	virtual ~RandomShockRealization() @+{}
-	void get(int n, Vector& out);
-	int numShocks() const
-		{@+ return factor.nrows();@+}
-protected:@;
-	void choleskyFactor(const TwoDMatrix& v);
-	void schurFactor(const TwoDMatrix& v);
-};
-
-@ This is just a matrix of finite numbers. It can be constructed from
-any |ShockRealization| with a given number of periods.
-
-@<|ExplicitShockRealization| class declaration@>=
-class ExplicitShockRealization : virtual public ShockRealization {
-	TwoDMatrix shocks;
-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();@+}
-};
-
-@ This represents a user given shock realization. The first matrix of
-the constructor is a covariance matrix of shocks, the second matrix is
-a rectangular matrix, where columns correspond to periods, rows to
-shocks. If an element of the matrix is {\tt NaN}, or {\tt Inf}, or
-{\tt -Inf}, then the random shock is taken instead of that element.
-
-In this way it is a generalization of both |RandomShockRealization|
-and |ExplicitShockRealization|.
-
-@<|GenShockRealization| class declaration@>=
-class GenShockRealization : public RandomShockRealization, public ExplicitShockRealization {
-public:@;
-	GenShockRealization(const TwoDMatrix& v, const TwoDMatrix& sh, int seed)
-		: RandomShockRealization(v, seed), ExplicitShockRealization(sh)@+
-		{
-			KORD_RAISE_IF(sh.nrows() != v.nrows() || v.nrows() != v.ncols(),
-						  "Wrong dimension of input matrix in GenShockRealization constructor");
-		}
-	void get(int n, Vector& out);
-	int numShocks() const
-		{@+ return RandomShockRealization::numShocks();@+}
-};
-
-@ End of {\tt decision\_rule.h} file.

dynare++/kord/dynamic_model.cweb

-@q $Id: dynamic_model.cweb 431 2005-08-16 15:41:01Z kamenik $ @>
-@q Copyright 2005, Ondra Kamenik @>
-
-@ Start of {\tt dynamic\_model.cpp} file.
-
-@c
-#include "dynamic_model.h"
-
-@<|NameList::print| code@>;
-@<|NameList::writeMat| code@>;
-@<|NameList::writeMatIndices| code@>;
-
-@ 
-@<|NameList::print| code@>=
-void NameList::print() const
-{
-	for (int i = 0; i < getNum(); i++)
-		printf("%s\n", getName(i));
-}
-
-@ 
-@<|NameList::writeMat| code@>=
-void NameList::writeMat(mat_t* fd, const char* vname) const
-{
-	int maxlen = 0;
-	for (int i = 0; i < getNum(); i++)
-		if (maxlen < (int)strlen(getName(i)))
-			maxlen = (int)strlen(getName(i));
-
-	if (maxlen == 0)
-		return;
-
-  char *m = new char[getNum()*maxlen];
-
-	for (int i = 0; i < getNum(); i++)
-		for (int j = 0; j < maxlen; j++)
-			if (j < (int)strlen(getName(i)))
-				m[j*getNum()+i] = getName(i)[j];
-			else
-				m[j*getNum()+i] = ' ';
-
-# if MATIO_MAJOR_VERSION > 1 || (MATIO_MAJOR_VERSION == 1 && MATIO_MINOR_VERSION >= 5)
-  size_t dims[2];
-  const matio_compression compression = MAT_COMPRESSION_NONE;
-# else
-  int dims[2];
-  const int compression = COMPRESSION_NONE;
-# endif
-  dims[0] = getNum();
-  dims[1] = maxlen;
-
-  matvar_t *v = Mat_VarCreate(vname, MAT_C_CHAR, MAT_T_UINT8, 2, dims, m, 0);
-
-  Mat_VarWrite(fd, v, compression);
-
-  Mat_VarFree(v);
-  delete[] m;
-}
-
-@ 
-@<|NameList::writeMatIndices| code@>=
-void NameList::writeMatIndices(mat_t* fd, const char* prefix) const
-{
-	char tmp[100];
-	TwoDMatrix aux(1,1);
-	for (int i = 0; i < getNum(); i++) {
-		sprintf(tmp, "%s_i_%s", prefix, getName(i));
-		aux.get(0,0) = i+1;
-		aux.writeMat(fd, tmp);
-	}
-}
-
-@ End of {\tt dynamic\_model.cpp} file.

dynare++/kord/dynamic_model.hweb

-@q $Id: dynamic_model.hweb 378 2005-07-21 15:50:20Z kamenik $ @>
-@q Copyright 2005, Ondra Kamenik @>
-
-@*2 Dynamic model abstraction. Start of {\tt dynamic\_model.h} file.
-
-This file only defines a generic interface to an SDGE model. The model
-takes the form:
-$$E_t\left[f(g^{**}(g^*(y,u_t),u_{t+1}),g(y,u),y,u_t)\right]=0$$
-The interface is defined via pure virtual class |DynamicModel|.
-
-@s NameList int
-@s DynamicModel int
-@c
-#ifndef DYNAMIC_MODEL_H
-#define DYNAMIC_MODEL_H
-
-#include "t_container.h"
-#include "sparse_tensor.h"
-
-#include "Vector.h"
-
-@<|NameList| class declaration@>;
-@<|DynamicModel| class declaration@>;
-
-#endif
-
-@ The class is a virtual pure class which provides an access to names
-of the variables.
-@<|NameList| class declaration@>=
-class NameList {
-public:@;
-	virtual ~NameList() {}
-	virtual int getNum() const =0;
-	virtual const char* getName(int i) const=0;
-	void print() const;
-	void writeMat(mat_t* fd, const char* vname) const;
-	void writeMatIndices(mat_t* fd, const char* prefix) const;
-};
-
-@ This is the interface to an information on a generic SDGE
-model. It is sufficient for calculations of policy rule Taylor
-approximations at some (not necessarily deterministic) steady state.
-
-We need to know a partitioning of endogenous variables $y$. We suppose
-that $y$ is partitioned as
-$$y=\left[\matrix{\hbox{static}\cr\hbox{pred}\cr\hbox{both}\cr\hbox{forward}}\right]$$
-of which we define
-$$y^*=\left[\matrix{\hbox{pred}\cr\hbox{both}}\right]\quad
-y^{**}=\left[\matrix{\hbox{both}\cr\hbox{forward}}\right]$$
-where ``static'' are meant those variables, which appear only at time
-$t$; ``pred'' are meant those variables, which appear only at $t$ and
-$t-1$; ``both'' are meant those variables, which appear at least at
-$t-1$ and $t+1$; and ``forward'' are meant those variables, which
-appear only at $t$ and $t+1$. This partitioning is given by methods
-|nstat()|, |npred()|, |nboth()|, and |nforw()|. The number of
-equations |numeq()| must be the same as a number of endogenous
-variables.
-
-In order to complete description, we need to know a number of
-exogenous variables, which is a size of $u$, hence |nexog()| method.
-
-The model contains an information about names of variables, the
-variance-covariance matrix of the shocks, the derivatives of equations
-of $f$ at some steady state, and the steady state. These can be
-retrieved by the corresponding methods.
-
-The derivatives of the system are calculated with respect to stacked
-variables, the stack looks as:
-$$\left[\matrix{y^{**}_{t+1}\cr y_t\cr y^*_{t-1}\cr u_t}\right].$$
-
-There are only three operations. The first
-|solveDeterministicSteady()| solves the deterministic steady steate
-which can be retrieved by |getSteady()| later. The method
-|evaluateSystem| calculates $f(y^{**},y,y^*,u)$, where $y$ and $u$ are
-passed, or $f(y^{**}_{t+1}, y_t, y^*_{t-1}, u)$, where $y^{**}_{t+1}$,
-$y_t$, $y^*_{t-1}$, $u$ are passed. Finally, the method
-|calcDerivativesAtSteady()| calculates derivatives of $f$ at the
-current steady state, and zero shocks. The derivatives can be
-retrieved with |getModelDerivatives()|. All the derivatives are done
-up to a given order in the model, which can be retrieved by |order()|.
-
-The model initialization is done in a constructor of the implementing
-class. The constructor usually calls a parser, which parses a given
-file (usually a text file), and retrieves all necessary information
-about the model, inluding variables, partitioning, variance-covariance
-matrix, information helpful for calculation of the deterministic
-steady state, and so on.
-
-@<|DynamicModel| class declaration@>=
-class DynamicModel {
-public:@;
-	virtual DynamicModel* clone() const =0;
-	virtual ~DynamicModel() {}
-
-	virtual int nstat() const =0;
-	virtual int nboth() const =0;
-	virtual int npred() const =0;
-	virtual int nforw() const =0;
-	virtual int nexog() const =0;
-	virtual int order() const =0;
-	int numeq() const
-		{@+ return nstat()+nboth()+npred()+nforw(); @+}
-
-	virtual const NameList& getAllEndoNames() const =0;
-	virtual const NameList& getStateNames() const =0;
-	virtual const NameList& getExogNames() const =0;
-	virtual const TwoDMatrix& getVcov() const =0;
-	virtual const TensorContainer<FSSparseTensor>& getModelDerivatives() const =0;
-	virtual const Vector& getSteady() const =0;
-	virtual Vector& getSteady() =0;
-
-	virtual void solveDeterministicSteady() =0;
-	virtual void evaluateSystem(Vector& out, const Vector& yy, const Vector& xx) =0;
-	virtual void evaluateSystem(Vector& out, const Vector& yym, const Vector& yy,
-								const Vector& yyp, const Vector& xx) =0;
-	virtual void calcDerivativesAtSteady() =0;
-};
-
-
-@ End of {\tt dynamic\_model.h} file.

dynare++/kord/faa_di_bruno.cweb

-@q $Id: faa_di_bruno.cweb 744 2006-05-09 13:16:07Z kamenik $ @>
-@q Copyright 2005, Ondra Kamenik @>
-
-@ Start of {\tt faa\_di\_bruno.cpp} file.
-
-@s FoldedFineContainer int
-@s UnfoldedFineContainer int
-
-@c
-#include "faa_di_bruno.h"
-#include "fine_container.h"
-
-#include <cmath>
-
-double FaaDiBruno::magic_mult = 1.5;
-@<|FaaDiBruno::calculate| folded sparse code@>;
-@<|FaaDiBruno::calculate| folded dense code@>;
-@<|FaaDiBruno::calculate| unfolded sparse code@>;
-@<|FaaDiBruno::calculate| unfolded dense code@>;
-@<|FaaDiBruno::estimRefinment| code@>;
-
-@ We take an opportunity to refine the stack container to avoid
-allocation of more memory than available.
-
-@<|FaaDiBruno::calculate| folded sparse code@>=
-void FaaDiBruno::calculate(const StackContainer<FGSTensor>& cont,
-						   const TensorContainer<FSSparseTensor>& f,
-						   FGSTensor& out)
-{
-	out.zeros();
-	for (int l = 1; l <= out.dimen(); l++) {
-		int mem_mb, p_size_mb;
-		int max = estimRefinment(out.getDims(), out.nrows(), l, mem_mb, p_size_mb);
-		FoldedFineContainer fine_cont(cont, max);
-		fine_cont.multAndAdd(l, f, out);
-		JournalRecord recc(journal);
-		recc << "dim=" << l << " avmem=" << mem_mb << " tmpmem=" << p_size_mb << " max=" << max
-			 << " stacks=" << cont.numStacks() << "->" << fine_cont.numStacks() << endrec;
-	}
-}
-
-@ Here we just simply evaluate |multAndAdd| for the dense
-container. There is no opportunity for tuning.
-
-@<|FaaDiBruno::calculate| folded dense code@>=
-void FaaDiBruno::calculate(const FoldedStackContainer& cont, const FGSContainer& g,
-						   FGSTensor& out)
-{
-	out.zeros();
-	for (int l = 1; l <= out.dimen(); l++) {
-		long int mem = SystemResources::availableMemory();
-		cont.multAndAdd(l, g, out);
-		JournalRecord rec(journal);
-		int mem_mb = mem/1024/1024;
-		rec << "dim=" << l << " avmem=" << mem_mb << endrec;
-	}
-}
-
-@ This is the same as |@<|FaaDiBruno::calculate| folded sparse
-code@>|. The only difference is that we construct unfolded fine
-container.
-
-@<|FaaDiBruno::calculate| unfolded sparse code@>=
-void FaaDiBruno::calculate(const StackContainer<UGSTensor>& cont,
-						   const TensorContainer<FSSparseTensor>& f,
-						   UGSTensor& out)
-{
-	out.zeros();
-	for (int l = 1; l <= out.dimen(); l++) {
-		int mem_mb, p_size_mb;
-		int max = estimRefinment(out.getDims(), out.nrows(), l, mem_mb, p_size_mb);
-		UnfoldedFineContainer fine_cont(cont, max);
-		fine_cont.multAndAdd(l, f, out);
-		JournalRecord recc(journal);
-		recc << "dim=" << l << " avmem=" << mem_mb << " tmpmem=" << p_size_mb << " max=" << max
-			 << " stacks=" << cont.numStacks() << "->" << fine_cont.numStacks() << endrec;
-	}
-}
-
-@ Again, no tuning opportunity here.
-@<|FaaDiBruno::calculate| unfolded dense code@>=
-void FaaDiBruno::calculate(const UnfoldedStackContainer& cont, const UGSContainer& g,
-					   UGSTensor& out)
-{
-	out.zeros();
-	for (int l = 1; l <= out.dimen(); l++) {
-		long int mem = SystemResources::availableMemory();
-		cont.multAndAdd(l, g, out);
-		JournalRecord rec(journal);
-		int mem_mb = mem/1024/1024;
-		rec << "dim=" << l << " avmem=" << mem_mb << endrec;
-	}
-}
-
-@ This function returns a number of maximum rows used for refinement of
-the stacked container. We want to set the maximum so that the expected
-memory consumption for the number of paralel threads would be less
-than available memory. On the other hand we do not want to be too
-pesimistic since a very fine refinement can be very slow.
-
-Besides memory needed for a dense unfolded slice of a tensor from
-|f|, each thread needs |magic_mult*per_size| bytes of memory. In the
-worst case, |magic_mult| will be equal to 2, this means memory
-|per_size| for target temporary (permuted symmetry) tensor plus one
-copy for intermediate result. However, this shows to be too
-pesimistic, so we set |magic_mult| to 1.5. The memory for permuted
-symmetry temporary tensor |per_size| is estimated as a weigthed
-average of unfolded memory of the |out| tensor and unfolded memory of
-a symetric tensor with the largest coordinate size. Some experiments
-showed that the best combination of the two is to take 100\% if the
-latter, so we set |lambda| to zero.
-
-The |max| number of rows in the refined |cont| must be such that each
-slice fits to remaining memory. Number of columns of the slice are
-never greater $max^l$. (This is not true, since stacks corresponing to
-unit/zero matrices cannot be further refined). We get en equation:
-
-$$nthreads\cdot max^l\cdot 8\cdot r = mem - 
-magic\_mult\cdot nthreads\cdot per\_size\cdot 8\cdot r,$$
-where |mem| is available memory in bytes, |nthreads| is a number of
-threads, $r$ is a number of rows, and $8$ is |sizeof(double)|.
-
-If the right hand side is less than zero, we set |max| to 10, just to
-let it do something.
-
-@<|FaaDiBruno::estimRefinment| code@>=
-int FaaDiBruno::estimRefinment(const TensorDimens& tdims, int nr, int l,
-							   int& avmem_mb, int& tmpmem_mb)
-{
-	int nthreads = THREAD_GROUP::max_parallel_threads;
-	long int per_size1 = tdims.calcUnfoldMaxOffset();
-	long int per_size2 = (long int)pow((double)tdims.getNVS().getMax(), l);
-	double lambda = 0.0;
-	long int per_size = sizeof(double)*nr
-		*(long int)(lambda*per_size1+(1-lambda)*per_size2);
-	long int mem = SystemResources::availableMemory();
-	int max = 0;
-	double num_cols = ((double)(mem-magic_mult*nthreads*per_size))
-		/nthreads/sizeof(double)/nr;
-	if (num_cols > 0) {
-		double maxd = pow(num_cols, ((double)1)/l);
-		max = (int)floor(maxd);
-	}
-	if (max == 0) {
-		max = 10;
-		JournalRecord rec(journal);
-		rec << "dim=" << l << " run out of memory, imposing max=" << max;
-		if (nthreads > 1)
-			rec << " (decrease number of threads)";
-		rec << endrec;
-	}
-	avmem_mb = mem/1024/1024;
-	tmpmem_mb = (nthreads*per_size)/1024/1024;
-	return max;
-}
-
-
-@ End of {\tt faa\_di\_bruno.cpp} file.

dynare++/kord/faa_di_bruno.hweb

-@q $Id: faa_di_bruno.hweb 744 2006-05-09 13:16:07Z kamenik $ @>
-@q Copyright 2005, Ondra Kamenik @>
-
-@*2 Faa Di Bruno evaluator. Start of {\tt faa\_di\_bruno.h} file.
-
-This defines a class which implements Faa Di Bruno Formula
-$$\left[B_{s^k}\right]_{\alpha_1\ldots\alpha_l}=\left[f_{z^l}\right]_{\beta_1\ldots\beta_l}
-\sum_{c\in M_{l,k}}\prod_{m=1}^l\left[z_{s^k(c_m)}\right]^{\beta_m}_{c_m(\alpha)}$$
-where $s^k$ is a general symmetry of dimension $k$ and $z$ is a stack of functions.
-
-@s FaaDiBruno int
-
-@c
-#ifndef FAA_DI_BRUNO_H
-#define FAA_DI_BRUNO_H
-
-#include "journal.h"
-#include "stack_container.h"
-#include "t_container.h"
-#include "sparse_tensor.h"
-#include "gs_tensor.h"
-
-@<|FaaDiBruno| class declaration@>;
-
-#endif
-
-@ Nothing special here. See |@<|FaaDiBruno::calculate| folded sparse
-code@>| for reason of having |magic_mult|.
-
-@<|FaaDiBruno| class declaration@>=
-class FaaDiBruno {
-	Journal& journal;
-public:@;
-	FaaDiBruno(Journal& jr)
-		: journal(jr)@+ {}
-	void calculate(const StackContainer<FGSTensor>& cont, const TensorContainer<FSSparseTensor>& f,
-				   FGSTensor& out);
-	void calculate(const FoldedStackContainer& cont, const FGSContainer& g,
-				   FGSTensor& out);
-	void calculate(const StackContainer<UGSTensor>& cont, const TensorContainer<FSSparseTensor>& f,
-				   UGSTensor& out);
-	void calculate(const UnfoldedStackContainer& cont, const UGSContainer& g,
-				   UGSTensor& out);
-protected:@;
-	int estimRefinment(const TensorDimens& tdims, int nr, int l, int& avmem_mb, int& tmpmem_mb);
-	static double magic_mult;
-};
-
-@ End of {\tt faa\_di\_bruno.h} file.

dynare++/kord/first_order.cweb

-@q $Id: first_order.cweb 2351 2009-09-03 14:58:03Z kamenik $ @>
-@q Copyright 2004, Ondra Kamenik @>
-
-@ Start of {\tt first\_order.cpp} file.
-
-@c
-
-#include "kord_exception.h"
-#include "first_order.h"
-
-#include <dynlapack.h>
-
-double qz_criterium = 1.000001;
-@<|order_eigs| function code@>;
-@<|FirstOrder::solve| code@>;
-@<|FirstOrder::journalEigs| code@>;
-
-@ This is a function which selects the eigenvalues pair used by
-|dgges|. See documentation to DGGES for details. Here we want
-to select (return true) the pairs for which $\alpha<\beta$.
-
-@<|order_eigs| function code@>=
-lapack_int order_eigs(const double* alphar, const double* alphai, const double* beta)
-{
-	return (*alphar * *alphar + *alphai * *alphai < *beta * *beta * qz_criterium);
-}
-
-
-@ Here we solve the linear approximation. The result are the matrices
-$g_{y^*}$ and $g_u$. The method solves the first derivatives of $g$ so
-that the following equation would be true:
-$$E_t[F(y^*_{t-1},u_t,u_{t+1},\sigma)] =
-E_t[f(g^{**}(g^*(y_{t-1}^*,u_t,\sigma), u_{t+1}, \sigma), g(y_{t-1}^*,u_t,\sigma),
-y^*_{t-1},u_t)]=0$$
-where $f$ is a given system of equations.
-
-It is known that $g_{y^*}$ is given by $F_{y^*}=0$, $g_u$ is given by
-$F_u=0$, and $g_\sigma$ is zero. The only input to the method are the
-derivatives |fd| of the system $f$, and partitioning of the vector $y$
-(from object).
-
-@<|FirstOrder::solve| code@>=
-void FirstOrder::solve(const TwoDMatrix& fd)
-{
-	JournalRecordPair pa(journal);
-	pa << "Recovering first order derivatives " << endrec;
-
-	::qz_criterium = FirstOrder::qz_criterium;
-
-	@<solve derivatives |gy|@>;
-	@<solve derivatives |gu|@>;
-	journalEigs();
-
-	if (! gy.isFinite() || ! gu.isFinite()) {
-		throw KordException(__FILE__, __LINE__,
-						  "NaN or Inf asserted in first order derivatives in FirstOrder::solve");
-	}
-}
-
-@ The derivatives $g_{y^*}$ are retrieved from the equation
-$F_{y^*}=0$. The calculation proceeds as follows:
-
-\orderedlist
-
-\li For each variable appearing at both $t-1$ and $t-1$ we add a dummy
-variable, so that the predetermined variables and forward looking would
-be disjoint. This is, the matrix of the first derivatives of the
-system written as:
-$$\left[\matrix{f_{y^{**}_+}&f_{ys}&f_{yp}&f_{yb}&f_{yf}&f_{y^*_-}}\right],$$
-where $f_{ys}$, $f_{yp}$, $f_{yb}$, and $f_{yf}$ are derivatives wrt
-static, predetermined, both, forward looking at time $t$, is rewritten
-to the matrix:
-$$\left[
-\matrix{f_{y^{**}_+}&f_{ys}&f_{yp}&f_{yb}&0&f_{yf}&f_{y^*_-}\cr
-        0           &0     &0     &I   &-I&0    &0}
- \right],$$
-where the second line has number of rows equal to the number of both variables.
-
-\li Next, provided that forward looking and predetermined are
-disjoint, the equation $F_{y^*}=0$ is written as:
-$$\left[f_+{y^{**}_+}\right]\left[g_{y^*}^{**}\right]\left[g_{y^*}^*\right]
-+\left[f_{ys}\right]\left[g^s_{y^*}\right]
-+\left[f_{y^*}\right]\left[g^*_{y^*}\right]
-+\left[f_{y^{**}}\right]\left[g^{**}_{y^*}\right]+\left[f_{y^*_-}\right]=0$$
-This is rewritten as
-$$\left[\matrix{f_{y^*}&0&f_{y^{**}_+}}\right]
-\left[\matrix{I\cr g^s_{y^*}\cr g^{**}_{y^*}}\right]\left[g_{y^*}^*\right]+ 
-\left[\matrix{f_{y^*_-}&f_{ys}&f_{y^{**}}}\right]
-\left[\matrix{I\cr g^s_{y^*}\cr g^{**}_{y^*}}\right]=0
-$$
-Now, in the above equation, there are the auxiliary variables standing
-for copies of both variables at time $t+1$. This equation is then
-rewritten as:
-$$
-\left[\matrix{f_{yp}&f_{yb}&0&f_{y^{**}_+}\cr 0&I&0&0}\right]
-\left[\matrix{I\cr g^s_{y^*}\cr g^{**}_{y^*}}\right]\left[g_{y^*}^*\right]+ 
-\left[\matrix{f_{y^*_-}&f_{ys}&0&f_{yf}\cr 0&0&-I&0}\right]
-\left[\matrix{I\cr g^s_{y^*}\cr g^{**}_{y^*}}\right]=0
-$$
-The two matrices are denoted as $D$ and $-E$, so the equation takes the form:
-$$D\left[\matrix{I\cr g^s_{y^*}\cr g^{**}_{y^*}}\right]\left[g_{y^*}^*\right]=
-E\left[\matrix{I\cr g^s_{y^*}\cr g^{**}_{y^*}}\right]$$
-
-\li Next we solve the equation by Generalized Schur decomposition:
-$$
-\left[\matrix{T_{11}&T_{12}\cr 0&T_{22}}\right]
-\left[\matrix{Z_{11}^T&Z_{21}^T\cr Z_{12}^T&Z_{22}^T}\right]
-\left[\matrix{I\cr X}\right]\left[g_{y^*}^*\right]=
-\left[\matrix{S_{11}&S_{12}\cr 0&S_{22}}\right]
-\left[\matrix{Z_{11}^T&Z_{21}^T\cr Z_{12}^T&Z_{22}^T}\right]
-\left[\matrix{I\cr X}\right]
-$$
-We reorder the eigenvalue pair so that $S_{ii}/T_{ii}$ with modulus
-less than one would be in the left-upper part.
-
-\li The Blanchard--Kahn stability argument implies that the pairs
-with modulus less that one will be in and only in $S_{11}/T_{11}$.
-The exploding paths will be then eliminated when
-$$
-\left[\matrix{Z_{11}^T&Z_{21}^T\cr Z_{12}^T&Z_{22}^T}\right]
-\left[\matrix{I\cr X}\right]=
-\left[\matrix{Y\cr 0}\right]
-$$
-From this we have, $Y=Z_{11}^{-1}$, and $X=Z_{21}Y$, or equivalently
-$X=-Z_{22}^{-T}Z_{12}^T$.  From the equation, we get
-$\left[g_{y^*}^*\right]=Y^{-1}T_{11}^{-1}S_{11}Y$, which is
-$Z_{11}T_{11}^{-1}S_{11}Z_{11}^{-1}$.
-
-\li We then copy the derivatives to storage |gy|. Note that the
-derivatives of both variables are in $X$ and in
-$\left[g_{y^*}^*\right]$, so we check whether the two submatrices are
-the same. The difference is only numerical error.
-
-\endorderedlist
-
-@<solve derivatives |gy|@>=
-	@<setup submatrices of |f|@>;
-	@<form matrix $D$@>;
-	@<form matrix $E$@>;
-	@<solve generalized Schur@>;
-	@<make submatrices of right space@>;
-	@<calculate derivatives of static and forward@>;
-	@<calculate derivatives of predetermined@>;
-	@<copy derivatives to |gy|@>;
-	@<check difference for derivatives of both@>;
-
-
-@ Here we setup submatrices of the derivatives |fd|.
-@<setup submatrices of |f|@>=
-	int off = 0;
-	ConstTwoDMatrix fyplus(fd, off, ypart.nyss());
-	off += ypart.nyss();
-	ConstTwoDMatrix fyszero(fd, off, ypart.nstat);
-	off += ypart.nstat;
-	ConstTwoDMatrix fypzero(fd, off, ypart.npred);
-	off += ypart.npred;
-	ConstTwoDMatrix fybzero(fd, off, ypart.nboth);
-	off += ypart.nboth;
-	ConstTwoDMatrix fyfzero(fd, off, ypart.nforw);
-	off += ypart.nforw;
-	ConstTwoDMatrix fymins(fd, off, ypart.nys());
-	off += ypart.nys();
-	ConstTwoDMatrix fuzero(fd, off, nu);
-	off += nu;
-
-@ 
-@<form matrix $D$@>=
-	lapack_int n = ypart.ny()+ypart.nboth;
-	TwoDMatrix matD(n, n);
-	matD.zeros();
-	matD.place(fypzero, 0, 0);
-	matD.place(fybzero, 0, ypart.npred);
-	matD.place(fyplus, 0, ypart.nys()+ypart.nstat);
-	for (int i = 0; i < ypart.nboth; i++)
-		matD.get(ypart.ny()+i, ypart.npred+i) = 1.0;
-
-@ 
-@<form matrix $E$@>=
-	TwoDMatrix matE(n, n);
-	matE.zeros();
-	matE.place(fymins, 0, 0);
-	matE.place(fyszero, 0, ypart.nys());
-	matE.place(fyfzero, 0, ypart.nys()+ypart.nstat+ypart.nboth);
-	for (int i = 0; i < ypart.nboth; i++)
-		matE.get(ypart.ny()+i, ypart.nys()+ypart.nstat+i) = -1.0;
-	matE.mult(-1.0);
-
-@ 
-@<solve generalized Schur@>=
-	TwoDMatrix vsl(n, n);
-	TwoDMatrix vsr(n, n);
-	lapack_int lwork = 100*n+16;
-	Vector work(lwork);
-	lapack_int *bwork = new lapack_int[n];
-	lapack_int info;
-	lapack_int sdim2 = sdim;
-	dgges("N", "V", "S", order_eigs, &n, matE.getData().base(), &n,
-				 matD.getData().base(), &n, &sdim2, alphar.base(), alphai.base(),
-				 beta.base(), vsl.getData().base(), &n, vsr.getData().base(), &n,
-				 work.base(), &lwork, bwork, &info);
-	sdim = sdim2;
-	bk_cond = (sdim == ypart.nys());
-	delete[] bwork;
-
-
-@ Here we setup submatrices of the matrix $Z$.
-@<make submatrices of right space@>=
-	ConstGeneralMatrix z11(vsr, 0, 0, ypart.nys(), ypart.nys());
-	ConstGeneralMatrix z12(vsr, 0, ypart.nys(), ypart.nys(), n-ypart.nys());
-	ConstGeneralMatrix z21(vsr, ypart.nys(), 0, n-ypart.nys(), ypart.nys());
-	ConstGeneralMatrix z22(vsr, ypart.nys(), ypart.nys(), n-ypart.nys(), n-ypart.nys());
-	
-@ Here we calculate $X=-Z_{22}^{-T}Z_{12}^T$, where $X$ is |sfder| in the code.
-@<calculate derivatives of static and forward@>=
-	GeneralMatrix sfder(z12, "transpose");
-	z22.multInvLeftTrans(sfder);
-	sfder.mult(-1);
-
-@ Here we calculate
-$g_{y^*}^*=Z_{11}T^{-1}_{11}S_{11}Z_{11}^{-1}
-=Z_{11}T^{-1}_{11}(Z_{11}^{-T}S^T_{11})^T$.
-
-@<calculate derivatives of predetermined@>=
-    ConstGeneralMatrix s11(matE, 0, 0, ypart.nys(), ypart.nys());
-	ConstGeneralMatrix t11(matD, 0, 0, ypart.nys(), ypart.nys());
-	GeneralMatrix dumm(s11, "transpose");
-	z11.multInvLeftTrans(dumm);
-	GeneralMatrix preder(dumm, "transpose");
-	t11.multInvLeft(preder);
-	preder.multLeft(z11);
-
-@ 
-@<copy derivatives to |gy|@>=
-	gy.place(preder, ypart.nstat, 0);
-	GeneralMatrix sder(sfder, 0, 0, ypart.nstat, ypart.nys());
-	gy.place(sder, 0, 0);
-	GeneralMatrix fder(sfder, ypart.nstat+ypart.nboth, 0, ypart.nforw, ypart.nys());
-	gy.place(fder, ypart.nstat+ypart.nys(), 0);
-
-@ 
-@<check difference for derivatives of both@>=
-	GeneralMatrix bder((const GeneralMatrix&)sfder, ypart.nstat, 0, ypart.nboth, ypart.nys());
-	GeneralMatrix bder2(preder, ypart.npred, 0, ypart.nboth, ypart.nys());
-	bder.add(-1, bder2);
-	b_error = bder.getData().getMax();
-
-@ The equation $F_u=0$ can be written as
-$$
-\left[f_{y^{**}_+}\right]\left[g^{**}_{y^*}\right]\left[g_u^*\right]+
-\left[f_y\right]\left[g_u\right]+\left[f_u\right]=0
-$$
-and rewritten as
-$$
-\left[f_y +
-\left[\matrix{0&f_{y^{**}_+}g^{**}_{y^*}&0}\right]\right]g_u=f_u
-$$
-This is exactly done here. The matrix
-$\left[f_y +\left[\matrix{0&f_{y^{**}_+}g^{**}_{y^*}&0}\right]\right]$ is |matA|
-in the code.
-
-@<solve derivatives |gu|@>=
-	GeneralMatrix matA(ypart.ny(), ypart.ny());
-	matA.zeros();
-	ConstGeneralMatrix gss(gy, ypart.nstat+ypart.npred, 0, ypart.nyss(), ypart.nys());
-	GeneralMatrix aux(fyplus, gss);
-	matA.place(aux, 0, ypart.nstat);
-	ConstGeneralMatrix fyzero(fd, 0, ypart.nyss(), ypart.ny(), ypart.ny());
-	matA.add(1.0, fyzero);
-	gu.zeros();
-	gu.add(-1.0, fuzero);
-	ConstGeneralMatrix(matA).multInvLeft(gu);
-
-@ 
-@<|FirstOrder::journalEigs| code@>=
-void FirstOrder::journalEigs()
-{
-	if (bk_cond) {
-		JournalRecord jr(journal);
-		jr << "Blanchard-Kahn conditition satisfied, model stable" << endrec;
-	} else {
-		JournalRecord jr(journal);
-		jr << "Blanchard-Kahn condition not satisfied, model not stable: sdim=" << sdim 
-		   << " " << "npred=" << ypart.nys() << endrec;
-	}
-	if (!bk_cond) {
-		for (int i = 0; i < alphar.length(); i++) {
-			if (i == sdim || i == ypart.nys()) {
-				JournalRecord jr(journal);
-				jr << "---------------------------------------------------- ";
-				if (i == sdim)
-					jr << "sdim";
-				else
-					jr << "npred";
-				jr << endrec;
-			}
-			JournalRecord jr(journal);
-			double mod = sqrt(alphar[i]*alphar[i]+alphai[i]*alphai[i]);
-			mod = mod/round(100000*std::abs(beta[i]))*100000;
-			jr << i << "\t(" << alphar[i] << "," << alphai[i] << ") / " << beta[i]
-			   << "  \t" << mod << endrec; 
-		}
-	}
-}
-
-
-@ End of {\tt first\_order.cpp} file.

dynare++/kord/first_order.hweb

-@q $Id: first_order.hweb 2345 2009-03-24 11:50:48Z kamenik $ @>
-@q Copyright 2004, Ondra Kamenik @>
-
-@*2 First order at deterministic steady. Start of {\tt first\_order.h} file.
-
-@s GeneralMatrix int
-@s ConstGeneralMatrix int
-@s FirstOrder int
-@s FirstOrderDerivs int
-@c
-
-#ifndef FIRST_ORDER_H
-#define FIRST_ORDER_H
-
-#include "korder.h"
-
-@<|FirstOrder| class declaration@>;
-@<|FirstOrderDerivs| class declaration@>;
-
-#endif
-
-@ 
-@<|FirstOrder| class declaration@>=
-template<int> class FirstOrderDerivs;
-class FirstOrder {
-	template <int> friend class FirstOrderDerivs;
-	PartitionY ypart;
-	int nu;
-	TwoDMatrix gy;
-	TwoDMatrix gu;
-	bool bk_cond;
-	double b_error;
-	int sdim;
-	Vector alphar;
-	Vector alphai;
-	Vector beta;
-	double qz_criterium;
-	Journal& journal;
-public:@;
-	FirstOrder(int num_stat, int num_pred, int num_both, int num_forw,
-			   int num_u, const FSSparseTensor& f, Journal& jr, double qz_crit)
-		: ypart(num_stat, num_pred, num_both, num_forw),
-		  nu(num_u),
-		  gy(ypart.ny(), ypart.nys()),
-		  gu(ypart.ny(), nu),
-		  alphar(ypart.ny()+ypart.nboth),
-		  alphai(ypart.ny()+ypart.nboth),
-		  beta(ypart.ny()+ypart.nboth),
-		  qz_criterium(qz_crit),
-		  journal(jr)
-		{@+ solve(FFSTensor(f)); @+}
-	bool isStable() const
-		{@+ return bk_cond;@+}
-	const TwoDMatrix& getGy() const
-		{@+ return gy;@+}
-	const TwoDMatrix& getGu() const
-		{@+ return gu;@+}
-protected:@;
-	void solve(const TwoDMatrix& f);
-	void journalEigs();
-};
-
-@ This class only converts the derivatives $g_{y^*}$ and $g_u$ to a
-folded or unfolded container.
-
-@<|FirstOrderDerivs| class declaration@>=
-template <int t>
-class FirstOrderDerivs : public ctraits<t>::Tg {
-public:@;
-	FirstOrderDerivs(const FirstOrder& fo)
-		: ctraits<t>::Tg(4)
-		{
-			IntSequence nvs(4);
-			nvs[0] = fo.ypart.nys(); nvs[1] = fo.nu; nvs[2] = fo.nu; nvs[3] = 1;
-			_Ttensor* ten = new _Ttensor(fo.ypart.ny(), TensorDimens(Symmetry(1,0,0,0),nvs));
-			ten->zeros(); ten->add(1.0, fo.gy);
-			this->insert(ten);
-			ten = new _Ttensor(fo.ypart.ny(), TensorDimens(Symmetry(0,1,0,0), nvs));
-			ten->zeros(); ten->add(1.0, fo.gu);
-			this->insert(ten);
-		}
-};
-
-
-@ End of {\tt first\_order.h} file.

dynare++/kord/global_check.cweb

-@q $Id: global_check.cweb 1830 2008-05-18 20:06:40Z kamenik $ @>
-@q Copyright 2005, Ondra Kamenik @>
-
-@ Start of {\tt global\_check.cpp} file.
-
-@c
-#include "SymSchurDecomp.h"
-
-#include "global_check.h"
-
-#include "smolyak.h"
-#include "product.h"
-#include "quasi_mcarlo.h"
-
-@<|ResidFunction| constructor code@>;
-@<|ResidFunction| copy constructor code@>;
-@<|ResidFunction| destructor code@>;
-@<|ResidFunction::setYU| code@>;
-@<|ResidFunction::eval| code@>;
-@<|GlobalChecker::check| vector code@>;
-@<|GlobalChecker::check| matrix code@>;
-@<|GlobalChecker::checkAlongShocksAndSave| code@>;
-@<|GlobalChecker::checkOnEllipseAndSave| code@>;
-@<|GlobalChecker::checkAlongSimulationAndSave| code@>;
-
-@ Here we just set a reference to the approximation, and create a new
-|DynamicModel|.
-
-@<|ResidFunction| constructor code@>=
-ResidFunction::ResidFunction(const Approximation& app)
-	: VectorFunction(app.getModel().nexog(), app.getModel().numeq()), approx(app),
-	  model(app.getModel().clone()),
-	  yplus(NULL), ystar(NULL), u(NULL), hss(NULL)
-{
-}
-
-@ 
-@<|ResidFunction| copy constructor code@>=
-ResidFunction::ResidFunction(const ResidFunction& rf)
-	: VectorFunction(rf), approx(rf.approx), model(rf.model->clone()),
-	  yplus(NULL), ystar(NULL), u(NULL), hss(NULL)
-{
-	if (rf.yplus)
-		yplus = new Vector(*(rf.yplus));
-	if (rf.ystar)
-		ystar = new Vector(*(rf.ystar));
-	if (rf.u)
-		u = new Vector(*(rf.u));
-	if (rf.hss)
-		hss = new FTensorPolynomial(*(rf.hss));
-}
-
-
-@ 
-@<|ResidFunction| destructor code@>=
-ResidFunction::~ResidFunction()
-{		
-	delete model;
-	@<delete |y| and |u| dependent data@>;
-}
-
-@ 
-@<delete |y| and |u| dependent data@>=
-	if (yplus)
-		delete yplus;
-	if (ystar)
-		delete ystar;
-	if (u)
-		delete u;
-	if (hss)
-		delete hss;
-
-
-@ This sets $y^*$ and $u$. We have to create |ystar|, |u|, |yplus| and
-|hss|.
-
-@<|ResidFunction::setYU| code@>=
-void ResidFunction::setYU(const Vector& ys, const Vector& xx)
-{
-	@<delete |y| and |u| dependent data@>;
-
-	ystar = new Vector(ys);
-	u = new Vector(xx);
-	yplus = new Vector(model->numeq());
-	approx.getFoldDecisionRule().evaluate(DecisionRule::horner,
-										  *yplus, *ystar, *u);
-
-	@<make a tensor polynomial of in-place subtensors from decision rule@>;
-	@<make |ytmp_star| be a difference of |yplus| from steady@>;
-	@<make |hss| and add steady to it@>;
-}
-
-@ Here we use a dirty tricky of converting |const| to non-|const| to
-obtain a polynomial of subtensor corresponding to non-predetermined
-variables. However, this new non-|const| polynomial will be used for a
-construction of |hss| and will be used in |const| context. So this
-dirty thing is safe.
-
-Note, that there is always a folded decision rule in |Approximation|.
-
-@<make a tensor polynomial of in-place subtensors from decision rule@>=
-	union {const FoldDecisionRule* c; FoldDecisionRule* n;} dr;
-	dr.c = &(approx.getFoldDecisionRule());
-	FTensorPolynomial dr_ss(model->nstat()+model->npred(), model->nboth()+model->nforw(),
-							*(dr.n));
-
-@ 
-@<make |ytmp_star| be a difference of |yplus| from steady@>=
-	Vector ytmp_star(ConstVector(*yplus, model->nstat(), model->npred()+model->nboth()));
-	ConstVector ysteady_star(dr.c->getSteady(), model->nstat(),
-							 model->npred()+model->nboth());
-	ytmp_star.add(-1.0, ysteady_star);
-
-@ Here is the |const| context of |dr_ss|.
-@<make |hss| and add steady to it@>=
-	hss = new FTensorPolynomial(dr_ss, ytmp_star);
-	ConstVector ysteady_ss(dr.c->getSteady(), model->nstat()+model->npred(),
-						   model->nboth()+model->nforw());
-	if (hss->check(Symmetry(0))) {
-		hss->get(Symmetry(0))->getData().add(1.0, ysteady_ss);
-	} else {
-		FFSTensor* ten = new FFSTensor(hss->nrows(), hss->nvars(), 0);
-		ten->getData() = ysteady_ss;
-		hss->insert(ten);
-	}
-
-@ Here we evaluate the residual $F(y^*,u,u')$. We have to evaluate |hss|
-for $u'=$|point| and then we evaluate the system $f$.
-
-@<|ResidFunction::eval| code@>=
-void ResidFunction::eval(const Vector& point, const ParameterSignal& sig, Vector& out)
-{
-	KORD_RAISE_IF(point.length() != hss->nvars(),
-				  "Wrong dimension of input vector in ResidFunction::eval");
-	KORD_RAISE_IF(out.length() != model->numeq(),
-				  "Wrong dimension of output vector in ResidFunction::eval");
-	Vector yss(hss->nrows());
-	hss->evalHorner(yss, point);
-	model->evaluateSystem(out, *ystar, *yplus, yss, *u);
-}
-
-@ This checks the $E[F(y^*,u,u')]$ for a given $y^*$ and $u$ by
-integrating with a given quadrature. Note that the input |ys| is $y^*$
-not whole $y$.
-
-@<|GlobalChecker::check| vector code@>=
-void GlobalChecker::check(const Quadrature& quad, int level,
-						  const ConstVector& ys, const ConstVector& x, Vector& out)
-{
-	for (int ifunc = 0; ifunc < vfs.getNum(); ifunc++)
-		((GResidFunction&)(vfs.getFunc(ifunc))).setYU(ys, x);
-	quad.integrate(vfs, level, out);
-}
-
-@ This method is a bulk version of |@<|GlobalChecker::check| vector
-code@>|. It decides between Smolyak and product quadrature according
-to |max_evals| constraint.
-
-Note that |y| can be either full (all endogenous variables including
-static and forward looking), or just $y^*$ (state variables). The
-method is able to recognize it.
-
-@<|GlobalChecker::check| matrix code@>=
-void GlobalChecker::check(int max_evals, const ConstTwoDMatrix& y,
-						  const ConstTwoDMatrix& x, TwoDMatrix& out)
-{
-	JournalRecordPair pa(journal);
-	pa << "Checking approximation error for " << y.ncols()
-	   << " states with at most " << max_evals << " evaluations" << endrec;
-
-	@<decide about type of quadrature@>;
-	Quadrature* quad;
-	int lev;
-	@<create the quadrature and report the decision@>;
-	@<check all column of |y| and |x|@>;
-	delete quad;
-}
-
-@ 
-@<decide about type of quadrature@>=
-	GaussHermite gh;
-
-	SmolyakQuadrature dummy_sq(model.nexog(), 1, gh);
-	int smol_evals;
-	int smol_level;
-	dummy_sq.designLevelForEvals(max_evals, smol_level, smol_evals);
-
-	ProductQuadrature dummy_pq(model.nexog(), gh);
-	int prod_evals;
-	int prod_level;
-	dummy_pq.designLevelForEvals(max_evals, prod_level, prod_evals);
-
-	bool take_smolyak = (smol_evals < prod_evals) && (smol_level >= prod_level-1);
-
-@ 
-@<create the quadrature and report the decision@>=
-	if (take_smolyak) {
-		quad = new SmolyakQuadrature(model.nexog(), smol_level, gh);
-		lev = smol_level;
-		JournalRecord rec(journal);
-		rec << "Selected Smolyak (level,evals)=(" << smol_level << ","
-			<< smol_evals << ") over product (" << prod_level << ","
-			<< prod_evals << ")" << endrec;
-	} else {
-		quad = new ProductQuadrature(model.nexog(), gh);
-		lev = prod_level;
-		JournalRecord rec(journal);
-		rec << "Selected product (level,evals)=(" << prod_level << ","
-			<< prod_evals << ") over Smolyak (" << smol_level << ","
-			<< smol_evals << ")" << endrec;
-	}
-
-@ 
-@<check all column of |y| and |x|@>=
-	int first_row = (y.nrows() == model.numeq())? model.nstat() : 0;
-	ConstTwoDMatrix ysmat(y, first_row, 0, model.npred()+model.nboth(), y.ncols());
-	for (int j = 0; j < y.ncols(); j++) {
-		ConstVector yj(ysmat, j);
-		ConstVector xj(x, j);
-		Vector outj(out, j);
-		check(*quad, lev, yj, xj, outj);
-	}
-
-
-
-@ This method checks an error of the approximation by evaluating
-residual $E[F(y^*,u,u')\vert y^*,u]$ for $y^*$ being the steady state, and
-changing $u$. We go through all elements of $u$ and vary them from
-$-mult\cdot\sigma$ to $mult\cdot\sigma$ in |m| steps.
-
-@<|GlobalChecker::checkAlongShocksAndSave| code@>=
-void GlobalChecker::checkAlongShocksAndSave(mat_t* fd, const char* prefix,
-											int m, double mult, int max_evals)
-{
-	JournalRecordPair pa(journal);
-	pa << "Calculating errors along shocks +/- "
-	   << mult << " std errors, granularity " << m << endrec;
-	@<setup |y_mat| of steady states for checking@>;
-	@<setup |exo_mat| for checking@>;
-
-	TwoDMatrix errors(model.numeq(), 2*m*model.nexog()+1);
-	check(max_evals, y_mat, exo_mat, errors);
-
-	@<report errors along shock and save them@>;
-}
-
-@ 
-@<setup |y_mat| of steady states for checking@>=
-	TwoDMatrix y_mat(model.numeq(), 2*m*model.nexog()+1);
-	for (int j = 0; j < 2*m*model.nexog()+1; j++) {
-		Vector yj(y_mat, j);
-		yj = (const Vector&)model.getSteady();
-	}
-
-@ 
-@<setup |exo_mat| for checking@>=
-	TwoDMatrix exo_mat(model.nexog(), 2*m*model.nexog()+1);
-	exo_mat.zeros();
-	for (int ishock = 0; ishock < model.nexog(); ishock++) {
-		double max_sigma = sqrt(model.getVcov().get(ishock,ishock));
-		for (int j = 0; j < 2*m; j++) {
-			int jmult = (j < m)? j-m: j-m+1;
-			exo_mat.get(ishock, 1+2*m*ishock+j) = 
-				mult*jmult*max_sigma/m;
-		}
-	}
-
-@ 
-@<report errors along shock and save them@>=
-	TwoDMatrix res(model.nexog(), 2*m+1);
-	JournalRecord rec(journal);
-	rec << "Shock    value         error" << endrec;
-	ConstVector err0(errors,0);
-	char shock[9];
-	char erbuf[17];
-	for (int ishock = 0; ishock < model.nexog(); ishock++) {
-		TwoDMatrix err_out(model.numeq(), 2*m+1);
-		sprintf(shock, "%-8s", model.getExogNames().getName(ishock));
-		for (int j = 0; j < 2*m+1; j++) {
-			int jj;
-			Vector error(err_out, j);
-			if (j != m) {
-				if (j < m)
-					jj = 1 + 2*m*ishock+j;
-				else
-					jj = 1 + 2*m*ishock+j-1;
-				ConstVector coljj(errors,jj);
-				error = coljj;
-			} else {
-				jj = 0;
-				error = err0;
-			}
-			JournalRecord rec1(journal);
-			sprintf(erbuf,"%12.7g    ", error.getMax());
-			rec1 << shock << " " << exo_mat.get(ishock,jj)
-				<< "\t" << erbuf << endrec;
-		}
-		char tmp[100];
-		sprintf(tmp, "%s_shock_%s_errors", prefix, model.getExogNames().getName(ishock));
-		err_out.writeMat(fd, tmp);
-	}
-
-
-@ This method checks errors on ellipse of endogenous states
-(predetermined variables). The ellipse is shaped according to
-covariance matrix of endogenous variables based on the first order
-approximation and scaled by |mult|. The points on the
-ellipse are chosen as polar images of the low discrepancy grid in a
-cube.
-
-The method works as follows: First we calculate symmetric Schur factor of
-covariance matrix of the states. Second we generate low discrepancy
-points on the unit sphere. Third we transform the sphere with the
-variance-covariance matrix factor and multiplier |mult| and initialize
-matrix of $u_t$ to zeros. Fourth we run the |check| method and save
-the results.
-
-@<|GlobalChecker::checkOnEllipseAndSave| code@>=
-void GlobalChecker::checkOnEllipseAndSave(mat_t* fd, const char* prefix,
-										  int m, double mult, int max_evals)
-{
-	JournalRecordPair pa(journal);
-	pa << "Calculating errors at " << m
-	   << " ellipse points scaled by " << mult << endrec; 
-	@<make factor of covariance of variables@>;
-	@<put low discrepancy sphere points to |ymat|@>;
-	@<transform sphere |ymat| and prepare |umat| for checking@>;
-	@<check on ellipse and save@>;
-}
-
-
-@ Here we set |ycovfac| to the symmetric Schur decomposition factor of
-a submatrix of covariances of all endogenous variables. The submatrix
-corresponds to state variables (predetermined plus both).
-
-@<make factor of covariance of variables@>=
-	TwoDMatrix* ycov = approx.calcYCov();
-	TwoDMatrix ycovpred((const TwoDMatrix&)*ycov, model.nstat(), model.nstat(),
-						model.npred()+model.nboth(), model.npred()+model.nboth());
-	delete ycov;
-	SymSchurDecomp ssd(ycovpred);
-	ssd.correctDefinitness(1.e-05);
-	TwoDMatrix ycovfac(ycovpred.nrows(), ycovpred.ncols());
-	KORD_RAISE_IF(! ssd.isPositiveSemidefinite(),
-				  "Covariance matrix of the states not positive \
-				  semidefinite in GlobalChecker::checkOnEllipseAndSave");
-	ssd.getFactor(ycovfac);
-
-
-@ Here we first calculate dimension |d| of the sphere, which is a
-number of state variables minus one. We go through the |d|-dimensional
-cube $\langle 0,1\rangle^d$ by |QMCarloCubeQuadrature| and make a
-polar transformation to the sphere. The polar transformation $f^i$ can
-be written recursively wrt. the dimension $i$ as:
-$$\eqalign{
-f^0() &= \left[1\right]\cr
-f^i(x_1,\ldots,x_i) &=
-\left[\matrix{cos(2\pi x_i)\cdot f^{i-1}(x_1,\ldots,x_{i-1})\cr sin(2\pi x_i)}\right]
-}$$
-
-@<put low discrepancy sphere points to |ymat|@>=
-   	int d = model.npred()+model.nboth()-1;
-	TwoDMatrix ymat(model.npred()+model.nboth(), (d==0)? 2:m);
-	if (d == 0) {
-		ymat.get(0,0) = 1;
-		ymat.get(0,1) = -1;
-	} else {
-		int icol = 0;
-		ReversePerScheme ps;
-		QMCarloCubeQuadrature qmc(d, m, ps);
-		qmcpit beg = qmc.start(m);
-		qmcpit end = qmc.end(m);
-		for (qmcpit run = beg; run != end; ++run, icol++) {
-			Vector ycol(ymat, icol);
-			Vector x(run.point());
-			x.mult(2*M_PI);
-			ycol[0] = 1;
-			for (int i = 0; i < d; i++) {
-				Vector subsphere(ycol, 0, i+1);
-				subsphere.mult(cos(x[i]));
-				ycol[i+1] = sin(x[i]);
-			}
-		}
-	}
-
-@ Here we multiply the sphere points in |ymat| with the Cholesky
-factor to obtain the ellipse, scale the ellipse by the given |mult|,
-and initialize matrix of shocks |umat| to zero.
-
-@<transform sphere |ymat| and prepare |umat| for checking@>=
-	TwoDMatrix umat(model.nexog(), ymat.ncols());
-	umat.zeros();
-	ymat.mult(mult);
-	ymat.multLeft(ycovfac);
-	ConstVector ys(model.getSteady(), model.nstat(),
-				   model.npred()+model.nboth());
-	for (int icol = 0; icol < ymat.ncols(); icol++) {
-		Vector ycol(ymat, icol);
-		ycol.add(1.0, ys);
-	}
-
-@ Here we check the points and save the results to MAT-4 file.
-@<check on ellipse and save@>=
-	TwoDMatrix out(model.numeq(), ymat.ncols());
-	check(max_evals, ymat, umat, out);
-
-	char tmp[100];
-	sprintf(tmp, "%s_ellipse_points", prefix);
-	ymat.writeMat(fd, tmp);
-	sprintf(tmp, "%s_ellipse_errors", prefix);
-	out.writeMat(fd, tmp);
-
-@ Here we check the errors along a simulation. We simulate, then set
-|x| to zeros, check and save results.
-
-@<|GlobalChecker::checkAlongSimulationAndSave| code@>=
-void GlobalChecker::checkAlongSimulationAndSave(mat_t* fd, const char* prefix,
-												int m, int max_evals)
-{
-	JournalRecordPair pa(journal);
-	pa << "Calculating errors at " << m
-	   << " simulated points" << endrec; 
-	RandomShockRealization sr(model.getVcov(), system_random_generator.int_uniform());
-	TwoDMatrix* y = approx.getFoldDecisionRule().simulate(DecisionRule::horner,
-														  m, model.getSteady(), sr);
-	TwoDMatrix x(model.nexog(), m);
-	x.zeros();
-	TwoDMatrix out(model.numeq(), m);
-	check(max_evals, *y, x, out);
-
-	char tmp[100];
-	sprintf(tmp, "%s_simul_points", prefix);
-	y->writeMat(fd, tmp);
-	sprintf(tmp, "%s_simul_errors", prefix);
-	out.writeMat(fd, tmp);
-
-	delete y;
-}
-
-
-@ End of {\tt global\_check.cpp} file.

dynare++/kord/global_check.hweb

-@q $Id: global_check.hweb 431 2005-08-16 15:41:01Z kamenik $ @>
-@q Copyright 2005, Ondra Kamenik @>
-
-@*2 Global check. Start of {\tt global\_check.h} file.
-
-The purpose of this file is to provide classes for checking error of
-approximation. If $y_t=g(y^*_{t-1},u)$ is an approximate solution,
-then we check for the error of residuals of the system equations. Let
-$F(y^*,u,u')=f(g^{**}(g^*(y^*,u'),u),g(y^*,u),y^*,u)$, then we
-calculate integral
-$$E[F(y^*,u,u')]$$@q'@>
-which we want to be zero for all $y^*$, and $u$.
-
-There are a few possibilities how and where the integral is
-evaluated. Currently we offer the following:
-
-\numberedlist
-\li Along shocks. The $y^*$ is set to steady state, and $u$ is set to
-zero but one element is going from minus through plus shocks in few
-steps. The user gives the scaling factor, for instance interval
-$\langle<-3\sigma,3\sigma\rangle$ (where $sigma$ is a standard error
-of the shock), and a number of steps. This is repeated for each shock
-(element of the $u$ vector).
-\li Along simulation. Some random simulation is run, and for each
-realization of $y^*$ and $u$ along the path we evaluate the residual.
-\li On ellipse. Let $V=AA^T$ be a covariance matrix of the
-predetermined variables $y^*$ based on linear approximation, then we
-calculate integral for points on the ellipse $\{Ax\vert\, \Vert
-x\Vert_2=1\}$. The points are selected by means of low discrepancy
-method and polar transformation. The shock $u$ are zeros.
-
-\li Unconditional distribution. 
-
-\endnumberedlist
-
-
-@s ResidFunction int
-@s GResidFunction int
-@s GlobalChecker int
-@s VectorFunction int
-@s ResidFunctionSig int
-@s GaussHermite int
-@s SmolyakQuadrature int
-@s ProductQuadrature int
-@s ParameterSignal int
-@s Quadrature int
-@s QMCarloCubeQuadrature int
-
-@c
-#ifndef GLOBAL_CHECK_H
-#define GLOBAL_CHECK_H
-
-#include <matio.h>
-
-#include "vector_function.h"
-#include "quadrature.h"
-
-#include "dynamic_model.h"
-#include "journal.h"
-#include "approximation.h"
-
-@<|ResidFunction| class declaration@>;
-@<|GResidFunction| class declaration@>;
-@<|GlobalChecker| class declaration@>;
-@<|ResidFunctionSig| class declaration@>;
-
-#endif
-
-@ This is a class for implementing |VectorFunction| interface
-evaluating the residual of equations, this is
-$$F(y^*,u,u')=f(g^{**}(g^*(y^*,u),u'),y^*,u)$$
-is written as a function of $u'$.
-
-When the object is constructed, one has to specify $(y^*,u)$, this is
-done by |setYU| method. The object has basically two states. One is
-after construction and before call to |setYU|. The second is after
-call |setYU|. We distinguish between the two states, an object in the
-second state contains |yplus|, |ystar|, |u|, and |hss|.
-
-The vector |yplus| is $g^*(y^*,u)$. |ystar| is $y^*$, and polynomial
-|hss| is partially evaluated $g^**(yplus, u)$.
-
-The pointer to |DynamicModel| is important, since the |DynamicModel|
-evaluates the function $f$. When copying the object, we have to make
-also a copy of |DynamicModel|.
-
-@<|ResidFunction| class declaration@>=
-class ResidFunction : public VectorFunction {
-protected:@;
-	const Approximation& approx;
-	DynamicModel* model;
-	Vector* yplus;
-	Vector* ystar;
-	Vector* u;
-	FTensorPolynomial* hss;
-public:@;
-	ResidFunction(const Approximation& app);
-	ResidFunction(const ResidFunction& rf);
-	virtual ~ResidFunction();
-	virtual VectorFunction* clone() const
-		{@+ return new ResidFunction(*this);@+}
-	virtual void eval(const Vector& point, const ParameterSignal& sig, Vector& out);
-	void setYU(const Vector& ys, const Vector& xx);
-};
-
-@ This is a |ResidFunction| wrapped with |GaussConverterFunction|.
-
-@<|GResidFunction| class declaration@>=
-class GResidFunction : public GaussConverterFunction {
-public:@;
-	GResidFunction(const Approximation& app)
-		: GaussConverterFunction(new ResidFunction(app), app.getModel().getVcov())@+ {}
-	GResidFunction(const GResidFunction& rf)
-		: GaussConverterFunction(rf)@+ {}
-	virtual ~GResidFunction()@+ {}
-	virtual VectorFunction* clone() const
-		{@+ return new GResidFunction(*this);@+}
-	void setYU(const Vector& ys, const Vector& xx)
-		{@+ ((ResidFunction*)func)->setYU(ys, xx);}
-};
-
-
-@ This is a class encapsulating checking algorithms. Its core routine
-is |check|, which calculates integral $E[F(y^*,u,u')\vert y^*,u]$ for
-given realizations of $y^*$ and $u$. The both are given in
-matrices. The methods checking along shocks, on ellipse and anlong a
-simulation path, just fill the matrices and call the core |check|.
-
-The method |checkUnconditionalAndSave| evaluates unconditional
-$E[F(y,u,u')]$.
-
-The object also maintains a set of |GResidFunction| functions |vfs| in
-order to save (possibly expensive) copying of |DynamicModel|s.
-
-@<|GlobalChecker| class declaration@>=
-class GlobalChecker {
-	const Approximation& approx;
-	const DynamicModel& model;
-	Journal& journal;
-	GResidFunction rf;
-	VectorFunctionSet vfs;
-public:@;
-	GlobalChecker(const Approximation& app, int n, Journal& jr)
-		: approx(app), model(approx.getModel()), journal(jr),
-		  rf(approx), vfs(rf, n)@+ {}
-	void check(int max_evals, const ConstTwoDMatrix& y,
-			   const ConstTwoDMatrix& x, TwoDMatrix& out);
-	void checkAlongShocksAndSave(mat_t* fd, const char* prefix,
-								 int m, double mult, int max_evals);
-	void checkOnEllipseAndSave(mat_t* fd, const char* prefix,
-							   int m, double mult, int max_evals);
-	void checkAlongSimulationAndSave(mat_t* fd, const char* prefix,
-									 int m, int max_evals);
-	void checkUnconditionalAndSave(mat_t* fd, const char* prefix,
-								   int m, int max_evals);
-protected:@;
-	void check(const Quadrature& quad, int level,
-			   const ConstVector& y, const ConstVector& x, Vector& out);
-};
-
-
-@ Signalled resid function. Not implemented yet. todo:
-@<|ResidFunctionSig| class declaration@>=
-class ResidFunctionSig : public ResidFunction {
-public:@;
-	ResidFunctionSig(const Approximation& app, const Vector& ys, const Vector& xx);
-};
-
-@ End of {\tt global\_check.h} file.