diff --git a/mex/sources/estimation/DsgeLikelihood.cpp b/mex/sources/estimation/DsgeLikelihood.cpp deleted file mode 100644 index ed0625d04..000000000 --- a/mex/sources/estimation/DsgeLikelihood.cpp +++ /dev/null @@ -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 . -*/ - - -#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& INorder_var, const vector& INmfys, const vector& INmf1, - Vector& INxparam1, const int INnum_dp, Vector& INdeepParams, - const Vector& INub, const Vector& INlb, const vector&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&INiv, vector&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* vll=new std::vector (nper); -// vll=new std::vector (numTimeObs);// vector of likelihoods - vll=new vector(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 *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 ⁣ -/******** deleting Vectors Crashes !************ - delete &SteadyState; - delete ¶m_ub; - delete ¶m_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;sslckk(0); - double qdelta=0; - for(i=0;i0) // 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;iparam_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&lgyidx2varobs= dynareParams.getIntVectorField(string("lgyidx2varobs")); - for (i=0;i 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. -*/ -#include "Estimation.h" -#include "k_ord_dynare.h" -#include "kalman.h" -#include "math.h" -#include "disclyap_fast.h" - -#include - -#include - - -class DsgeLikelihood - { - double likelihood; // sum of vector of KF step log likelihoods - vector* 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& order_var; - const vector& mfys; - const vector& 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&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&iv; //= dynareParams.getIntVectorField(string("restrict_var_list")); - vector⁣ //= 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¶ms); - // 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&iv,vector&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¶ms, 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& INorder_var, const vector& INmfys, const vector& INmf, - Vector& INxparam1, const int INnum_dp, Vector& INdeepParams, - const Vector& INub, const Vector& INlb, const vector&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&iv, vector&ic - ,const int jcols, const char *dfExt); //, KordpDynare& inModel, Approximation& INapprox ); - DsgeLikelihood( const Vector¶ms,const GeneralMatrix&data, const vector& 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& getLikVector() {return *vll;} // vector of log likelihoods for each Kalman step - //GeneralMatrix&lyapunov_symm(const GeneralMatrix &G, const GeneralMatrix & V); - }; - - diff --git a/mex/sources/estimation/EstimatedParameter.hh b/mex/sources/estimation/EstimatedParameter.hh new file mode 100644 index 000000000..596e4c18f --- /dev/null +++ b/mex/sources/estimation/EstimatedParameter.hh @@ -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 . + */ +/////////////////////////////////////////////////////////// +// 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 +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 subSampleIDs; +}; + +#endif // !defined(D879C8AE_5B69_4fc3_83BD_FA5A99030ECF__INCLUDED_) diff --git a/mex/sources/estimation/EstimatedParametersDescription.hh b/mex/sources/estimation/EstimatedParametersDescription.hh new file mode 100644 index 000000000..bf1a0bb6b --- /dev/null +++ b/mex/sources/estimation/EstimatedParametersDescription.hh @@ -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 . + */ +/////////////////////////////////////////////////////////// +// 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 &estSubsamples, std::vector &estParams); +//Vector getAllParams(int period); +//int getNumPeriods(); + EstimationSubsample&getSubsample(int period); + std::vector&getAllSubsamples(); +//EstimatedParameter& getParameter(size_t ID); + +//private: + + std::vector estSubsamples; + std::vector estParams; +//vector prior; +}; + +#endif // !defined(E8F2C096_A301_42e8_80FF_A643291BF995__INCLUDED_) diff --git a/mex/sources/estimation/EstimationSubsample.hh b/mex/sources/estimation/EstimationSubsample.hh new file mode 100644 index 000000000..47b3d5c63 --- /dev/null +++ b/mex/sources/estimation/EstimationSubsample.hh @@ -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 . + */ + +/////////////////////////////////////////////////////////// +// 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 +#include + +/** + * 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 subSampleParamIDs; +}; + +#endif // !defined(34FB1351_4066_4993_B384_56AA65BC32C4__INCLUDED_) diff --git a/mex/sources/estimation/InitializeKalmanFilter.cc b/mex/sources/estimation/InitializeKalmanFilter.cc index 76e5d7330..da16e46f5 100644 --- a/mex/sources/estimation/InitializeKalmanFilter.cc +++ b/mex/sources/estimation/InitializeKalmanFilter.cc @@ -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 &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &zeta_static_arg, - const Matrix &llincidence, double qz_criterium, const std::vector &order_var_arg, const std::vector &inv_order_var_arg, + double qz_criterium, const std::vector &order_var_arg, const std::vector &inv_order_var_arg, const std::vector &rivIN, const std::vector &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), diff --git a/mex/sources/estimation/InitializeKalmanFilter.hh b/mex/sources/estimation/InitializeKalmanFilter.hh index ab227b6c1..27af85c8d 100644 --- a/mex/sources/estimation/InitializeKalmanFilter.hh +++ b/mex/sources/estimation/InitializeKalmanFilter.hh @@ -42,7 +42,7 @@ class InitializeKalmanFilter public: InitializeKalmanFilter(const std::string& modName, size_t n_endo, size_t n_exo, const std::vector &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &zeta_static_arg, - const Matrix &llincidence, double qz_criterium, const std::vector &order_var_arg, const std::vector &inv_order_var_arg, + double qz_criterium, const std::vector &order_var_arg, const std::vector &inv_order_var_arg, const std::vector &riv, const std::vector &ric, double lyapunov_tol, int &info); virtual ~InitializeKalmanFilter(); // initialise all KF matrices diff --git a/mex/sources/estimation/KalmanFilter.cc b/mex/sources/estimation/KalmanFilter.cc index 9b4510f80..2ed8b40f6 100644 --- a/mex/sources/estimation/KalmanFilter.cc +++ b/mex/sources/estimation/KalmanFilter.cc @@ -33,7 +33,7 @@ KalmanFilter::~KalmanFilter() KalmanFilter::KalmanFilter(const std::string &modName, size_t n_endo, size_t n_exo, const std::vector &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &zeta_static_arg, - const Matrix &ll_incidence, double qz_criterium, const std::vector &order_var, + double qz_criterium, const std::vector &order_var, const std::vector &inv_order_var, const std::vector &varobs_arg, const std::vector &riv, const std::vector &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); diff --git a/mex/sources/estimation/KalmanFilter.hh b/mex/sources/estimation/KalmanFilter.hh index c65e665a0..4edaa193d 100644 --- a/mex/sources/estimation/KalmanFilter.hh +++ b/mex/sources/estimation/KalmanFilter.hh @@ -50,7 +50,7 @@ public: virtual ~KalmanFilter(); KalmanFilter(const std::string &modName, size_t n_endo, size_t n_exo, const std::vector &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &zeta_static_arg, - const Matrix &llincidence, double qz_criterium, const std::vector &order_var_arg, const std::vector &inv_order_var, + double qz_criterium, const std::vector &order_var_arg, const std::vector &inv_order_var, const std::vector &varobs_arg, const std::vector &riv, const std::vector &ric, double riccati_tol, double lyapunov_tol, int &info); diff --git a/mex/sources/estimation/LogLikelihoodMain.hh b/mex/sources/estimation/LogLikelihoodMain.hh new file mode 100644 index 000000000..4f6582333 --- /dev/null +++ b/mex/sources/estimation/LogLikelihoodMain.hh @@ -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 . + */ +/////////////////////////////////////////////////////////// +// 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 &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 &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, + const std::vector &zeta_static_arg, const double qz_criterium, const std::vector &order_var, + const std::vector &inv_order_var, const std::vector &varobs, const std::vector &riv, + const std::vector &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_) diff --git a/mex/sources/estimation/LogLikelihoodSubSample.cc b/mex/sources/estimation/LogLikelihoodSubSample.cc index ab66589db..4193aae9d 100644 --- a/mex/sources/estimation/LogLikelihoodSubSample.cc +++ b/mex/sources/estimation/LogLikelihoodSubSample.cc @@ -31,10 +31,10 @@ LogLikelihoodSubSample::~LogLikelihoodSubSample() LogLikelihoodSubSample::LogLikelihoodSubSample(const std::string &modName, size_t n_endo, size_t n_exo, const std::vector &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, - const std::vector &zeta_static_arg, const Matrix &ll_incidence, const double qz_criterium, const std::vector &order_var, + const std::vector &zeta_static_arg, const double qz_criterium, const std::vector &order_var, const std::vector &inv_order_var, const std::vector &varobs, const std::vector &riv, const std::vector &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) { }; diff --git a/mex/sources/estimation/LogLikelihoodSubSample.hh b/mex/sources/estimation/LogLikelihoodSubSample.hh index 06d115eab..23b9146c2 100644 --- a/mex/sources/estimation/LogLikelihoodSubSample.hh +++ b/mex/sources/estimation/LogLikelihoodSubSample.hh @@ -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 &zeta_fwrd_arg, const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, - const std::vector &zeta_static_arg, const Matrix &ll_incidence, const double qz_criterium, const std::vector &order_var, + const std::vector &zeta_static_arg, const double qz_criterium, const std::vector &order_var, const std::vector &inv_order_var, const std::vector &varobs, const std::vector &riv, const std::vector &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); }; diff --git a/mex/sources/estimation/MexUtils.cpp b/mex/sources/estimation/MexUtils.cpp deleted file mode 100644 index aaaa12ec3..000000000 --- a/mex/sources/estimation/MexUtils.cpp +++ /dev/null @@ -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 . -*/ - -#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::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& -MexStruct::getStringVectorField(const string& field) - { - charArraySt * cap = getCharArrayField(field); - vector *sv=(new vector(cap->len)); - for (int i= 0; ilen;++i) - (*sv)[i]=string(cap->charArrayPtr[i]); - return *sv; - } - -vector& -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 *vars = (new vector(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 &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& -MexStructParam::getStringVectorField(const string& field) - { - charArraySt * cap = getCharArrayField( field); - vector *sv= (new vector(cap->len)); - for (int i= 0; ilen;++i) - (*sv)[i]=string(cap->charArrayPtr[i]); - return *sv; - } - -vector& -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 *vars = (new vector(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 *vOrder) -MexStructParam::ReorderCols(GeneralMatrix &tdx, const vector &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 *vOrder) -MexStructParam::ReorderCols(GeneralMatrix *tdx, const vector &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; - } diff --git a/mex/sources/estimation/MexUtils.h b/mex/sources/estimation/MexUtils.h deleted file mode 100644 index 7a1ccdf37..000000000 --- a/mex/sources/estimation/MexUtils.h +++ /dev/null @@ -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 . -*/ -#ifndef MEXUTILS_H -#define MEXUTILS_H -#include -#include -#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 parStructBase; - vector parStruct; // array of struct pointers - map parNamStructMap; // which struct par belongs - const int numParStruct; - const string structName; - - mxArray* getMxField( const string& field) - { - map::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::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*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& getStringVectorField(const string& field); - vector& 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&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*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& getStringVectorField(const string& field); - vector& 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&vOrder) - {ReorderCols(tdx, vOrder, "stat");}; - static void ReorderCols(GeneralMatrix& tdx, const vector&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 \ No newline at end of file diff --git a/mex/sources/estimation/ModelSolution.cc b/mex/sources/estimation/ModelSolution.cc index 1d0a79f50..67ef6a305 100644 --- a/mex/sources/estimation/ModelSolution.cc +++ b/mex/sources/estimation/ModelSolution.cc @@ -32,15 +32,22 @@ */ ModelSolution::ModelSolution(const std::string& modName, size_t n_endo_arg, size_t n_exo_arg, const std::vector& zeta_fwrd_arg, const std::vector& zeta_back_arg, const std::vector& zeta_mixed_arg, - const std::vector& zeta_static_arg, const Matrix& llincidence, double INqz_criterium) + const std::vector& 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); diff --git a/mex/sources/estimation/ModelSolution.hh b/mex/sources/estimation/ModelSolution.hh index 7ae4cacc4..69ef4d737 100644 --- a/mex/sources/estimation/ModelSolution.hh +++ b/mex/sources/estimation/ModelSolution.hh @@ -40,7 +40,7 @@ class ModelSolution{ public: ModelSolution(const std::string& modName, size_t n_endo, size_t n_exo, const std::vector& zeta_fwrd_arg, const std::vector& zeta_back_arg, const std::vector& zeta_mixed_arg, - const std::vector& zeta_static_arg, const Matrix& llincidence, double qz_criterium); + const std::vector& 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 zeta_fwrd_mixed, zeta_back_mixed; Matrix jacobian; Vector residual; Matrix Mx; diff --git a/mex/sources/estimation/Prior.hh b/mex/sources/estimation/Prior.hh new file mode 100644 index 000000000..e04ab7595 --- /dev/null +++ b/mex/sources/estimation/Prior.hh @@ -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 . + */ + +/////////////////////////////////////////////////////////// +// 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_) diff --git a/mex/sources/estimation/dr1.cpp b/mex/sources/estimation/dr1.cpp deleted file mode 100644 index f3a7f81d4..000000000 --- a/mex/sources/estimation/dr1.cpp +++ /dev/null @@ -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 . - */ - -/*********************************************** -% 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 - vectorspan2nx(sss-exo_nbr); - for (i=0;ispannx2end(exo_nbr); - for (i=0;i1 - { - 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); - vectorkk(0); - int maximum_endo_lag= (int)dynareParams.getDoubleField(string("maximum_endo_lag")); - for(i=0;ikk2(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; - vectork1klag(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;igetParams()=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 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::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::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() - diff --git a/mex/sources/estimation/dynare_resolve.cpp b/mex/sources/estimation/dynare_resolve.cpp deleted file mode 100644 index 5bbfb62d8..000000000 --- a/mex/sources/estimation/dynare_resolve.cpp +++ /dev/null @@ -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 . -*/ -/***************************************************** -% 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&iv,vector&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); - //aux(:,2) = aux(:,2) + nstatic; - vectork(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(); - vectori_n_iv(n_iv); - for (int i=0;i 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;ispan2end(A.numRows()-n_iv); - for (int i=0;i. -*/ -#ifndef IOUTILS_H -#define IOUTILS_H -#include "GeneralMatrix.h" -#include "SylvException.h" -#include -#include -//using namespace std; - - -struct charArraySt - { - char ** charArrayPtr; - int len; - }; - -class GeneralParams - { -// map 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& - getStringVectorField(const string& field)=0; - virtual vector& - 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&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 diff --git a/mex/sources/estimation/libmat/VDVEigDecomposition.cc b/mex/sources/estimation/libmat/VDVEigDecomposition.cc new file mode 100644 index 000000000..488e0f400 --- /dev/null +++ b/mex/sources/estimation/libmat/VDVEigDecomposition.cc @@ -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 . + */ + +#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; +} diff --git a/mex/sources/estimation/libmat/VDVEigDecomposition.hh b/mex/sources/estimation/libmat/VDVEigDecomposition.hh new file mode 100644 index 000000000..54e08e63d --- /dev/null +++ b/mex/sources/estimation/libmat/VDVEigDecomposition.hh @@ -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 . + */ + +/** + * Calculates vector of Eigen values D and matrix of Eigen vectors V + * of input matrix m + */ + +#if !defined(VDVEigDec_INCLUDE) +#define VDVEigDec_INCLUDE + +#include +#include "Vector.hh" +#include "Matrix.hh" +#include + +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 + 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 +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 diff --git a/mex/sources/estimation/tests/DsgeLikelihood_mexf.cpp b/mex/sources/estimation/tests/DsgeLikelihood_mexf.cpp deleted file mode 100644 index d2c67ce74..000000000 --- a/mex/sources/estimation/tests/DsgeLikelihood_mexf.cpp +++ /dev/null @@ -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 . - */ - - -// 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&mfys=dynareParams.getIntVectorField("mfys"); - vector&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&order_var = dr.getIntVectorField(string("order_var")); -#ifdef DEBUG - for ( gg=0;gg&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&iv= dynareParams.getIntVectorField(string("restrict_var_list")); - vector&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= 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& vll=dl.getLikVector(); - plhs[5] = mxCreateDoubleMatrix(nper,1, mxREAL); - double * mxll= mxGetPr(plhs[5]); - // assign likelihood array - for (int j=0;j