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

Merge branch 'master' of kirikou.dynare.org:/srv/d_kirikou/git/dynare

time-shift
Stéphane Adjemian (Scylla) 2010-05-06 17:41:01 +02:00
parent 79184b2ded ba0503628e
commit 7f8afdd8b0
27 changed files with 543 additions and 2759 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -1,797 +0,0 @@
/*
* 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/>.
*/
#include "mexutils.h"
#include "DsgeLikelihood.h"
#include "utils.h"
#include "disclyap_fast.h"
#include "ts_exception.h"
#include "dynamic_dll.h"
DsgeLikelihood::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,
const int INorder, const int INendo_nbr, const int INexo_nbr, const int INnstatic,
const int INnpred, const int INnfwrd, const int INnum_of_observations, const bool INno_more_missing_observations,
const vector<int>& INorder_var, const vector<int>& INmfys, const vector<int>& INmf1,
Vector& INxparam1, const int INnum_dp, Vector& INdeepParams,
const Vector& INub, const Vector& INlb, const vector<int>&INpshape,
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,
GeneralMatrix& INaux, vector<int>&INiv, vector<int>&INic
,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),
nfwrd(INnfwrd), number_of_observations(INnum_of_observations), no_more_missing_observations(INno_more_missing_observations),
order_var(INorder_var), mfys(INmfys), mf(INmf1), xparam1(INxparam1),
num_dp(INnum_dp), deepParams(INdeepParams), //num_dp((int)dynareParams.getDoubleField(string("np),// no of deep params
param_ub(INub), param_lb(INlb),pshape(INpshape), //pshape(dynareParams.getIntVectorField(string("pshape),
p6(INp6), p7(INp7), p3(INp3), p4(INp4), SteadyState(INSteadyState),
constant(INconstant), dynareParams(INdynareParams),
dr(INdr), kstate(INkstate), ghx(INghx),ghu(INghu),
aux(INaux), iv(INiv), ic(INic)
//, model(kOrdModel), approx(INapprox )
{
/*****
bayestOptions_("caller","bayestopt_");
options_("caller","options_");
M_Options_("caller","M_");
dr_("caller","dr");
oo_("caller","oo_");
*********/
// setting some frequently used common variables that do not need updating
//std::vector<double>* vll=new std::vector<double> (nper);
// vll=new std::vector<double> (numTimeObs);// vector of likelihoods
vll=new vector<double>(numTimeObs);// vector of likelihoods
kalman_algo=(int)dynareParams.getDoubleField(string("kalman_algo"));
presampleStart=1+(int)dynareParams.getDoubleField(string("presample"));
mode_compute=(int)dynareParams.getDoubleField(string("mode_compute"));
#ifdef DEBUG
mexPrintf("mode_compute=%d presampleStart=%d \n", mode_compute,presampleStart);
#endif
// Pepare data for Constructing k-order-perturbation classes
//const char *
string& fname=dynareParams.getStringField(string("fname"));
fName = (char *)fname.c_str();
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 = npred-nBoth; // correct nPred for nBoth.
//vector<int> *var_order_vp = &order_var;
#ifdef DEBUG
mexPrintf("fName %s qz_criterium=%f nMax_lag=%d nPred=%d :Construction of vCov\n", fName,qz_criterium,nMax_lag,nPred);
#endif
vCov = new TwoDMatrix(Q);
// the lag, current and lead blocks of the jacobian respectively
llincidence = new TwoDMatrix (dynareParams.getMatrixField(string("lead_lag_incidence")));
#ifdef DEBUG
mexPrintf("Construction of casOrdEndoNames\n");
#endif
charArraySt * casOrdEndoNames=dynareParams.getCharArrayField(string("var_order_endo_names"));
#ifdef DEBUG
mexPrintf("Construction of endoNamesMX\n");
#endif
const char **endoNamesMX=(const char ** )casOrdEndoNames->charArrayPtr;
#ifdef DEBUG
for (int i = 0; i < endo_nbr; i++)
{
mexPrintf("k_ord_perturbation: EndoNameList[%d][0]= %s.\n", i, endoNamesMX[i]);
}
#endif
charArraySt * casExoNames=dynareParams.getCharArrayField(string("exo_names"));
const char **exoNamesMX=(const char ** )casExoNames->charArrayPtr;
Vector &NNZD =dynareParams.getDoubleVectorField(string("NNZDerivatives"));
const int nSteps = 0; // Dynare++ solving steps, for time being default to 0 = deterministic steady state
const double sstol = 1.e-13; //NL solver tolerance from
// Construct k-order-perturbation classes
THREAD_GROUP::max_parallel_threads = 2; //params.num_threads;
try
{
// make journal name and journal
std::string jName(fName); //params.basename);
jName += ".jnl";
journal= new Journal(jName.c_str());
#ifdef DEBUG
mexPrintf("k_order_perturbation: Calling dynamicDLL constructor.\n");
#endif
dynamicDLLp=new DynamicModelDLL (fName, endo_nbr, jcols, nMax_lag, exo_nbr, dfExt);
// intiate tensor library
#ifdef DEBUG
mexPrintf("k_order_perturbation: Call tls init order:%d, size: %d\n", order, nstatic+2*nPred+3*nBoth+2*nfwrd+exo_nbr);
#endif
tls.init(order, nstatic+2*nPred+3*nBoth+2*nfwrd+exo_nbr);
#ifdef DEBUG
mexPrintf("estimation: Calling dynare model constructor .\n");
#endif
// 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,
sstol, &order_var /*var_order_vp*/, llincidence, qz_criterium);
// construct main K-order approximation class
#ifdef DEBUG
mexPrintf("estimation: Call Approximation constructor with qz_criterium=%f \n", qz_criterium);
#endif
approx= new Approximation(*model, *journal, nSteps, false, qz_criterium);
// run stochastic steady
#ifdef DEBUG
mexPrintf("estimation:k_order_perturbation and Approximation created.\n");
#endif
}
catch (const KordException &e)
{
printf("Caugth Kord exception: ");
e.print();
mexPrintf("Caugth Kord exception: %s", e.get_message());
std::string errfile(fName); //(params.basename);
errfile += "_error.log";
FILE *errfd = NULL;
if (NULL == (errfd = fopen(errfile.c_str(), "wb")))
{
fprintf(stderr, "Couldn't open %s for writing.\n", errfile.c_str());
return; // e.code();
}
fprintf(errfd, "Caugth Kord exception: %s", e.get_message());
fclose(errfd);
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());
std::string errfile(fName); //(params.basename);
errfile += "_error.log";
FILE *errfd = NULL;
if (NULL == (errfd = fopen(errfile.c_str(), "wb")))
{
fprintf(stderr, "Couldn't open %s for writing.\n", errfile.c_str());
return; // e.code();
}
fprintf(errfd, "Caugth KordDynare exception: %s", e.message());
fclose(errfd);
return; // 255;
}
catch (const ogu::Exception &e)
{
printf("Caught ogu::Exception: ");
e.print();
mexPrintf("Caugth general exception: %s", e.message());
std::string errfile(fName); //(params.basename);
errfile += "_error.log";
FILE *errfd = NULL;
if (NULL == (errfd = fopen(errfile.c_str(), "wb")))
{
fprintf(stderr, "Couldn't open %s for writing.\n", errfile.c_str());
return; // e.code();
}
e.print(errfd);
fclose(errfd);
return; // 255;
} //catch
};
DsgeLikelihood::~DsgeLikelihood()
{
// delete llincidence;
// delete vCov;
delete approx;
delete model;
delete dynamicDLLp;
delete journal;
delete vll;
delete &H;
delete &Q;
delete &kstate;
delete &pshape;
delete &ghx;
delete &ghu;
delete &dynareParams;
delete &dr;
delete &aux;
delete &iv;
delete &ic;
/******** deleting Vectors Crashes !************
delete &SteadyState;
delete &param_ub;
delete &param_lb;
delete &p6;
delete &p7;
delete &p3;
delete &p4;
delete &xparam1;
delete &deepParams;
**********/
}
double
DsgeLikelihood::CalcLikelihood(Vector& xparams)
// runs all routines needed to calculate likelihood
{
likelihood=0.0;
info=0;
xparam1=xparams;
/*******************************
* loop for each sub-sample period
********************************/
for (int sslc=0;sslc<numPeriods;++sslc)
{
/*****************************************************************************--
% 1. Get the structural parameters & define penalties
******************************************************************************-*/
cost_flag = 1;
int i;
if (mode_compute != 1)
{
// vector<int>kk(0);
double qdelta=0;
for(i=0;i<xparam1.length();++i)
{
if(xparam1[i]<param_lb[i]) // kk.push_back[i+1];
qdelta+=(xparam1[i]-param_lb[i])*(xparam1[i]-param_lb[i]);
}
if ( qdelta>0) // i.e. kk.size()>0)
{
// fval = bayestopt_.penalty+sum((bayestopt_.lb(k)-xparam1(k)).^2);
likelihood = penalty+qdelta;
cost_flag = 0;
info = 41;
return likelihood;
}
qdelta=0;
// kk.clear();
for(i=0;i<xparam1.length();++i)
{
if(xparam1[i]>param_ub[i]) // kk.push_back[i+1];
qdelta+=(xparam1[i]-param_ub[i])*(xparam1[i]-param_ub[i]);
}
if ( qdelta>0) // i.e. kk.size()>0)
{
//fval = bayestopt_.penalty+sum((xparam1(k)-bayestopt_.ub(k)).^2);
likelihood = penalty+qdelta;
cost_flag = 0;
info = 42;
return likelihood;
}
} // mode compute
#ifdef DEBUG
mexPrintf("Calling of updataeHQparams\n");
#endif
if(info=updateQHparams()) // updates Q and H matrices and deep parameters
{
#ifdef DEBUG
mexPrintf("Failing of updataeHQparams info =%d\n", info);
#endif
return likelihood;
}
/*****************************************************************************--
% 2. call model setup & reduction program and pre-filter data
// dynare_resolve(() // ... comes here doing:
// resol:
// check if ys is steady state and calculate one i not
// dr
// kalman_transition_matrix(out: A,B, in dr)
// IN: bayestopt_.restrict_var_list, bayestopt_.restrict_columns, bayestopt_.restrict_aux, )
***************************************************************/
#ifdef DEBUG
mexPrintf(" *********** Calling dynareResolveDR *********** \n");
#endif
if (info = dynareResolveDR (iv,ic, aux)) //OUT: [T,R,SteadyState],
{
#ifdef DEBUG
mexPrintf("Failing of dynareResolveDR info =%d\n", info);
#endif
return likelihood=penalty;
}
/*****************************************************************************--
% 2.b pre-filter and detrend data
******************************************************************************-*/
#ifdef DEBUG
mexPrintf("*********** pre-filter and detrend data *********** \n");
#endif
//if options_.noconstant
if ((int)dynareParams.getDoubleField(string("noconstant")))
constant.zeros();
else
{
//if options_.loglinear
if ((int)dynareParams.getDoubleField(string("loglinear")))
{
for (i =0;i<numVarobs;++i)
constant[i] = log(SteadyState[mfys[i]-1]);
}
else
{
for (i =0;i<numVarobs;++i)
constant[i] = SteadyState[mfys[i]-1];
}
}
Vector trend_coeff(numVarobs);
//trend = repmat(constant,1,gend);
GeneralMatrix constMx(constant.base(),numVarobs,1);
#ifdef DEBUG
mexPrintf("Calling constMx.repmat numTimeObs=%d\n",numTimeObs);
#endif
GeneralMatrix&trend = constMx.repmat(1,numTimeObs);
//if bayestopt_.with_trend
/************************************
if ((int)dynareParams.getDoubleField(string("with_trend")))
{
trend_coeff.zeros();
// GeneralMatrix& mt = dynareParams.getMatrixField(string("trend_coeffs"));
// Vector vt(mt.base, MAX(mt.numCols(), mt.numRows()));
Vector& vtc = dynareParams.getDoubleVectorField(string("trend_coeffs"));
for (i=0;i<vtc.length();++i)
{
if (vtc[i]!=0.0)
trend_coeff[i] = vtc[i];
}
//trend = repmat(constant,1,gend)+trend_coeff*[1:gend];
GeneralMatrix trend_coefMx(numVarobs, numTimeObs);
for (i=1;i<=numTimeObs;++i)
for (int j=0;j<numVarobs;++j)
trend_coefMx.get(j,i)=trend_coeff[j]*i;
trend.add(1,trend_coefMx);
}
*************************************/
presampleStart =(int) dynareParams.getDoubleField(string("presample"))+1;
int no_missing_data_flag = (number_of_observations==numTimeObs*numVarobs);
//Y =data-trend;
Y=data;
#ifdef DEBUG
mexPrintf("Calling Y.add( trend) in GeneralMatrices\n");
#endif
Y.add(-1,trend);
/*****************************************************************************
% 3. Initial condition of the Kalman filter
*******************************************************************************/
#ifdef DEBUG
mexPrintf("Calling InitiateKalmanMatrices\n");
#endif
if( InitiateKalmanMatrices())
return likelihood=penalty;
/*****************************************************************************
// 4. Likelihood evaluation
// choose and run KF to get likelihood fval
*****************************************************************************/
likelihood+=KalmanFilter(0.000001, false);// calls Kalman
/****************************************************************************
// Adds prior if necessary
****************************************************************************/
//likelihood-= priordens(xparam1,pshape,p6,p7,p3,p4);//fval = (likelihood-lnprior);
//options_.kalman_algo = kalman_algo;
} // end sub-sample loop
#ifdef DEBUG
mexPrintf("End of CallcLlikelihood returning likelihood=%f\n", likelihood);
#endif
return likelihood;
}
/**************************************************
* lower level, private Member functions definitions
***************************************************/
/*****************************************************************************--
% 1. Get the structural parameters & define penalties
******************************************************************************-*/
int
DsgeLikelihood::updateQHparams()// updates Q and H matrices and deep parameters
{
int i=0, offset=0, nv=0, k, k1, k2, info=0;
#ifdef DEBUG
mexPrintf("Setting of Q \n");
#endif
delete &Q;
Q = dynareParams.getMatrixField(string("Sigma_e"));
#ifdef DEBUG
mexPrintf("Setting of H \n");
#endif
delete &H;
H = dynareParams.getMatrixField(string("H"));
nv=(int)dynareParams.getDoubleField(string("nvx"));
if(nv)
{
#ifdef DEBUG
mexPrintf("Setting of Q var_exo\n");
#endif
// if(&estvx) delete &estvx;
GeneralMatrix&estvx=dynareParams.getMatrixField(string("var_exo"));
for (i=0;i<nv;++i)
{
k =(int)estvx.get(i,0)-1;
Q.get(k,k) = xparam1[i]*xparam1[i];
}
offset = nv;
delete &estvx;
}
nv=(int)dynareParams.getDoubleField(string("nvn"));
if(nv)
{
#ifdef DEBUG
mexPrintf("Setting of H var_endo\n");
#endif
GeneralMatrix&estvn=dynareParams.getMatrixField(string("var_endo"));
for (i=0;i<nv;++i)
{
k =(int)estvn.get(i,0)-1;
H.get(k,k) = xparam1[i+offset]*xparam1[i+offset];
}
offset += nv;
delete &estvn;
}
//if estim_params_.ncx
//for i=1:estim_params_.ncx
nv=(int)dynareParams.getDoubleField(string("ncx"));
if(nv)
{
#ifdef DEBUG
mexPrintf("Setting of Q corrx\n");
#endif
GeneralMatrix&corrx=dynareParams.getMatrixField(string("corrx"));
for (i=0;i<nv;++i)
{
k1 =(int)corrx.get(i,0)-1;
k2 =(int)corrx.get(i,1)-1;
Q.get(k1,k2) = xparam1[i+offset]*sqrt(Q.get(k1,k1)*Q.get(k2,k2));
Q.get(k2,k1) = Q.get(k1,k2);
}
// [CholQ,testQ] = chol(Q);
delete &corrx;
int testQ=0;
try
{
NormCholesky chol(Q);
}
catch(const TSException &e)
{
// if (string(e.getMessage())==sting("The matrix is not positive definite in NormCholesky constructor"))
if (0==strncmp(e.getMessage(),"The matrix is not positive definite in NormCholesky constructor",35))
testQ=1;
else
{
printf("Caugth unhandled TS exception with Q matrix: ");
likelihood=penalty;
TS_RAISE(e.getMessage());
}
}
if (testQ)
{
// The variance-covariance matrix of the structural innovations is not definite positive.
// We have to compute the eigenvalues of this matrix in order to build the penalty.
double delta=0;
VDVFact eigQ(Q); // get eigenvalues
//k = find(a < 0);
if(eigQ.hasConverged())
{
const Vector& evQ=eigQ.getD();
for (i=0;i<evQ.length();++i)
if (evQ[i]<0)
delta-=evQ[i];
}
likelihood = penalty+delta;// +sum(-a(k));
cost_flag = 0;
info = 43;
return info;
}
//offset = offset+estim_params_.ncx;
offset += nv;
}//end
//if estim_params_.ncn
//for i=1:estim_params_.ncn
nv=(int)dynareParams.getDoubleField(string("ncn"));
if(nv)
{
#ifdef DEBUG
mexPrintf("Setting of H corrn\n");
#endif
GeneralMatrix&corrn=dynareParams.getMatrixField(string("corrn"));
vector<int>&lgyidx2varobs= dynareParams.getIntVectorField(string("lgyidx2varobs"));
for (i=0;i<nv;++i)
{
// k1 = options_.lgyidx2varobs(estim_params_.corrn(i,1));
// k2 = options_.lgyidx2varobs(estim_params_.corrn(i,2));
k1 = lgyidx2varobs[(int)corrn.get(i,0)-1];
k2 = lgyidx2varobs[(int)corrn.get(i,1)-1];
// H(k1,k2) = xparam1(i+offset)*sqrt(H(k1,k1)*H(k2,k2));
// H(k2,k1) = H(k1,k2);
H.get(k1,k2) = xparam1[i+offset]*sqrt(H.get(k1,k1)*H.get(k2,k2));
H.get(k2,k1) = H.get(k1,k2);
}
delete &corrn;
//[CholH,testH] = chol(H);
int testH=0;
try
{
NormCholesky chol(H);
}
catch(const TSException &e)
{
// if (string(e.getMessage)==sting("The matrix is not positive definite in NormCholesky constructor"))
if (0==strncmp(e.getMessage(),"The matrix is not positive definite in NormCholesky constructor",35))
testH=1;
else
{
printf("Caugth unhandled TS exception with H matrix: ");
likelihood=penalty;
TS_RAISE((const char*)e.getMessage());
}
}
if (testH)
{
//a = diag(eig(H));
double delta=0;
VDVFact eigH(H); // get eigenvalues
//k = find(a < 0);
if(eigH.hasConverged())
{
const Vector& evH=eigH.getD();
for (i=0;i<evH.length();++i)
if (evH[i]<0)
delta-=evH[i];
}
likelihood = penalty+delta; // +sum(-a(k));
cost_flag = 0;
info = 44;
return info;
}; // end if
//offset = offset+estim_params_.ncn;
offset += nv;
}
//if estim_params_.np > 0 // i.e. num of deep parameters >0
// M_.params(estim_params_.param_vals(:,1)) = xparam1(offset+1:end);
if(num_dp > 0)
{
// if(xparam1.length()>=offset+num_dp)
// memcpy(deepParams.base(), xparam1.base()+offset*sizeof(double),num_dp*sizeof(double));
if(xparam1.length()>=offset+num_dp)
deepParams=Vector(xparam1,offset,num_dp);
else
TS_RAISE("Inssuficient length of the xparam1 parameters vector");
}
#ifdef DEBUG
mexPrintf("End of Setting of HQ params\n");
#endif
/**********
M_.Sigma_e = Q;
M_.H = H;
******************/
return info;
};
/*****************************************************************************--
% 3. Initial condition of the Kalman filter
******************************************************************************-*/
int
DsgeLikelihood::InitiateKalmanMatrices()
{
int np = T.numRows();// size(T,1);
double lyapunov_tol=dynareParams.getDoubleField(string("lyapunov_tol"));
// int lik_init=(int)dynareParams.getDoubleField(string("lik_init"));
//if options_.lik_init == 1 % Kalman filter
// GeneralMatrix RQRt(R,Q); // R*Q
GeneralMatrix RQ(R,Q); // R*Q
#ifdef DEBUG
mexPrintf("Calling RQRt.multRightTrans(R)\n");
#endif
// GeneralMatrix::md_length=RQRt.numRows();
//RQRt.md_length=RQRt.numRows();
// RQRt.multRightTrans(R); // R*Q*Rt
// GeneralMatrix RQRt(np,np);
// RQRt.zeros();
// RQRt.multAndAdd(RQ,R, "T", 1.0);
GeneralMatrix RQRt(RQ,R, "T");
#ifdef DEBUG
for (int i=0;i<RQRt.numRows();++i)
{
for (int j=0;j<RQRt.numCols();++j)
mexPrintf(" %f ", RQRt.get(i,j));
mexPrintf("\n");
}
#endif
GeneralMatrix Ptmp(np,np);
Ptmp.zeros();
//Pstar = lyapunov_symm (T,R*Q*R',options_.qz_criterium,options_.lyapunov_complex_threshold)
#ifdef DEBUG
mexPrintf("Calling disclyap_fast to initialise Pstar:\n");
#endif
try
{
disclyap_fast(T,RQRt,Ptmp, lyapunov_tol, 0); // 1 to check chol
Pstar=Ptmp;
#ifdef DEBUG
Pstar.print();
mexPrintf("Initialise Pinf\n");
#endif
//Pinf=[]
Pinf.zeros();
}
catch(const TSException &e)
{
if (0==strncmp(e.getMessage(),"The matrix is not positive definite in NormCholesky constructor",35))
{
printf(e.getMessage());
#ifdef MATLAB
mexPrintf(e.getMessage());
#endif
likelihood=penalty;
return 1;
}
else
{
printf("Caugth unhandled TS exception with H matrix: ");
#ifdef MATLAB
mexPrintf("Caugth unhandled TS exception with H matrix: ");
#endif
likelihood=penalty;
TS_RAISE((const char*)e.getMessage());
}
}
//a=zeros(size(T,1),1);
a_init.zeros();
//if (lik_init == 2)// Old Diffuse Kalman filter
// Pstar = options_.Harvey_scale_factor*eye(np);
//Pinf = [];
//else if (lik_init == 3) // Diffuse Kalman filter
// else ...
return 0;
}
/*****************************************************************************
// 4. Likelihood evaluation
// choose and run KF to get likelihood fval
*****************************************************************************/
double
DsgeLikelihood::KalmanFilter(double riccatiTol=0.000001,bool uni = false)
{
bool diffuse=false;
double loglik;
try
{
// make input matrices
int start = presampleStart;
int nper=Y.numCols();
#ifdef DEBUG
mexPrintf("kalman_filter: periods=%d start=%d, a.length=%d, uni=%d diffuse=%d\n", nper, start,a_init.length(), uni, diffuse);
Pstar.print(); Z.print(); H.print(); T.print(); R.print(); Q.print(); Y.print();
#endif
// make storage for output
int per;
int d;
// create state init
StateInit* init = NULL;
if (diffuse||uni)
{
if (diffuse)
{
init = new StateInit(Pstar, Pinf, a_init);
}
else
{
init = new StateInit(Pstar, a_init);
}
// fork, create objects and do filtering
KalmanTask kt(Y, Z, H, T, R, Q, *init);
if (uni)
{
KalmanUniTask kut(kt);
loglik = kut.filter(per, d, (start-1), vll);
per = per / Y.numRows();
d = d / Y.numRows();
}
else
{
#ifdef KF_TIMING_LOOP
mexPrintf("kalman_filter: starting 1000 loops\n");
for (int tt=0;tt<1000;++tt)
{
#endif
loglik = kt.filter(per, d, (start-1), vll);
#ifdef KF_TIMING_LOOP
}
mexPrintf("kalman_filter: finished 1000 loops\n");
#endif
}
}
else // basic Kalman
{
init = new StateInit(Pstar, a_init);
BasicKalmanTask bkt(Y, Z, H, T, R, Q, *init, riccatiTol);
#ifdef KF_TIMING_LOOP
mexPrintf("kalman_filter: starting 1000 loops\n");
for (int tt=0;tt<1000;++tt)
{
#endif
loglik = bkt.filter( per, d, (start-1), vll);
#ifdef DEBUG
mexPrintf("Basickalman_filter: loglik=%f \n", loglik);
#endif
#ifdef KF_TIMING_LOOP
}
mexPrintf("Basickalman_filter: finished 1000 loops\n");
#endif
}
// destroy init
delete init;
}
catch (const TSException& e)
{
mexErrMsgTxt(e.getMessage());
}
catch (SylvException& e)
{
char mes[300];
e.printMessage(mes, 299);
mexErrMsgTxt(mes);
}
return loglik;
}

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

