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

Amend Estimation C++ ModelSolution, some utils (new Matrix version of dynamic_dll) and a test routine

time-shift
unknown 2010-02-18 10:34:59 +00:00 committed by George Peredia
parent c1e7ade63e
commit 9a77efe440
6 changed files with 618 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

109
mex/sources/estimation/ModelSolution.cc Normal file
View File

@ -0,0 +1,109 @@
/*
* Copyright (C) 2010 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
///////////////////////////////////////////////////////////
// ModelSolution.cpp
// Implementation of the Class ModelSolution
// Created on: 02-Feb-2010 13:06:35
///////////////////////////////////////////////////////////
#include "ModelSolution.hh"
/**
* compute the steady state (2nd stage), and computes first order approximation
*/
ModelSolution::ModelSolution(const std::string& modName, size_t n_endo_arg, size_t n_exo_arg, const std::vector<size_t>& zeta_fwrd_arg,
const std::vector<size_t>& zeta_back_arg, const std::vector<size_t>& zeta_mixed_arg,
const std::vector<size_t>& zeta_static_arg, const Matrix& llincidence, double INqz_criterium)
: n_endo(n_endo_arg), n_exo(n_exo_arg), // n_jcols = Num of Jacobian columns = nStat+2*nPred+3*nBoth+2*nForw+nExog
n_jcols (n_exo+n_endo+ zeta_back_arg.size() /*nsPred*/ + zeta_fwrd_arg.size() /*nsForw*/ +2*zeta_mixed_arg.size()),
ll_incidence(llincidence), jacobian (n_endo,n_jcols), residual(n_endo), Mx(2,n_exo),
decisionRules ( n_endo_arg, n_exo_arg, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, INqz_criterium)
{
std::string sExt(""); // use a pre-constructed model_dynamic.ext file name
dynamicDLLp = new DynamicModelDLL(modName, n_endo, n_jcols, /* nMax_lag= */ 1, n_exo, sExt) ;
}
ModelSolution::~ModelSolution()
{
delete dynamicDLLp;
}
void
ModelSolution::compute(Vector& steadyState, const Vector& deepParams, Matrix& ghx, Matrix& ghu)
{
// compute Steady State
ComputeSteadyState(steadyState, deepParams);
// then get jacobian and
ComputeModelSolution( steadyState, deepParams, ghx, ghu);
}
void
ModelSolution::ComputeModelSolution(Vector& steadyState, const Vector& deepParams, Matrix& ghx, Matrix& ghu)
{
// set extended Steady State
Vector llXsteadyState(n_jcols-n_exo);
try
{
for (int ll_row = 0; ll_row < ll_incidence.getRows(); ll_row++)
{
// populate (non-sparse) vector with ysteady values
for (int i = 0; i < n_endo; i++)
{
if (ll_incidence(ll_row, i))
llXsteadyState(((int) ll_incidence(ll_row, i))-1) = steadyState(i);
}
}
#ifdef DEBUG
// std::cout << "Vector llXsteadyState: " << std::endl << llXsteadyState << std::endl;
mexPrintf(" get jacobian \n");
#endif
//get jacobian
dynamicDLLp->eval(llXsteadyState, Mx, &deepParams, 1, residual, &jacobian, NULL, NULL);
#ifdef DEBUG
std::cout << "jacobian: " << std::endl << jacobian << std::endl;
mexPrintf(" compute rules \n");
#endif
//compute rules
decisionRules.compute(jacobian,ghx, ghu);
}
catch (const ModelSolutionException &e)
{
mexPrintf("Caught ModelSolution exception in LLxSteady: ");
e.print();
}
catch (...)
{
mexPrintf("Caught unknown error in ModelSolution::ComputeModelSolution\n");
}
}
void
ModelSolution::ComputeSteadyState(Vector& steadyState, const Vector& deepParams)
{
// does nothig for time being.
}

74
mex/sources/estimation/ModelSolution.hh Normal file
View File

