diff --git a/matlab/get_the_name.m b/matlab/get_the_name.m
index e625d1bd4..b3cf31b9a 100644
--- a/matlab/get_the_name.m
+++ b/matlab/get_the_name.m
@@ -53,7 +53,7 @@ elseif k <= (nvx+nvn)
vname = deblank(options_.varobs(estim_params_.var_endo(k-estim_params_.nvx,1),:));
nam=['SE_EOBS_',vname];
if TeX
- tname = deblank(options_.TeX_varobs(estim_params_.var_endo(k-estim_params_.nvx,1),:));
+ tname = deblank(M_.endo_names_tex(estim_params_.var_endo(k-estim_params_.nvx,1),:));
texnam = ['$ SE_{' tname '} $'];
end
elseif k <= (nvx+nvn+ncx)
diff --git a/matlab/mult_elimination.m b/matlab/mult_elimination.m
index f5f11d009..8c9fbf53b 100644
--- a/matlab/mult_elimination.m
+++ b/matlab/mult_elimination.m
@@ -1,120 +1,124 @@
-function dr=mult_elimination(varlist,M_, options_, oo_)
-% function mult_elimination()
-% replaces Lagrange multipliers in Ramsey policy by lagged value of state
-% and shock variables
-%
-% INPUT
-% none
-%
-% OUTPUT
-% dr: a structure with the new decision rule
-%
-% SPECIAL REQUIREMENTS
-% none
-
-% Copyright (C) 2003-2008 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 .
-
-dr = oo_.dr;
-
-nstatic = dr.nstatic;
-npred = dr.npred;
-order_var = dr.order_var;
-nstates = M_.endo_names(order_var(nstatic+(1:npred)),:);
-
-il = strmatch('mult_',nstates);
-nil = setdiff(1:dr.npred,il);
-m_nbr = length(il);
-nm_nbr = length(nil);
-
-AA1 = dr.ghx(:,nil);
-AA2 = dr.ghx(:,il);
-A1 = dr.ghx(nstatic+(1:npred),nil);
-A2 = dr.ghx(nstatic+(1:npred),il);
-B = dr.ghu(nstatic+(1:npred),:);
-A11 = A1(nil,:);
-A21 = A1(il,:);
-A12 = A2(nil,:);
-A22 = A2(il,:);
-
-[Q1,R1,E1] = qr(A2);
-n1 = sum(abs(diag(R1)) > 1e-8);
-
-Q1_12 = Q1(1:nm_nbr,n1+1:end);
-Q1_22 = Q1(nm_nbr+1:end,n1+1:end);
-[Q2,R2,E2] = qr(Q1_22');
-n2 = sum(abs(diag(R2)) > 1e-8);
-
-R2_1 = inv(R2(1:n2,1:n2));
-
-M1(order_var,:) = AA1 - AA2*E2*[R2_1*Q2(:,1:n2)'*Q1_12'; zeros(m_nbr-n2,nm_nbr)];
-M2(order_var,:) = AA2*E2*[R2_1*Q2(:,1:n2)'*[Q1_12' Q1_22']*A1; zeros(m_nbr-n2,length(nil))];
-M3(order_var,:) = dr.ghu;
-M4(order_var,:) = AA2*E2*[R2_1*Q2(:,1:n2)'*[Q1_12' Q1_22']*B; zeros(m_nbr-n2,size(B,2))];
-
-endo_nbr = M_.orig_model.endo_nbr;
-exo_nbr = M_.exo_nbr;
-
-lead_lag_incidence = M_.lead_lag_incidence(:,1:endo_nbr+exo_nbr);
-lead_lag_incidence1 = lead_lag_incidence(1,:) > 0;
-maximum_lag = M_.maximum_lag;
-for i=1:maximum_lag-1
- lead_lag_incidence1 = [lead_lag_incidence1; lead_lag_incidence(i,:)| ...
- lead_lag_incidence(i+1,:)];
-end
-lead_lag_incidence1 = [lead_lag_incidence1; ...
- lead_lag_incidence(M_.maximum_lag,:) > 0];
-k = find(lead_lag_incidence1');
-lead_lag_incidence1 = zeros(size(lead_lag_incidence1'));
-lead_lag_incidence1(k) = 1:length(k);
-lead_lag_incidence1 = lead_lag_incidence1';
-
-kstate = zeros(0,2);
-for i=maximum_lag:-1:1
- k = find(lead_lag_incidence(i,:));
- kstate = [kstate; [k' repmat(i+1,length(k),1)]];
-end
-
-dr.M1 = M1;
-dr.M2 = M2;
-dr.M3 = M3;
-dr.M4 = M4;
-
-nvar = length(varlist);
-nspred = dr.nspred;
-nspred = 6;
-if nvar > 0
- res_table = zeros(2*(nspred+M_.exo_nbr),nvar);
- headers = 'Variables';
- for i=1:length(varlist)
- k = strmatch(varlist{i},M_.endo_names(dr.order_var,:),'exact');
- headers = strvcat(headers,varlist{i});
-
- res_table(1:nspred,i) = M1(k,:)';
- res_table(nspred+(1:nspred),i) = M2(k,:)';
- res_table(2*nspred+(1:M_.exo_nbr),i) = M3(k,:)';
- res_table(2*nspred+M_.exo_nbr+(1:M_.exo_nbr),i) = M4(k,:)';
- end
-
- my_title='ELIMINATION OF THE MULTIPLIERS';
- lab1 = M_.endo_names(dr.order_var(dr.nstatic+[ 1 2 5:8]),:);
- labels = strvcat(lab1,lab1,M_.exo_names,M_.exo_names);
- lh = size(labels,2)+2;
- dyntable(my_title,headers,labels,res_table,lh,10,6);
- disp(' ')
-end
-
+function dr=mult_elimination(varlist,M_, options_, oo_)
+% function mult_elimination()
+% replaces Lagrange multipliers in Ramsey policy by lagged value of state
+% and shock variables
+%
+% INPUT
+% none
+%
+% OUTPUT
+% dr: a structure with the new decision rule
+%
+% SPECIAL REQUIREMENTS
+% none
+
+% Copyright (C) 2003-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 .
+
+ dr = oo_.dr;
+
+ nstatic = dr.nstatic;
+ npred = dr.npred;
+ order_var = dr.order_var;
+ nstates = M_.endo_names(order_var(nstatic+(1:npred)),:);
+
+ il = strmatch('mult_',nstates);
+ nil = setdiff(1:dr.npred,il);
+ m_nbr = length(il);
+ nm_nbr = length(nil);
+
+ AA1 = dr.ghx(:,nil);
+ AA2 = dr.ghx(:,il);
+ A1 = dr.ghx(nstatic+(1:npred),nil);
+ A2 = dr.ghx(nstatic+(1:npred),il);
+ B = dr.ghu(nstatic+(1:npred),:);
+ A11 = A1(nil,:);
+ A21 = A1(il,:);
+ A12 = A2(nil,:);
+ A22 = A2(il,:);
+ B1 = B(nil,:);
+ B2 = B(il,:);
+
+ [Q1,R1,E1] = qr([A12; A22]);
+ n1 = sum(abs(diag(R1)) > 1e-8);
+
+ Q1_12 = Q1(1:nm_nbr,n1+1:end);
+ Q1_22 = Q1(nm_nbr+(1:m_nbr),n1+1:end);
+ [Q2,R2,E2] = qr(Q1_22');
+ n2 = sum(abs(diag(R2)) > 1e-8);
+
+ R2_1 = inv(R2(1:n2,1:n2));
+
+ M1 = AA1 - AA2*E2*[R2_1*Q2(:,1:n2)'*Q1_12'; zeros(m_nbr-n2,nm_nbr)];
+ M2 = AA2*E2*[R2_1*Q2(:,1:n2)'*[Q1_12' Q1_22']*[A11;A21]; zeros(m_nbr-n2,length(nil))];
+ M3 = dr.ghu;
+ M4 = AA2*E2*[R2_1*Q2(:,1:n2)'*[Q1_12' Q1_22']*[B1;B2]; zeros(m_nbr-n2,size(B,2))];
+
+ k1 = nstatic+(1:npred);
+ k1 = k1(nil);
+
+ endo_nbr = M_.orig_model.endo_nbr;
+ exo_nbr = M_.exo_nbr;
+
+ lead_lag_incidence = M_.lead_lag_incidence(:,1:endo_nbr+exo_nbr);
+ lead_lag_incidence1 = lead_lag_incidence(1,:) > 0;
+ maximum_lag = M_.maximum_lag;
+ for i=1:maximum_lag-1
+ lead_lag_incidence1 = [lead_lag_incidence1; lead_lag_incidence(i,:)| ...
+ lead_lag_incidence(i+1,:)];
+ end
+ lead_lag_incidence1 = [lead_lag_incidence1; ...
+ lead_lag_incidence(M_.maximum_lag,:) > 0];
+ k = find(lead_lag_incidence1');
+ lead_lag_incidence1 = zeros(size(lead_lag_incidence1'));
+ lead_lag_incidence1(k) = 1:length(k);
+ lead_lag_incidence1 = lead_lag_incidence1';
+
+ kstate = zeros(0,2);
+ for i=maximum_lag:-1:1
+ k = find(lead_lag_incidence(i,:));
+ kstate = [kstate; [k' repmat(i+1,length(k),1)]];
+ end
+
+ dr.M1 = M1;
+ dr.M2 = M2;
+ dr.M3 = M3;
+ dr.M4 = M4;
+
+ nvar = length(varlist);
+ nspred = dr.nspred;
+ nspred = 6;
+ if nvar > 0
+ res_table = zeros(2*(nspred+M_.exo_nbr),nvar);
+ headers = 'Variables';
+ for i=1:length(varlist)
+ k = strmatch(varlist{i},M_.endo_names(dr.order_var,:),'exact');
+ headers = strvcat(headers,varlist{i});
+
+ res_table(1:nspred,i) = M1(k,:)';
+ res_table(nspred+(1:nspred),i) = M2(k,:)';
+ res_table(2*nspred+(1:M_.exo_nbr),i) = M3(k,:)';
+ res_table(2*nspred+M_.exo_nbr+(1:M_.exo_nbr),i) = M4(k,:)';
+ end
+
+ my_title='ELIMINATION OF THE MULTIPLIERS';
+ lab1 = M_.endo_names(dr.order_var(dr.nstatic+[ 1 2 5:8]),:);
+ labels = strvcat(lab1,lab1,M_.exo_names,M_.exo_names);
+ lh = size(labels,2)+2;
+ dyntable(my_title,headers,labels,res_table,lh,10,6);
+ disp(' ')
+ end
diff --git a/mex/sources/estimation/DetrendData.cc b/mex/sources/estimation/DetrendData.cc
index ebd1d282b..6e045f5f6 100644
--- a/mex/sources/estimation/DetrendData.cc
+++ b/mex/sources/estimation/DetrendData.cc
@@ -31,7 +31,7 @@ DetrendData::DetrendData(const bool INlogLinear) //, Vector& INtrendCoeff)
};
void
-DetrendData::detrend(const Vector &SteadyState, const MatrixView &dataView, Matrix &Y)
+DetrendData::detrend(const VectorView &SteadyState, const MatrixConstView &dataView, Matrix &Y)
{
// dummy
Y=dataView;
diff --git a/mex/sources/estimation/DetrendData.hh b/mex/sources/estimation/DetrendData.hh
index cf1ee8681..674e2ec81 100644
--- a/mex/sources/estimation/DetrendData.hh
+++ b/mex/sources/estimation/DetrendData.hh
@@ -33,7 +33,7 @@ class DetrendData
public:
virtual ~DetrendData(){};
DetrendData(const bool logLinear); // add later Vector& trendCoeff);
- void detrend(const Vector &SteadyState, const MatrixView &dataView, Matrix &Y);
+ void detrend(const VectorView &SteadyState, const MatrixConstView &dataView, Matrix &Y);
private:
const bool logLinear;
diff --git a/mex/sources/estimation/EstimatedParameter.cc b/mex/sources/estimation/EstimatedParameter.cc
new file mode 100644
index 000000000..b66b576f6
--- /dev/null
+++ b/mex/sources/estimation/EstimatedParameter.cc
@@ -0,0 +1,40 @@
+/*
+ * 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.cpp
+// Implementation of the Class EstimatedParameter
+// Created on: 02-Feb-2010 13:06:35
+///////////////////////////////////////////////////////////
+
+#include "EstimatedParameter.hh"
+
+EstimatedParameter::EstimatedParameter(const EstimatedParameter::pType type_arg,
+ size_t ID1_arg, size_t ID2_arg, const std::vector &subSampleIDs_arg,
+ double lower_bound_arg, double upper_bound_arg, Prior prior_arg) :
+ ptype(type_arg), ID1(ID1_arg), ID2(ID2_arg),
+ lower_bound(lower_bound_arg), upper_bound(upper_bound_arg), prior(prior_arg),
+ subSampleIDs(subSampleIDs_arg)
+{
+}
+
+EstimatedParameter::~EstimatedParameter()
+{
+}
+
diff --git a/mex/sources/estimation/EstimatedParameter.hh b/mex/sources/estimation/EstimatedParameter.hh
index 596e4c18f..3a83b6270 100644
--- a/mex/sources/estimation/EstimatedParameter.hh
+++ b/mex/sources/estimation/EstimatedParameter.hh
@@ -27,18 +27,11 @@
#include "Prior.hh"
#include
+#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;
+ // parameter types
enum pType
{
shock_SD = 1, // standard deviation of a structural shock
@@ -47,9 +40,18 @@ public:
measureErr_Corr = 4, // correlation between two measurement errors
deepPar = 5 // deep parameter
};
- const pType ptype;
- const double lower_boound;
- const double upper_boound;
+
+ EstimatedParameter(const EstimatedParameter::pType type,
+ size_t ID1, size_t ID2, const std::vector &subSampleIDs,
+ double lower_bound, double upper_bound, Prior prior
+ );
+ virtual ~EstimatedParameter();
+
+ const enum pType ptype;
+ const size_t ID1;
+ const size_t ID2;
+ const double lower_bound;
+ const double upper_bound;
const Prior prior;
const std::vector subSampleIDs;
};
diff --git a/mex/sources/estimation/EstimatedParametersDescription.cc b/mex/sources/estimation/EstimatedParametersDescription.cc
new file mode 100644
index 000000000..6141f6b46
--- /dev/null
+++ b/mex/sources/estimation/EstimatedParametersDescription.cc
@@ -0,0 +1,36 @@
+/*
+ * 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 .
+ */
+
+///////////////////////////////////////////////////////////
+// EstimatedParametersDescription.cpp
+// Implementation of the Class EstimatedParametersDescription
+// Created on: 02-Feb-2010 13:06:47
+///////////////////////////////////////////////////////////
+
+#include "EstimatedParametersDescription.hh"
+
+EstimatedParametersDescription::~EstimatedParametersDescription()
+{
+}
+
+EstimatedParametersDescription::EstimatedParametersDescription(std::vector &INestSubsamples, std::vector &INestParams) :
+ estSubsamples(INestSubsamples), estParams(INestParams)
+{
+}
+
diff --git a/mex/sources/estimation/EstimatedParametersDescription.hh b/mex/sources/estimation/EstimatedParametersDescription.hh
index bf1a0bb6b..8d6aaf009 100644
--- a/mex/sources/estimation/EstimatedParametersDescription.hh
+++ b/mex/sources/estimation/EstimatedParametersDescription.hh
@@ -46,21 +46,9 @@ 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:
-
+ EstimatedParametersDescription(std::vector &estSubsamples, std::vector &estParams);
std::vector estSubsamples;
std::vector estParams;
-//vector prior;
};
#endif // !defined(E8F2C096_A301_42e8_80FF_A643291BF995__INCLUDED_)
diff --git a/mex/sources/estimation/Estimation.h b/mex/sources/estimation/EstimationSubsample.cc
similarity index 55%
rename from mex/sources/estimation/Estimation.h
rename to mex/sources/estimation/EstimationSubsample.cc
index 508383fa6..f8e937377 100644
--- a/mex/sources/estimation/Estimation.h
+++ b/mex/sources/estimation/EstimationSubsample.cc
@@ -1,44 +1,36 @@
-/*
- * 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 ESTIMATION_H
-#define ESTIMATION_H
-//#include "k_ord_dynare.h"
-#include "kalman.h"
-#include "math.h"
-#include "disclyap_fast.h"
-#include "Mexutils.h"
-#include "estutils.h"
-
-#include
-
-#include
-
-class DsgeLikelihood;//GeneralParams;
-
-class Estimation
- {
- double logLikelihood;
- //instance MexStruct
- GeneralParams& gParams;
- DsgeLikelihood& DsgeLik;
-public:
- double getLikelihood()
- {return logLikelihood;};
- };
-#endif
\ No newline at end of file
+/*
+ * 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.cpp
+// Implementation of the Class EstimationSubsample
+// Created on: 02-Feb-2010 13:06:01
+///////////////////////////////////////////////////////////
+
+#include "EstimationSubsample.hh"
+
+EstimationSubsample::EstimationSubsample(size_t INstartPeriod, size_t INendPeriod) :
+ startPeriod(INstartPeriod), endPeriod(INendPeriod)
+{
+}
+
+EstimationSubsample::~EstimationSubsample()
+{
+}
+
diff --git a/mex/sources/estimation/InitializeKalmanFilter.cc b/mex/sources/estimation/InitializeKalmanFilter.cc
index da16e46f5..d55064fce 100644
--- a/mex/sources/estimation/InitializeKalmanFilter.cc
+++ b/mex/sources/estimation/InitializeKalmanFilter.cc
@@ -33,46 +33,48 @@ 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,
- 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),
+ const std::vector &zeta_varobs_back_mixed_arg,
+ double qz_criterium_arg,
+ double lyapunov_tol_arg, int &info) :
+ lyapunov_tol(lyapunov_tol_arg),
+ zeta_varobs_back_mixed(zeta_varobs_back_mixed_arg),
detrendData(false), modelSolution(modName, n_endo_arg, n_exo_arg, zeta_fwrd_arg, zeta_back_arg,
- 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),
- Rt(n_exo_arg, riv.size()), RQ(riv.size(), n_exo_arg)
+ zeta_mixed_arg, zeta_static_arg, qz_criterium_arg),
+ discLyapFast(zeta_varobs_back_mixed.size()),
+ g_x(n_endo_arg, zeta_back_arg.size() + zeta_mixed_arg.size()),
+ g_u(n_endo_arg, n_exo_arg),
+ Rt(n_exo_arg, zeta_varobs_back_mixed.size()),
+ RQ(zeta_varobs_back_mixed.size(), n_exo_arg)
{
- size_t n_pred = ghx.getCols();
- size_t n_static = zeta_static_arg.size();
- size_t j = 0;
- for (size_t i = 0; i < n_endo_arg; ++i)
- if (inv_order_var_arg[i]-n_static < n_pred && inv_order_var_arg[i]-n_static >= 0)
- inv_ric[j++] = ricIN[inv_order_var_arg[i]-n_static];
+ std::vector zeta_back_mixed;
+ set_union(zeta_back_arg.begin(), zeta_back_arg.end(),
+ zeta_mixed_arg.begin(), zeta_mixed_arg.end(),
+ back_inserter(zeta_back_mixed));
+ for (size_t i = 0; i < zeta_back_mixed.size(); i++)
+ pi_bm_vbm.push_back(find(zeta_varobs_back_mixed.begin(), zeta_varobs_back_mixed.end(),
+ zeta_back_mixed[i]) - zeta_varobs_back_mixed.begin());
+
}
// initialise parameter dependent KF matrices only but not Ps
void
-InitializeKalmanFilter::initialize(Vector &steadyState, const Vector &deepParams, Matrix &R, const Matrix &Z,
+InitializeKalmanFilter::initialize(VectorView &steadyState, const Vector &deepParams, Matrix &R,
const Matrix &Q, Matrix &RQRt, Matrix &T,
- double &penalty, const MatrixView &dataView, Matrix &Y, int &info)
+ double &penalty, const MatrixConstView &dataView, Matrix &Y, int &info)
{
- modelSolution.compute(steadyState, deepParams, ghx_raw, ghu_raw);
+ modelSolution.compute(steadyState, deepParams, g_x, g_u);
detrendData.detrend(steadyState, dataView, Y);
- mat::reorderRowsByVectors(ghx, mat::nullVec, ghx_raw, order_var);
- mat::reorderRowsByVectors(ghu, mat::nullVec, ghu_raw, order_var);
-
setT(T, info);
setRQR(R, Q, RQRt, info);
}
// initialise all KF matrices
void
-InitializeKalmanFilter::initialize(Vector &steadyState, const Vector &deepParams, Matrix &R, const Matrix &Z,
+InitializeKalmanFilter::initialize(VectorView &steadyState, const Vector &deepParams, Matrix &R,
const Matrix &Q, Matrix &RQRt, Matrix &T, Matrix &Pstar, Matrix &Pinf,
- double &penalty, const MatrixView &dataView, Matrix &Y, int &info)
+ double &penalty, const MatrixConstView &dataView, Matrix &Y, int &info)
{
- initialize(steadyState, deepParams, R, Z, Q, RQRt, T, penalty, dataView, Y, info);
+ initialize(steadyState, deepParams, R, Q, RQRt, T, penalty, dataView, Y, info);
setPstar(Pstar, Pinf, T, RQRt, info);
}
@@ -80,26 +82,19 @@ InitializeKalmanFilter::initialize(Vector &steadyState, const Vector &deepParams
void
InitializeKalmanFilter::setT(Matrix &T, int &info)
{
- // here is the content of [T R]=[A,B] = kalman_transition_matrix(oo_.dr,iv,ic,aux,M_.exo_nbr);
-
+ // Initialize the empty columns of T to zero
T.setAll(0.0);
-
- //T(i_n_iv,ic) = dr.ghx(iv,:);
- mat::assignByVectors(T, mat::nullVec, inv_ric, ghx, riv, mat::nullVec);
+ mat::assignByVectors(T, mat::nullVec, pi_bm_vbm, g_x, zeta_varobs_back_mixed, mat::nullVec);
}
void
InitializeKalmanFilter::setRQR(Matrix &R, const Matrix &Q, Matrix &RQRt, int &info)
{
- R.setAll(0.0);
- //B(i_n_iv,:) = dr.ghu(iv,:); and R=B;
- mat::assignByVectors(R, mat::nullVec, mat::nullVec, ghu, riv, mat::nullVec);
+ mat::assignByVectors(R, mat::nullVec, mat::nullVec, g_u, zeta_varobs_back_mixed, mat::nullVec);
// Matrix RQRt=R*Q*R'
- RQ.setAll(0.0);
- blas::gemm("N", "N", 1.0, R, Q, 1.0, RQ); // R*Q
- RQRt.setAll(0.0);
- blas::gemm("N", "T", 1.0, RQ, R, 1.0, RQRt); // R*Q*R'
+ blas::gemm("N", "N", 1.0, R, Q, 0.0, RQ); // R*Q
+ blas::gemm("N", "T", 1.0, RQ, R, 0.0, RQRt); // R*Q*R'
}
void
diff --git a/mex/sources/estimation/InitializeKalmanFilter.hh b/mex/sources/estimation/InitializeKalmanFilter.hh
index 27af85c8d..0383fc7c2 100644
--- a/mex/sources/estimation/InitializeKalmanFilter.hh
+++ b/mex/sources/estimation/InitializeKalmanFilter.hh
@@ -40,29 +40,32 @@ class InitializeKalmanFilter
{
public:
+ /*!
+ \param[in] zeta_varobs_back_mixed_arg The union of indices of observed, backward and mixed variables
+ */
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,
- 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);
+ const std::vector &zeta_varobs_back_mixed_arg,
+ double qz_criterium_arg, double lyapunov_tol_arg, int &info);
virtual ~InitializeKalmanFilter();
// initialise all KF matrices
- void initialize(Vector &steadyState, const Vector &deepParams, Matrix &R, const Matrix &Z, const Matrix &Q, Matrix &RQRt,
- Matrix &T, Matrix &Pstar, Matrix &Pinf, double &penalty, const MatrixView &dataView, Matrix &Y, int &info);
+ void initialize(VectorView &steadyState, const Vector &deepParams, Matrix &R, const Matrix &Q, Matrix &RQRt,
+ Matrix &T, Matrix &Pstar, Matrix &Pinf, double &penalty, const MatrixConstView &dataView, Matrix &Y, int &info);
// initialise parameter dependent KF matrices only but not Ps
- void initialize(Vector &steadyState, const Vector &deepParams, Matrix &R, const Matrix &Z, const Matrix &Q, Matrix &RQRt,
- Matrix &T, double &penalty, const MatrixView &dataView, Matrix &Y, int &info);
+ void initialize(VectorView &steadyState, const Vector &deepParams, Matrix &R, const Matrix &Q, Matrix &RQRt,
+ Matrix &T, double &penalty, const MatrixConstView &dataView, Matrix &Y, int &info);
private:
- const std::vector riv; // restrict_var_list
- std::vector inv_ric; // inverse restrict_columns
- const std::vector order_var;
const double lyapunov_tol;
+ const std::vector zeta_varobs_back_mixed;
+ //! Indices of back+mixed zetas inside varobs+back+mixed zetas
+ std::vector pi_bm_vbm;
DetrendData detrendData;
ModelSolution modelSolution;
DiscLyapFast discLyapFast; //Lyapunov solver
- Matrix ghx, ghx_raw;
- Matrix ghu, ghu_raw;
+ Matrix g_x;
+ Matrix g_u;
Matrix Rt, RQ;
void setT(Matrix &T, int &info);
void setRQR(Matrix &R, const Matrix &Q, Matrix &RQRt, int &info);
diff --git a/mex/sources/estimation/KalmanFilter.cc b/mex/sources/estimation/KalmanFilter.cc
index 2ed8b40f6..47b082c83 100644
--- a/mex/sources/estimation/KalmanFilter.cc
+++ b/mex/sources/estimation/KalmanFilter.cc
@@ -33,42 +33,58 @@ 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,
- 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) :
- Z(varobs_arg.size(), riv.size()), T(riv.size()), R(riv.size(), n_exo),
- Pstar(riv.size(), riv.size()), Pinf(riv.size(), riv.size()), RQRt(riv.size(), riv.size()),
- Ptmp(riv.size(), riv.size()), F(varobs_arg.size(), varobs_arg.size()),
- Finv(varobs_arg.size(), varobs_arg.size()), K(riv.size(), varobs_arg.size()), KFinv(riv.size(), varobs_arg.size()),
- 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),
+ double qz_criterium_arg, const std::vector &varobs_arg,
+ double riccati_tol_arg, double lyapunov_tol_arg, int &info) :
+ zeta_varobs_back_mixed(compute_zeta_varobs_back_mixed(zeta_back_arg, zeta_mixed_arg, varobs_arg)),
+ Z(varobs_arg.size(), zeta_varobs_back_mixed.size()), T(zeta_varobs_back_mixed.size()), R(zeta_varobs_back_mixed.size(), n_exo),
+ Pstar(zeta_varobs_back_mixed.size(), zeta_varobs_back_mixed.size()), Pinf(zeta_varobs_back_mixed.size(), zeta_varobs_back_mixed.size()), RQRt(zeta_varobs_back_mixed.size(), zeta_varobs_back_mixed.size()),
+ Ptmp(zeta_varobs_back_mixed.size(), zeta_varobs_back_mixed.size()), F(varobs_arg.size(), varobs_arg.size()),
+ Finv(varobs_arg.size(), varobs_arg.size()), K(zeta_varobs_back_mixed.size(), varobs_arg.size()), KFinv(zeta_varobs_back_mixed.size(), varobs_arg.size()),
+ oldKFinv(zeta_varobs_back_mixed.size(), varobs_arg.size()), Finverter(varobs_arg.size()), a_init(zeta_varobs_back_mixed.size(), 1),
+ a_new(zeta_varobs_back_mixed.size(), 1), vt(varobs_arg.size(), 1), vtFinv(1, varobs_arg.size()), vtFinvVt(1), riccati_tol(riccati_tol_arg),
initKalmanFilter(modName, n_endo, n_exo, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg,
- zeta_static_arg, qz_criterium, order_var, inv_order_var, riv, ric, lyapunov_tol, info)
+ zeta_static_arg, zeta_varobs_back_mixed, qz_criterium_arg, lyapunov_tol_arg, info)
{
a_init.setAll(0.0);
Z.setAll(0.0);
for (size_t i = 0; i < varobs_arg.size(); ++i)
- Z(i, varobs_arg[i]) = 1.0;
+ {
+ size_t j = find(zeta_varobs_back_mixed.begin(), zeta_varobs_back_mixed.end(),
+ varobs_arg[i]) - zeta_varobs_back_mixed.begin();
+ Z(i, j) = 1.0;
+ }
+}
+std::vector
+KalmanFilter::compute_zeta_varobs_back_mixed(const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &varobs_arg)
+{
+ std::vector varobs_sorted = varobs_arg;
+ sort(varobs_sorted.begin(), varobs_sorted.end());
+ std::vector zeta_back_mixed, zeta_varobs_back_mixed;
+ set_union(zeta_back_arg.begin(), zeta_back_arg.end(),
+ zeta_mixed_arg.begin(), zeta_mixed_arg.end(),
+ back_inserter(zeta_back_mixed));
+ set_union(zeta_back_mixed.begin(), zeta_back_mixed.end(),
+ varobs_arg.begin(), varobs_arg.end(),
+ back_inserter(zeta_varobs_back_mixed));
+ return zeta_varobs_back_mixed;
}
double
-KalmanFilter::compute(const MatrixView &dataView, Vector &steadyState,
+KalmanFilter::compute(const MatrixConstView &dataView, VectorView &steadyState,
const Matrix &Q, const Matrix &H, const Vector &deepParams,
VectorView &vll, size_t start, size_t period, double &penalty, int &info)
{
- double ll;
Matrix Y(dataView.getRows(), dataView.getCols()); // data
if(period==0) // initialise all KF matrices
- initKalmanFilter.initialize(steadyState, deepParams, R, Z, Q, RQRt, T, Pstar, Pinf,
+ initKalmanFilter.initialize(steadyState, deepParams, R, Q, RQRt, T, Pstar, Pinf,
penalty, dataView, Y, info);
else // initialise parameter dependent KF matrices only but not Ps
- initKalmanFilter.initialize(steadyState, deepParams, R, Z, Q, RQRt, T,
+ initKalmanFilter.initialize(steadyState, deepParams, R, Q, RQRt, T,
penalty, dataView, Y, info);
- return ll = filter(Y, H, vll, start, info);
+ return filter(Y, H, vll, start, info);
};
@@ -79,7 +95,7 @@ double
KalmanFilter::filter(const Matrix &dataView, const Matrix &H, VectorView &vll, size_t start, int &info)
{
double loglik=0.0, ll, logFdet, Fdet;
- int p = Finv.getRows();
+ size_t p = Finv.getRows();
bool nonstationary = true;
for (size_t t = 0; t < dataView.getCols(); ++t)
@@ -87,8 +103,7 @@ KalmanFilter::filter(const Matrix &dataView, const Matrix &H, VectorView &vll,
if (nonstationary)
{
// K=PZ'
- K.setAll(0.0);
- blas::gemm("N", "T", 1.0, Pstar, Z, 1.0, K);
+ blas::gemm("N", "T", 1.0, Pstar, Z, 0.0, K);
//F=ZPZ' +H = ZK+H
F = H;
@@ -99,11 +114,10 @@ KalmanFilter::filter(const Matrix &dataView, const Matrix &H, VectorView &vll,
mat::set_identity(Finv);
Finverter.invMult("N", F, Finv); // F now contains its LU decomposition!
// KFinv gain matrix
- KFinv.setAll(0.0);
- blas::gemm("N", "N", 1.0, K, Finv, 1.0, KFinv);
+ blas::gemm("N", "N", 1.0, K, Finv, 0.0, KFinv);
// deteminant of F:
Fdet = 1;
- for (int d = 0; d < p; ++d)
+ for (size_t d = 0; d < p; ++d)
Fdet *= (-F(d, d));
logFdet=log(fabs(Fdet));
@@ -114,8 +128,7 @@ KalmanFilter::filter(const Matrix &dataView, const Matrix &H, VectorView &vll,
blas::gemm("N", "T", -1.0, KFinv, K, 1.0, Ptmp);
// 2) Ptmp= T*Ptmp
Pstar = Ptmp;
- Ptmp.setAll(0.0);
- blas::gemm("N", "N", 1.0, T, Pstar, 1.0, Ptmp);
+ blas::gemm("N", "N", 1.0, T, Pstar, 0.0, Ptmp);
// 3) Pt+1= Ptmp*T' +RQR'
Pstar = RQRt;
blas::gemm("N", "T", 1.0, Ptmp, T, 1.0, Pstar);
@@ -125,7 +138,6 @@ KalmanFilter::filter(const Matrix &dataView, const Matrix &H, VectorView &vll,
}
// err= Yt - Za
- //VectorView yt(dataView.getData()+t*dataView.getRows(),dataView.getRows(),1); // current observation vector
MatrixConstView yt(dataView, 0, t, dataView.getRows(), 1); // current observation vector
vt = yt;
blas::gemm("N", "N", -1.0, Z, a_init, 1.0, vt);
diff --git a/mex/sources/estimation/KalmanFilter.hh b/mex/sources/estimation/KalmanFilter.hh
index 4edaa193d..70e295fa5 100644
--- a/mex/sources/estimation/KalmanFilter.hh
+++ b/mex/sources/estimation/KalmanFilter.hh
@@ -29,7 +29,7 @@
#include "InitializeKalmanFilter.hh"
/**
- * •vanilla Kalman filter without constant and with measurement error (use scalar
+ * Vanilla Kalman filter without constant and with measurement error (use scalar
* 0 when no measurement error).
* If multivariate filter is faster, do as in Matlab: start with multivariate
* filter and switch to univariate filter only in case of singularity
@@ -50,18 +50,16 @@ 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,
- 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);
+ double qz_criterium_arg, const std::vector &varobs_arg,
+ double riccati_tol_arg, double lyapunov_tol_arg, int &info);
- double compute(const MatrixView &dataView, Vector &steadyState,
+ double compute(const MatrixConstView &dataView, VectorView &steadyState,
const Matrix &Q, const Matrix &H, const Vector &deepParams,
VectorView &vll, size_t start, size_t period, double &penalty, int &info);
-protected:
-// static double calcStepLogLik(const PLUFact &Finv, const Vector &v);
-
private:
+ const std::vector zeta_varobs_back_mixed;
+ static std::vector compute_zeta_varobs_back_mixed(const std::vector &zeta_back_arg, const std::vector &zeta_mixed_arg, const std::vector &varobs_arg);
Matrix Z; //nob*mm matrix mapping endogeneous variables and observations
Matrix T; //mm*mm transition matrix of the state equation.
Matrix R; //mm*rr matrix, mapping structural innovations to state variables.
@@ -77,7 +75,9 @@ private:
Matrix vtFinv, vtFinvVt; // intermeiate observation error *Finv vector
double riccati_tol;
InitializeKalmanFilter initKalmanFilter; //Initialise KF matrices
- double filter(const Matrix &dataView, const Matrix &H, VectorView &vll, size_t start, int &info);
+
+ // Method
+ double filter(const Matrix &data, const Matrix &H, VectorView &vll, size_t start, int &info);
};
diff --git a/mex/sources/estimation/LogLikelihoodMain.cc b/mex/sources/estimation/LogLikelihoodMain.cc
new file mode 100644
index 000000000..c96af71e5
--- /dev/null
+++ b/mex/sources/estimation/LogLikelihoodMain.cc
@@ -0,0 +1,64 @@
+/*
+ * 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.cpp
+// Implementation of the Class LogLikelihoodMain
+// Created on: 02-Feb-2010 12:57:09
+///////////////////////////////////////////////////////////
+
+#include "LogLikelihoodMain.hh"
+
+LogLikelihoodMain::LogLikelihoodMain( //const Matrix &data_arg, Vector &deepParams_arg,
+ 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 &varobs, double riccati_tol, double lyapunov_tol, int &info_arg)
+
+ : estSubsamples(estiParDesc.estSubsamples),
+ logLikelihoodSubSample(modName, estiParDesc, n_endo, n_exo, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, qz_criterium,
+ varobs, riccati_tol, lyapunov_tol, info_arg),
+ //deepParams(deepParams_arg),
+ vll(estSubsamples[estSubsamples.size()-1].endPeriod) // time dimension size of data
+
+{
+
+}
+
+LogLikelihoodMain::~LogLikelihoodMain()
+{
+
+}
+
+double
+LogLikelihoodMain::compute(Matrix &steadyState, const Vector &estParams, Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H, size_t start, int &info)
+{
+ for (size_t i = 1; i <= (size_t) estSubsamples.size(); ++i)
+ {
+ VectorView vSteadyState = mat::get_col(steadyState, i);
+
+ MatrixConstView dataView(data, 0, estSubsamples[i].startPeriod,
+ dataView.getRows(), estSubsamples[i].endPeriod-estSubsamples[i].startPeriod+1);
+ VectorView vllView(vll, estSubsamples[i].startPeriod, estSubsamples[i].endPeriod-estSubsamples[i].startPeriod+1);
+ logLikelihood += logLikelihoodSubSample.compute(vSteadyState, dataView, estParams, deepParams,
+ Q, H, vllView, info, start, i);
+ }
+ return logLikelihood;
+};
+
diff --git a/mex/sources/estimation/LogLikelihoodMain.hh b/mex/sources/estimation/LogLikelihoodMain.hh
index 4f6582333..2f428b4e6 100644
--- a/mex/sources/estimation/LogLikelihoodMain.hh
+++ b/mex/sources/estimation/LogLikelihoodMain.hh
@@ -25,7 +25,6 @@
#if !defined(E126AEF5_AC28_400a_821A_3BCFD1BC4C22__INCLUDED_)
#define E126AEF5_AC28_400a_821A_3BCFD1BC4C22__INCLUDED_
-//#include "EstimatedParametersDescription.hh"
#include "LogLikelihoodSubSample.hh"
class LogLikelihoodMain {
@@ -33,26 +32,25 @@ 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;
+ Vector vll; // vector of all KF step likelihoods
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);
+ LogLikelihoodMain( //const Matrix &data, Vector &deepParams_arg, //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_arg, const std::vector &varobs_arg,
+ double riccati_tol_arg, double lyapunov_tol_arg, 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
+ /**
+ * Compute method Inputs:
+ * Matrix &steadyState; Matrix of sub-sample periods column-vectors of steady states, one column vectro for each sub-sample period
+ * vectors of deep deepParams and estimated estParams
+ * Matrix &data input data reference
+ * Q and H KF matrices of shock and measurement error varinaces and covariances
+ * KF logLikelihood calculation start period.
+ */
+ double compute(Matrix &steadyState, const Vector &estParams, Vector &deepParams, const MatrixConstView &data, Matrix &Q, Matrix &H, size_t presampleStart, int &info); // for calls from estimation and to set Steady State
Vector &
getVll()
diff --git a/mex/sources/estimation/LogLikelihoodSubSample.cc b/mex/sources/estimation/LogLikelihoodSubSample.cc
index 4193aae9d..1a29c38ed 100644
--- a/mex/sources/estimation/LogLikelihoodSubSample.cc
+++ b/mex/sources/estimation/LogLikelihoodSubSample.cc
@@ -23,28 +23,162 @@
// Created on: 14-Jan-2010 22:39:14
///////////////////////////////////////////////////////////
-#include "LogLikelihoodSubSample.hh"
+//#include "LogLikelihoodSubSample.hh"
+#include "LogLikelihoodMain.hh" // use ...Main.hh for testing only
+#include
+#include "LapackBindings.hh"
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 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) :
+LogLikelihoodSubSample::LogLikelihoodSubSample(const std::string &modName, EstimatedParametersDescription &INestiParDesc, 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 &varobs, double riccati_tol, double lyapunov_tol, int &INinfo) :
+ estiParDesc(INestiParDesc),
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)
+ varobs, riccati_tol, lyapunov_tol, INinfo), eigQ(n_exo), eigH(varobs.size()), info(INinfo)
{
};
double
-LogLikelihoodSubSample::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)
+LogLikelihoodSubSample::compute(VectorView &steadyState, const MatrixConstView &dataView, const Vector &estParams, Vector &deepParams,
+ Matrix &Q, Matrix &H, VectorView &vll, int &info, size_t start, size_t period)
{
- logLikelihood = kalmanFilter.compute(dataView, steadyState, Q, H, deepParams, vll, start, period, penalty, info);
+
+ updateParams(estParams, deepParams, Q, H, period);
+ if (info == 0)
+ logLikelihood = kalmanFilter.compute(dataView, steadyState, Q, H, deepParams, vll, start, period, penalty, info);
+ // else
+ // logLikelihood+=penalty;
+
return logLikelihood;
};
+void
+LogLikelihoodSubSample::updateParams(const Vector &estParams, Vector &deepParams,
+ Matrix &Q, Matrix &H, size_t period)
+{
+ size_t i, j, k, k1, k2;
+ int test;
+ bool found;
+ std::vector::const_iterator it;
+ info = 0;
+
+ for (i = 0; i < estParams.getSize(); ++i)
+ {
+ found = false;
+ it = find(estiParDesc.estParams[i].subSampleIDs.begin(),
+ estiParDesc.estParams[i].subSampleIDs.end(), period);
+ if (it != estiParDesc.estParams[i].subSampleIDs.end())
+ found = true;
+ if (found)
+ {
+ switch (estiParDesc.estParams[i].ptype)
+ {
+ case EstimatedParameter::shock_SD:
+ k = (int) estiParDesc.estParams[i].ID1;
+ Q(k, k) = estParams(i)*estParams(i);
+ break;
+
+ case EstimatedParameter::measureErr_SD:
+#ifdef DEBUG
+ mexPrintf("Setting of H var_endo\n");
+#endif
+ k = (int) estiParDesc.estParams[i].ID1;
+ H(k, k) = estParams(i)*estParams(i);
+ break;
+
+ case EstimatedParameter::shock_Corr:
+#ifdef DEBUG
+ mexPrintf("Setting of Q corrx\n");
+#endif
+ k1 = (int) estiParDesc.estParams[i].ID1;
+ k2 = (int) estiParDesc.estParams[i].ID2;
+ Q(k1, k2) = estParams(i)*sqrt(Q(k1, k1)*Q(k2, k2));
+ Q(k2, k1) = Q(k1, k2);
+ // [CholQ,testQ] = chol(Q);
+ test = lapack::choleskyDecomp(Q, "L");
+ if (test > 0)
+ {
+ mexPrintf("Caugth unhandled exception with cholesky of Q matrix: ");
+ logLikelihood = penalty;
+ info = 1;
+ }
+ else if (test < 0)
+ {
+ // 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;
+ eigQ.calculate(Q); // get eigenvalues
+ //k = find(a < 0);
+ if (eigQ.hasConverged())
+ {
+ const Vector &evQ = eigQ.getD();
+ for (i = 0; i < evQ.getSize(); ++i)
+ if (evQ(i) < 0)
+ delta -= evQ(i);
+ }
+
+ logLikelihood = penalty+delta;
+ info = 43;
+ } // if
+ break;
+
+ case EstimatedParameter::measureErr_Corr:
+#ifdef DEBUG
+ mexPrintf("Setting of H corrn\n");
+#endif
+ k1 = (int) estiParDesc.estParams[i].ID1;
+ k2 = (int) estiParDesc.estParams[i].ID2;
+ // H(k1,k2) = xparam1(i)*sqrt(H(k1,k1)*H(k2,k2));
+ // H(k2,k1) = H(k1,k2);
+ H(k1, k2) = estParams(i)*sqrt(H(k1, k1)*H(k2, k2));
+ H(k2, k1) = H(k1, k2);
+
+ //[CholH,testH] = chol(H);
+ test = lapack::choleskyDecomp(H, "L");
+ if (test > 0)
+ {
+ mexPrintf("Caugth unhandled exception with cholesky of Q matrix: ");
+ logLikelihood = penalty;
+ info = 1;
+ }
+ else if (test < 0)
+ {
+ // The variance-covariance matrix of the measurement errors is not definite positive.
+ // We have to compute the eigenvalues of this matrix in order to build the penalty.
+ //a = diag(eig(H));
+ double delta = 0;
+ eigH.calculate(H); // get eigenvalues
+ //k = find(a < 0);
+ if (eigH.hasConverged())
+ {
+ const Vector &evH = eigH.getD();
+ for (i = 0; i < evH.getSize(); ++i)
+ if (evH(i) < 0)
+ delta -= evH(i);
+ }
+ logLikelihood = penalty+delta;
+ info = 44;
+ } // end if
+ break;
+
+ //if estim_params_.np > 0 // i.e. num of deep parameters >0
+ case EstimatedParameter::deepPar:
+ k = (int) estiParDesc.estParams[i].ID1;
+ deepParams(k) = estParams(i);
+ break;
+ default:
+ logLikelihood = penalty;
+ info = 1;
+ } // end switch
+ } // end found
+#ifdef DEBUG
+ mexPrintf("End of Setting of HQ params\n");
+#endif
+ } //end for
+};
+
diff --git a/mex/sources/estimation/LogLikelihoodSubSample.hh b/mex/sources/estimation/LogLikelihoodSubSample.hh
index 23b9146c2..29a28f67a 100644
--- a/mex/sources/estimation/LogLikelihoodSubSample.hh
+++ b/mex/sources/estimation/LogLikelihoodSubSample.hh
@@ -34,12 +34,11 @@ class LogLikelihoodSubSample {
public:
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 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);
+ 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 &varobs_arg, double riccati_tol_in, double lyapunov_tol, int &info);
- double compute(Vector &steadyState, const MatrixView &dataView, Vector &deepParams,
+ double compute(VectorView &steadyState, const MatrixConstView &dataView, const Vector &estParams, Vector &deepParams,
Matrix &Q, Matrix &H, VectorView &vll, int &info, size_t start, size_t period);
virtual ~LogLikelihoodSubSample();
@@ -50,6 +49,7 @@ private:
KalmanFilter kalmanFilter;
VDVEigDecomposition eigQ;
VDVEigDecomposition eigH;
+ int &info;
// methods
void updateParams(const Vector &estParams, Vector &deepParams,
Matrix &Q, Matrix &H, size_t period);
diff --git a/mex/sources/estimation/Makefile.am b/mex/sources/estimation/Makefile.am
index bcfa414d9..e820c961e 100644
--- a/mex/sources/estimation/Makefile.am
+++ b/mex/sources/estimation/Makefile.am
@@ -9,11 +9,21 @@ endif
EXTRA_DIST = \
DecisionRules.cc \
DecisionRules.hh \
- dr1.cpp \
- DsgeLikelihood.cpp \
- DsgeLikelihood.h \
- dynare_resolve.cpp \
- Estimation.h \
- ioutils.h \
- MexUtils.cpp \
- MexUtils.h
+ DetrendData.cc \
+ DetrendData.hh \
+ EstimatedParameter.hh \
+ EstimatedParametersDescription.hh \
+ EstimationSubsample.hh \
+ InitializeKalmanFilter.cc \
+ InitializeKalmanFilter.hh \
+ KalmanFilter.cc \
+ KalmanFilter.hh \
+ LogLikelihoodMain.hh \
+ LogLikelihoodSubSample.cc \
+ LogLikelihoodSubSample.hh \
+ ModelSolution.cc \
+ ModelSolution.hh \
+ Prior.hh \
+ utils/dynamic_dll.cc \
+ utils/dynamic_dll.hh \
+ utils/ts_exception.h
diff --git a/mex/sources/estimation/ModelSolution.cc b/mex/sources/estimation/ModelSolution.cc
index 67ef6a305..7bb25ea40 100644
--- a/mex/sources/estimation/ModelSolution.cc
+++ b/mex/sources/estimation/ModelSolution.cc
@@ -51,7 +51,7 @@ ModelSolution::ModelSolution(const std::string& modName, size_t n_endo_arg, siz
}
void
-ModelSolution::compute(Vector& steadyState, const Vector& deepParams, Matrix& ghx, Matrix& ghu) throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException)
+ModelSolution::compute(VectorView& steadyState, const Vector& deepParams, Matrix& ghx, Matrix& ghu) throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException)
{
// compute Steady State
ComputeSteadyState(steadyState, deepParams);
@@ -63,7 +63,7 @@ ModelSolution::compute(Vector& steadyState, const Vector& deepParams, Matrix& gh
}
void
-ModelSolution::ComputeModelSolution(Vector& steadyState, const Vector& deepParams, Matrix& ghx, Matrix& ghu) throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException)
+ModelSolution::ComputeModelSolution(VectorView &steadyState, const Vector& deepParams, Matrix& ghx, Matrix& ghu) throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException)
{
// set extended Steady State
@@ -85,7 +85,7 @@ ModelSolution::ComputeModelSolution(Vector& steadyState, const Vector& deepParam
decisionRules.compute(jacobian,ghx, ghu);
}
void
-ModelSolution::ComputeSteadyState(Vector& steadyState, const Vector& deepParams)
+ModelSolution::ComputeSteadyState(VectorView& steadyState, const Vector& deepParams)
{
// does nothig for time being.
}
diff --git a/mex/sources/estimation/ModelSolution.hh b/mex/sources/estimation/ModelSolution.hh
index 69ef4d737..2870ba9ab 100644
--- a/mex/sources/estimation/ModelSolution.hh
+++ b/mex/sources/estimation/ModelSolution.hh
@@ -42,7 +42,7 @@ public:
const std::vector& zeta_back_arg, const std::vector& zeta_mixed_arg,
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);
+ void compute( VectorView& steadyState, const Vector& deepParams, Matrix& ghx, Matrix& ghu ) throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException);
private:
const size_t n_endo;
@@ -55,8 +55,8 @@ private:
DecisionRules decisionRules;
DynamicModelDLL dynamicDLLp;
//Matrix jacobian;
- void ComputeModelSolution( Vector& steadyState, const Vector& deepParams, Matrix& ghx, Matrix& ghu ) throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException);
- void ComputeSteadyState( Vector& steadyState, const Vector& deepParams);
+ void ComputeModelSolution( VectorView& steadyState, const Vector& deepParams, Matrix& ghx, Matrix& ghu ) throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException);
+ void ComputeSteadyState( VectorView& steadyState, const Vector& deepParams);
};
diff --git a/mex/sources/estimation/Prior.cc b/mex/sources/estimation/Prior.cc
new file mode 100644
index 000000000..32e1a78c6
--- /dev/null
+++ b/mex/sources/estimation/Prior.cc
@@ -0,0 +1,38 @@
+/*
+ * 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.cpp
+// Implementation of the Class Prior
+// Created on: 02-Feb-2010 13:06:20
+///////////////////////////////////////////////////////////
+
+#include "Prior.hh"
+
+Prior::Prior(Prior::pShape shape_arg, double mean_arg, double mode_arg, double lower_bound_arg, double upper_bound_arg, double fhp_arg, double shp_arg) :
+ shape(shape_arg), mean(mean_arg), mode(mode_arg), lower_bound(lower_bound_arg), upper_bound(upper_bound_arg), fhp(fhp_arg), shp(shp_arg)
+{
+
+}
+
+Prior::~Prior()
+{
+
+}
+
diff --git a/mex/sources/estimation/Prior.hh b/mex/sources/estimation/Prior.hh
index e04ab7595..ace5a442f 100644
--- a/mex/sources/estimation/Prior.hh
+++ b/mex/sources/estimation/Prior.hh
@@ -30,9 +30,6 @@ 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
{
@@ -43,11 +40,15 @@ public:
Uniform = 5,
Inv_gamma_2 = 6 //Inverse gamma (type 2) density
};
+
+ Prior(Prior::pShape shape, double mean, double mode, double lower_bound, double upper_bound, double fhp, double shp);
+ virtual ~Prior();
+
const pShape shape;
const double mean;
const double mode;
- const double lower_boound;
- const double upper_boound;
+ const double lower_bound;
+ const double upper_bound;
/**
* first shape parameter
*/
diff --git a/mex/sources/estimation/libmat/Makefile.am b/mex/sources/estimation/libmat/Makefile.am
index b462ac588..ef4a5a450 100644
--- a/mex/sources/estimation/libmat/Makefile.am
+++ b/mex/sources/estimation/libmat/Makefile.am
@@ -5,7 +5,8 @@ endif
endif
EXTRA_DIST = \
- BlasBindings.hh
+ BlasBindings.hh \
+ DiscLyapFast.hh \
GeneralizedSchurDecomposition.cc \
GeneralizedSchurDecomposition.hh \
LUSolver.cc \
diff --git a/mex/sources/estimation/libmat/Matrix.cc b/mex/sources/estimation/libmat/Matrix.cc
index bc43da795..20eecdd86 100644
--- a/mex/sources/estimation/libmat/Matrix.cc
+++ b/mex/sources/estimation/libmat/Matrix.cc
@@ -63,16 +63,6 @@ MatrixView::MatrixView(double *data_arg, size_t rows_arg, size_t cols_arg, size_
{
}
-MatrixView::MatrixView(Matrix &arg, size_t row_offset, size_t col_offset,
- size_t rows_arg, size_t cols_arg) :
- data(arg.getData() + row_offset + col_offset*arg.getLd()), rows(rows_arg), cols(cols_arg), ld(arg.getLd())
-{
- assert(row_offset < arg.getRows()
- && row_offset + rows_arg <= arg.getRows()
- && col_offset < arg.getCols()
- && col_offset + cols_arg <= arg.getCols());
-}
-
std::ostream &
operator<<(std::ostream &out, const MatrixView &M)
{
@@ -94,16 +84,6 @@ MatrixConstView::MatrixConstView(const double *data_arg, size_t rows_arg, size_t
{
}
-MatrixConstView::MatrixConstView(const Matrix &arg, size_t row_offset, size_t col_offset,
- size_t rows_arg, size_t cols_arg) :
- data(arg.getData() + row_offset + col_offset*arg.getLd()), rows(rows_arg), cols(cols_arg), ld(arg.getLd())
-{
- assert(row_offset < arg.getRows()
- && row_offset + rows_arg <= arg.getRows()
- && col_offset < arg.getCols()
- && col_offset + cols_arg <= arg.getCols());
-}
-
std::ostream &
operator<<(std::ostream &out, const MatrixConstView &M)
{
diff --git a/mex/sources/estimation/libmat/Matrix.hh b/mex/sources/estimation/libmat/Matrix.hh
index b2ffdb20f..3ba5b94c9 100644
--- a/mex/sources/estimation/libmat/Matrix.hh
+++ b/mex/sources/estimation/libmat/Matrix.hh
@@ -102,8 +102,16 @@ class MatrixView
const size_t rows, cols, ld;
public:
MatrixView(double *data_arg, size_t rows_arg, size_t cols_arg, size_t ld_arg);
- MatrixView(Matrix &arg, size_t row_offset, size_t col_offset,
- size_t rows_arg, size_t cols_arg);
+ template
+ MatrixView(Mat &arg, size_t row_offset, size_t col_offset,
+ size_t rows_arg, size_t cols_arg) :
+ data(arg.getData() + row_offset + col_offset*arg.getLd()), rows(rows_arg), cols(cols_arg), ld(arg.getLd())
+ {
+ assert(row_offset < arg.getRows()
+ && row_offset + rows_arg <= arg.getRows()
+ && col_offset < arg.getCols()
+ && col_offset + cols_arg <= arg.getCols());
+ }
virtual ~MatrixView(){};
inline size_t getRows() const { return rows; }
inline size_t getCols() const { return cols; }
@@ -140,8 +148,16 @@ class MatrixConstView
const size_t rows, cols, ld;
public:
MatrixConstView(const double *data_arg, size_t rows_arg, size_t cols_arg, size_t ld_arg);
- MatrixConstView(const Matrix &arg, size_t row_offset, size_t col_offset,
- size_t rows_arg, size_t cols_arg);
+ template
+ MatrixConstView(const Mat &arg, size_t row_offset, size_t col_offset,
+ size_t rows_arg, size_t cols_arg) :
+ data(arg.getData() + row_offset + col_offset*arg.getLd()), rows(rows_arg), cols(cols_arg), ld(arg.getLd())
+ {
+ assert(row_offset < arg.getRows()
+ && row_offset + rows_arg <= arg.getRows()
+ && col_offset < arg.getCols()
+ && col_offset + cols_arg <= arg.getCols());
+ }
virtual ~MatrixConstView(){};
inline size_t getRows() const { return rows; }
inline size_t getCols() const { return cols; }
diff --git a/mex/sources/estimation/libmat/VDVEigDecomposition.hh b/mex/sources/estimation/libmat/VDVEigDecomposition.hh
index 54e08e63d..b22c4c6d0 100644
--- a/mex/sources/estimation/libmat/VDVEigDecomposition.hh
+++ b/mex/sources/estimation/libmat/VDVEigDecomposition.hh
@@ -33,7 +33,7 @@
class VDVEigDecomposition
{
lapack_int lda, n;
- int lwork, info;
+ lapack_int lwork, info;
double tmpwork;
double *work;
bool converged;
@@ -45,9 +45,9 @@ public:
class VDVEigException
{
public:
- const int info;
+ const lapack_int info;
std::string message;
- VDVEigException(int info_arg, std::string message_arg) :
+ VDVEigException(lapack_int info_arg, std::string message_arg) :
info(info_arg), message(message_arg) {
};
};
@@ -101,12 +101,12 @@ VDVEigDecomposition::calculate(const Mat &m) throw(VDVEigException)
if (m.getCols() != (size_t) n || m.getLd() != (size_t) lda)
throw(VDVEigException(info, "Matrix not matching VDVEigDecomposition class"));
- int tmplwork = -1;
+ lapack_int tmplwork = -1;
V = m;
dsyev("V", "U", &n, V.getData(), &lda, D.getData(), &tmpwork, &tmplwork, &info);
- if (lwork < (int) tmpwork)
+ if (lwork < tmpwork)
{
- lwork = (int) tmpwork;
+ lwork = tmpwork;
delete[] work;
work = new double[lwork];
}
diff --git a/mex/sources/estimation/libmat/Vector.cc b/mex/sources/estimation/libmat/Vector.cc
index f22993e7f..ad012c8a5 100644
--- a/mex/sources/estimation/libmat/Vector.cc
+++ b/mex/sources/estimation/libmat/Vector.cc
@@ -52,16 +52,23 @@ operator<<(std::ostream &out, const Vector &V)
return out;
}
-VectorView::VectorView(Vector &arg, size_t offset, size_t size_arg)
- : data(arg.getData() + offset), size(size_arg), stride(1)
-{
-}
-
VectorView::VectorView(double *data_arg, size_t size_arg, size_t stride_arg)
: data(data_arg), size(size_arg), stride(stride_arg)
{
}
+VectorView::VectorView(Vector &arg, size_t offset, size_t size_arg)
+ : data(arg.getData() + offset*arg.getStride()), size(size_arg), stride(arg.getStride())
+{
+ assert(offset < arg.getSize() && offset + size <= arg.getSize());
+}
+
+VectorView::VectorView(VectorView &arg, size_t offset, size_t size_arg)
+ : data(arg.getData() + offset*arg.getStride()), size(size_arg), stride(arg.getStride())
+{
+ assert(offset < arg.getSize() && offset + size <= arg.getSize());
+}
+
VectorView &
VectorView::operator= (const VectorView &arg)
{
@@ -73,6 +80,29 @@ VectorView::operator= (const VectorView &arg)
return *this;
}
+VectorConstView::VectorConstView(const double *data_arg, size_t size_arg, size_t stride_arg)
+ : data(data_arg), size(size_arg), stride(stride_arg)
+{
+}
+
+VectorConstView::VectorConstView(const Vector &arg, size_t offset, size_t size_arg)
+ : data(arg.getData() + offset*arg.getStride()), size(size_arg), stride(arg.getStride())
+{
+ assert(offset < arg.getSize() && offset + size <= arg.getSize());
+}
+
+VectorConstView::VectorConstView(const VectorView &arg, size_t offset, size_t size_arg)
+ : data(arg.getData() + offset*arg.getStride()), size(size_arg), stride(arg.getStride())
+{
+ assert(offset < arg.getSize() && offset + size <= arg.getSize());
+}
+
+VectorConstView::VectorConstView(const VectorConstView &arg, size_t offset, size_t size_arg)
+ : data(arg.getData() + offset*arg.getStride()), size(size_arg), stride(arg.getStride())
+{
+ assert(offset < arg.getSize() && offset + size <= arg.getSize());
+}
+
std::ostream &
operator<<(std::ostream &out, const VectorView &V)
{
@@ -80,16 +110,6 @@ operator<<(std::ostream &out, const VectorView &V)
return out;
}
-VectorConstView::VectorConstView(const Vector &arg, size_t offset, size_t size_arg)
- : data(arg.getData() + offset), size(size_arg), stride(1)
-{
-}
-
-VectorConstView::VectorConstView(const double *data_arg, size_t size_arg, size_t stride_arg)
- : data(data_arg), size(size_arg), stride(stride_arg)
-{
-}
-
std::ostream &
operator<<(std::ostream &out, const VectorConstView &V)
{
diff --git a/mex/sources/estimation/libmat/Vector.hh b/mex/sources/estimation/libmat/Vector.hh
index 7b77772c5..cf42438ce 100644
--- a/mex/sources/estimation/libmat/Vector.hh
+++ b/mex/sources/estimation/libmat/Vector.hh
@@ -97,8 +97,11 @@ private:
double *const data;
const size_t size, stride;
public:
- VectorView(Vector &arg, size_t offset, size_t size_arg);
VectorView(double *data_arg, size_t size_arg, size_t stride_arg);
+ /* Can't use a template for the 2 constructors below: this would override the
+ constructor which uses a pointer, because the argument list is the same */
+ VectorView(Vector &arg, size_t offset, size_t size_arg);
+ VectorView(VectorView &arg, size_t offset, size_t size_arg);
virtual ~VectorView(){};
inline size_t getSize() const { return size; }
inline size_t getStride() const { return stride; }
@@ -135,8 +138,13 @@ private:
const double *const data;
const size_t size, stride;
public:
- VectorConstView(const Vector &arg, size_t offset, size_t size_arg);
VectorConstView(const double *data_arg, size_t size_arg, size_t stride_arg);
+ /* Can't use a template for the 3 constructors below: this would override the
+ constructor which uses a pointer, because the argument list is the same */
+ VectorConstView(const Vector &arg, size_t offset, size_t size_arg);
+ VectorConstView(const VectorView &arg, size_t offset, size_t size_arg);
+ VectorConstView(const VectorConstView &arg, size_t offset, size_t size_arg);
+
virtual ~VectorConstView(){};
inline size_t getSize() const { return size; }
inline size_t getStride() const { return stride; }
diff --git a/mex/sources/estimation/tests/testInitKalman.cc b/mex/sources/estimation/tests/testInitKalman.cc
index ef7981f6d..c261f9b71 100644
--- a/mex/sources/estimation/tests/testInitKalman.cc
+++ b/mex/sources/estimation/tests/testInitKalman.cc
@@ -100,47 +100,32 @@ main(int argc, char **argv)
for (int i = 0; i < 2; ++i)
zeta_fwrd_arg.push_back(fwd[i]-1);
- size_t nriv = 6, nric = 4;
- size_t sriv[] = {7, 8, 10, 11, 12, 13};
- size_t sric[] = {3, 4, 5, 6};
+ size_t varobs[] = {12, 11};
+ std::vector varobs_arg;
+ for (size_t i = 0; i < 2; ++i)
+ varobs_arg.push_back(varobs[i]-1);
- std::vector riv;
- for (size_t i = 0; i < nriv; ++i)
- riv.push_back(sriv[i]-1);
- std::vector ric;
- for (size_t i = 0; i < nric; ++i)
- ric.push_back(sric[i]-1);
+ // Compute zeta_varobs_back_mixed
+ sort(varobs_arg.begin(), varobs_arg.end());
+ std::vector zeta_back_mixed, zeta_varobs_back_mixed;
+ set_union(zeta_back_arg.begin(), zeta_back_arg.end(),
+ zeta_mixed_arg.begin(), zeta_mixed_arg.end(),
+ back_inserter(zeta_back_mixed));
+ set_union(zeta_back_mixed.begin(), zeta_back_mixed.end(),
+ varobs_arg.begin(), varobs_arg.end(),
+ back_inserter(zeta_varobs_back_mixed));
- size_t nobs = 2;
- size_t svarobs[] = {2, 1};
- std::vector varobs;
- for (size_t i = 0; i < nobs; ++i)
- varobs.push_back(svarobs[i]-1);
+ size_t n_vbm = zeta_varobs_back_mixed.size();
Matrix
- T(riv.size(), riv.size()), R(riv.size(), n_exo),
- RQRt(riv.size(), riv.size()), Pstar(riv.size(), riv.size()),
- Pinf(riv.size(), riv.size()), Z(varobs.size(), riv.size()), Q(n_exo);
- Z.setAll(0.0);
- for (size_t i = 0; i < varobs.size(); ++i)
- Z(i, varobs[i]) = 1.0;
+ T(n_vbm, n_vbm), R(n_vbm, n_exo),
+ RQRt(n_vbm, n_vbm), Pstar(n_vbm, n_vbm),
+ Pinf(n_vbm, n_vbm), Q(n_exo);
MatrixView
vCovVW(vcov, n_exo, n_exo, n_exo);
Q = vCovVW;
- size_t sorder_var[] =
- { 4, 5, 6, 8, 9, 10, 11, 12, 14, 1, 7, 13, 2, 3, 15};
- std::vector order_var;
- for (size_t ii = 0; ii < n_endo; ++ii)
- order_var.push_back(sorder_var[ii]-1); //= (1:endo_nbr)';
-
- size_t sinv_order_var[] =
- { 10, 13, 14, 1, 2, 3, 11, 4, 5, 6, 7, 8, 12, 9, 15};
- std::vector inv_order_var;
- for (size_t ii = 0; ii < n_endo; ++ii)
- inv_order_var.push_back(sinv_order_var[ii]-1); //= (1:endo_nbr)';
-
double lyapunov_tol = 1e-16;
int info = 0;
const MatrixView
@@ -153,13 +138,13 @@ main(int argc, char **argv)
InitializeKalmanFilter
initializeKalmanFilter(modName, n_endo, n_exo,
- zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg, qz_criterium,
- order_var, inv_order_var, riv, ric, lyapunov_tol, info);
+ zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg, zeta_static_arg,
+ zeta_varobs_back_mixed, qz_criterium,
+ lyapunov_tol, info);
- std::cout << "Initilise KF with Q: " << std::endl << Q << std::endl;
- std::cout << "and Z" << std::endl << Z << std::endl;
+ std::cout << "Initialize KF with Q: " << std::endl << Q << std::endl;
- initializeKalmanFilter.initialize(steadyState, deepParams, R, Z, Q, RQRt, T, Pstar, Pinf,
+ initializeKalmanFilter.initialize(steadyState, deepParams, R, Q, RQRt, T, Pstar, Pinf,
penalty, dataView, yView, info);
std::cout << "Matrix T: " << std::endl << T << std::endl;
diff --git a/mex/sources/estimation/tests/testKalman.cc b/mex/sources/estimation/tests/testKalman.cc
index 947690fdb..b240a1bce 100644
--- a/mex/sources/estimation/tests/testKalman.cc
+++ b/mex/sources/estimation/tests/testKalman.cc
@@ -100,51 +100,20 @@ main(int argc, char **argv)
for (int i = 0; i < 2; ++i)
zeta_fwrd_arg.push_back(fwd[i]-1);
- size_t nriv = 6, nric = 4;
- size_t sriv[] = {7, 8, 10, 11, 12, 13};
- size_t sric[] = {3, 4, 5, 6};
-
- std::vector riv;
- for (size_t i = 0; i < nriv; ++i)
- riv.push_back(sriv[i]-1);
- std::vector ric;
- for (size_t i = 0; i < nric; ++i)
- ric.push_back(sric[i]-1);
-
size_t nobs = 2;
- size_t svarobs[] = {2, 1};
- std::vector varobs;
+ size_t varobs[] = {12, 11};
+ std::vector varobs_arg;
for (size_t i = 0; i < nobs; ++i)
- varobs.push_back(svarobs[i]-1);
+ varobs_arg.push_back(varobs[i]-1);
Matrix Q(n_exo), H(nobs);
-/*
- T(riv.size(), riv.size()), R(riv.size(), n_exo),
- RQRt(riv.size(), riv.size()), Pstar(riv.size(), riv.size()),
- Pinf(riv.size(), riv.size()), Z(varobs.size(), riv.size()),
- */
-
H.setAll(0.0);
-/* for (size_t i = 0; i < varobs.size(); ++i)
- Z(i, varobs[i]) = 1.0;
-*/
+
MatrixView
vCovVW(vcov, n_exo, n_exo, n_exo);
Q = vCovVW;
std::cout << "Matrix Q: " << std::endl << Q << std::endl;
- size_t sorder_var[] =
- { 4, 5, 6, 8, 9, 10, 11, 12, 14, 1, 7, 13, 2, 3, 15};
- std::vector order_var;
- for (size_t ii = 0; ii < n_endo; ++ii)
- order_var.push_back(sorder_var[ii]-1); //= (1:endo_nbr)';
-
- size_t sinv_order_var[] =
- { 10, 13, 14, 1, 2, 3, 11, 4, 5, 6, 7, 8, 12, 9, 15};
- std::vector inv_order_var;
- for (size_t ii = 0; ii < n_endo; ++ii)
- inv_order_var.push_back(sinv_order_var[ii]-1); //= (1:endo_nbr)';
-
double lyapunov_tol = 1e-16;
double riccati_tol = 1e-16;
int info = 0;
@@ -158,10 +127,8 @@ main(int argc, char **argv)
KalmanFilter kalman(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);
+ varobs_arg, riccati_tol, lyapunov_tol, info);
- std::cout << "Initilise KF with Q: " << std::endl << Q << std::endl;
- // std::cout << "and Z" << std::endl << Z << std::endl;
size_t start=0, period=0;
double ll=kalman.compute(dataView, steadyState, Q, H, deepParams,
vwll, start, period, penalty, info);
diff --git a/tests/optimal_policy/mult_elimination_test.mod b/tests/optimal_policy/mult_elimination_test.mod
new file mode 100644
index 000000000..0e87973b6
--- /dev/null
+++ b/tests/optimal_policy/mult_elimination_test.mod
@@ -0,0 +1,67 @@
+// Test of mult_elimination function
+
+// parameters value is set to posterior mean as computed by ls1.mod
+// dR in objective function
+
+var A de dq dR pie pie_obs pie_s R R_obs y y_obs y_s ;
+varexo e_A e_pies e_q e_ys ;
+
+parameters psi1 psi2 psi3 rho_R tau alpha rr k rho_q rho_A rho_ys rho_pies ls1_r s1 s2 s3 s4 s5;
+
+psi1 = 1.54;
+psi2 = 0.25;
+psi3 = 0.25;
+rho_R = 0.5;
+alpha = 0.3;
+rr = 2.51;
+k = 0.5;
+tau = 0.5;
+rho_q = 0.4;
+rho_A = 0.2;
+rho_ys = 0.9;
+rho_pies = 0.7;
+
+model(linear);
+y = y(+1) - (tau +alpha*(2-alpha)*(1-tau))*(R-pie(+1))-alpha*(tau +alpha*(2-alpha)*(1-tau))*dq(+1) + alpha*(2-alpha)*((1-tau)/tau)*(y_s-y_s(+1))-A(+1);
+pie = exp(-rr/400)*pie(+1)+alpha*exp(-rr/400)*dq(+1)-alpha*dq+(k/(tau+alpha*(2-alpha)*(1-tau)))*y+alpha*(2-alpha)*(1-tau)/(tau*(tau+alpha*(2-alpha)*(1-tau)))*y_s;
+pie = de+(1-alpha)*dq+pie_s;
+//R = rho_R*R(-1)+(1-rho_R)*(psi1*pie+psi2*(y+alpha*(2-alpha)*((1-tau)/tau)*y_s)+psi3*de)+e_R;
+dq = rho_q*dq(-1)+e_q;
+y_s = rho_ys*y_s(-1)+e_ys;
+pie_s = rho_pies*pie_s(-1)+e_pies;
+A = rho_A*A(-1)+e_A;
+y_obs = y-y(-1)+A;
+pie_obs = 4*pie;
+R_obs = 4*R;
+dR = R-R(-1);
+end;
+
+shocks;
+var e_q = 2.5^2;
+var e_A = 1.89;
+var e_ys = 1.89;
+var e_pies = 1.89;
+end;
+
+planner_objective 0.25*pie_obs^2+y^2+0.1*dR^2;
+
+ramsey_policy(order=1,irf=0,planner_discount=0.95);
+
+dr2 = mult_elimination({'R'},M_,options_,oo_);
+
+k1 = oo_.dr.nstatic+(1:oo_.dr.npred);
+k2 = strmatch('mult_',M_.endo_names(oo_.dr.order_var(k1),:));
+k3 = k1(setdiff(1:oo_.dr.npred,k2));
+k4 = oo_.dr.order_var(k3);
+
+V0 = oo_.var(k4,k4);
+
+
+Atest = [dr2.M1(k3,:) dr2.M2(k3,:) dr2.M4(k3,:); eye(6) zeros(6,10);zeros(4,16)];
+Btest = [dr2.M3(k3,:); zeros(6,4); eye(4)];
+
+V1=lyapunov_symm(Atest,Btest*M_.Sigma_e*Btest',options_.qz_criterium,options_.lyapunov_complex_threshold);
+
+if max(max(abs(V1(1:6,1:6)-V0)))
+ disp('Test OK')
+end
\ No newline at end of file