@ -1,144 +0,0 @@
/*
* 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/>.
*/
#include "Estimation.h"
#include "k_ord_dynare.h"
#include "kalman.h"
#include "math.h"
#include "disclyap_fast.h"
#include <cstring>
#include <cctype>
class DsgeLikelihood
{
double likelihood; // sum of vector of KF step log likelihoods
vector<double>* vll; // vector of KF step log likelihoods
Vector& a_init;//initial value of the state, usually set to 0.
GeneralMatrix& Q;// Kalman Matrices
GeneralMatrix& R;
GeneralMatrix& T;
GeneralMatrix& Z;
GeneralMatrix& Pstar;
GeneralMatrix& Pinf;
GeneralMatrix& H;
const GeneralMatrix& data;
GeneralMatrix& Y;
//GeneralMatrix& currentDataSubSample;
//int periodStart, periodEnd; // start and end of current sub sample
const int numPeriods;//=1; number of structural change periods
const int numVarobs; // number of observed variables in the observation vector at time t.
const int numTimeObs; // number of obsevations (vectors) in the time series
const int order;
const int endo_nbr;
const int exo_nbr;
const int nstatic;
const int npred;
const int nfwrd;
char *fName;
int presampleStart;
int kalman_algo; // type of kalman algorithm: multi- or uni-variate
int mode_compute;
int info;
double penalty;
int cost_flag;
const int number_of_observations;
const bool no_more_missing_observations;
const vector<int>& order_var;
const vector<int>& mfys;
const vector<int>& mf; // positions of observed variables in restricted state vector for likelihood computation.
Vector& xparam1; // all estimated parameters incl sderrs
const int num_dp; // number of deep parameters
Vector& deepParams; // estimated deep parameters subset of xparam1 only
const Vector& param_ub; // upper and lower bounds
const Vector& param_lb;
const vector<int>&pshape;
const Vector& p6;
const Vector& p7;
const Vector& p3;
const Vector& p4;
Vector& SteadyState;
Vector& constant;
GeneralParams& dynareParams;
//GeneralParams& parameterDescription;
GeneralParams& dr;
GeneralMatrix& aux; //= dynareParams.getMatrixField(string("restrict_aux"));
vector<int>&iv; //= dynareParams.getIntVectorField(string("restrict_var_list"));
vector<int>&ic; //= dynareParams.getIntVectorField(string("restrict_columns"));
GeneralMatrix& kstate;
GeneralMatrix& ghx;
GeneralMatrix& ghu;
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;
//KalmanUniTask ukt;// univariate
// member functions
MexStructParam& SetDRModel(MexStructParam&params);
// void disclyap_fast(const GeneralMatrix &G, const GeneralMatrix & V, GeneralMatrix &X, double tol = 1e-16, int flag_ch=0);
GeneralMatrix& SolveDRkOrderPert();//calls k-order pert or whatever;
int dynareResolveDR(vector<int>&iv,vector<int>&ic,GeneralMatrix& aux); // returns int info, ys, and TT and RR Decision Rule
int SolveDRModel(const int endo_nbr, const int exo_nbr, const int nstatic, const int npred, int nfwrd);//int dr1(); // returns int info and updated dr
int updateQHparams();// updates Q and H matrices and deep parameters
int InitiateKalmanMatrices();
void DataPreparation(MexStructParam&params, const GeneralMatrix &data);
double KalmanFilter(double riccatiTol,bool uni);// calls Kalman
public:
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,
const int INorder, const int INendo_nbr, const int INexo_nbr, const int INnstatic,
const int INnpred, const int INnfwrd, const int INnum_of_observations, const bool INno_more_missing_observations,
const vector<int>& INorder_var, const vector<int>& INmfys, const vector<int>& INmf,
Vector& INxparam1, const int INnum_dp, Vector& INdeepParams,
const Vector& INub, const Vector& INlb, const vector<int>&INpshape,
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,
GeneralMatrix& aux, vector<int>&iv, vector<int>&ic
,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,
KordpDynare& model); // constructor
~DsgeLikelihood();// destructor
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);
};