@ -0,0 +1,74 @@
/*
* Copyright (C) 2010 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
///////////////////////////////////////////////////////////
// ComputeModelSolution.h
// Implementation of the Class ModelSolution
// Created on: 15-Jan-2010 07:37:47
///////////////////////////////////////////////////////////
#if !defined(ModelSolution_5ADFF920_9C74_46f5_9FE9_88AD4D4BBF19__INCLUDED_)
#define ModelSolution_5ADFF920_9C74_46f5_9FE9_88AD4D4BBF19__INCLUDED_
#include "DecisionRules.hh"
#include "dynamic_dll.hh"
/**
* compute the steady state (2nd stage), and
* computes first order approximation
*
*/
class ModelSolution{
public:
class ModelSolutionException: public TSException
{
public:
const lapack_int info, n;
ModelSolutionException(lapack_int info_arg, lapack_int n_arg, std::string& mes)
: TSException(__FILE__, __LINE__, mes), info(info_arg), n(n_arg) {};
};
virtual ~ModelSolution();
ModelSolution(const std::string& modName, size_t n_endo, size_t n_exo, const std::vector<size_t>& zeta_fwrd_arg,
const std::vector<size_t>& zeta_back_arg, const std::vector<size_t>& zeta_mixed_arg,
const std::vector<size_t>& zeta_static_arg, const Matrix& llincidence, double qz_criterium);
void compute( Vector& steadyState, const Vector& deepParams, Matrix& ghx, Matrix& ghu );
private:
const int n_endo;
const int n_exo;
const int n_jcols; // Num of Jacobian columns
const Matrix ll_incidence; // leads and lags indices
Matrix jacobian;
Vector residual;
Matrix Mx;
DecisionRules decisionRules;
DynamicModelDLL* dynamicDLLp;
//Matrix jacobian;
void ComputeModelSolution( Vector& steadyState, const Vector& deepParams, Matrix& ghx, Matrix& ghu );
void ComputeSteadyState( Vector& steadyState, const Vector& deepParams);
};
#endif // !defined(5ADFF920_9C74_46f5_9FE9_88AD4D4BBF19__INCLUDED_)

128
mex/sources/estimation/tests/testModelSolution.cc Normal file
View File

