From a11817cfa49a936080e0664da54772faa9eaa145 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Fri, 17 Dec 2010 18:34:23 +0100 Subject: [PATCH] k-order: added support for m-files, added tests and modified manual --- doc/manual.xml | 6 +- mex/build/k_order_perturbation.am | 6 +- .../dynamic_abstract_class.cc | 56 +++++++++++ .../dynamic_abstract_class.hh | 33 +++++++ .../k_order_perturbation/dynamic_dll.cc | 5 +- .../k_order_perturbation/dynamic_dll.hh | 12 ++- mex/sources/k_order_perturbation/dynamic_m.cc | 61 ++++++++++++ mex/sources/k_order_perturbation/dynamic_m.hh | 44 +++++++++ .../k_order_perturbation/k_ord_dynare.cc | 13 ++- .../k_order_perturbation/k_ord_dynare.hh | 7 +- .../k_order_perturbation.cc | 15 ++- preprocessor/ModFile.cc | 4 +- tests/Makefile.am | 9 +- tests/k_order_perturbation/fs2000k2_m.mod | 87 ++++++++++++++++ .../{fs2000k2.mod => fs2000k2_use_dll.mod} | 0 tests/k_order_perturbation/fs2000k3_m.mod | 80 +++++++++++++++ .../{fs2000k3.mod => fs2000k3_use_dll.mod} | 0 tests/k_order_perturbation/fs2000k_1_m.mod | 98 +++++++++++++++++++ .../{fs2000k_1.mod => fs2000k_1_use_dll.mod} | 4 +- 19 files changed, 510 insertions(+), 30 deletions(-) create mode 100644 mex/sources/k_order_perturbation/dynamic_abstract_class.cc create mode 100644 mex/sources/k_order_perturbation/dynamic_abstract_class.hh create mode 100644 mex/sources/k_order_perturbation/dynamic_m.cc create mode 100644 mex/sources/k_order_perturbation/dynamic_m.hh create mode 100644 tests/k_order_perturbation/fs2000k2_m.mod rename tests/k_order_perturbation/{fs2000k2.mod => fs2000k2_use_dll.mod} (100%) create mode 100644 tests/k_order_perturbation/fs2000k3_m.mod rename tests/k_order_perturbation/{fs2000k3.mod => fs2000k3_use_dll.mod} (100%) create mode 100644 tests/k_order_perturbation/fs2000k_1_m.mod rename tests/k_order_perturbation/{fs2000k_1.mod => fs2000k_1_use_dll.mod} (95%) diff --git a/doc/manual.xml b/doc/manual.xml index c6993dce1..a47611901 100644 --- a/doc/manual.xml +++ b/doc/manual.xml @@ -122,7 +122,7 @@ Some installation instructions for GNU Octave can be found on the Dynare Wiki. - If you plan to use the option (in particular when computing third order approximations with ), you will need to install the necessary requirements for compiling MEX files on your machine. If you are using MATLAB under Windows, install a C++ compiler on your machine and configure it with MATLAB: see instructions on the Dynare wiki. Users of Octave under Linux should install the package for MEX file compilation (under Debian or Ubuntu, it is called octave3.2-headers or octave3.0-headers). If you are using Octave or MATLAB under Mac OS X, you should install the latest version of XCode: see instructions on the Dynare wiki. Mac OS X Octave users will also need to install gnuplot if they want graphing capabilities. Users of MATLAB under Linux and Mac OS X, and users of Octave under Windows, normally need to do nothing, since a working compilation environment is available by default. + If you plan to use the option, you will need to install the necessary requirements for compiling MEX files on your machine. If you are using MATLAB under Windows, install a C++ compiler on your machine and configure it with MATLAB: see instructions on the Dynare wiki. Users of Octave under Linux should install the package for MEX file compilation (under Debian or Ubuntu, it is called octave3.2-headers or octave3.0-headers). If you are using Octave or MATLAB under Mac OS X, you should install the latest version of XCode: see instructions on the Dynare wiki. Mac OS X Octave users will also need to install gnuplot if they want graphing capabilities. Users of MATLAB under Linux and Mac OS X, and users of Octave under Windows, normally need to do nothing, since a working compilation environment is available by default. @@ -2243,11 +2243,11 @@ steady; - Order of Taylor approximation. Acceptable values are 1, 2 and 3. Note that for third order, option is implied (in particular, you must specify the option on the block and you need a working compilation environment, see for more details), and only empirical moments are available (you must provide a value for option). Default: 2 + Order of Taylor approximation. Acceptable values are 1, 2 and 3. Note that for third order, option is implied and only empirical moments are available (you must provide a value for option). Default: 2 - Use a k-order solver, implemented in C++, instead of the default Dynare solver. When using this option, you must specify the option on the block, and you need a working compilation environment, i.e. a working mex command (see for more details). Default: disabled for order 1 and 2, enabled otherwise + Use a k-order solver (implemented in C++) instead of the default Dynare solver. This option is not yet compatible with the option. Default: disabled for order 1 and 2, enabled otherwise = INTEGER diff --git a/mex/build/k_order_perturbation.am b/mex/build/k_order_perturbation.am index 57e8a6131..28deafc96 100644 --- a/mex/build/k_order_perturbation.am +++ b/mex/build/k_order_perturbation.am @@ -15,4 +15,8 @@ nodist_k_order_perturbation_SOURCES = \ k_ord_dynare.cc \ k_ord_dynare.hh \ dynamic_dll.cc \ - dynamic_dll.hh + dynamic_dll.hh \ + dynamic_abstract_class.cc \ + dynamic_abstract_class.hh \ + dynamic_m.cc \ + dynamic_m.hh diff --git a/mex/sources/k_order_perturbation/dynamic_abstract_class.cc b/mex/sources/k_order_perturbation/dynamic_abstract_class.cc new file mode 100644 index 000000000..9825c1344 --- /dev/null +++ b/mex/sources/k_order_perturbation/dynamic_abstract_class.cc @@ -0,0 +1,56 @@ +/* + * 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 "dynamic_abstract_class.hh" + +void +DynamicModelAC::copyDoubleIntoTwoDMatData(double *dm, TwoDMatrix *tdm, int rows, int cols) +{ + int dmIdx = 0; + for (int j=0; jget(i,j) = dm[dmIdx++]; +} + +double * +DynamicModelAC::unpackSparseMatrix(mxArray *sparseMat) +{ + int totalCols = mxGetN(sparseMat); + mwIndex *rowIdxVector = mxGetIr(sparseMat); + mwSize sizeRowIdxVector = mxGetNzmax(sparseMat); + mwIndex *colIdxVector = mxGetJc(sparseMat); + + double *ptr = mxGetPr(sparseMat); + double *newMat = (double *)malloc(sizeRowIdxVector*3*sizeof(double)); + + int rind = 0; + int retvalind0 = 0; + int retvalind1 = sizeRowIdxVector; + int retvalind2 = sizeRowIdxVector*2; + + for (int i=0; i. + */ + +#ifndef _DYNAMICMODELAC_HH +#define _DYNAMICMODELAC_HH + +#include "k_ord_dynare.hh" + +class DynamicModelAC +{ +public: + double *unpackSparseMatrix(mxArray *sparseMatrix); + void copyDoubleIntoTwoDMatData(double *dm, TwoDMatrix *tdm, int rows, int cols); + virtual void eval(const Vector &y, const Vector &x, const Vector ¶ms, const Vector &ySteady, + Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) throw (DynareException) = 0; +}; +#endif diff --git a/mex/sources/k_order_perturbation/dynamic_dll.cc b/mex/sources/k_order_perturbation/dynamic_dll.cc index a7943db9e..64f3b0aba 100644 --- a/mex/sources/k_order_perturbation/dynamic_dll.cc +++ b/mex/sources/k_order_perturbation/dynamic_dll.cc @@ -17,7 +17,6 @@ * along with Dynare. If not, see . */ -#include "k_ord_dynare.hh" #include "dynamic_dll.hh" #include @@ -36,7 +35,7 @@ DynamicModelDLL::DynamicModelDLL(const string &modName, const string &sExt) thro dynamicHinstance = LoadLibrary(fName.c_str()); if (dynamicHinstance == NULL) throw 1; - Dynamic = (DynamicFn) GetProcAddress(dynamicHinstance, "Dynamic"); + Dynamic = (DynamicDLLFn) GetProcAddress(dynamicHinstance, "Dynamic"); if (Dynamic == NULL) { FreeLibrary(dynamicHinstance); // Free the library @@ -49,7 +48,7 @@ DynamicModelDLL::DynamicModelDLL(const string &modName, const string &sExt) thro cerr << dlerror() << endl; throw 1; } - Dynamic = (DynamicFn) dlsym(dynamicHinstance, "Dynamic"); + Dynamic = (DynamicDLLFn) dlsym(dynamicHinstance, "Dynamic"); if ((Dynamic == NULL) || dlerror()) { dlclose(dynamicHinstance); // Free the library diff --git a/mex/sources/k_order_perturbation/dynamic_dll.hh b/mex/sources/k_order_perturbation/dynamic_dll.hh index e9ce127fc..937e7483d 100644 --- a/mex/sources/k_order_perturbation/dynamic_dll.hh +++ b/mex/sources/k_order_perturbation/dynamic_dll.hh @@ -17,6 +17,9 @@ * along with Dynare. If not, see . */ +#ifndef _DYNAMIC_DLL_HH +#define _DYNAMIC_DLL_HH + #if defined(_WIN32) || defined(__CYGWIN32__) # define NOMINMAX // Do not define "min" and "max" macros # include @@ -25,12 +28,12 @@ #endif #include -#include +#include "dynamic_abstract_class.hh" #include "dynare_exception.h" // _Dynamic DLL pointer -typedef void (*DynamicFn) +typedef void (*DynamicDLLFn) (const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, double *residual, double *g1, double *g2, double *g3); @@ -38,10 +41,10 @@ typedef void (*DynamicFn) * creates pointer to Dynamic function inside _dynamic.dll * and handles calls to it. **/ -class DynamicModelDLL +class DynamicModelDLL : public DynamicModelAC { private: - DynamicFn Dynamic; // pointer to the Dynamic function in DLL + DynamicDLLFn Dynamic; // pointer to the Dynamic function in DLL #if defined(_WIN32) || defined(__CYGWIN32__) HINSTANCE dynamicHinstance; // DLL instance pointer in Windows #else @@ -56,3 +59,4 @@ public: void eval(const Vector &y, const Vector &x, const Vector ¶ms, const Vector &ySteady, Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) throw (DynareException); }; +#endif diff --git a/mex/sources/k_order_perturbation/dynamic_m.cc b/mex/sources/k_order_perturbation/dynamic_m.cc new file mode 100644 index 000000000..4096c4c36 --- /dev/null +++ b/mex/sources/k_order_perturbation/dynamic_m.cc @@ -0,0 +1,61 @@ +/* + * 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 "dynamic_m.hh" + +DynamicModelMFile::DynamicModelMFile(const string &modName) throw (DynareException) : + DynamicMFilename(modName + "_dynamic") +{ +} + +DynamicModelMFile::~DynamicModelMFile() +{ +} + +void +DynamicModelMFile::eval(const Vector &y, const Vector &x, const Vector &modParams, const Vector &ySteady, + Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) throw (DynareException) +{ + mxArray *prhs[nrhs_dynamic], *plhs[nlhs_dynamic]; + + prhs[0] = mxCreateDoubleMatrix(y.length(), 1, mxREAL); + prhs[1] = mxCreateDoubleMatrix(1, x.length(), mxREAL); + prhs[2] = mxCreateDoubleMatrix(modParams.length(), 1, mxREAL); + prhs[3] = mxCreateDoubleScalar(1.0); + + memcpy((void *)(mxGetPr(prhs[0])), (void *)y.base(), y.length()*sizeof(double)); + memcpy((void *)(mxGetPr(prhs[1])), (void *)x.base(), x.length()*sizeof(double)); + memcpy((void *)(mxGetPr(prhs[2])), (void *)modParams.base(), modParams.length()*sizeof(double)); + + int retVal = mexCallMATLAB(nlhs_dynamic, plhs, nrhs_dynamic, prhs, DynamicMFilename.c_str()); + if (retVal != 0) + throw DynareException(__FILE__, __LINE__, "Trouble calling " + DynamicMFilename); + + residual = Vector(mxGetPr(plhs[0]), residual.skip(), (int)mxGetM(plhs[0])); + copyDoubleIntoTwoDMatData(mxGetPr(plhs[1]), g1, (int)mxGetM(plhs[1]), (int)mxGetN(plhs[1])); + if (g2 != NULL) + copyDoubleIntoTwoDMatData(unpackSparseMatrix(plhs[2]), g2, (int)mxGetNzmax(plhs[2]), 3); + if (g3 != NULL) + copyDoubleIntoTwoDMatData(unpackSparseMatrix(plhs[3]), g3, (int)mxGetNzmax(plhs[3]), 3); + + for (int i=0; i. + */ + +#ifndef _DYNAMIC_M_HH +#define _DYNAMIC_M_HH + +#include "dynamic_abstract_class.hh" + +#include "mex.h" +#include + +/** + * handles calls to _dynamic.m + * + **/ +class DynamicModelMFile : public DynamicModelAC +{ +private: + const string DynamicMFilename; + const static int nlhs_dynamic = 4; + const static int nrhs_dynamic = 4; +public: + DynamicModelMFile(const string &modName) throw (DynareException); + virtual ~DynamicModelMFile(); + void eval(const Vector &y, const Vector &x, const Vector ¶ms, const Vector &ySteady, + Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) throw (DynareException); +}; +#endif diff --git a/mex/sources/k_order_perturbation/k_ord_dynare.cc b/mex/sources/k_order_perturbation/k_ord_dynare.cc index 907ed7381..e24de39fe 100644 --- a/mex/sources/k_order_perturbation/k_ord_dynare.cc +++ b/mex/sources/k_order_perturbation/k_ord_dynare.cc @@ -21,14 +21,17 @@ #include #include "first_order.h" -#include "k_ord_dynare.hh" -#include "dynamic_dll.hh" +#include "dynamic_abstract_class.hh" #include #include #include "memory_file.h" + +#include +#include + /**************************************************************************************/ /* Dynare DynamicModel class */ /**************************************************************************************/ @@ -38,13 +41,13 @@ KordpDynare::KordpDynare(const vector &endo, int num_endo, Vector &ysteady, TwoDMatrix &vcov, Vector &inParams, int nstat, int npred, int nforw, int nboth, const int jcols, const Vector &nnzd, const int nsteps, int norder, - Journal &jr, DynamicModelDLL &dynamicDLL, double sstol, + Journal &jr, DynamicModelAC *dynamicModelFile_arg, double sstol, const vector &var_order, const TwoDMatrix &llincidence, double criterium) throw (TLException) : nStat(nstat), nBoth(nboth), nPred(npred), nForw(nforw), nExog(nexog), nPar(npar), nYs(npred + nboth), nYss(nboth + nforw), nY(num_endo), nJcols(jcols), NNZD(nnzd), nSteps(nsteps), nOrder(norder), journal(jr), ySteady(ysteady), params(inParams), vCov(vcov), md(1), dnl(*this, endo), denl(*this, exo), dsnl(*this, dnl, denl), ss_tol(sstol), varOrder(var_order), - ll_Incidence(llincidence), qz_criterium(criterium), dynamicDLL(dynamicDLL) + ll_Incidence(llincidence), qz_criterium(criterium), dynamicModelFile(dynamicModelFile_arg) { ReorderDynareJacobianIndices(); @@ -113,7 +116,7 @@ KordpDynare::calcDerivativesAtSteady() throw (DynareException) Vector llxSteady(nJcols-nExog); LLxSteady(ySteady, llxSteady); - dynamicDLL.eval(llxSteady, xx, params, ySteady, out, &g1, g2p, g3p); + dynamicModelFile->eval(llxSteady, xx, params, ySteady, out, &g1, g2p, g3p); populateDerivativesContainer(g1, 1, JacobianIndices); diff --git a/mex/sources/k_order_perturbation/k_ord_dynare.hh b/mex/sources/k_order_perturbation/k_ord_dynare.hh index f7998a610..d275ed900 100644 --- a/mex/sources/k_order_perturbation/k_ord_dynare.hh +++ b/mex/sources/k_order_perturbation/k_ord_dynare.hh @@ -80,13 +80,16 @@ public: /*********************************************/ // The following only implements DynamicModel with help of ogdyn::DynareModel // instantiation of pure abstract DynamicModel decl. in dynamic_model.h +class DynamicModelAC; class DynamicModelDLL; +class DynamicModelMFile; class KordpDynare : public DynamicModel { friend class DynareNameList; friend class DynareStateNameList; friend class DynamicModelDLL; + friend class DynamicModelMFile; const int nStat; const int nBoth; @@ -120,7 +123,7 @@ public: Vector &ySteady, TwoDMatrix &vCov, Vector ¶ms, int nstat, int nPred, int nforw, int nboth, const int nJcols, const Vector &NNZD, const int nSteps, const int ord, - Journal &jr, DynamicModelDLL &dynamicDLL, double sstol, + Journal &jr, DynamicModelAC *dynamicModelFile_arg, double sstol, const vector &varOrder, const TwoDMatrix &ll_Incidence, double qz_criterium) throw (TLException); @@ -222,7 +225,7 @@ public: void evaluateSystem(Vector &out, const Vector &yym, const Vector &yy, const Vector &yyp, const Vector &xx) throw (DynareException); void calcDerivativesAtSteady() throw (DynareException); - DynamicModelDLL &dynamicDLL; + DynamicModelAC *dynamicModelFile; DynamicModel * clone() const { diff --git a/mex/sources/k_order_perturbation/k_order_perturbation.cc b/mex/sources/k_order_perturbation/k_order_perturbation.cc index 530c8929c..699ee0db5 100644 --- a/mex/sources/k_order_perturbation/k_order_perturbation.cc +++ b/mex/sources/k_order_perturbation/k_order_perturbation.cc @@ -24,8 +24,7 @@ 1) dr 2) M_ 3) options - 4) oo_ - 5) string containing the MEX extension (with a dot at the beginning) + 4) string containing the MEX extension (with a dot at the beginning) Outputs: - if order == 1: only g_1 @@ -33,7 +32,7 @@ - if order == 3: g_0, g_1, g_2, g_3 */ -#include "k_ord_dynare.hh" +#include "dynamic_m.hh" #include "dynamic_dll.hh" #include @@ -76,6 +75,7 @@ extern "C" { const mxArray *dr = prhs[0]; const mxArray *M_ = prhs[1]; const mxArray *options_ = prhs[2]; + int use_dll = (int)mxGetScalar(mxGetField(options_, 0, "use_dll")); mxArray *mFname = mxGetField(M_, 0, "fname"); if (!mxIsChar(mFname)) @@ -199,7 +199,12 @@ extern "C" { std::string jName(fName); //params.basename); jName += ".jnl"; Journal journal(jName.c_str()); - DynamicModelDLL dynamicDLL(fName, dfExt); + + DynamicModelAC *dynamicModelFile; + if (use_dll == 1) + dynamicModelFile = new DynamicModelDLL(fName, dfExt); + else + dynamicModelFile = new DynamicModelMFile(fName); // intiate tensor library tls.init(kOrder, nStat+2*nPred+3*nBoth+2*nForw+nExog); @@ -207,7 +212,7 @@ extern "C" { // make KordpDynare object KordpDynare dynare(endoNames, nEndo, exoNames, nExog, nPar, ySteady, vCov, modParams, nStat, nPred, nForw, nBoth, - jcols, NNZD, nSteps, kOrder, journal, dynamicDLL, + jcols, NNZD, nSteps, kOrder, journal, dynamicModelFile, sstol, var_order_vp, llincidence, qz_criterium); // construct main K-order approximation class diff --git a/preprocessor/ModFile.cc b/preprocessor/ModFile.cc index 3c603959a..bbdb1c4d4 100644 --- a/preprocessor/ModFile.cc +++ b/preprocessor/ModFile.cc @@ -129,9 +129,9 @@ ModFile::checkPass() exit(EXIT_FAILURE); } - if (mod_file_struct.k_order_solver && !use_dll) + if (mod_file_struct.k_order_solver && byte_code) { - cerr << "ERROR: When using option 'k_order_solver' (which is implicit if order >= 3), you must specify option 'use_dll' on the 'model' block" << endl; + cerr << "ERROR: 'k_order_solver' (which is implicit if order >= 3), is not yet compatible with 'bytecode'." << endl; exit(EXIT_FAILURE); } diff --git a/tests/Makefile.am b/tests/Makefile.am index 0f74f47d5..17e3376a0 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -27,9 +27,12 @@ MODS = \ block_bytecode/ireland.mod \ block_bytecode/ramst_normcdf_and_friends.mod \ k_order_perturbation/fs2000k2a.mod \ - k_order_perturbation/fs2000k2.mod \ - k_order_perturbation/fs2000k_1.mod \ - k_order_perturbation/fs2000k3.mod \ + k_order_perturbation/fs2000k2_use_dll.mod \ + k_order_perturbation/fs2000k_1_use_dll.mod \ + k_order_perturbation/fs2000k3_use_dll.mod \ + k_order_perturbation/fs2000k2_m.mod \ + k_order_perturbation/fs2000k_1_m.mod \ + k_order_perturbation/fs2000k3_m.mod \ partial_information/PItest3aHc0PCLsimModPiYrVarobsAll.mod \ partial_information/PItest3aHc0PCLsimModPiYrVarobsCNR.mod \ arima/mod1.mod \ diff --git a/tests/k_order_perturbation/fs2000k2_m.mod b/tests/k_order_perturbation/fs2000k2_m.mod new file mode 100644 index 000000000..8d7a606ae --- /dev/null +++ b/tests/k_order_perturbation/fs2000k2_m.mod @@ -0,0 +1,87 @@ +/* Checks that, for order = 2, k_order_solver = 0 (fs2000k2a) + and k_order_solver = 1 (this file) give the same results */ + +var m P c e W R k d n l gy_obs gp_obs y dA ; +varexo e_a e_m; + +parameters alp bet gam mst rho psi del; + +alp = 0.33; +bet = 0.99; +gam = 0.003; +mst = 1.011; +rho = 0.7; +psi = 0.787; +del = 0.02; + +model; +dA = exp(gam+e_a); +log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m; +-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0; +W = l/n; +-(psi/(1-psi))*(c*P/(1-n))+l/n = 0; +R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W; +1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0; +c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1); +P*c = m; +m-1+d = l; +e = exp(e_a); +y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a)); +gy_obs = dA*y/y(-1); +gp_obs = (P/P(-1))*m(-1)/dA; +end; + +initval; +m = mst; +P = 2.25; +c = 0.45; +e = 1; +W = 4; +R = 1.02; +k = 6; +d = 0.85; +n = 0.19; +l = 0.86; +y = 0.6; +gy_obs = exp(gam); +gp_obs = exp(-gam); +dA = exp(gam); +end; + +shocks; +var e_a; stderr 0.014; +var e_m; stderr 0.005; +end; + +steady; + +stoch_simul(order=2,k_order_solver,periods=1000); + +if ~exist('fs2000k2a_results.mat','file'); + error('fs2000k2a must be run first'); +end; + +oo1 = load('fs2000k2a_results','oo_'); + +dr0 = oo1.oo_.dr; +dr = oo_.dr; + +if max(max(abs(dr0.ghx - dr.ghx))) > 1e-12; + error('error in ghx'); +end; +if max(max(abs(dr0.ghu - dr.ghu))) > 1e-12; + error('error in ghu'); +end; +if max(max(abs(dr0.ghxx - dr.ghxx))) > 1e-12; + error('error in ghxx'); +end; +if max(max(abs(dr0.ghuu - dr.ghuu))) > 1e-12; + error('error in ghuu'); +end; +if max(max(abs(dr0.ghxu - dr.ghxu))) > 1e-12; + error('error in ghxu'); +end; +if max(max(abs(dr0.ghs2 - dr.ghs2))) > 1e-12; + error('error in ghs2'); +end; + diff --git a/tests/k_order_perturbation/fs2000k2.mod b/tests/k_order_perturbation/fs2000k2_use_dll.mod similarity index 100% rename from tests/k_order_perturbation/fs2000k2.mod rename to tests/k_order_perturbation/fs2000k2_use_dll.mod diff --git a/tests/k_order_perturbation/fs2000k3_m.mod b/tests/k_order_perturbation/fs2000k3_m.mod new file mode 100644 index 000000000..a3130b141 --- /dev/null +++ b/tests/k_order_perturbation/fs2000k3_m.mod @@ -0,0 +1,80 @@ +// checks whether second order coefficients are the same with order=2 and order=3 with k_order_solver=1 + +var m P c e W R k d n l gy_obs gp_obs y dA ; +varexo e_a e_m; + +parameters alp bet gam mst rho psi del; + +alp = 0.33; +bet = 0.99; +gam = 0.003; +mst = 1.011; +rho = 0.7; +psi = 0.787; +del = 0.02; + +model; +dA = exp(gam+e_a); +log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m; +-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0; +W = l/n; +-(psi/(1-psi))*(c*P/(1-n))+l/n = 0; +R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W; +1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0; +c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1); +P*c = m; +m-1+d = l; +e = exp(e_a); +y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a)); +gy_obs = dA*y/y(-1); +gp_obs = (P/P(-1))*m(-1)/dA; +end; + +initval; +m = mst; +P = 2.25; +c = 0.45; +e = 1; +W = 4; +R = 1.02; +k = 6; +d = 0.85; +n = 0.19; +l = 0.86; +y = 0.6; +gy_obs = exp(gam); +gp_obs = exp(-gam); +dA = exp(gam); +end; + +shocks; +var e_a; stderr 0.014; +var e_m; stderr 0.005; +end; + +steady; + +stoch_simul(order=3,periods=1000); + +if ~exist('fs2000k2a_results.mat','file'); + error('fs2000k2a must be run first'); +end; + +oo1 = load('fs2000k2a_results','oo_'); + +dr0 = oo1.oo_.dr; +dr = oo_.dr; + +if max(max(abs(dr0.ghxx - dr.ghxx))) > 1e-12; + error('error in ghxx'); +end; +if max(max(abs(dr0.ghuu - dr.ghuu))) > 1e-12; + error('error in ghuu'); +end; +if max(max(abs(dr0.ghxu - dr.ghxu))) > 1e-12; + error('error in ghxu'); +end; +if max(max(abs(dr0.ghs2 - dr.ghs2))) > 1e-12; + error('error in ghs2'); +end; + diff --git a/tests/k_order_perturbation/fs2000k3.mod b/tests/k_order_perturbation/fs2000k3_use_dll.mod similarity index 100% rename from tests/k_order_perturbation/fs2000k3.mod rename to tests/k_order_perturbation/fs2000k3_use_dll.mod diff --git a/tests/k_order_perturbation/fs2000k_1_m.mod b/tests/k_order_perturbation/fs2000k_1_m.mod new file mode 100644 index 000000000..bb19c46fe --- /dev/null +++ b/tests/k_order_perturbation/fs2000k_1_m.mod @@ -0,0 +1,98 @@ +/* Checks that, for order = 2 and k_order_solver = 1, a model with 2 leads + and the same model with one lead (using auxiliary vars) give the same result */ + +var m m_1 P P_1 c e W R k d n l gy_obs gp_obs y dA AUXv; +varexo e_a e_m; + +parameters alp bet gam mst rho psi del; + +alp = 0.33; +bet = 0.99; +gam = 0.003; +mst = 1.011; +rho = 0.7; +psi = 0.787; +del = 0.02; + +model; +dA = exp(gam+e_a); +log(m) = (1-rho)*log(mst) + rho*log(m_1(-1))+e_m; +-P/(c(+1)*P(+1)*m)+AUXv(+1)=0; +W = l/n; +-(psi/(1-psi))*(c*P/(1-n))+l/n = 0; +R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W; +1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0; +c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1); +P*c = m; +m-1+d = l; +e = exp(e_a); +y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a)); +gy_obs = dA*y/y(-1); +gp_obs = (P/P_1(-1))*m_1(-1)/dA; +m_1 = m; +P_1 = P; +AUXv = bet*P*(alp*exp(-alp*(gam+log(e)))*k(-1)^(alp-1)*n^(1-alp)+(1-del)*exp(-(gam+log(e))))/(c(+1)*P(+1)*m); +end; + +initval; +m = mst; +m_1=mst; +P = 2.25; +P_1 = 2.25; +c = 0.45; +e = 1; +W = 4; +R = 1.02; +k = 6; +d = 0.85; +n = 0.19; +l = 0.86; +y = 0.6; +gy_obs = exp(gam); +gp_obs = exp(-gam); +dA = exp(gam); +AUXv = 1; +end; + +shocks; +var e_a; stderr 0.014; +var e_m; stderr 0.005; +end; + +steady; + +stoch_simul(order=2,k_order_solver,irf=0); + +if ~exist('fs2000k2_m_results.mat','file'); + error('fs2000k2 must be run first'); +end; + +oo1 = load('fs2000k2_m_results','oo_'); + +dr0 = oo1.oo_.dr; +dr = oo_.dr; + +ikr = [2:10 1 13:17]; +ikc = [1 3 4 2]; +ikc2 = [1 3 4 2 9 11 12 10 13 15 16 14 5 7 8 6]; +ikc2u = [1 2 5 6 7 8 3 4]; + +if max(max(abs(dr0.ghx - dr.ghx(ikr,ikc)))) > 1e-12; + error('error in ghx'); +end; +if max(max(abs(dr0.ghu - dr.ghu(ikr,:)))) > 1e-12; + error('error in ghu'); +end; +if max(max(abs(dr0.ghxx - dr.ghxx(ikr,ikc2)))) > 1e-12; + error('error in ghxx'); +end; +if max(max(abs(dr0.ghuu - dr.ghuu(ikr,:)))) > 1e-12; + error('error in ghuu'); +end; +if max(max(abs(dr0.ghxu - dr.ghxu(ikr,ikc2u)))) > 1e-12; + error('error in ghxu'); +end; +if max(max(abs(dr0.ghs2 - dr.ghs2(ikr,:)))) > 1e-12; + error('error in ghs2'); +end; + diff --git a/tests/k_order_perturbation/fs2000k_1.mod b/tests/k_order_perturbation/fs2000k_1_use_dll.mod similarity index 95% rename from tests/k_order_perturbation/fs2000k_1.mod rename to tests/k_order_perturbation/fs2000k_1_use_dll.mod index 40fcd5ce8..5fcf5b4cd 100644 --- a/tests/k_order_perturbation/fs2000k_1.mod +++ b/tests/k_order_perturbation/fs2000k_1_use_dll.mod @@ -63,11 +63,11 @@ steady; stoch_simul(order=2,k_order_solver,irf=0); -if ~exist('fs2000k2_results.mat','file'); +if ~exist('fs2000k2_use_dll_results.mat','file'); error('fs2000k2 must be run first'); end; -oo1 = load('fs2000k2_results','oo_'); +oo1 = load('fs2000k2_use_dll_results','oo_'); dr0 = oo1.oo_.dr; dr = oo_.dr;