57
mex/sources/estimation/EstimatedParameter.hh Normal file
View File

@ -0,0 +1,57 @@
/*
* Copyright (C) 2009-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/>.
*/
///////////////////////////////////////////////////////////
// EstimatedParameter.h
// Implementation of the Class EstimatedParameter
// Created on: 02-Feb-2010 13:06:35
///////////////////////////////////////////////////////////
#if !defined(D879C8AE_5B69_4fc3_83BD_FA5A99030ECF__INCLUDED_)
#define D879C8AE_5B69_4fc3_83BD_FA5A99030ECF__INCLUDED_
#include "Prior.hh"
#include <cstdlib>
struct EstimatedParameter
{
public:
EstimatedParameter(
size_t ID1, size_t ID2, size_t type, size_t startPeriod, size_t endPeriod,
double lower_boound, double upper_boound, Prior prior
);
virtual ~EstimatedParameter();
public:
const size_t ID1;
const size_t ID2;
enum pType
{
shock_SD = 1, // standard deviation of a structural shock
measureErr_SD = 2, // standard deviation of a measurement error
shock_Corr = 3, // correlation betwwen two structural shocks
measureErr_Corr = 4, // correlation between two measurement errors
deepPar = 5 // deep parameter
};
const pType ptype;
const double lower_boound;
const double upper_boound;
const Prior prior;
const std::vector<size_t> subSampleIDs;
};
#endif // !defined(D879C8AE_5B69_4fc3_83BD_FA5A99030ECF__INCLUDED_)

66
mex/sources/estimation/EstimatedParametersDescription.hh Normal file
View File

@ -0,0 +1,66 @@
/*
* Copyright (C) 2009-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/>.
*/
///////////////////////////////////////////////////////////
// EstimatedParametersDescription.h
// Implementation of the Class EstimatedParametersDescription
// Created on: 02-Feb-2010 13:06:46
///////////////////////////////////////////////////////////
#if !defined(E8F2C096_A301_42e8_80FF_A643291BF995__INCLUDED_)
#define E8F2C096_A301_42e8_80FF_A643291BF995__INCLUDED_
#include "EstimationSubsample.hh"
#include "EstimatedParameter.hh"
#include "Prior.hh"
/**
* all estimation periods held in the vectors of EstPeriod-s are held in the
* EstimatedParametersDescription together with the time-invariant attributes (e.g.
* IDs..) for integration..
*
* The EstimatedParametersDescription. structure (or class with protected
* elements) integrates all parametr and time period information and makes it
* available to (friend ?) class LogPosteriorDensity on as needed basis and
* updates all parameters (inc. H and Q) when time slot xparam1 is supplied.
*
*
*
*/
class EstimatedParametersDescription
{
public:
virtual ~EstimatedParametersDescription();
public:
EstimatedParametersDescription(size_t ncx, size_t ncn, size_t nvx, size_t nvn, size_t np,
std::vector<EstimationSubsample> &estSubsamples, std::vector<EstimatedParameter> &estParams);
//Vector getAllParams(int period);
//int getNumPeriods();
EstimationSubsample&getSubsample(int period);
std::vector<EstimationSubsample>&getAllSubsamples();
//EstimatedParameter& getParameter(size_t ID);
//private:
std::vector<EstimationSubsample> estSubsamples;
std::vector<EstimatedParameter> estParams;
//vector<Prior> prior;
};
#endif // !defined(E8F2C096_A301_42e8_80FF_A643291BF995__INCLUDED_)

72
mex/sources/estimation/EstimationSubsample.hh Normal file
View File

@ -0,0 +1,72 @@
/*
* Copyright (C) 2009-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/>.
*/
///////////////////////////////////////////////////////////
// EstimationSubsample.h
// Implementation of the Class EstimationSubsample
// Created on: 02-Feb-2010 13:06:01
///////////////////////////////////////////////////////////
#if !defined(EstSub34FB1351_4066_4993_B384_56AA65BC32C4__INCLUDED_)
#define EstSub34FB1351_4066_4993_B384_56AA65BC32C4__INCLUDED_
#include <cstdlib>
#include <vector>
/**
* Contains longest common periods for different parameters with different start
* and end periods:
*
* It defines start and end of that period.which has start higher/equal than any
* individual parameter start and its end lower|equal than any individual
* parameter end in that perod window.) integrated xparam1 and priors for that
* period
*
* All those EstPeriod-s need to be constructed at some point once at a start of
* estimation by EstPeriod costructor from the information such as invidual
* parameter start-end periods and the relevant associated priors for those
* periods. That initial info-set is held in the vector of vectors of
* EstParamSubSampleDescription
*
* it constructs and contains indices of all estimated xparam1 parameters valid
* for this period in the larger, extended xparam1x vector passed in. from the
* gradient function.
*
* The EstimatedParametersDescription. structure (or class with protected
* elements) integrates all that information and makes it available to (friend ?)
* class LogPosteriorDensity on as needed basis and updates all parameters (inc. H
* and Q) when time slot xparam1 is supplied.
*
*
*/
class EstimationSubsample {
public:
EstimationSubsample(size_t startPeriod, size_t endPeriod);
virtual ~EstimationSubsample();
const size_t startPeriod;
const size_t endPeriod;
/**
* indices of all estimated xparam1 parameters valid for this period in the larger,
* extended xparam1x vector passed in. from the gradient function will be added at a later stage of development
*/
// const std::vector<size_t> subSampleParamIDs;
};
#endif // !defined(34FB1351_4066_4993_B384_56AA65BC32C4__INCLUDED_)

4
mex/sources/estimation/InitializeKalmanFilter.cc
View File

@ -33,11 +33,11 @@ InitializeKalmanFilter::~InitializeKalmanFilter()
InitializeKalmanFilter::InitializeKalmanFilter(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 qz_criterium, const std::vector<size_t> &order_var_arg, const std::vector<size_t> &inv_order_var_arg,
double qz_criterium, const std::vector<size_t> &order_var_arg, const std::vector<size_t> &inv_order_var_arg,
const std::vector<size_t> &rivIN, const std::vector<size_t> &ricIN, double lyapunov_tolIN, int &info) :
riv(rivIN), inv_ric(ricIN.size()), order_var(order_var_arg), lyapunov_tol(lyapunov_tolIN),
detrendData(false), modelSolution(modName, n_endo_arg, n_exo_arg, zeta_fwrd_arg, zeta_back_arg,
zeta_mixed_arg, zeta_static_arg, llincidence, qz_criterium), discLyapFast(riv.size()),
zeta_mixed_arg, zeta_static_arg, qz_criterium), discLyapFast(riv.size()),
ghx(n_endo_arg, zeta_back_arg.size() + zeta_mixed_arg.size()),
ghx_raw(n_endo_arg, zeta_back_arg.size() + zeta_mixed_arg.size()),
ghu(n_endo_arg, n_exo_arg), ghu_raw(n_endo_arg, n_exo_arg),

2
mex/sources/estimation/InitializeKalmanFilter.hh
View File

@ -42,7 +42,7 @@ class InitializeKalmanFilter
public:
InitializeKalmanFilter(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, const std::vector<size_t> &order_var_arg, const std::vector<size_t> &inv_order_var_arg,
double qz_criterium, const std::vector<size_t> &order_var_arg, const std::vector<size_t> &inv_order_var_arg,
const std::vector<size_t> &riv, const std::vector<size_t> &ric, double lyapunov_tol, int &info);
virtual ~InitializeKalmanFilter();
// initialise all KF matrices

4
mex/sources/estimation/KalmanFilter.cc
View File