@ -0,0 +1,128 @@
/*
* Copyright (C) 2010 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
// Test for ModelSolution
// Uses fs2000k2.mod and its ..._dynamic.mexw32
#include "ModelSolution.hh"
int
main (int argc, char** argv)
{
std::string modName("fs2000k2_dynamic.mexw32");
const int npar = 7; //(int)mxGetM(mxFldp);
const size_t n_endo=15, n_exo=2;
std::vector<size_t> zeta_fwrd_arg;
std::vector<size_t> zeta_back_arg;
std::vector<size_t> zeta_mixed_arg;
std::vector<size_t> zeta_static_arg;
//std::vector<size_t>
double qz_criterium=1.0e-9;
Vector steadyState(n_endo), deepParams(npar);
double dYSparams [] = { // 27 mxGetData(mxFldp);
1.0110, 2.2582, 0.4477, 1.0000,
4.5959, 1.0212, 5.8012, 0.8494,
0.1872, 0.8604, 1.0030, 1.0080,
0.5808, 1.0030, 2.2093 //2.2582, 0.4477
};
double vcov[] = { //(double *) mxGetData(mxFldp);
0.1960e-3, 0.0,
0.0, 0.0250e-3
};
int nVCVpar = 2; //(int)mxGetN(mxFldp);
MatrixView vCovVW(vcov,nVCVpar,nVCVpar,nVCVpar);
Matrix vCov (nVCVpar, nVCVpar);
vCov = vCovVW;
Matrix ll_incidence(3,n_endo); // leads and lags indices
double inllincidence[]={
// 1, 2, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 4, 0, 0,
// 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
// 0, 20, 21, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 22
1, 5, 0,
2, 6, 20,
0, 7, 21,
0, 8, 0,
0, 9, 0,
0, 10, 0,
3, 11, 0,
0, 12, 0,
0, 13, 0,
0, 14, 0,
0, 15, 0,
0, 16, 0,
4, 17, 0,
0, 18, 0,
0, 19, 22,
};
MatrixView llincidence(inllincidence, 3, n_endo,3); // leads and lags indices
ll_incidence= llincidence;
double dparams[] = { 0.3300,
0.9900,
0.0030,
1.0110,
0.7000,
0.7870,
0.0200};
VectorView modParamsVW (dparams, npar,1);
deepParams=modParamsVW;
VectorView steadyStateVW(dYSparams,n_endo,1);
steadyState=steadyStateVW;
//std::cout << "VectorView deepParamsVW: " << std::endl << modParamsVW << std::endl;
std::cout << "Vector deepParams: " << std::endl << deepParams << std::endl;
std::cout << "Matrix vCov: " << std::endl << vCov << std::endl;
std::cout << "MatrixVw llincidence: " << std::endl << llincidence << std::endl;
std::cout << "Matrix ll_incidence: " << std::endl << ll_incidence << std::endl;
std::cout << "Vector steadyState: " << std::endl << steadyState << std::endl;
//order_var = [ stat_var(:); pred_var(:); both_var(:); fwrd_var(:)];
size_t statc[]={ 4, 5, 6, 8, 9, 10, 11, 12, 14};
size_t back[]={1,7,13};
size_t both[]={2};
size_t fwd[]={ 3,15};
for (int i=0;i<9;++i)
zeta_static_arg.push_back(statc[i]);
for (int i=0;i<3;++i)
zeta_back_arg.push_back(back[i]);
for (int i=0;i<1;++i)
zeta_mixed_arg.push_back(both[i]);
for (int i=0;i<2;++i)
zeta_fwrd_arg.push_back(fwd[i]);
Matrix ghx(n_endo, zeta_back_arg.size() + zeta_mixed_arg.size());
Matrix ghu(n_endo,n_exo);
// exit(0);
ModelSolution modelSolution( modName, n_endo, n_exo
, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, ll_incidence, qz_criterium);
// exit(0);
modelSolution.compute(steadyState, deepParams, ghx, ghu);
std::cout << "Matrix ghx: " << std::endl << ghx << std::endl;
std::cout << "Matrix ghu: " << std::endl << ghu << std::endl;
}

159
mex/sources/estimation/utils/dynamic_dll.cc Normal file
View File

@ -0,0 +1,159 @@
/*
* Copyright (C) 2010 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
//#include "k_ord_dynare.hh"
#include "dynamic_dll.hh"
#include <sstream>
/***********************************
* Members of DynamicModelDLL for handling loading and calling
* <model>_dynamic () function
**************************************/
DynamicModelDLL::DynamicModelDLL(const std::string &modName, const int y_length, const int j_cols,
const int n_max_lag, const int n_exog, const std::string &sExt) throw (TSException) :
length(y_length), jcols(j_cols), nMax_lag(n_max_lag), nExog(n_exog)
{
std::string fName;
#if !defined(__CYGWIN32__) && !defined(_WIN32)
fName = "./";
#endif
if (sExt.size()>0) //construct modelNmae_dynamic file name with the given extension
fName += modName + "_dynamic" + sExt;
else // i.e. an already pre-constructed modelNmae_dynamic file name together with the appropriate extension
fName += modName;
// end if
try
{
#if defined(__CYGWIN32__) || defined(_WIN32)
dynamicHinstance = LoadLibrary(fName.c_str());
if (dynamicHinstance == NULL)
throw 1;
Dynamic = (DynamicFn) GetProcAddress(dynamicHinstance, "Dynamic");
if (Dynamic == NULL)
{
FreeLibrary(dynamicHinstance); // Free the library
throw 2;
}
#else // Linux or Mac
dynamicHinstance = dlopen(fName.c_str(), RTLD_NOW);
if ((dynamicHinstance == NULL) || dlerror())
{
cerr << dlerror() << endl;
throw 1;
}
Dynamic = (DynamicFn) dlsym(dynamicHinstance, "Dynamic");
if ((Dynamic == NULL) || dlerror())
{
dlclose(dynamicHinstance); // Free the library
cerr << dlerror() << endl;
throw 2;
}
#endif
}
catch (int i)
{
std::ostringstream msg;
msg << "Error when loading " << fName << " (";
if (i == 1)
msg << "can't dynamically load the file";
if (i == 2)
msg << "can't locate the 'Dynamic' symbol";
msg << ")";
throw TSException(__FILE__, __LINE__, msg.str());
}
catch (...)
{
throw TSException(__FILE__, __LINE__, std::string("Can't find Dynamic function in ") + fName);
}
}
DynamicModelDLL::~DynamicModelDLL()
{
#if defined(__CYGWIN32__) || defined(_WIN32)
bool result = FreeLibrary(dynamicHinstance);
if (result == 0)
throw TSException(__FILE__, __LINE__, std::string("Can't free the *_dynamic DLL"));
#else
dlclose(dynamicHinstance);
#endif
}
void
DynamicModelDLL::eval(double *y, double *x, int nb_row_x, double *params,
int it_, double *residual, double *g1, double *g2, double *g3)
{
Dynamic(y, x, nb_row_x, params, it_, residual, g1, g2, g3);
}
void
DynamicModelDLL::eval(const Vector &y, const Matrix &x, const Vector *modParams,
int it_, Vector &residual, Matrix *g1, Matrix *g2, Matrix *g3) throw (TSException)
{
double *dresidual, *dg1 = NULL, *dg2 = NULL, *dg3 = NULL;
if ((jcols-nExog) != y.getSize())
throw TSException(__FILE__, __LINE__, "DLL Error: (jcols-nExog)!=ys.length()");
if (g1 != NULL)
{
if (g1->getRows() != length) // dummy
{
delete g1;
g1 = new Matrix(length, jcols); // and get a new one
g1->setAll(0.0); //zero
}
dg1 = const_cast<double *>(g1->getData());
}
if (g2 != NULL)
dg2 = const_cast<double *>(g2->getData());
//dresidual = const_cast<double *>(residual.getData());
if (g3 != NULL)
dg3 = const_cast<double *>(g3->getData());
dresidual = const_cast<double *>(residual.getData());
double *dy = const_cast<double *>(y.getData());
double *dx = const_cast<double *>(x.getData());
double *dbParams = const_cast<double *>(modParams->getData());
Dynamic(dy, dx, nExog, dbParams, it_, dresidual, dg1, dg2, dg3);
}
void
DynamicModelDLL::eval(const Vector &y, const Matrix &x, const Vector *modParams,
Vector &residual, Matrix *g1, Matrix *g2, Matrix *g3) throw (TSException)
{
eval(y, x, modParams, nMax_lag, residual, g1, g2, g3);
}
void
DynamicModelDLL::eval(const Vector &y, const Vector &x, const Vector *modParams,
Vector &residual, Matrix *g1, Matrix *g2, Matrix *g3) throw (TSException)
{
/** ignore given exogens and create new 2D x matrix since
* when calling <model>_dynamic(z,x,params,it_) x must be equal to
* zeros(M_.maximum_lag+1,M_.exo_nbr)
**/
Matrix &mx = *(new Matrix(nMax_lag+1, nExog));
mx.setAll(0); // initialise shocks to 0s
eval(y, mx, modParams, nMax_lag, residual, g1, g2, g3);
}

