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

Adding new test mexFunction + some minor updates and corrections 

git-svn-id: https://www.dynare.org/svn/dynare/trunk@3103 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
george 2009-10-30 13:21:03 +00:00
parent c6ae95e084
commit a479888c1b
4 changed files with 274 additions and 45 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

61
mex/sources/estimation/DsgeLikelihood.cpp
View File

@ -37,7 +37,7 @@ DsgeLikelihood::DsgeLikelihood( Vector& inA_init, GeneralMatrix& inQ, GeneralMa
const Vector&INp6, const Vector&INp7, const Vector&INp3, const Vector&INp4,
Vector& INSteadyState, Vector& INconstant, GeneralParams& INdynareParams, //GeneralParams& parameterDescription,
GeneralParams& INdr, GeneralMatrix& INkstate, GeneralMatrix& INghx, GeneralMatrix& INghu
,char *dfExt)//, KordpDynare& kOrdModel, Approximation& INapprox )
,const int jcols, const char *dfExt)//, KordpDynare& kOrdModel, Approximation& INapprox )
:a_init(inA_init), Q(inQ), R(inR), T(inT), Z(inZ), Pstar(inPstar), Pinf(inPinf), H(inH), data(inData), Y(inY),
numPeriods(INnumPeriods), numVarobs(inData.numRows()), numTimeObs(inData.numCols()),
order(INorder), endo_nbr(INendo_nbr ), exo_nbr(INexo_nbr), nstatic(INnstatic), npred(INnpred),
@ -59,10 +59,11 @@ DsgeLikelihood::DsgeLikelihood( Vector& inA_init, GeneralMatrix& inQ, GeneralMa
oo_("caller","oo_");
*********/
// setting some frequently used common variables that do not need updating
// bayestopt_.mf = bayestopt_.mf1
//std::vector<double>* vll=new std::vector<double> (nper);
vll=new std::vector<double> (numTimeObs);// vector of likelihoods
kalman_algo=(int)dynareParams.getDoubleField(string("kalman_algo"));
presampleStart=1+(int)dynareParams.getDoubleField(string("presample"));
mode_compute=(int)dr.getDoubleField(string("mode_compute"));
// Pepare data for Constructing k-order-perturbation classes
//const char *
@ -71,20 +72,11 @@ DsgeLikelihood::DsgeLikelihood( Vector& inA_init, GeneralMatrix& inQ, GeneralMa
double qz_criterium = dynareParams.getDoubleField(string("qz_criterium"));//qz_criterium = 1+1e-6;
int nMax_lag =(int)dynareParams.getDoubleField(string("maximum_lag"));
const int nBoth=(int)dr.getDoubleField(string("nboth"));
const int nPred = endo_nbr-nBoth; // correct nPred for nBoth.
const int nPred = npred-nBoth; // correct nPred for nBoth.
//vector<int> *var_order_vp = &order_var;
TwoDMatrix *vCov = new TwoDMatrix(Q);
vCov = new TwoDMatrix(Q);
// the lag, current and lead blocks of the jacobian respectively
TwoDMatrix *llincidence // = new TwoDMatrix(nrows, npar, dparams);
=(TwoDMatrix *)&(dynareParams.getMatrixField(string("lead_lag_incidence")));
//const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns
int nsPred=(int)dr.getDoubleField(string("nspred"));
int nsForw=(int)dr.getDoubleField(string("nsfwrd"));
const int jcols = exo_nbr+endo_nbr+nsPred+nsForw;
mexPrintf("k_order_perturbation: jcols= %d .\n", jcols);
llincidence = new TwoDMatrix (dynareParams.getMatrixField(string("lead_lag_incidence")));
charArraySt * casOrdEndoNames=dynareParams.getCharArrayField(string("var_order_endo_names"));
const char **endoNamesMX=(const char ** )casOrdEndoNames->charArrayPtr;
@ -113,14 +105,11 @@ DsgeLikelihood::DsgeLikelihood( Vector& inA_init, GeneralMatrix& inQ, GeneralMa
// make journal name and journal
std::string jName(fName); //params.basename);
jName += ".jnl";
Journal *journal= new Journal(jName.c_str());
journal= new Journal(jName.c_str());
#ifdef DEBUG
mexPrintf("k_order_perturbation: Calling dynamicDLL constructor.\n");
#endif
// DynamicFn * pDynamicFn = loadModelDynamicDLL (fname);
// DynamicModelDLL dynamicDLL(fName, nEndo, jcols, nMax_lag, nExog, dfExt);
//DynamicModelDLL* dynamicDLLp=new DynamicModelDLL (fName, endo_nbr, jcols, nMax_lag, exo_nbr, dfExt);
DynamicModelDLL dynamicDLLp(fName, endo_nbr, jcols, nMax_lag, exo_nbr, dfExt);
dynamicDLLp=new DynamicModelDLL (fName, endo_nbr, jcols, nMax_lag, exo_nbr, dfExt);
// intiate tensor library
#ifdef DEBUG
@ -134,7 +123,7 @@ DsgeLikelihood::DsgeLikelihood( Vector& inA_init, GeneralMatrix& inQ, GeneralMa
// make KordpDynare object
model=new KordpDynare (endoNamesMX, endo_nbr, exoNamesMX, exo_nbr, num_dp,//nPar, // paramNames,
&SteadyState, vCov, &deepParams /*modParams*/, nstatic, nPred, nfwrd, nBoth,
jcols, &NNZD, nSteps, order, *journal, dynamicDLLp,
jcols, &NNZD, nSteps, order, *journal, *dynamicDLLp,
sstol, &order_var /*var_order_vp*/, llincidence, qz_criterium);
// construct main K-order approximation class
@ -212,6 +201,32 @@ DsgeLikelihood::DsgeLikelihood( Vector& inA_init, GeneralMatrix& inQ, GeneralMa
};
DsgeLikelihood::~DsgeLikelihood()
{
delete approx;
delete model;
delete dynamicDLLp;
delete journal;
delete llincidence;
delete vCov;
delete vll;
delete &H;
delete &Q;
delete &SteadyState;
delete &kstate;
delete &param_ub;
delete &param_lb;
delete &pshape;
delete &p6;
delete &p7;
delete &p3;
delete &p4;
delete &xparam1;
delete &deepParams;
delete &ghx;
delete &ghu;
}
double
DsgeLikelihood::CalcLikelihood(Vector& xparams)
// runs all routines needed to calculate likelihood
@ -542,7 +557,7 @@ DsgeLikelihood::InitiateKalmanMatrices()
disclyap_fast(T,RQRt,Ptmp, lyapunov_tol, 1); // 1 to check chol
Pstar=Ptmp;
//Pinf=[]
Pinf = *(new GeneralMatrix(np,np));
//Pinf = *(new GeneralMatrix(np,np));
Pinf.zeros();
//Z= zeros(size(data,1), size(T,2))
@ -563,7 +578,7 @@ DsgeLikelihood::InitiateKalmanMatrices()
for (int i = 0;i<numVarobs;++i)
Ztmp.get(i,mf[i])=1.0;
Z=Ztmp;
delete &Ztmp;
}

20
mex/sources/estimation/DsgeLikelihood.h
View File

@ -31,13 +31,7 @@ class DsgeLikelihood
{
double likelihood; // sum of vector of KF step log likelihoods
vector<double>* vll; // vector of KF step log likelihoods
/*******
MexStruct &options_;
MexStruct &M_;
MexStruct &bayestOptions_;
MexStruct &dr_;
MexStruct &oo_;
***************/
Vector& a_init;//initial value of the state, usually set to 0.
GeneralMatrix& Q;// Kalman Matrices
GeneralMatrix& R;
@ -91,9 +85,12 @@ class DsgeLikelihood
GeneralMatrix& ghx;
GeneralMatrix& ghu;
//DynamicModelDLL* dynamicDLLp;
DynamicModelDLL* dynamicDLLp;
Journal *journal;
KordpDynare* model;// to be initialised by high level calling function
Approximation* approx;
TwoDMatrix *llincidence;
TwoDMatrix *vCov;
//friend class BasicKalmanTask;
//BasicKalmanTask bkt;
//friend class KalmanUniTask;
@ -110,7 +107,7 @@ class DsgeLikelihood
double KalmanFilter(double riccatiTol,bool uni);// calls Kalman
public:
DsgeLikelihood( Vector& inA_init, GeneralMatrix& inQ, GeneralMatrix& R,
DsgeLikelihood( Vector& inA_init, GeneralMatrix& inQ, GeneralMatrix& inR,
GeneralMatrix& inT, GeneralMatrix& inZ, GeneralMatrix& inPstar, GeneralMatrix& inPinf,
GeneralMatrix& inH, const GeneralMatrix&inData, GeneralMatrix&inY,
const int INnumPeriods, // const int INnumVarobs, // const int INnumTimeObs,
@ -123,7 +120,7 @@ class DsgeLikelihood
Vector& INSteadyState, Vector& INconstant, GeneralParams& INdynareParams,
//GeneralParams& parameterDescription,
GeneralParams& INdr, GeneralMatrix& INkstate, GeneralMatrix& INghx, GeneralMatrix& INghu
,char *dfExt); //, KordpDynare& inModel, Approximation& INapprox );
,const int jcols, const char *dfExt); //, KordpDynare& inModel, Approximation& INapprox );
DsgeLikelihood( const Vector&params,const GeneralMatrix&data, const vector<int>& data_index, const int gend,
const int number_of_observations, const bool no_more_missing_observations);//, KordpDynare& model ); // constructor, and
DsgeLikelihood( GeneralParams& options_,GeneralParams& M_,GeneralParams& bayestopt_, GeneralMatrix& inData,
@ -132,6 +129,9 @@ class DsgeLikelihood
double CalcLikelihood(Vector& xparams);// runs all routines needed to calculate likelihood
double getLik(){return likelihood;}
int getInfo(){return info;}
int getCostFlag(){return cost_flag;}
Vector& getSteadyState(){ return SteadyState;}
vector<double>& getLikVector() {return *vll;} // vector of log likelihoods for each Kalman step
//GeneralMatrix&lyapunov_symm(const GeneralMatrix &G, const GeneralMatrix & V);
};

15
mex/sources/estimation/dr1.cpp
View File

@ -28,22 +28,13 @@
********************************************************************/
#include "mexutils.h"
#include "DsgeLikelihood.h"
#include "dynamic_dll.h"
int
DsgeLikelihood::SolveDRModel(const int endo_nbr, const int exo_nbr, const int nstatic, const int npred, int nfwrd)//dr1()
{
int info = 0;
int i;
/*
int xlen = (int)dynareParams.getDoubleField(string("maximum_endo_lead")) +
(int)dynareParams.getDoubleField(string("M_.maximum_endo_lag)) + 1;
int klen = (int)dynareParams.getDoubleField(string("maximum_endo_lag")) +
(int)dynareParams.getDoubleField(string("M_.maximum_endo_lead")) + 1;
*/
// % expanding system for Optimal Linear Regulator
if ((int)dynareParams.getDoubleField(string("ramsey_policy")))
throw SYLV_MES_EXCEPTION("K-order-perturbation can not soleve for Ramsey policy");
@ -59,11 +50,11 @@ DsgeLikelihood::SolveDRModel(const int endo_nbr, const int exo_nbr, const int ns
vector<int>span2nx(sss-exo_nbr);
for (i=0;i<sss-exo_nbr;++i)
span2nx[i]=i+1;
ghx= *(new GeneralMatrix(ghx_u, nullVec,span2nx));//ghx_u(:,1:sss-M_.exo_nbr);
ghx= ( GeneralMatrix(ghx_u, nullVec,span2nx));//ghx_u(:,1:sss-M_.exo_nbr);
vector<int>spannx2end(exo_nbr);
for (i=0;i<exo_nbr;++i)
spannx2end[i]=sss-exo_nbr+i+1;
ghu= *(new GeneralMatrix(ghx_u, nullVec,spannx2end)); //ghx_u(:,sss-M_.exo_nbr+1:end);
ghu= ( GeneralMatrix(ghx_u, nullVec,spannx2end)); //ghx_u(:,sss-M_.exo_nbr+1:end);
}
catch(int e)
{

223
mex/sources/estimation/tests/DsgeLikelihood_mexf.cpp Normal file
View File

@ -0,0 +1,223 @@
/*
* Copyright (C) 2008-2009 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/>.
*/
// set_state_space and kstate preamble
// should be performed before calling DageLikelihood, not repeatedly withing dr1.
/******************************************
* mexFunction: Matlab Inerface point and the main application driver
* for DsgeLikelihood
*****************************************************************%function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations)
% function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations)
% Evaluates the posterior kernel of a dsge model.
%
% INPUTS
% xparam1 [double] vector of model parameters.
% gend [integer] scalar specifying the number of observations.
% data [double] matrix of data
% data_index [cell] cell of column vectors
% number_of_observations [integer]
% no_more_missing_observations [integer]
%
% OUTPUTS
% fval : value of the posterior kernel at xparam1.
% cost_flag : zero if the function returns a penalty, one otherwise.
% ys : steady state of original endogenous variables
% trend_coeff :
% info : vector of informations about the penalty:
% 41: one (many) parameter(s) do(es) not satisfied the lower bound
% 42: one (many) parameter(s) do(es) not satisfied the upper bound
% vll : vector of time-step log-likelihoods at xparam1.
%
*****************************************************************/
#include "DsgeLikelihood.h"
#ifdef MATLAB
#include "mexutils.h"
#endif
extern "C" {
void
mexFunction(int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[])
{
if (nrhs < 6)
mexErrMsgTxt("Must have at least 6 input parameters.\n");
if (nlhs == 0)
mexErrMsgTxt("Must have at least 1 output parameter.\n");
GeneralMatrix params(mxGetPr(prhs[0]), mxGetM(prhs[0]), mxGetN(prhs[0]));
if(1!= MIN(params.numCols(),params.numRows()))
throw SYLV_MES_EXCEPTION("Vextor is 2D Matrix!");
Vector xparam1(params.getData());//(params.base(), MAX(params.numCols(),params.numRows()));
const int nper = (const int)mxGetScalar(prhs[1]); //gend
GeneralMatrix data(mxGetPr(prhs[2]), mxGetM(prhs[2]), mxGetN(prhs[2]));
const int num_of_observations = (const int)mxGetScalar(prhs[4]);
const bool no_more_missing_observations= (const bool)mxGetScalar(prhs[5]);
const char *dfExt = NULL; //Dyanamic file extension, e.g.".dll" or .mexw32;
if (prhs[5] != NULL)
{
const mxArray *mexExt = prhs[6];
dfExt = mxArrayToString(mexExt);
}
#ifdef DEBUG
mexPrintf("k_order_perturbation: mexExt=%s.\n", dfExt);
#endif
/***********
***************/
int numPeriods=1;
//MexStruct dynareParams();
GeneralParams& dynareParams=*(new MexStruct());
//MexStructParam& dr=dynareParams.getStructField("dr");
GeneralParams& dr=dynareParams.getStructField("dr");
vector<int>&mfys=dynareParams.getIntVectorField("mfys");
vector<int>&mf=dynareParams.getIntVectorField("mf1");
int numVarobs=data.numRows();
Vector& SteadyState=dr.getDoubleVectorField(string("ys"));
Vector constant(numVarobs);//=*(new Vector(numVarobs));// = *(new Vector(nobs));
GeneralMatrix&kstate = dr.getMatrixField(string("kstate"));
vector<int>&order_var = dr.getIntVectorField(string("order_var"));
int order=(int)dynareParams.getDoubleField(string("order"));
int endo_nbr = (int)dynareParams.getDoubleField(string("endo_nbr"));
int exo_nbr = (int)dynareParams.getDoubleField(string("exo_nbr"));
int nstatic = (int)dr.getDoubleField(string("nstatic"));
int npred = (int)dr.getDoubleField(string("npred"));
int nfwrd = (int)dr.getDoubleField(string("nfwrd"));
Vector& ub=dynareParams.getDoubleVectorField(string("ub"));
Vector& lb=dynareParams.getDoubleVectorField(string("lb"));
int num_dp=(int)dynareParams.getDoubleField(string("np"));// no of deep params
Vector& deepParams=*(new Vector(num_dp));
vector<int>&pshape=dynareParams.getIntVectorField(string("pshape"));
Vector& p6= dynareParams.getDoubleVectorField(string("p6"));
Vector& p7= dynareParams.getDoubleVectorField(string("p7"));
Vector& p3= dynareParams.getDoubleVectorField(string("p3"));
Vector& p4= dynareParams.getDoubleVectorField(string("p4"));
//const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns
int nsPred=(int)dr.getDoubleField(string("nspred"));
int nsForw=(int)dr.getDoubleField(string("nsfwrd"));
const int jcols = exo_nbr+endo_nbr+nsPred+nsForw;
GeneralMatrix& aux= dr.getMatrixField(string("transition_auxiliary_variables"));
int nr=endo_nbr+aux.numRows();
Vector& a_init=*(new Vector(numVarobs));
GeneralMatrix& Q = dynareParams.getMatrixField(string("Sigma_e"));
GeneralMatrix& H = dynareParams.getMatrixField(string("H"));
GeneralMatrix Y(data.numRows(),data.numCols());
GeneralMatrix T(nr,nr);
GeneralMatrix Z(numVarobs,nr);
GeneralMatrix Pstar(nr,nr);
GeneralMatrix R(nr,exo_nbr);
GeneralMatrix ghx(endo_nbr,jcols-exo_nbr);
GeneralMatrix ghu(endo_nbr,exo_nbr);
//Pinf=[]
GeneralMatrix Pinf (nr,nr);
Pinf.zeros();
try
{
DsgeLikelihood dl( a_init, Q, R,T, Z, Pstar, Pinf, H,data,Y,
numPeriods, // const int INnumVarobs, // const int INnumTimeObs,
order, endo_nbr, exo_nbr, nstatic, npred,nfwrd, num_of_observations,
no_more_missing_observations, order_var, mfys, mf, xparam1,
num_dp, deepParams, ub, lb, pshape, p6, p7, p3, p4, SteadyState, constant,
dynareParams, dr, kstate, ghx, ghu, jcols, dfExt);
double loglikelihood=dl.CalcLikelihood(xparam1);
/*****************************************************************
% OUTPUTS
% fval : value of the posterior kernel at xparam1.
% cost_flag : zero if the function returns a penalty, one otherwise.
% ys : steady state of original endogenous variables
% trend_coeff :
% info : vector of informations about the penalty:
% vll : vector of time-step log-likelihoods at xparam1.
*****************************************************************/
if (nlhs >= 1)
plhs[0] = mxCreateDoubleScalar(loglikelihood);
if (nlhs >= 2)
plhs[1] = mxCreateDoubleScalar((double)dl.getCostFlag());
if (nlhs >= 3)
{
plhs[2] = mxCreateDoubleMatrix(endo_nbr, 1, mxREAL);
Vector vss(mxGetPr(plhs[2]),endo_nbr);
vss= dl.getSteadyState();
}
if (nlhs >= 4)
plhs[3] = mxCreateDoubleScalar((double)dl.getCostFlag());
if (nlhs >= 5)
plhs[4] = mxCreateDoubleMatrix(numVarobs,1, mxREAL);//dummy trend_coeff
if (nlhs >= 6)
{
// output full log-likelihood array
/* Set the output pointer to the array of log likelihood. */
std::vector<double>& vll=dl.getLikVector();
plhs[5] = mxCreateDoubleMatrix(nper,1, mxREAL);
double * mxll= mxGetPr(plhs[5]);
// assign likelihood array
for (int j=0;j<nper;++j)
mxll[j]=vll[j];
}
}
catch (const KordException &e)
{
printf("Caugth Kord exception: ");
e.print();
mexPrintf("Caugth Kord exception: %s", e.get_message());
return; // e.code();
}
catch (const TLException &e)
{
printf("Caugth TL exception: ");
e.print();
return; // 255;
}
catch (SylvException &e)
{
printf("Caught Sylv exception: ");
e.printMessage();
return; // 255;
}
catch (const DynareException &e)
{
printf("Caught KordpDynare exception: %s\n", e.message());
mexPrintf("Caugth Dynare exception: %s", e.message());
return; // 255;
}
catch (const ogu::Exception &e)
{
printf("Caught ogu::Exception: ");
e.print();
mexPrintf("Caugth general exception: %s", e.message());
return; // 255;
} //catch
}; // end of mexFunction()
}; // end of extern C