@ -33,7 +33,7 @@ KalmanFilter::~KalmanFilter()
KalmanFilter::KalmanFilter(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 &ll_incidence, double qz_criterium, const std::vector<size_t> &order_var,
double qz_criterium, const std::vector<size_t> &order_var,
const std::vector<size_t> &inv_order_var, const std::vector<size_t> &varobs_arg,
const std::vector<size_t> &riv, const std::vector<size_t> &ric,
double riccati_tol_in, double lyapunov_tol, int &info) :
@ -44,7 +44,7 @@ KalmanFilter::KalmanFilter(const std::string &modName, size_t n_endo, size_t n_e
oldKFinv(riv.size(), varobs_arg.size()), Finverter(varobs_arg.size()), a_init(riv.size(), 1),
a_new(riv.size(), 1), vt(varobs_arg.size(), 1), vtFinv(1, varobs_arg.size()), vtFinvVt(1), riccati_tol(riccati_tol_in),
initKalmanFilter(modName, n_endo, n_exo, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg,
zeta_static_arg, ll_incidence, qz_criterium, order_var, inv_order_var, riv, ric, lyapunov_tol, info)
zeta_static_arg, qz_criterium, order_var, inv_order_var, riv, ric, lyapunov_tol, info)
{
a_init.setAll(0.0);
Z.setAll(0.0);

2
mex/sources/estimation/KalmanFilter.hh
View File

@ -50,7 +50,7 @@ public:
virtual ~KalmanFilter();
KalmanFilter(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, const std::vector<size_t> &order_var_arg, const std::vector<size_t> &inv_order_var,
double qz_criterium, const std::vector<size_t> &order_var_arg, const std::vector<size_t> &inv_order_var,
const std::vector<size_t> &varobs_arg, const std::vector<size_t> &riv, const std::vector<size_t> &ric,
double riccati_tol, double lyapunov_tol, int &info);

70
mex/sources/estimation/LogLikelihoodMain.hh Normal file
View File

@ -0,0 +1,70 @@
/*
* Copyright (C) 2009-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/>.
*/
///////////////////////////////////////////////////////////
// LogLikelihoodMain.h
// Implementation of the Class LogLikelihoodMain
// Created on: 02-Feb-2010 12:57:09
///////////////////////////////////////////////////////////
#if !defined(E126AEF5_AC28_400a_821A_3BCFD1BC4C22__INCLUDED_)
#define E126AEF5_AC28_400a_821A_3BCFD1BC4C22__INCLUDED_
//#include "EstimatedParametersDescription.hh"
#include "LogLikelihoodSubSample.hh"
class LogLikelihoodMain {
private:
double logLikelihood;
std::vector<EstimationSubsample> &estSubsamples; // reference to member of EstimatedParametersDescription
LogLikelihoodSubSample logLikelihoodSubSample;
Vector deepParams;
Vector &vll; // vector of all KF step likelihoods
Matrix &data; // input data
Matrix &steadyState;
//GeneralParams& estimOptions;
int presampleStart;
Matrix Q;
Matrix H;
public:
virtual ~LogLikelihoodMain();
LogLikelihoodMain(const Matrix &data, //GeneralParams& estimOptions,
const std::string &modName, EstimatedParametersDescription &estiParDesc, 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 double qz_criterium, const std::vector<size_t> &order_var,
const std::vector<size_t> &inv_order_var, const std::vector<size_t> &varobs, const std::vector<size_t> &riv,
const std::vector<size_t> &ric, double riccati_tol_in, double lyapunov_tol, int &info);
double compute(Matrix &steadyState, Vector &estParams, int &info); // for calls from estimation and to set Steady State
double compute(Vector &estParams); // for calls in loop from optimiser
Vector &
getVll()
{
return vll;
};
double
getLogLikelihood()
{
return logLikelihood;
};
};
#endif // !defined(E126AEF5_AC28_400a_821A_3BCFD1BC4C22__INCLUDED_)

4
mex/sources/estimation/LogLikelihoodSubSample.cc
View File

@ -31,10 +31,10 @@ LogLikelihoodSubSample::~LogLikelihoodSubSample()
LogLikelihoodSubSample::LogLikelihoodSubSample(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 &ll_incidence, const double qz_criterium, const std::vector<size_t> &order_var,
const std::vector<size_t> &zeta_static_arg, const double qz_criterium, const std::vector<size_t> &order_var,
const std::vector<size_t> &inv_order_var, const std::vector<size_t> &varobs, const std::vector<size_t> &riv,
const std::vector<size_t> &ric, double riccati_tol, double lyapunov_tol, int &info) :
kalmanFilter(modName, n_endo, n_exo, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, ll_incidence, qz_criterium,
kalmanFilter(modName, n_endo, n_exo, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, qz_criterium,
order_var, inv_order_var, varobs, riv, ric, riccati_tol, lyapunov_tol, info)
{
};

22
mex/sources/estimation/LogLikelihoodSubSample.hh
View File

@ -27,26 +27,32 @@
#define DF8B7AF5_8169_4587_9037_2CD2C82E2DDF__INCLUDED_
#include "KalmanFilter.hh"
#include "EstimatedParametersDescription.hh"
#include "VDVEigDecomposition.hh"
class LogLikelihoodSubSample {
public:
LogLikelihoodSubSample(const std::string &modName, size_t n_endo, size_t n_exo,
LogLikelihoodSubSample(const std::string &modName, EstimatedParametersDescription &estiParDesc, 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 &ll_incidence, const double qz_criterium, const std::vector<size_t> &order_var,
const std::vector<size_t> &zeta_static_arg, const double qz_criterium, const std::vector<size_t> &order_var,
const std::vector<size_t> &inv_order_var, const std::vector<size_t> &varobs, const std::vector<size_t> &riv,
const std::vector<size_t> &ric, double riccati_tol_in, double lyapunov_tol, int &info);
double compute(Vector &steadyState, const MatrixView &dataView, const Vector &deepParams, //const estPeriod &estPeriod,
VectorView &vll, int &info, size_t start, size_t period, const Matrix &Q, const Matrix &H);
double compute(Vector &steadyState, const MatrixView &dataView, Vector &deepParams,
Matrix &Q, Matrix &H, VectorView &vll, int &info, size_t start, size_t period);
virtual ~LogLikelihoodSubSample();
public:
KalmanFilter kalmanFilter;
private:
double penalty;
double logLikelihood;
double logLikelihood;
EstimatedParametersDescription &estiParDesc;
KalmanFilter kalmanFilter;
VDVEigDecomposition eigQ;
VDVEigDecomposition eigH;
// methods
void updateParams(const Vector &estParams, Vector &deepParams,
Matrix &Q, Matrix &H, size_t period);
};

573
mex/sources/estimation/MexUtils.cpp
View File

@ -1,573 +0,0 @@
/*
* 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/>.
*/
#include "mexutils.h"
MexStruct::MexStruct( const int numparstruct):
numParStruct(numparstruct), structName(string(""))
{
// get Dynare mexSturcture pointers and store them locally
#ifdef DEBUG
mexPrintf("MexStruct reserve=%d \n", numParStruct);
#endif
parStruct.reserve(numParStruct);
parStructBase.reserve(numParStruct);
for (int i=0;i<numParStruct;++i)
{
parStruct[i]=mexGetVariablePtr("caller", DynareParamStructsNm[i]);
parStructBase[i]=caller;
#ifdef DEBUG
mexPrintf("MexStruct to insert i=%d parStructNm[i]=%s using base[i]=%d %s\n", i, DynareParamStructsNm[i], parStructBase[i], mexBase[parStructBase[i]]);
#endif
// get field names into the map:
pair <map<string, int>::iterator, bool> ret;
int j=0;
const char* field;
while ((field=mxGetFieldNameByNumber(parStruct[i],j))!=NULL)
{
#ifdef DEBUG
mexPrintf("MexStruct insert field= %s\n", field);
#endif
ret=parNamStructMap.insert(make_pair(string(field),i));
if (!ret.second)
mexPrintf("MexStruct failed to insert field %s from struct %d %s using base %d %s\n"
, field, i, DynareParamStructsNm[i], parStructBase[i], mexBase[parStructBase[i]]);
// mexErrMsgTxt("MexStruct: Failed to insert param \n");
j++;
}
}
#ifdef DEBUG
mexPrintf("MexStruct insert finished \n");
#endif
}
MexStructParam&
MexStruct::getMexStructField( const string& field)
{
mxArray* mxf = getMxField( field);
return *(new MexStructParam(mxf, this, field));
}
void *
MexStruct::getField( const string& field)
{
mxArray* mxf = getMxField( field);
return mxGetData(mxf);
}
GeneralMatrix &
MexStruct::getMatrixField( const string& field)
{
mxArray* mxf = getMxField( field);
return *(new GeneralMatrix(mxGetPr(mxf), mxGetM(mxf), mxGetN(mxf))) ;
}
double&
MexStruct::getDoubleField(const string& field)
{
const mxArray* mxf = getMxField( field);
if (!mxIsDouble(mxf))
mexErrMsgTxt("Input must be of type double.");
// double ret=mxGetScalar(mxf);
// return ret;
double *ret=(new double);
*ret=mxGetScalar(mxf);
return *ret;
}
string&
MexStruct::getStringField(const string& field)
{
mxArray* mxf = getMxField( field);
if (!mxIsChar(mxf))
mexErrMsgTxt("Input must be of type char.");
return *(new string(mxArrayToString(mxf)));
}
char *
MexStruct::getCharField( const string& field)
{
mxArray* mxf = getMxField( field);
if (!mxIsChar(mxf))
mexErrMsgTxt("Input must be of type char.");
return mxArrayToString(mxf);
}
vector<string>&
MexStruct::getStringVectorField(const string& field)
{
charArraySt * cap = getCharArrayField(field);
vector <string>*sv=(new vector<string>(cap->len));
for (int i= 0; i<cap->len;++i)
(*sv)[i]=string(cap->charArrayPtr[i]);
return *sv;
}
vector<int>&
MexStruct::getIntVectorField(const string& field)
{
mxArray* mxfp = getMxField( field);
if (MIN(mxGetM(mxfp),mxGetN(mxfp))!=1)
throw SYLV_MES_EXCEPTION("Int vector is a 2D Matrix .");
double* dparams = (double *) mxGetData(mxfp);
int npar = (int) MAX(mxGetM(mxfp), mxGetN(mxfp));
vector<int> *vars = (new vector<int>(npar));
for (int v = 0; v < npar; v++)
{
(*vars)[v] = (int)(*(dparams++)); //[v];
#ifdef DEBUG
mexPrintf("%s[%d]=%d.\n", field.c_str() , v, (*vars)[v]);
#endif
}
return *vars;
};
Vector&
MexStruct::getDoubleVectorField(const string& field)
{
mxArray* mxfp = getMxField( field);
if (MIN(mxGetM(mxfp), mxGetN(mxfp))!=1)
throw SYLV_MES_EXCEPTION("Double Vector is a 2D Matrix .");
double* dparams = (double *) mxGetData(mxfp);
int npar = (int) MAX(mxGetM(mxfp), mxGetN(mxfp));
Vector *vars = (new Vector(dparams,npar));
return *vars;
};
charArraySt *
MexStruct::getCharArrayField( const string& field)
{
mxArray* mxfp = getMxField( field);
const int len = (int) mxGetM(mxfp);
const int width = (int) mxGetN(mxfp);
if (!mxIsChar(mxfp))
mexErrMsgTxt("Input must be of type char.");
charArraySt * cap=new charArraySt;
cap->charArrayPtr = (char**)MxArrayToStringArray(mxfp, len, width);
cap->len=len;
return cap;
}
void
MexStruct::ReorderCols(GeneralMatrix &tdx, const vector<int> &vOrder)
{MexStructParam::ReorderCols(tdx,vOrder, "static");};
void
MexStruct::setField( const string& field, const string& val)
{
mxArray *newVal=mxCreateString(val.c_str());
UpdateMexField( field, newVal);
}
/*******
{
mxArray *fp = mxGetField(mexStruct, 0, field);
mxDestroy(fp)
mxSetField(mexStruct,0,field, newVal);
}
***************/
void
MexStruct::setField( const string& field, const double val)
{
mxArray *newVal=mxCreateDoubleScalar(val);
UpdateMexField( field, newVal);
}
void
MexStruct::setField( const string& field, const GeneralMatrix& gmval)
{
mxArray *newVal=mxCreateDoubleMatrix(gmval.numRows(),gmval.numCols(),mxREAL);
memcpy(mxGetPr(newVal),gmval.base(),gmval.numRows()*gmval.numCols()*sizeof(double));
UpdateMexField( field, newVal);
}
void
MexStruct::UpdateMexField( const string& field, mxArray *newVal)
{
StructBaseNameMap* dpsm=getMxFieldStruct(field);
mxArray *sp =mexGetVariable(dpsm->baseName, dpsm->structName);
if (sp == NULL )
{
mexPrintf("Variable not found : base: %s, structName %s\n", dpsm->baseName, dpsm->structName);
mexErrMsgTxt(" \n");
}
mxArray *fp = mxGetField(sp, 0, field.c_str());
mxDestroyArray(fp);
mxSetField(sp,0,field.c_str(), newVal);
mexPutVariable(dpsm->baseName, dpsm->structName, sp);
//mxDestroyArray(newVal);
}
/***********************************************
* MexStructParam
* holds single Matlab structure passed as parameter
***********************************************/
MexStructParam&
MexStructParam::getMexStructField( const string& field)
{
mxArray* mxf = getMxField( field);
return *(new MexStructParam(mxf, this, field));
}
void *
MexStructParam::getField( const string& field)
{
mxArray* mxf = getMxField( field);
return mxGetData(mxf);
}
double&
MexStructParam::getDoubleField(const string& field)
{
const mxArray* mxf = getMxField( field);
if (!mxIsDouble(mxf))
mexErrMsgTxt("Input must be of type double.");
double *ret=(new double);
*ret=mxGetScalar(mxf);
return *ret;
}
GeneralMatrix&
MexStructParam::getMatrixField( const string& field)
{
mxArray* mxf = getMxField( field);
return *(new GeneralMatrix(mxGetPr(mxf), mxGetM(mxf), mxGetN(mxf)) );
}
string&
MexStructParam::getStringField(const string& field)
{
mxArray* mxf = getMxField( field);
if (!mxIsChar(mxf))
mexErrMsgTxt("Input must be of type char.");
return *(new string(mxArrayToString(mxf)));
}
char *
MexStructParam::getCharField( const string& field)
{
mxArray* mxf = getMxField( field);
if (!mxIsChar(mxf))
mexErrMsgTxt("Input must be of type char.");
return mxArrayToString(mxf);
}
vector<string>&
MexStructParam::getStringVectorField(const string& field)
{
charArraySt * cap = getCharArrayField( field);
vector <string>*sv= (new vector<string>(cap->len));
for (int i= 0; i<cap->len;++i)
(*sv)[i]=string(cap->charArrayPtr[i]);
return *sv;
}
vector<int>&
MexStructParam::getIntVectorField(const string& field)
{
mxArray* mxfp = getMxField( field);
if (MIN(mxGetM(mxfp), mxGetN(mxfp))!=1)
throw SYLV_MES_EXCEPTION("Int vector is a 2D Matrix .");
double* dparams = (double *) mxGetData(mxfp);
int npar = (int) MAX(mxGetM(mxfp), mxGetN(mxfp));
vector<int> *vars = (new vector<int>(npar));
for (int v = 0; v < npar; v++)
{
(*vars)[v] = (int)(*(dparams++)); //[v];
#ifdef DEBUG
mexPrintf("%s[%d]=%d.\n", field.c_str(), v, (*vars)[v]);
#endif
}
return *vars;
};
Vector&
MexStructParam::getDoubleVectorField(const string& field)
{
mxArray* mxfp = getMxField( field);
if (MIN(mxGetM(mxfp), mxGetN(mxfp))!=1)
throw SYLV_MES_EXCEPTION("Double Vector is a 2D Matrix .");
double* dparams = (double *) mxGetData(mxfp);
int npar = (int) MAX(mxGetM(mxfp), mxGetN(mxfp));
Vector *vars = (new Vector(dparams,npar));
return *vars;
};
charArraySt *
MexStructParam::getCharArrayField( const string& field)
{
mxArray* mxfp = getMxField( field);
const int len = (int) mxGetM(mxfp);
const int width = (int) mxGetN(mxfp);
if (!mxIsChar(mxfp))
mexErrMsgTxt("Input must be of type char.");
charArraySt * cap=new charArraySt;
cap->charArrayPtr = (char**)MxArrayToStringArray(mxfp, len, width);
cap->len=len;
return cap;
}
void
MexStructParam::setField( const string& field, const string& val)
{
mxArray *newVal=mxCreateString(val.c_str());
mxSetField((mxArray*)parStruct,0,field.c_str(), newVal);
//UpdateMexField( field, newVal);
}
/*******
{
mxArray *fp = mxGetField(mexStruct, 0, field);
mxDestroy(fp)
mxSetField(mexStruct,0,field, newVal);
}
***************/
void
MexStructParam::setField( const string& field, const double val)
{
mxArray *newVal=mxCreateDoubleScalar(val);
mxSetField((mxArray*)parStruct,0,field.c_str(), newVal);
//UpdateMexField( field, newVal);
}
void
MexStructParam::setField( const string& field, const GeneralMatrix& gmval)
{
mxArray *newVal=mxCreateDoubleMatrix(gmval.numRows(),gmval.numCols(),mxREAL);
memcpy(mxGetPr(newVal),gmval.base(),gmval.numRows()*gmval.numCols()*sizeof(double));
mxSetField((mxArray*)parStruct,0,field.c_str(), newVal);
//UpdateMexField( field, newVal);
}
void
MexStructParam::UpdateMexField( const string& field, mxArray *newVal)
{
mxSetField((mxArray*)parStruct,0,field.c_str(), newVal);
/************
if (parStructParent!=NULL)
parStructParent->setField(field,newVal);
else
{
mxArray *sp =mexGetVariable("caller",structName);
if (sp == NULL )
{
mexPrintf("Variable not found : base: %s, structName %s\n", dpsm->baseName, dpsm->structName);
mexErrMsgTxt(" \n");
}
mxArray *fp = mxGetField(sp, 0, field.c_str());
mxDestroyArray(fp);
mxSetField(sp,0,field.c_str(), newVal);
mexPutVariable(dpsm->baseName, dpsm->structName, sp);
//mxDestroyArray(newVal);
}
************/
}
/************************************
* Reorder first variables in a vector
* according to order given in varsOrder
************************************/
void
//MexStructParam::ReorderCols(GeneralMatrix &tdx, const vector<int> *vOrder)
MexStructParam::ReorderCols(GeneralMatrix &tdx, const vector<int> &vOrder, char* stat)
{
if (tdx.numCols() > vOrder.size())
{
mexPrintf(" Error in ReorderColumns - size of order var is too small");
return;
}
// GeneralMatrix tmp(*tdx); // temporary 2D matrix
GeneralMatrix tmpR(tdx); // temporary 2D matrix
// GeneralMatrix &tmpR = tmp;
tdx.zeros(); // empty original matrix
// reorder the columns
try
{
for (int i = 0; i < tdx.numCols(); i++)
tdx.copyColumns(tmpR, (vOrder)[i],vOrder[i], i);
// tdx->copyColumn(tmpR, (*vOrder)[i], i);
}
catch (const SylvException &e)
{
printf("Caugth exception in ReorderColumns: ");
e.printMessage();
return; // 255;
}
catch (...)
{
mexPrintf(" Error in ReorderColumns - wrong index?");
}
}
void
MexStructParam::ReorderCols(GeneralMatrix &tdx, const int *vOrder)
{
// GeneralMatrix tmp(*tdx); // temporary 2D matrix
GeneralMatrix tmpR(tdx); // temporary 2D matrix
// GeneralMatrix &tmpR = tdztmp;
tdx.zeros(); // empty original matrix
// reorder the columns
try
{
for (int i = 0; i < tdx.numCols(); i++)
tdx.copyColumns(tmpR, vOrder[i],vOrder[i], i);
}
catch (const SylvException &e)
{
printf("Caugth SYLV_EXCEPTION in ReorderColumns: ");
e.printMessage();
return; // 255;
}
catch (...)
{
mexPrintf(" Error in ReorderColumns - wrong index?");
}
}
/**************
void
//MexStructParam::ReorderCols(GeneralMatrix *tdx, const vector<int> *vOrder)
MexStructParam::ReorderCols(GeneralMatrix *tdx, const vector<int> &vOrder)
{
if (tdx->ncols() > vOrder->size())
{
mexPrintf(" Error in ReorderColumns - size of order var is too small");
return;
}
GeneralMatrix tmp(*tdx); // temporary 2D matrix
// GeneralMatrix& tmpR(tdx); // temporary 2D matrix
GeneralMatrix &tmpR = tmp;
tdx->zeros(); // empty original matrix
// reorder the columns
try
{
for (int i = 0; i < tdx->ncols(); i++)
tdx->copyColumn(tmpR, (vOrder)[i], i);
// tdx->copyColumn(tmpR, (*vOrder)[i], i);
}
catch (const TLException &e)
{
printf("Caugth TL exception in ReorderColumns: ");
e.print();
return; // 255;
}
catch (...)
{
mexPrintf(" Error in ReorderColumns - wrong index?");
}
}
void
MexStructParam::ReorderCols(GeneralMatrix *tdx, const int *vOrder)
{
GeneralMatrix tmp(*tdx); // temporary 2D matrix
// GeneralMatrix tmpR(tdx); // temporary 2D matrix
GeneralMatrix &tmpR = tdztmp;
tdx->zeros(); // empty original matrix
// reorder the columns
try
{
for (int i = 0; i < tdx->ncols(); i++)
tdx->copyColumn(tmpR, vOrder[i], i);
}
catch (const TLException &e)
{
printf("Caugth TL exception in ReorderColumns: ");
e.print();
return; // 255;
}
catch (...)
{
mexPrintf(" Error in ReorderColumns - wrong index?");
}
}
*******************/
//////////////////////////////////////////////////////
// Convert Matlab string array to C type array of string pointers
// Poblem is that Matlab mx function returns a long string concatenated by columns rather than rows
// hence a rather low level approach is needed
///////////////////////////////////////////////////////
const char **
MxArrayToStringArray(const mxArray *mxFldp, const int len, const int width)
{
char *cNamesCharStr = mxArrayToString(mxFldp);
const char **ret = MxArrayToStringArray(cNamesCharStr, len, width);
return ret;
}
const char **
MxArrayToStringArray(const char *cNamesCharStr, const int len=1, const int width=1)
{
char cNamesMX[len][width+1]; //
#ifdef DEBUG
mexPrintf("loop MxArrayToStringArray cNamesCharStr = %s \n", cNamesCharStr);
#endif
for (int i = 0; i < width; i++)
{
for (int j = 0; j < len; j++)
{
// Allow alphanumeric and underscores "_" only:
if (isalnum(cNamesCharStr[j+i*len]) || ('_' == cNamesCharStr[j+i*len]))
{
cNamesMX[j][i] = cNamesCharStr[j+i*len];
}
else cNamesMX[j][i] = '\0';
}
}
const char **ret = (const char **) mxCalloc(len, sizeof(char *));
for (int j = 0; j < len; j++)
{
cNamesMX[j][width] = '\0';
#ifdef DEBUG
// mexPrintf("String [%d]= %s \n", j, cNamesMX[j]);
#endif
char *token = (char *) mxCalloc(strlen(cNamesMX[j])+1, sizeof(char));
strcpy(token, cNamesMX[j]);
ret[j] = token;
#ifdef DEBUG
mexPrintf("ret [%d]= %s \n", j, ret[j]);
#endif
}
return ret;
}

229
mex/sources/estimation/MexUtils.h
View File

@ -1,229 +0,0 @@
/*
* 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/>.
*/
#ifndef MEXUTILS_H
#define MEXUTILS_H
#include <utility>
#include <vector>
#include "ioutils.h"
#include "mex.h"
#include "matrix.h"
//#include "SylvException.h"
#define MIN(X,Y) ((X) < (Y) ? (X) : (Y))
#define MAX(X,Y) ((X) > (Y) ? (X) : (Y))
//#include "k_ord_dynare.h"
//using namespace std;
enum {base, caller, global};
enum {M_, oo_, options_,bayestopt_, estim_params_, dr};
extern const char *DynareParamStructsNm [];
extern const char* mexBase[];
const int numParStructs=5;
//#define numParStructsM 6
/**
struct DynareParamsStructMap
{ char * paramName;
int dynParStruct;
}
struct DynareParamsStructBasePair
{
int dynParStruct;
int mexBase;
}
**/
struct StructBaseNameMap
{ char* baseName;
char* structName;
};
class MexStructParam;
class MexStruct :public virtual GeneralParams
{
vector <int> parStructBase;
vector <const mxArray*> parStruct; // array of struct pointers
map <string, int> parNamStructMap; // which struct par belongs
const int numParStruct;
const string structName;
mxArray* getMxField( const string& field)
{
map<string, int>::iterator it=parNamStructMap.find(field);
if (it==parNamStructMap.end())
{
mexPrintf("getMxField:no parameter with such name");
throw(SYLV_MES_EXCEPTION("getMxField:no parameter with such name"));
}
return mxGetField(parStruct[it->second], 0, field.c_str() );
}
StructBaseNameMap* getMxFieldStruct( const string& field)
{
map<string, int>::iterator it=parNamStructMap.find(field);
if (it==parNamStructMap.end())
throw(SYLV_MES_EXCEPTION("no parameter with such name"));
StructBaseNameMap* dpsm=new StructBaseNameMap;
dpsm->baseName=(char*)mexBase[parStructBase[it->second]];
dpsm->structName=(char*)DynareParamStructsNm[it->second];
return dpsm;
}
// void ReorderCols(GeneralMatrix* tdx, const vector<int>*vOrder);
// {KordpDynare::ReorderCols((TwoDMatrix*) tdx, *vOrder)};
// void ReorderCols(GeneralMatrix* tdx, const char*vOrder);
// {KordpDynare::ReorderCols((TwoDMatrix*) tdx, *vOrder)};
public:
MexStruct(int numParStr=1);
/**
MexStruct( char *sBases, char * mexStructNames, int numParStr=1)
: base(sBases), structName(mexStructNames), numParStruct(numParStr)
{
mexStruct=mexGetVariable(base, structName);
}
**/
~MexStruct(){};
void *
getMxArrayField(const string& field)
{
return getMxField( field);
}
string& name(){return *(new string(structName));};
MexStructParam& getMexStructField( const string& field);
GeneralParams& getStructField( const string& field)
{return (GeneralParams&) getMexStructField( field);};
void * getField( const string& field);
double& getDoubleField(const string& field);
string& getStringField(const string& field);
vector<string>& getStringVectorField(const string& field);
vector<int>& getIntVectorField(const string& field);
Vector& getDoubleVectorField(const string& field);
char * getCharField( const string& field);
charArraySt * getCharArrayField( const string& field);
GeneralMatrix & getMatrixField( const string& field);
void ReorderCols(GeneralMatrix& tdx, const vector<int>&vOrder);
// (MexStructParam::ReorderCols(tdx, vOrder););
// static void ReorderCols(GeneralMatrix& tdx, const int*vOrder);
// (MexStructParam::ReorderCols(tdx, vOrder););
void setField(const string& field, const string& newVal);
void setField(const string& field, const double val);
void setField(const string& field, const GeneralMatrix& val);
void UpdateMexField( const string& field, mxArray *newVal);
//void putMexStruct();//{mexPutVariable(base, structName, mexStruct);};
};
/***********************************************
* MexStructParam
* holds single Matlab structure passed as parameter
***********************************************/
class MexStructParam :public virtual GeneralParams
{
const mxArray* parStruct; // struct pointer
const GeneralParams* parStructParent; // if any
const string structName; // if any, param name of the structure in its parent.
mxArray* getMxField( const string& field)
{
return mxGetField(parStruct, 0, field.c_str());
}
// void ReorderCols(GeneralMatrix* tdx, const vector<int>*vOrder);
// void ReorderCols(GeneralMatrix* tdx, const char*vOrder);
public:
MexStructParam(const mxArray* paramStruct, const GeneralParams* parent, const string& name):
parStruct(paramStruct), parStructParent(parent), structName(name) {};
MexStructParam( const mxArray* paramStruct, const GeneralParams* parent):
parStruct(paramStruct), parStructParent(parent), structName(string("")){};
MexStructParam(const mxArray* paramStruct, const string& name):
parStruct(paramStruct), parStructParent(NULL), structName(name){};
~MexStructParam(){};
void *
getMxArrayField(const string& field)
{
return getMxField( field);
}
string& name(){return *(new string(structName));};
MexStructParam& getMexStructField( const string& field);
GeneralParams& getStructField( const string& field)
{return (GeneralParams&) getMexStructField( field);};
void * getField( const string& field);
double& getDoubleField(const string& field);
string& getStringField(const string& field);
vector<string>& getStringVectorField(const string& field);
vector<int>& getIntVectorField(const string& field);
Vector& getDoubleVectorField(const string& field);
char * getCharField( const string& field);
charArraySt * getCharArrayField( const string& field);
GeneralMatrix & getMatrixField( const string& field);
void ReorderCols(GeneralMatrix& tdx, const vector<int>&vOrder)
{ReorderCols(tdx, vOrder, "stat");};
static void ReorderCols(GeneralMatrix& tdx, const vector<int>&vOrder, char*stat);
static void ReorderCols(GeneralMatrix& tdx, const int*vOrder);
void setField(const string& field, const string& newVal);
void setField(const string& field, const double val);
void setField(const string& field, const GeneralMatrix& val);
void UpdateMexField( const string& field, mxArray *newVal);
//void putMexStruct();//{mexPutVariable(base, structName, mexStruct);};
};
/***************
class ConstMexStruct :public GeneralParams
{
const mxArray* mexStruct;
DynareParamsStructMap dynpsm;
const char *base; // one of: base, caller, global
const char *structName;
MexStruct( char *sBase, char * mexStructName)
: base(sBase), structName(mexStructName)
{
mexStruct=mexGetVariablePtr(base, structName);
}
virtual ~MexStruct(){};
mxArray *
getField(string& field)
{
return mxGetField(mexStruct, 0, field);
}
};
*******************/
//////////////////////////////////////////////////////
// Convert Matlab Dynare string array to C type array of string pointers
// Poblem is that Matlab mx function returns a long string concatenated by columns rather than rows
// hence a rather low level approach is needed
///////////////////////////////////////////////////////
const char **
MxArrayToStringArray(const mxArray *mxFldp, const int len, const int width);
const char **
MxArrayToStringArray(const char *cNamesCharStr, const int len, const int width);
#endif

28
mex/sources/estimation/ModelSolution.cc
View File

@ -32,15 +32,22 @@
*/
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)
const std::vector<size_t>& zeta_static_arg, 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),
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),
dynamicDLLp(modName, n_endo, n_jcols, /* nMax_lag= */ 1, n_exo, std::string(""))
{
Mx.setAll(0.0);
jacobian.setAll(0.0);
set_union(zeta_fwrd_arg.begin(), zeta_fwrd_arg.end(),
zeta_mixed_arg.begin(), zeta_mixed_arg.end(),
back_inserter(zeta_fwrd_mixed));
set_union(zeta_back_arg.begin(), zeta_back_arg.end(),
zeta_mixed_arg.begin(), zeta_mixed_arg.end(),
back_inserter(zeta_back_mixed));
}
void
@ -62,15 +69,14 @@ ModelSolution::ComputeModelSolution(Vector& steadyState, const Vector& deepParam
Vector llXsteadyState(n_jcols-n_exo);
for (size_t ll_row = 0; ll_row < ll_incidence.getRows(); ll_row++)
{
// populate (non-sparse) vector with ysteady values
for (size_t i = 0; i < n_endo; i++)
{
if (ll_incidence(ll_row, i))
llXsteadyState(((size_t) ll_incidence(ll_row, i))-1) = steadyState(i);
}
}
for (size_t i = 0; i < zeta_back_mixed.size(); i++)
llXsteadyState(i) = steadyState(zeta_back_mixed[i]);
for (size_t i = 0; i < n_endo; i++)
llXsteadyState(zeta_back_mixed.size() + i) = steadyState(i);
for (size_t i = 0; i < zeta_fwrd_mixed.size(); i++)
llXsteadyState(zeta_back_mixed.size() + n_endo + i) = steadyState(zeta_fwrd_mixed[i]);
//get jacobian
dynamicDLLp.eval(llXsteadyState, Mx, &deepParams, 1, residual, &jacobian, NULL, NULL);

4
mex/sources/estimation/ModelSolution.hh
View File

@ -40,7 +40,7 @@ class ModelSolution{
public:
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);
const std::vector<size_t>& zeta_static_arg, double qz_criterium);
virtual ~ModelSolution(){};
void compute( Vector& steadyState, const Vector& deepParams, Matrix& ghx, Matrix& ghu ) throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException);
@ -48,7 +48,7 @@ private:
const size_t n_endo;
const size_t n_exo;
const size_t n_jcols; // Num of Jacobian columns
const Matrix ll_incidence; // leads and lags indices
std::vector<size_t> zeta_fwrd_mixed, zeta_back_mixed;
Matrix jacobian;
Vector residual;
Matrix Mx;

