diff --git a/mex/sources/estimation/ModelSolution.cc b/mex/sources/estimation/ModelSolution.cc new file mode 100644 index 000000000..8207dc479 --- /dev/null +++ b/mex/sources/estimation/ModelSolution.cc @@ -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 . + */ + +/////////////////////////////////////////////////////////// +// 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& zeta_fwrd_arg, + const std::vector& zeta_back_arg, const std::vector& zeta_mixed_arg, + const std::vector& 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. +} + diff --git a/mex/sources/estimation/ModelSolution.hh b/mex/sources/estimation/ModelSolution.hh new file mode 100644 index 000000000..51e9a1cb9 --- /dev/null +++ b/mex/sources/estimation/ModelSolution.hh @@ -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 . + */ + +/////////////////////////////////////////////////////////// +// 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& zeta_fwrd_arg, + const std::vector& zeta_back_arg, const std::vector& zeta_mixed_arg, + const std::vector& 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_) diff --git a/mex/sources/estimation/tests/testModelSolution.cc b/mex/sources/estimation/tests/testModelSolution.cc new file mode 100644 index 000000000..8222361ab --- /dev/null +++ b/mex/sources/estimation/tests/testModelSolution.cc @@ -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 . + */ + +// 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 zeta_fwrd_arg; + std::vector zeta_back_arg; + std::vector zeta_mixed_arg; + std::vector zeta_static_arg; + //std::vector + 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; + + +} + diff --git a/mex/sources/estimation/utils/dynamic_dll.cc b/mex/sources/estimation/utils/dynamic_dll.cc new file mode 100644 index 000000000..86b337c09 --- /dev/null +++ b/mex/sources/estimation/utils/dynamic_dll.cc @@ -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 . + */ + +//#include "k_ord_dynare.hh" +#include "dynamic_dll.hh" + +#include + + +/*********************************** +* Members of DynamicModelDLL for handling loading and calling +* _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(g1->getData()); + } + if (g2 != NULL) + dg2 = const_cast(g2->getData()); + //dresidual = const_cast(residual.getData()); + if (g3 != NULL) + dg3 = const_cast(g3->getData()); + dresidual = const_cast(residual.getData()); + double *dy = const_cast(y.getData()); + double *dx = const_cast(x.getData()); + double *dbParams = const_cast(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 _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); +} diff --git a/mex/sources/estimation/utils/dynamic_dll.hh b/mex/sources/estimation/utils/dynamic_dll.hh new file mode 100644 index 000000000..a9df5871f --- /dev/null +++ b/mex/sources/estimation/utils/dynamic_dll.hh @@ -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 . + */ + +// 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 +# ifdef _MSC_VER +# define K_ORDER_PERTURBATION_API __declspec(dllexport) +# endif +#else +# include // unix/linux DLL (.so) handling routines +#endif + +#include +#include +#include "Matrix.hh" + +#include "ts_exception.h" + +// _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 _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); +}; diff --git a/mex/sources/estimation/utils/ts_exception.h b/mex/sources/estimation/utils/ts_exception.h new file mode 100644 index 000000000..381a4c55e --- /dev/null +++ b/mex/sources/estimation/utils/ts_exception.h @@ -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 . + */ + +/* derived from c++kalman_filter library by O. Kamenik */ + +#ifndef TS_EXCEPTION_H +#define TS_EXCEPTION_H + +#include +#include + +#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 +