81
mex/sources/estimation/utils/dynamic_dll.hh Normal file
View File

@ -0,0 +1,81 @@
/*
* Copyright (C) 2010 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
// The following ifdef block is the standard way of creating macros which make exporting
// from a DLL simpler. All files within this DLL are compiled with the K_ORDER_PERTURBATION_EXPORTS
// symbol defined on the command line. this symbol should not be defined on any project
// that uses this DLL. This way any other project whose source files include this file see
// K_ORDER_PERTURBATION_API functions as being imported from a DLL, wheras this DLL sees symbols
// defined with this macro as being exported.
#if defined(_WIN32) || defined(__CYGWIN32__)
# include <windows.h>
# ifdef _MSC_VER
# define K_ORDER_PERTURBATION_API __declspec(dllexport)
# endif
#else
# include <dlfcn.h> // unix/linux DLL (.so) handling routines
#endif
#include <string>
#include <dynmex.h>
#include "Matrix.hh"
#include "ts_exception.h"
// <model>_Dynamic DLL pointer
typedef void (*DynamicFn)
(double *y, double *x, int nb_row_x, double *params,
int it_, double *residual, double *g1, double *g2, double *g3);
/**
* creates pointer to Dynamic function inside <model>_dynamic.dll
* and handles calls to it.
**/
class DynamicModelDLL
{
private:
DynamicFn Dynamic; // pointer to the Dynamic function in DLL
const int length; // tot num vars = Num of Jacobian rows
const int jcols; // tot num var t-1, t and t+1 instances + exogs = Num of Jacobian columns
const int nMax_lag; // no of lags
const int nExog; // no of exogenous
#if defined(_WIN32) || defined(__CYGWIN32__)
HINSTANCE dynamicHinstance; // DLL instance pointer in Windows
#else
void *dynamicHinstance; // and in Linux or Mac
#endif
public:
// construct and load Dynamic model DLL
DynamicModelDLL(const std::string &modName, const int length, const int jcols,
const int nMax_lag, const int nExog, const std::string &sExt) throw (TSException);
virtual ~DynamicModelDLL();
// evaluate Dynamic model DLL
void eval(double *y, double *x, int nb_row_x, double *params,
int it_, double *residual, double *g1, double *g2, double *g3);
void eval(const Vector &y, const Vector &x, const Vector *params,
Vector &residual, Matrix *g1, Matrix *g2, Matrix *g3) throw (TSException);
void eval(const Vector &y, const Matrix &x, const Vector *params,
int it_, Vector &residual, Matrix *g1, Matrix *g2, Matrix *g3) throw (TSException);
void eval(const Vector &y, const Matrix &x, const Vector *params,
Vector &residual, Matrix *g1, Matrix *g2, Matrix *g3) throw (TSException);
};