62
mex/sources/estimation/Prior.hh Normal file
View File

@ -0,0 +1,62 @@
/*
* Copyright (C) 2009-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/>.
*/
///////////////////////////////////////////////////////////
// Prior.h
// Implementation of the Class Prior
// Created on: 02-Feb-2010 13:06:20
///////////////////////////////////////////////////////////
#if !defined(Prior_8D5F562F_C831_43f3_B390_5C4EF4433756__INCLUDED_)
#define Prior_8D5F562F_C831_43f3_B390_5C4EF4433756__INCLUDED_
struct Prior
{
public:
Prior(int shape, double mean, double mode, double lower_boound, double upper_boound, double fhp, double shp);
virtual ~Prior();
//! probablity density functions
enum pShape
{
Beta = 1,
Gamma = 2,
Gaussian = 3, // i.e. Normal density
Inv_gamma_1 = 4, // Inverse gamma (type 1) density
Uniform = 5,
Inv_gamma_2 = 6 //Inverse gamma (type 2) density
};
const pShape shape;
const double mean;
const double mode;
const double lower_boound;
const double upper_boound;
/**
* first shape parameter
*/
const double fhp;
/**
* second shape parameter
*/
const double shp;
};
#endif // !defined(8D5F562F_C831_43f3_B390_5C4EF4433756__INCLUDED_)

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

