Commit a11817cf authored by Houtan Bastani's avatar Houtan Bastani Committed by Sébastien Villemot

k-order: added support for m-files, added tests and modified manual

parent d0853e16
......@@ -122,7 +122,7 @@
<para>Some installation instructions for GNU Octave can be found on the <ulink url="http://www.dynare.org/DynareWiki/DynareOctave">Dynare Wiki</ulink>.</para>
<para>If you plan to use the <xref linkend="use_dll"/> option (in particular when computing third order approximations with <xref linkend="k_order_solver"/>), 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 <ulink url="http://www.dynare.org/DynareWiki/ConfigureMatlabWindowsForMexCompilation">instructions on the Dynare wiki</ulink>. Users of Octave under Linux should install the package for MEX file compilation (under Debian or Ubuntu, it is called <filename>octave3.2-headers</filename> or <filename>octave3.0-headers</filename>). If you are using Octave or MATLAB under Mac OS X, you should install the latest version of XCode: see <ulink url="http://www.dynare.org/DynareWiki/InstallOnMacOSX">instructions on the Dynare wiki</ulink>. 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.</para>
<para>If you plan to use the <xref linkend="use_dll"/> 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 <ulink url="http://www.dynare.org/DynareWiki/ConfigureMatlabWindowsForMexCompilation">instructions on the Dynare wiki</ulink>. Users of Octave under Linux should install the package for MEX file compilation (under Debian or Ubuntu, it is called <filename>octave3.2-headers</filename> or <filename>octave3.0-headers</filename>). If you are using Octave or MATLAB under Mac OS X, you should install the latest version of XCode: see <ulink url="http://www.dynare.org/DynareWiki/InstallOnMacOSX">instructions on the Dynare wiki</ulink>. 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.</para>
</sect1>
......@@ -2243,11 +2243,11 @@ steady;
</varlistentry>
<varlistentry id="order" xreflabel="order">
<term><option>order = <replaceable>INTEGER</replaceable></option></term>
<listitem><para>Order of Taylor approximation. Acceptable values are <literal>1</literal>, <literal>2</literal> and <literal>3</literal>. Note that for third order, <link linkend="k_order_solver"><option>k_order_solver</option></link> option is implied (in particular, you must specify the <xref linkend="use_dll"/> option on the <xref linkend="model"/> block and you need a working compilation environment, see <xref linkend="software-requirements"/> for more details), and only empirical moments are available (you must provide a value for <option>periods</option> option). Default: <literal>2</literal></para></listitem>
<listitem><para>Order of Taylor approximation. Acceptable values are <literal>1</literal>, <literal>2</literal> and <literal>3</literal>. Note that for third order, <link linkend="k_order_solver"><option>k_order_solver</option></link> option is implied and only empirical moments are available (you must provide a value for <option>periods</option> option). Default: <literal>2</literal></para></listitem>
</varlistentry>
<varlistentry id="k_order_solver">
<term><option>k_order_solver</option></term>
<listitem><para>Use a k-order solver, implemented in C++, instead of the default Dynare solver. When using this option, you must specify the <xref linkend="use_dll"/> option on the <xref linkend="model"/> block, and you need a working compilation environment, <foreignphrase>i.e.</foreignphrase> a working <literal>mex</literal> command (see <xref linkend="software-requirements"/> for more details). Default: disabled for order 1 and 2, enabled otherwise</para></listitem>
<listitem><para>Use a k-order solver (implemented in C++) instead of the default Dynare solver. This option is not yet compatible with the <xref linkend="bytecode"/> option. Default: disabled for order 1 and 2, enabled otherwise</para></listitem>
</varlistentry>
<varlistentry>
<term><option>periods</option> = <replaceable>INTEGER</replaceable></term>
......
......@@ -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
/*
* 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 "dynamic_abstract_class.hh"
void
DynamicModelAC::copyDoubleIntoTwoDMatData(double *dm, TwoDMatrix *tdm, int rows, int cols)
{
int dmIdx = 0;
for (int j=0; j<cols; j++)
for (int i=0; i<rows; i++)
tdm->get(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<totalCols; i++)
for (int j=0;j<(int)(colIdxVector[i+1]-colIdxVector[i]); j++, rind++)
{
newMat[retvalind0++] = rowIdxVector[rind] + 1;
newMat[retvalind1++] = i + 1;
newMat[retvalind2++] = ptr[rind];
}
return newMat;
}
/*
* 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/>.
*/
#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 &params, const Vector &ySteady,
Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) throw (DynareException) = 0;
};
#endif
......@@ -17,7 +17,6 @@
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#include "k_ord_dynare.hh"
#include "dynamic_dll.hh"
#include <sstream>
......@@ -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
......
......@@ -17,6 +17,9 @@
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef _DYNAMIC_DLL_HH
#define _DYNAMIC_DLL_HH
#if defined(_WIN32) || defined(__CYGWIN32__)
# define NOMINMAX // Do not define "min" and "max" macros
# include <windows.h>
......@@ -25,12 +28,12 @@
#endif
#include <string>
#include <dynmex.h>
#include "dynamic_abstract_class.hh"
#include "dynare_exception.h"
// <model>_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 <model>_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 &params, const Vector &ySteady,
Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) throw (DynareException);
};
#endif
/*
* 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 "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<nrhs_dynamic; i++)
mxDestroyArray(prhs[i]);
for (int i=0; i<nlhs_dynamic; i++)
mxDestroyArray(plhs[i]);
}
/*
* 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/>.
*/
#ifndef _DYNAMIC_M_HH
#define _DYNAMIC_M_HH
#include "dynamic_abstract_class.hh"
#include "mex.h"
#include <dynmex.h>
/**
* handles calls to <model>_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 &params, const Vector &ySteady,
Vector &residual, TwoDMatrix *g1, TwoDMatrix *g2, TwoDMatrix *g3) throw (DynareException);
};
#endif
......@@ -21,14 +21,17 @@
#include <vector>
#include "first_order.h"
#include "k_ord_dynare.hh"
#include "dynamic_dll.hh"
#include "dynamic_abstract_class.hh"
#include <cmath>
#include <sstream>
#include "memory_file.h"
#include <iostream>
#include <fstream>
/**************************************************************************************/
/* Dynare DynamicModel class */
/**************************************************************************************/
......@@ -38,13 +41,13 @@ KordpDynare::KordpDynare(const vector<string> &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<int> &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);
......
......@@ -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 &params, 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<int> &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
{
......
......@@ -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 <cmath>
......@@ -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
......
......@@ -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);
}
......
......@@ -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 \
......
/* 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;
// 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;
/* 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');