67
mex/sources/estimation/utils/ts_exception.h Normal file
View File

@ -0,0 +1,67 @@
/*
* Copyright (C) 2010 Dynare Team
*
* This file is part of Dynare.
*
* Dynare is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Dynare is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
/* derived from c++kalman_filter library by O. Kamenik */
#ifndef TS_EXCEPTION_H
#define TS_EXCEPTION_H
#include <stdio.h>
#include <string>
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
#include "mex.h"
#endif
#define TS_RAISE(mes) \
throw TSException(__FILE__, __LINE__, mes);
#define TS_RAISE_IF(expr, mes) \
if (expr) throw TSException(__FILE__, __LINE__, mes);
class TSException
{
std::string fname; //char fname[50];
int lnum;
std::string message;// char message[500];
public:
TSException(const char*f,int l,const std::string &mes)
{
fname =std::string(f); // strncpy(fname,f,50);fname[49]= '\0';
message=mes;//strncpy(message,mes,500);message[499]= '\0';
lnum= l;
}
virtual ~TSException(){};
virtual void print()const
{
printf("At %s:%d:%s\n",fname.c_str(),lnum,message.c_str());
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
mexPrintf("At %s:%d:%s\n",fname.c_str(),lnum,message.c_str());
#endif
}
virtual const std::string getMessage()const
{return message;}
};
;
#endif