@ -1,303 +0,0 @@
/*
* 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/>.
*/
/***********************************************
% Based on Dynare Matlab
% function [dr,info,M_,options_,oo_] = SetDRModel(dr,task,M_,options_,oo_)
% formerly known as dr1 in Matalb Dynare
% Computes the reduced form solution of a rational expectation model (first or second order
% approximation of the stochastic model around the deterministic steady state).
%
%
********************************************************************/
#include "mexutils.h"
#include "DsgeLikelihood.h"
int
DsgeLikelihood::SolveDRModel(const int endo_nbr, const int exo_nbr, const int nstatic, const int npred, int nfwrd)//dr1()
{
int infoDR = 0;
int i;
// % 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");
else // dr=set_state_space(dr,M_); - to be done prior to calling estaimte!!
{
if (order ==1)
{
try
{
GeneralMatrix& ghx_u=SolveDRkOrderPert();//(dr,task,M_,options_, oo_ , ['.' mexext]);
//SteadyState=ysteady;
if (ghx_u.isZero())
{
mexPrintf("******** ghx_u is Zero ! *******\n");
//throw(1);
}
/**********/
int sss= ghx_u.numCols();
#ifdef DEBUG
mexPrintf("*********GHX_U colos %d then Allocate GHX and GHU *********\n", sss);
ghx_u.print();
#endif
vector<int>span2nx(sss-exo_nbr);
for (i=0;i<sss-exo_nbr;++i)
span2nx[i]=i+1;
//ghx= ( (ghx_u, nullVec,span2nx));//ghx_u(:,1:sss-M_.exo_nbr);
GeneralMatrix gh_x (ghx_u, nullVec,span2nx);//ghx_u(:,1:sss-M_.exo_nbr);
#ifdef DEBUG
mexPrintf("*********GH_X*********\n");
gh_x.print();
#endif
ghx= gh_x;
#ifdef DEBUG
mexPrintf("*********GHX*********\n");
ghx.print();
#endif
vector<int>spannx2end(exo_nbr);
for (i=0;i<exo_nbr;++i)
spannx2end[i]=sss-exo_nbr+i+1;
ghu= ( GeneralMatrix(ghx_u, nullVec,spannx2end)); //ghx_u(:,sss-M_.exo_nbr+1:end);
/**********
//Test only:
ghu=dr.getMatrixField(string("ghu"));
ghx=dr.getMatrixField(string("ghx"));
********/
#ifdef DEBUG
mexPrintf("*********GHU*********\n");
ghu.print();
// ghx.print();
#endif
// end test
delete &ghx_u;
}
catch(int e)
{
throw SYLV_MES_EXCEPTION("Problem with using k_order perturbation solver ");
info = 4;
penalty = 1000; // info(2)
infoDR=info;
return infoDR;
}//end
}
else //if options_.order >1
{
throw SYLV_MES_EXCEPTION("can not use order != 1 for estimation yet!");
info = 4;
penalty = 1000;//info(2) = 1000;
infoDR=info;
return infoDR;
};// end if
if ((int)dynareParams.getDoubleField(string("loglinear")) == 1)
{
#ifdef DEBUG
mexPrintf("setting SolveDRModel loglinear results\n");
#endif
//k = find(dr.kstate(:,2) <= M_.maximum_endo_lag+1);
vector<int>kk(0);
int maximum_endo_lag= (int)dynareParams.getDoubleField(string("maximum_endo_lag"));
for(i=0;i<kstate.numRows();++i)
if ( kstate.get(i,1)<=maximum_endo_lag+1)
kk.push_back(i+1);
//klag = dr.kstate(k,[1 2]);
vector<int>kk2(2);
kk2[0]=1;
kk2[1]=2;
#ifdef DEBUG
mexPrintf("setting klag for loglinear results\n");
#endif
GeneralMatrix klag (kstate, kk,kk2);
//k1 = dr.order_var;
vector<int>k1klag(0);
#ifdef DEBUG
mexPrintf("setting k1klag for loglinear results\n");
#endif
for (i=0; i< klag.numRows();++i)
if ((int) klag.get(i, 0)>0)
k1klag.push_back(order_var[(int) klag.get(i, 0)-1]);
//dr.ghx = repmat(1./dr.ys(k1),1,size(dr.ghx,2)).*dr.ghx.* ...
// ...repmat(dr.ys(k1(klag(:,1)))',size(dr.ghx,1),1);
Vector invOrdSS(endo_nbr);//SteadyState.length());
for (i=0;i<endo_nbr;++i)
invOrdSS[i]=1/SteadyState[order_var[i]-1];
#ifdef DEBUG
mexPrintf("setting mInvOrdSS for loglinear results\n");
#endif
GeneralMatrix mInvOrdSS(invOrdSS.base(),endo_nbr,1);
#ifdef DEBUG
mInvOrdSS.print();
mexPrintf("SolveDRModel Call repmat 1 for loglinear ghx results\n");
#endif
GeneralMatrix&repInvSSx=mInvOrdSS.repmat(1,ghx.numCols());
Vector k1klagSS(k1klag.size());
for (i=0;i<k1klag.size();++i)
k1klagSS[i]=SteadyState[k1klag[i]-1];
// GeneralMatrix mSSt(SteadyState.base(),1,endo_nbr);
// GeneralMatrix mk1klagSSt(mSSt, k1klag,nullVec);
GeneralMatrix mk1klagSSt(k1klagSS.base(), 1,k1klag.size());
#ifdef DEBUG
mk1klagSSt.print();
repInvSSx.print();
mexPrintf("SolveDRModel Call repmat 2 for loglinear ghx results\n");
#endif
GeneralMatrix&repk1klagSSt=mk1klagSSt.repmat(ghx.numRows(),1);
#ifdef DEBUG
mexPrintf("Final setting SolveDRModel loglinear ghx results\n");
#endif
ghx.multElements(repInvSSx);
ghx.multElements(repk1klagSSt);
//dr.ghu = repmat(1./dr.ys(k1),1,size(dr.ghu,2)).*dr.ghu;
#ifdef DEBUG
mexPrintf("SolveDRModel Call repmat 1 for loglinear ghu results\n");
#endif
GeneralMatrix&repInvSSu=mInvOrdSS.repmat(1,ghu.numCols());
#ifdef DEBUG
mexPrintf("Final setting SolveDRModel loglinear ghu results\n");
#endif
ghu.multElements(repInvSSu);
delete &repk1klagSSt;
delete &repInvSSu;
delete &repInvSSx;
#ifdef DEBUG
mexPrintf("Final loglinear ghu and ghx results\n");
ghu.print();
ghx.print();
#endif
};//end if
}
return infoDR;
}
/********************************************************************
* Solve the reduced DR k-order-model
*********************************************************************/
GeneralMatrix&
DsgeLikelihood::SolveDRkOrderPert() //kOrderPerturbation
{
// GeneralMatrix nullMat(0,0);
model->getParams()=deepParams;
model->getSteady()=SteadyState;
GeneralMatrix* dgyu=new GeneralMatrix (ghu.numRows(), ghx.numCols()+ghu.numCols());
dgyu->zeros();
try
{
#ifdef DEBUG
mexPrintf(" DeepParams:.\n");
deepParams.print();
mexPrintf(" Calling walkStochSteady with Params:.\n");
model->getParams().print();
#endif
approx->walkStochSteady();
#ifdef DEBUG
mexPrintf("End of walkStochSteady - write map.\n");
#endif
// Write derivative outputs into memory map
map<string, ConstTwoDMatrix> mm;
approx->getFoldDecisionRule().writeMMap(&mm);
#ifdef DEBUG
mexPrintf("k_order_perturbation: Map print: \n");
approx->getFoldDecisionRule().print();
#endif
#ifdef DEBUG
mexPrintf("k_order_perturbation: Map print: \n");
for (map<string, ConstTwoDMatrix>::const_iterator cit = mm.begin();
cit != mm.end(); ++cit)
{
mexPrintf("k_order_perturbation: Map print: string: %s , g:\n", (*cit).first.c_str());
(*cit).second.print();
}
#endif
// get latest ysteady
// SteadyState=model->getSteady();
#ifdef DEBUG
// mexPrintf("Steady State\n");
// SteadyState.print();
#endif
// developement of the output.
#ifdef DEBUG
mexPrintf("k_order_perturbation: Filling outputs.\n");
#endif
int ii=1;
// Set the output pointer to the combined output matrix gyu = [gy gu].
for (map<string, ConstTwoDMatrix>::const_iterator cit = mm.begin();
((cit != mm.end()) && (ii < 4)); ++cit)
{
if ((*cit).first!="g_0" && ii==2)
{
dgyu->getData() = (*cit).second.getData();
#ifdef DEBUG
mexPrintf("k_order_perturbation: cit %d numRows %d numCols %d print: \n", ii, (*cit).second.numRows(), (*cit).second.numCols());
(*cit).second.print();
mexPrintf("k_order_perturbation: dguy output %d print: \n", ii);
dgyu->print(); //!! This print Crashes???
#endif
return *dgyu;
}
++ii;
}
return *dgyu;
}
catch (const KordException &e)
{
printf("Caugth Kord exception in SolveDRkOrderPert: ");
e.print();
mexPrintf("Caugth Kord exception: %s", e.get_message());
}
catch (const TLException &e)
{
printf("Caugth TL exception in SolveDRkOrderPert: ");
mexPrintf("Caugth TL exception in SolveDRkOrderPert: ");
e.print();
}
catch (SylvException &e)
{
printf("Caught Sylv exception in SolveDRkOrderPert: ");
mexPrintf("Caught Sylv exception in SolveDRkOrderPert: ");
e.printMessage();
}
catch (const DynareException &e)
{
printf("Caught KordpDynare exception in SolveDRkOrderPert: %s\n", e.message());
mexPrintf("Caugth Dynare exception in SolveDRkOrderPert: %s", e.message());
}
catch (const ogu::Exception &e)
{
printf("Caught ogu::Exception in SolveDRkOrderPert: ");
e.print();
mexPrintf("Caugth general exception inSolveDRkOrderPert: %s", e.message());
} //catch
}; // end of mexFunction()

199
mex/sources/estimation/dynare_resolve.cpp
View File

@ -1,199 +0,0 @@
/*
* 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/>.
*/
/*****************************************************
% Based on Matlab Dynare
% function [int info] = dynare_resolve(ys,iv,ic,aux)
%
% function [A,B,ys,info] = dynare_resolve(iv,ic,aux)
% Computes the linear approximation and the matrices A and B of the
% transition equation and doing:
% check if ys is steady state
% dr
% kalman_transition_matrix
%
% INPUTS
% iv: selected variables (observed and state variables)
% ic: state variables position in the transition matrix columns
% aux: indices for auxiliary equations
%
% MODIFIES
% A: matrix of predetermined variables effects in linear solution (ghx)
% B: matrix of shocks effects in linear solution (ghu)
% ys: steady state of original endogenous variables
%
% OUTPUTS
% info=1: the model doesn't determine the current variables '...' uniquely
% info=2: MJDGGES returns the following error code'
% info=3: Blanchard Kahn conditions are not satisfied: no stable '...' equilibrium
% info=4: Blanchard Kahn conditions are not satisfied:'...' indeterminacy
% info=5: Blanchard Kahn conditions are not satisfied:'...' indeterminacy due to rank failure
% info=20: can't find steady state info(2) contains sum of sqare residuals
% info=30: variance can't be computed
%
% SPECIAL REQUIREMENTS
% none
********************************************************************/
#include "dsgeLikelihood.h"
int
DsgeLikelihood::dynareResolveDR(vector<int>&iv,vector<int>&ic,GeneralMatrix& aux) // i.e. dynare_resolve()
{
//here comes a subset of [oo_.dr,info] = resol(oo_.steady_state,0);
// check if ys is steady state and calcluate new one if not
// testing for steadystate file: To Be Icluded at a later stage
#ifdef DEBUG
mexPrintf("Calling SolveDRModel\n");
#endif
// [dr,info,M_,options_,oo_] = dr1(dr,check_flag,M_,options_,oo_);
int DRinfo = SolveDRModel(endo_nbr,exo_nbr,nstatic, npred, nfwrd); //formerly known as dr1 in Matalb Dynare, i.e.
if (DRinfo)
{
info = DRinfo;
return DRinfo;
}
// End of resol: now rest of dynare_resolve:
// if nargin == 0
// if (iv.size()==0)
#ifdef DEBUG
mexPrintf(" if iv==NULL\n");
#endif
if (&iv==NULL)
{
#ifdef DEBUG
mexPrintf(" construct iv\n");
#endif
//iv = (1:endo_nbr)';
for (int i=1;i<=endo_nbr;++i)
iv.push_back(i);//= (1:endo_nbr)';
}
// if (ic.size()==0)
#ifdef DEBUG
mexPrintf(" if ic==NULL\n");
#endif
if (&ic==NULL)
{
#ifdef DEBUG
mexPrintf(" construct ic\n");
#endif
//ic = [ nstatic+(1:npred) endo_nbr+(1:size(oo_.dr.ghx,2)-npred) ]';
//ic(npred+ghx.numCols());
for(int i=0;i<npred;++i)
{
ic.push_back(i+nstatic+1);
}
for(int j=0;j<ghx.numCols()-npred;++j)
ic.push_back(j+endo_nbr+1);
}
#ifdef DEBUG
mexPrintf(" if aux==NULL\n");
#endif
if (&aux==NULL)
{
#ifdef DEBUG
mexPrintf(" construct aux\n");
#endif
int i;
aux =dr.getMatrixField(string("transition_auxiliary_variables"));
//k = find(aux(:,2) > npred);
//aux(:,2) = aux(:,2) + nstatic;
vector<int>k(0);
for (i=0; i< aux.numRows();++i)
{
if (aux.get(i,1)>npred)
k.push_back(i+1);
aux.get(i, 1)+=nstatic;
}
//aux(k,2) = aux(k,2) + oo_.dr.nfwrd;
for ( i=0; i< k.size();++i)
aux.get(k[i]-1,1) +=nfwrd;
}//end if
#ifdef DEBUG
mexPrintf("[A,B] = kalman_transition_matrix\n");
#endif
// here is the content of [T R]=[A,B] = kalman_transition_matrix(oo_.dr,iv,ic,aux,M_.exo_nbr);
int n_iv = iv.size();//length(iv);
int n_ir1 = aux.numRows();// size(aux,1);
int nr = n_iv + n_ir1;
// GeneralMatrix A=*(new GeneralMatrix (nr,nr));
// A.zeros();
// GeneralMatrix B=*(new GeneralMatrix(nr,exo_nbr));
// B.zeros();
T.zeros();
R.zeros();
vector<int>i_n_iv(n_iv);
for (int i=0;i<n_iv;++i)
i_n_iv[i]=i+1;//= (1:n_iv)';
#ifdef DEBUG
mexPrintf("T=A assign by vec\n");
#endif
//A(i_n_iv,ic) = dr.ghx(iv,:);
//A.AssignByVectors (i_n_iv,ic, ghx, iv, nullVec);//dr.ghx(iv,:);
T.AssignByVectors (i_n_iv,ic, ghx, iv, nullVec);//dr.ghx(iv,:);
#ifdef DEBUG
mexPrintf("Completed T=A assign by vec\n");
T.print();
#endif
if (n_ir1 > 0)
{
//A(n_iv+1:end,:) = sparse(aux(:,1),aux(:,2),ones(n_ir1,1),n_ir1,nr);
#ifdef DEBUG
mexPrintf("Create sparse\n");
#endif
GeneralMatrix sparse(n_ir1,nr);
for (int i=0;i<n_ir1;++i)
sparse.get((int)aux.get(i,0)-1,(int)aux.get(i,1)-1)=1;
/* vector<int>span2end(A.numRows()-n_iv);
for (int i=0;i<A.numRows()-n_iv;++i)
span2end[i]=i+n_iv+1;
A.place(sparse,n_iv,0); and T=A;
*/
#ifdef DEBUG
mexPrintf("T.place (sparse,n_iv,0)\n");
#endif
T.place(sparse,n_iv,0);
}
#ifdef DEBUG
mexPrintf("R=B assign by vec R rows %d cols %d \n", R.numRows(), R.numCols());
#endif
//B(i_n_iv,:) = dr.ghu(iv,:); and R=B;
R.AssignByVectors(i_n_iv, nullVec, ghu, iv, nullVec);
#ifdef DEBUG
mexPrintf("Completed R=B assign by vec\n");
R.print();
#endif
// ys = oo_.dr.ys;
// GeneralMatrix&ysmx = dr.getMatrixField(string("ys"));
// SteadyState=ysmx.getData();
return DRinfo;
}

75
mex/sources/estimation/ioutils.h
View File

@ -1,75 +0,0 @@
/*
* 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/>.
*/
#ifndef IOUTILS_H
#define IOUTILS_H
#include "GeneralMatrix.h"
#include "SylvException.h"
#include <map>
#include <string>
//using namespace std;
struct charArraySt
{
char ** charArrayPtr;
int len;
};
class GeneralParams
{
// map <string, int> params;
// const char *structName;
public:
GeneralParams(){};
virtual ~GeneralParams(){};
virtual string& name()=0;
virtual void *
getField(const string& field)=0;
virtual double&
getDoubleField(const string& field)=0;
virtual string &
getStringField(const string& field)=0;
virtual vector<string>&
getStringVectorField(const string& field)=0;
virtual vector<int>&
getIntVectorField(const string& field)=0;
virtual Vector&
getDoubleVectorField(const string& field)=0;
virtual GeneralParams&
getStructField( const string& field)=0;
virtual char *
getCharField(const string& field)=0;
virtual charArraySt *
getCharArrayField( const string& field)=0;
// uses General Matrix from sylv library
virtual GeneralMatrix&
getMatrixField(const string& field)=0;
virtual void ReorderCols(GeneralMatrix& tdx, const vector<int>&vOrder)=0;
// virtual void ReorderCols(GeneralMatrix& tdx, const int*vOrder)=0;
virtual void
// overloaded pramater update "set" methods:
setField(const string& field, const string&newVal)=0;
virtual void
setField(const string& field, const double val)=0;
virtual void
setField(const string& field, const GeneralMatrix& val)=0;
};
#endif

48
mex/sources/estimation/libmat/VDVEigDecomposition.cc Normal file
View File

@ -0,0 +1,48 @@
/*
* 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 "VDVEigDecomposition.hh"
VDVEigDecomposition::VDVEigDecomposition(const Matrix &m) throw(VDVEigException) :
lda(m.getLd()), n(m.getCols()), lwork(-1),
info(0), converged(false), V(m), D(n)
{
if (m.getRows() != m.getCols())
throw(VDVEigException(info, "Matrix is not square in VDVEigDecomposition constructor"));
dsyev("V", "U", &n, V.getData(), &lda, D.getData(), &tmpwork, &lwork, &info);
lwork = (int) tmpwork;
work = new double[lwork];
if (info < 0)
throw(VDVEigException(info, "Internal error in VDVEigDecomposition constructor"));
}
VDVEigDecomposition::VDVEigDecomposition(size_t inn) throw(VDVEigException) :
lda(inn), n(inn), lwork(3*inn-1),
info(0), converged(false), V(inn), D(inn)
{
work = new double[lwork];
};
std::ostream &
operator<<(std::ostream &out, const VDVEigDecomposition::VDVEigException &e)
{
out << " info " << e.info << ": " << e.message << "\n";
return out;
}

124
mex/sources/estimation/libmat/VDVEigDecomposition.hh Normal file
View File

@ -0,0 +1,124 @@
/*
* 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/>.
*/
/**
* Calculates vector of Eigen values D and matrix of Eigen vectors V
* of input matrix m
*/
#if !defined(VDVEigDec_INCLUDE)
#define VDVEigDec_INCLUDE
#include <cstdlib>
#include "Vector.hh"
#include "Matrix.hh"
#include <dynlapack.h>
class VDVEigDecomposition
{
lapack_int lda, n;
int lwork, info;
double tmpwork;
double *work;
bool converged;
Matrix V;
Vector D;
public:
public:
class VDVEigException
{
public:
const int info;
std::string message;
VDVEigException(int info_arg, std::string message_arg) :
info(info_arg), message(message_arg) {
};
};
/**
* This constructor only creates optimal workspace using
* input matrix m
*/
VDVEigDecomposition(const Matrix &m) throw(VDVEigException);
/**
* This constructoro only crates workspace using the size of
* the input matrix m
*/
VDVEigDecomposition(size_t n) throw(VDVEigException);
virtual ~VDVEigDecomposition()
{
if (work)
delete[] work;
};
template <class Mat>
void calculate(const Mat &H) throw(VDVEigException);
// get eigenvalues
Vector &
getD()
{
return D;
};
// get eigen vector
Matrix &
getV()
{
return V;
};
// check if converged
bool
hasConverged()
{
return converged;
};
};
template <class Mat>
void
VDVEigDecomposition::calculate(const Mat &m) throw(VDVEigException)
{
info = 0;
if (m.getRows() != m.getCols())
throw(VDVEigException(info, "Matrix is not square in VDVEigDecomposition calculate"));
if (m.getCols() != (size_t) n || m.getLd() != (size_t) lda)
throw(VDVEigException(info, "Matrix not matching VDVEigDecomposition class"));
int tmplwork = -1;
V = m;
dsyev("V", "U", &n, V.getData(), &lda, D.getData(), &tmpwork, &tmplwork, &info);
if (lwork < (int) tmpwork)
{
lwork = (int) tmpwork;
delete[] work;
work = new double[lwork];
}
dsyev("V", "U", &n, V.getData(), &lda, D.getData(), work, &lwork, &info);
if (info < 0)
throw(VDVEigException(info, "Internal error in VDVEigDecomposition calculation"));
converged = true;
if (info)
converged = false;
}
std::ostream &operator<<(std::ostream &out, const VDVEigDecomposition::VDVEigException &e);
#endif

302
mex/sources/estimation/tests/DsgeLikelihood_mexf.cpp
View File

@ -1,302 +0,0 @@
/*
* 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)
% 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"
#include "mexutils.h"
extern const char *DynareParamStructsNm []={"M_", "oo_", "options_", "bayestopt_", "estim_params_", "dr"};
extern const char* mexBase[]={"base", "caller", "global"};
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("estimation: mexExt=%s.\n", dfExt);
#endif
/***********
***************/
int numPeriods=1;
//MexStruct dynareParams();
MexStruct& dynareParams=*(new MexStruct(numParStructs));
#ifdef DEBUG
mexPrintf("getting dr\n");
#endif
MexStructParam& dr=dynareParams.getMexStructField("dr");
vector<int>&mfys=dynareParams.getIntVectorField("mfys");
vector<int>&mf=dynareParams.getIntVectorField("mf1");
#ifdef DEBUG
mexPrintf("getting SS\n");
#endif
Vector& SteadyState=dr.getDoubleVectorField(string("ys"));
#ifdef DEBUG
int gg;
for ( gg=0;gg<SteadyState.length();++gg)
mexPrintf("SteadyState %d = %f\n", gg, SteadyState[gg]);
#endif
int numVarobs=data.numRows();
Vector constant(numVarobs);//=*(new Vector(numVarobs));// = *(new Vector(nobs));
GeneralMatrix&kstate = dr.getMatrixField(string("kstate"));
vector<int>&order_var = dr.getIntVectorField(string("order_var"));
#ifdef DEBUG
for ( gg=0;gg<order_var.size();++gg)
mexPrintf("order_var %d = %d\n", gg, order_var[gg]);
#endif
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"));
#ifdef DEBUG
for ( gg=0;gg<lb.length();++gg)
mexPrintf("lb %d = %f\n", gg, lb[gg]);
#endif
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;
#ifdef DEBUG
mexPrintf("jcols = %d, exo_nbr=%d\n", jcols, exo_nbr);
#endif
GeneralMatrix& aux = dynareParams.getMatrixField(string("restrict_aux"));
vector<int>&iv= dynareParams.getIntVectorField(string("restrict_var_list"));
vector<int>&ic= dynareParams.getIntVectorField(string("restrict_columns"));
int nr=iv.size()+aux.numRows(); // Size of T matrix
Vector& a_init=*(new Vector(nr));
a_init.zeros();
GeneralMatrix& Q = dynareParams.getMatrixField(string("Sigma_e"));
#ifdef DEBUG
Q.print();
#endif
GeneralMatrix& Hrtmp = dynareParams.getMatrixField(string("H"));
GeneralMatrix * Hp;
if (Hrtmp.numCols()==0 || Hrtmp.numRows()==0)
{
delete &Hrtmp;
Hp = new GeneralMatrix(numVarobs,numVarobs);
Hp->zeros();
#ifdef DEBUG
mexPrintf("finished local initialising of H \n");
#endif
}
else
Hp=&Hrtmp;
GeneralMatrix& H=*Hp;
GeneralMatrix Y(data.numRows(),data.numCols());
GeneralMatrix T(nr,nr);
GeneralMatrix Z(numVarobs,nr);
Z.zeros();
for (int i = 0;i<numVarobs;++i)
Z.get(i,mf[i]-1)=1;
GeneralMatrix Pstar(nr,nr);
GeneralMatrix R(nr,exo_nbr);
GeneralMatrix ghx(endo_nbr,nsPred);
GeneralMatrix ghu(endo_nbr,exo_nbr);
//Pinf=[]
GeneralMatrix Pinf (nr,nr);
Pinf.zeros();
double loglikelihood;
try
{
#ifdef DEBUG
mexPrintf("Try construction of DsgeLikelihood\n");
#endif
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, aux, iv, ic, jcols, dfExt);
#ifdef DEBUG
mexPrintf("Try CalcLikelihood\n");
#endif
#ifdef LL_TIMING_LOOP
mexPrintf("DsgeLikelihood::CalcLikelihood: starting 1000 loops\n");
for (int tt=0;tt<1000;++tt)
{
#endif
loglikelihood=dl.CalcLikelihood(xparam1);
#ifdef LL_TIMING_LOOP
}
mexPrintf("DsgeLikelihood::CalcLikelihood: finished 1000 loops\n");
#endif
/*****************************************************************
% 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.
*****************************************************************/
#ifdef DEBUG
mexPrintf("Try Outputs with nper=%d, loglikelihood = %f\n",nper,loglikelihood);
#endif
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);
#ifdef DEBUG
mexPrintf("SteadyState size %d \n", dl.getSteadyState().length());
dl.getSteadyState().print() ;
mexPrintf("Try getSteadyState into vss size %d \n", vss.length());
#endif
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];
#ifdef DEBUG
mexPrintf("mxll[%d]=%f vll[%d]=%f\n",j, mxll[j], i, vll[j]);
#endif
}
}
*********************/
}
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

32
mex/sources/estimation/tests/global_access_test.cpp
View File

@ -1,32 +0,0 @@
#include "mex.h"
#include "matrix.h"
extern "C" {
// mexFunction: Matlab Inerface point and the main application driver
void
mexFunction(int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[])
{
// const char *base="base";
mxArray* options_ = mexGetVariable("caller", "options_");
int kOrder;
mxArray *mxFldp = mxGetField(options_, 0, "order");
if (mxIsNumeric(mxFldp))
kOrder = (int) mxGetScalar(mxFldp);
else
kOrder = 1;
mexPrintf("order: %d \n",kOrder);
int new_order=23;
mxArray *mxNewOrder=mxCreateDoubleScalar((double) new_order);
mxSetField(options_,0,"order", mxNewOrder);
mexPutVariable("caller", "options_", options_);
}
}

27
mex/sources/estimation/tests/testInitKalman.cc
View File

@ -66,29 +66,6 @@ main(int argc, char **argv)
0.0, 0.000078535044
};
Matrix
ll_incidence(3, n_endo); // leads and lags indices
double inllincidence[] = {
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.3560,
0.9930,
@ -106,8 +83,6 @@ main(int argc, char **argv)
steadyStateVW(dYSparams, n_endo, 1);
steadyState = steadyStateVW;
std::cout << "Vector deepParams: " << std::endl << deepParams << 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;
// Set zeta vectors [0:(n-1)] from Matlab indices [1:n] so that:
@ -178,7 +153,7 @@ main(int argc, char **argv)
InitializeKalmanFilter
initializeKalmanFilter(modName, n_endo, n_exo,
zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, ll_incidence, qz_criterium,
zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, qz_criterium,
order_var, inv_order_var, riv, ric, lyapunov_tol, info);
std::cout << "Initilise KF with Q: " << std::endl << Q << std::endl;

27
mex/sources/estimation/tests/testKalman.cc
View File

@ -66,29 +66,6 @@ main(int argc, char **argv)
0.0, 0.000078535044
};
Matrix
ll_incidence(3, n_endo); // leads and lags indices
double inllincidence[] = {
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.3560,
0.9930,
@ -106,8 +83,6 @@ main(int argc, char **argv)
steadyStateVW(dYSparams, n_endo, 1);
steadyState = steadyStateVW;
std::cout << "Vector deepParams: " << std::endl << deepParams << 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;
// Set zeta vectors [0:(n-1)] from Matlab indices [1:n] so that:
@ -182,7 +157,7 @@ main(int argc, char **argv)
double penalty = 1e8;
KalmanFilter kalman(modName, n_endo, n_exo,
zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, ll_incidence, qz_criterium,
zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, qz_criterium,
order_var, inv_order_var, varobs, riv, ric, riccati_tol, lyapunov_tol, info);
std::cout << "Initilise KF with Q: " << std::endl << Q << std::endl;

25
mex/sources/estimation/tests/testModelSolution.cc
View File

@ -56,27 +56,6 @@ main (int argc, char** argv)
Matrix vCov (nVCVpar, nVCVpar);
vCov = vCovVW;
Matrix ll_incidence(3,n_endo); // leads and lags indices
double inllincidence[]={
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,
@ -91,8 +70,6 @@ main (int argc, char** argv)
steadyState=steadyStateVW;
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;
// Set zeta vectors [0:(n-1)] from Matlab indices [1:n]
@ -114,7 +91,7 @@ main (int argc, char** argv)
Matrix ghu(n_endo,n_exo);
ModelSolution modelSolution( modName, n_endo, n_exo
, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, ll_incidence, qz_criterium);
, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, qz_criterium);
modelSolution.compute(steadyState, deepParams, ghx, ghu);