Adding InitializeKalmanFilter.cc/hh, its supporting files: DetrendData.cc/hh (a basic one), libmat/DiscLyapFast.cc, a test driver routine and a minor amendment to ModelSolution.hh.
parent
882ea0ed75
commit
881b1168fe
|
@ -0,0 +1,37 @@
|
|||
/*
|
||||
* Copyright (C) 2009-2010 Dynare Team
|
||||
*
|
||||
* This file is part of Dynare.
|
||||
*
|
||||
* Dynare is free software: you can redistribute it and/or modify
|
||||
* it under the terms of the GNU General Public License as published by
|
||||
* the Free Software Foundation, either version 3 of the License, or
|
||||
* (at your option) any later version.
|
||||
*
|
||||
* Dynare is distributed in the hope that it will be useful,
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
* GNU General Public License for more details.
|
||||
*
|
||||
* You should have received a copy of the GNU General Public License
|
||||
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
///////////////////////////////////////////////////////////
|
||||
// DetrendData.cpp
|
||||
// Implementation of the Class DetrendData
|
||||
// Created on: 02-Feb-2010 13:01:15
|
||||
///////////////////////////////////////////////////////////
|
||||
|
||||
#include "DetrendData.hh"
|
||||
|
||||
DetrendData::DetrendData(const bool INlogLinear) //, Vector& INtrendCoeff)
|
||||
: logLinear(INlogLinear) //trendCoeff(INtrendCoeff)
|
||||
{
|
||||
};
|
||||
|
||||
void
|
||||
DetrendData::detrend(const Vector &SteadyState, const MatrixView &dataView, Matrix &Yiew)
|
||||
{
|
||||
};
|
||||
|
|
@ -0,0 +1,44 @@
|
|||
/*
|
||||
* Copyright (C) 2009-2010 Dynare Team
|
||||
*
|
||||
* This file is part of Dynare.
|
||||
*
|
||||
* Dynare is free software: you can redistribute it and/or modify
|
||||
* it under the terms of the GNU General Public License as published by
|
||||
* the Free Software Foundation, either version 3 of the License, or
|
||||
* (at your option) any later version.
|
||||
*
|
||||
* Dynare is distributed in the hope that it will be useful,
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
* GNU General Public License for more details.
|
||||
*
|
||||
* You should have received a copy of the GNU General Public License
|
||||
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
///////////////////////////////////////////////////////////
|
||||
// DetrendData.h
|
||||
// Implementation of the Class DetrendData
|
||||
// Created on: 02-Feb-2010 13:01:15
|
||||
///////////////////////////////////////////////////////////
|
||||
|
||||
#if !defined(DetrendData_312823A1_6248_4af0_B204_DB22F1237E9B__INCLUDED_)
|
||||
#define DetrendData_312823A1_6248_4af0_B204_DB22F1237E9B__INCLUDED_
|
||||
|
||||
#include "Matrix.hh"
|
||||
|
||||
class DetrendData
|
||||
{
|
||||
|
||||
public:
|
||||
virtual ~DetrendData(){};
|
||||
DetrendData(const bool logLinear); // add later Vector& trendCoeff);
|
||||
void detrend(const Vector &SteadyState, const MatrixView &dataView, Matrix &Yiew);
|
||||
|
||||
private:
|
||||
const bool logLinear;
|
||||
//Vector trendCoeff;
|
||||
|
||||
};
|
||||
|
||||
#endif // !defined(312823A1_6248_4af0_B204_DB22F1237E9B__INCLUDED_)
|
|
@ -0,0 +1,125 @@
|
|||
/*
|
||||
* Copyright (C) 2010 Dynare Team
|
||||
*
|
||||
* This file is part of Dynare.
|
||||
*
|
||||
* Dynare is free software: you can redistribute it and/or modify
|
||||
* it under the terms of the GNU General Public License as published by
|
||||
* the Free Software Foundation, either version 3 of the License, or
|
||||
* (at your option) any later version.
|
||||
*
|
||||
* Dynare is distributed in the hope that it will be useful,
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
* GNU General Public License for more details.
|
||||
*
|
||||
* You should have received a copy of the GNU General Public License
|
||||
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
///////////////////////////////////////////////////////////
|
||||
// InitializeKalmanFilter.cpp
|
||||
// Implementation of the Class InitializeKalmanFilter
|
||||
// Created on: 02-Feb-2010 12:25:28
|
||||
///////////////////////////////////////////////////////////
|
||||
|
||||
#include "InitializeKalmanFilter.hh"
|
||||
#include "BlasBindings.hh"
|
||||
#include "disclyap_fast.hh"
|
||||
|
||||
InitializeKalmanFilter::~InitializeKalmanFilter()
|
||||
{
|
||||
}
|
||||
|
||||
InitializeKalmanFilter::InitializeKalmanFilter(const std::string &modName, size_t n_endo_arg, size_t n_exo_arg,
|
||||
const std::vector<size_t> &zeta_fwrd_arg, const std::vector<size_t> &zeta_back_arg,
|
||||
const std::vector<size_t> &zeta_mixed_arg, const std::vector<size_t> &zeta_static_arg,
|
||||
const Matrix &llincidence, double qz_criterium, const std::vector<size_t> &order_var_arg, const std::vector<size_t> &inv_order_var_arg,
|
||||
const std::vector<size_t> &rivIN, const std::vector<size_t> &ricIN, double lyapunov_tolIN, int &info) :
|
||||
riv(rivIN), inv_ric(ricIN.size()), order_var(order_var_arg), lyapunov_tol(lyapunov_tolIN),
|
||||
detrendData(false), modelSolution(modName, n_endo_arg, n_exo_arg, zeta_fwrd_arg, zeta_back_arg,
|
||||
zeta_mixed_arg, zeta_static_arg, llincidence, qz_criterium), discLyapFast(riv.size()),
|
||||
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), RQR(riv.size(), riv.size())
|
||||
{
|
||||
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];
|
||||
}
|
||||
|
||||
void
|
||||
InitializeKalmanFilter::initialize(Vector &steadyState, Vector &deepParams, const Vector &xparams1,
|
||||
Matrix &R, Matrix &Z, Matrix &Q, Matrix &T, Matrix &Pstar, Matrix &Pinf,
|
||||
double &penalty, const MatrixView &dataView, MatrixView &yView, int &info)
|
||||
{
|
||||
modelSolution.compute(steadyState, deepParams, ghx_raw, ghu_raw);
|
||||
|
||||
mat::reorderRowsByVectors(ghx, mat::nullVec, ghx_raw, order_var);
|
||||
mat::reorderRowsByVectors(ghu, mat::nullVec, ghu_raw, order_var);
|
||||
|
||||
setT(T, info);
|
||||
setR(R, info);
|
||||
setPstar(Pstar, Pinf, T, R, Q, info);
|
||||
}
|
||||
|
||||
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);
|
||||
|
||||
T.setAll(0.0);
|
||||
|
||||
//T(i_n_iv,ic) = dr.ghx(iv,:);
|
||||
mat::assignByVectors(T, mat::nullVec, inv_ric, ghx, riv, mat::nullVec);
|
||||
}
|
||||
|
||||
void
|
||||
InitializeKalmanFilter::setR(Matrix &R, 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);
|
||||
}
|
||||
|
||||
void
|
||||
InitializeKalmanFilter::setPstar(Matrix &Pstar, Matrix &Pinf, Matrix &T, Matrix &R, Matrix &Q, int &info)
|
||||
{
|
||||
|
||||
// Matrix RQR=R*Q*R'
|
||||
RQ.setAll(0.0);
|
||||
blas::gemm("N", "N", 1.0, R, Q, 0.0, RQ); // R*Q
|
||||
RQR.setAll(0.0);
|
||||
//mat::transpose(Rt, R);
|
||||
blas::gemm("N", "T", 1.0, RQ, R, 0.0, RQR); // R*Q*R'
|
||||
//mat::transpose(RQR);
|
||||
|
||||
try
|
||||
{
|
||||
// disclyap_fast(T, RQR, Pstar, lyapunov_tol, 0); // 1 to check chol
|
||||
discLyapFast.solve_lyap(T, RQR, Pstar, lyapunov_tol, 0); // 1 to check chol
|
||||
|
||||
Pinf.setAll(0.0);
|
||||
}
|
||||
catch(const DiscLyapFast::DLPException &e)
|
||||
{
|
||||
if (e.info > 0) // The matrix is not positive definite in NormCholesky calculator
|
||||
{
|
||||
printf(e.message.c_str());
|
||||
info = -1; //likelihood = penalty;
|
||||
return;
|
||||
}
|
||||
else if (e.info < 0)
|
||||
{
|
||||
printf("Caugth unhandled TS exception with Pstar matrix: ");
|
||||
printf(e.message.c_str());
|
||||
info = -1; //likelihood = penalty;
|
||||
throw;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,69 @@
|
|||
/*
|
||||
* Copyright (C) 2010 Dynare Team
|
||||
*
|
||||
* This file is part of Dynare.
|
||||
*
|
||||
* Dynare is free software: you can redistribute it and/or modify
|
||||
* it under the terms of the GNU General Public License as published by
|
||||
* the Free Software Foundation, either version 3 of the License, or
|
||||
* (at your option) any later version.
|
||||
*
|
||||
* Dynare is distributed in the hope that it will be useful,
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
* GNU General Public License for more details.
|
||||
*
|
||||
* You should have received a copy of the GNU General Public License
|
||||
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
///////////////////////////////////////////////////////////
|
||||
// InitializeKalmanFilter.h
|
||||
// Implementation of the Class InitializeKalmanFilter
|
||||
// Created on: 02-Feb-2010 12:25:28
|
||||
///////////////////////////////////////////////////////////
|
||||
|
||||
#if !defined(C3D996B8_22AB_4b77_B693_BA4777AFB091__INCLUDED_)
|
||||
#define C3D996B8_22AB_4b77_B693_BA4777AFB091__INCLUDED_
|
||||
|
||||
#include "DetrendData.hh"
|
||||
#include "ModelSolution.hh"
|
||||
#include "DiscLyapFast.hh"
|
||||
#include <string>
|
||||
|
||||
/**
|
||||
* if model is declared stationary ?compute covariance matrix of endogenous
|
||||
* variables () by doubling algorithm
|
||||
*
|
||||
*/
|
||||
class InitializeKalmanFilter
|
||||
{
|
||||
|
||||
public:
|
||||
InitializeKalmanFilter(const std::string& modName, size_t n_endo, size_t n_exo, const std::vector<size_t> &zeta_fwrd_arg,
|
||||
const std::vector<size_t> &zeta_back_arg, const std::vector<size_t> &zeta_mixed_arg, const std::vector<size_t> &zeta_static_arg,
|
||||
const Matrix &llincidence, double qz_criterium, const std::vector<size_t> &order_var_arg, const std::vector<size_t> &inv_order_var_arg,
|
||||
const std::vector<size_t> &riv, const std::vector<size_t> &ric, double lyapunov_tol, int &info);
|
||||
virtual ~InitializeKalmanFilter();
|
||||
void initialize(Vector &steadyState, Vector &deepParams, const Vector &xparams1, Matrix &R, Matrix &Z, Matrix &Q,
|
||||
Matrix &T, Matrix &Pstar, Matrix &Pinf, double &penalty, const MatrixView &dataView, MatrixView &yView, int &info);
|
||||
|
||||
private:
|
||||
const std::vector<size_t> riv; // restrict_var_list
|
||||
std::vector<size_t> inv_ric; // inverse restrict_columns
|
||||
const std::vector<size_t> order_var;
|
||||
const double lyapunov_tol;
|
||||
|
||||
DetrendData detrendData;
|
||||
ModelSolution modelSolution;
|
||||
DiscLyapFast discLyapFast; //Lyapunov solver
|
||||
Matrix ghx, ghx_raw;
|
||||
Matrix ghu, ghu_raw;
|
||||
Matrix Rt, RQ, RQR;
|
||||
void setT(Matrix &T, int &info);
|
||||
void setR(Matrix &R, int &info);
|
||||
void setPstar(Matrix &Pstar, Matrix &Pinf, Matrix &T, Matrix &R, Matrix &Q, int &info);
|
||||
|
||||
};
|
||||
|
||||
#endif // !defined(C3D996B8_22AB_4b77_B693_BA4777AFB091__INCLUDED_)
|
|
@ -41,12 +41,13 @@ public:
|
|||
ModelSolution(const std::string& modName, size_t n_endo, size_t n_exo, const std::vector<size_t>& zeta_fwrd_arg,
|
||||
const std::vector<size_t>& zeta_back_arg, const std::vector<size_t>& zeta_mixed_arg,
|
||||
const std::vector<size_t>& zeta_static_arg, const Matrix& llincidence, double qz_criterium);
|
||||
virtual ~ModelSolution(){};
|
||||
void compute( Vector& steadyState, const Vector& deepParams, Matrix& ghx, Matrix& ghu ) throw (DecisionRules::BlanchardKahnException, GeneralizedSchurDecomposition::GSDException);
|
||||
|
||||
private:
|
||||
const int n_endo;
|
||||
const int n_exo;
|
||||
const int n_jcols; // Num of Jacobian columns
|
||||
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
|
||||
Matrix jacobian;
|
||||
Vector residual;
|
||||
|
|
|
@ -0,0 +1,121 @@
|
|||
/*
|
||||
* Copyright (C) 2010 Dynare Team
|
||||
*
|
||||
* This file is part of Dynare.
|
||||
*
|
||||
* Dynare is free software: you can redistribute it and/or modify
|
||||
* it under the terms of the GNU General Public License as published by
|
||||
* the Free Software Foundation, either version 3 of the License, or
|
||||
* (at your option) any later version.
|
||||
*
|
||||
* Dynare is distributed in the hope that it will be useful,
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
* GNU General Public License for more details.
|
||||
*
|
||||
* You should have received a copy of the GNU General Public License
|
||||
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
/****************************************************************
|
||||
% a wrapper class for function X=disclyap_fast(G,V,ch)
|
||||
%
|
||||
% Solve the discrete Lyapunov Equation
|
||||
% X=G*X*G'+V
|
||||
% Using the Doubling Algorithm
|
||||
%
|
||||
% If ch is defined then the code will check if the resulting X
|
||||
% is positive definite and generate an error message if it is not
|
||||
%
|
||||
% based on work of Joe Pearlman and Alejandro Justiniano
|
||||
% 3/5/2005
|
||||
% C++ version 28/07/09 by Dynare team
|
||||
****************************************************************/
|
||||
|
||||
#if !defined(DiscLyapFast_INCLUDE)
|
||||
#define DiscLyapFast_INCLUDE
|
||||
|
||||
#include "dynlapack.h"
|
||||
#include "Matrix.hh"
|
||||
#include "BlasBindings.hh"
|
||||
|
||||
class DiscLyapFast
|
||||
{
|
||||
Matrix A0, A1, Ptmp, P0, P1, I;
|
||||
|
||||
public:
|
||||
class DLPException
|
||||
{
|
||||
public:
|
||||
const int info;
|
||||
std::string message;
|
||||
DLPException(int info_arg, std::string message_arg) :
|
||||
info(info_arg), message(message_arg) {
|
||||
};
|
||||
};
|
||||
|
||||
DiscLyapFast(size_t n) :
|
||||
A0(n), A1(n), Ptmp(n), P0(n), P1(n), I(n)
|
||||
{
|
||||
mat::set_identity(I);
|
||||
};
|
||||
virtual ~DiscLyapFast(){};
|
||||
template <class MatG, class MatV, class MatX >
|
||||
void solve_lyap(const MatG &G, const MatV &V, MatX &X, double tol, size_t flag_ch) throw(DLPException);
|
||||
|
||||
};
|
||||
|
||||
template <class MatG, class MatV, class MatX >
|
||||
void
|
||||
DiscLyapFast::solve_lyap(const MatG &G, const MatV &V, MatX &X, double tol = 1e-16, size_t flag_ch = 0) throw(DLPException)
|
||||
{
|
||||
P0 = V;
|
||||
P1 = V;
|
||||
A0 = G;
|
||||
|
||||
const double alpha = 1.0;
|
||||
const double half = 0.5;
|
||||
const double omega = 0.0;
|
||||
|
||||
bool matd = true;
|
||||
while (matd) // matrix diff > tol
|
||||
{
|
||||
//P1=P0+A0*P0*A0';
|
||||
// first step Ptmp=P0*A0';
|
||||
// DGEMM: C := alpha*op( A )*op( B ) + beta*C,
|
||||
blas::gemm("N", "T", alpha, P0, A0, omega, Ptmp);
|
||||
// P1=P0+A0*Ptmp;
|
||||
blas::gemm("N", "N", alpha, A0, Ptmp, alpha, P1);
|
||||
// A1=A0*A0;
|
||||
blas::gemm("N", "N", alpha, A0, A0, omega, A1);
|
||||
|
||||
// ensure symmetry of P1=(P1+P1')/2;
|
||||
Ptmp = P1;
|
||||
blas::gemm("T", "N", half, Ptmp, I, half, P1);
|
||||
|
||||
// check if max( max( abs( P1 - P0 ) ) )>tol
|
||||
matd = mat::isDiffSym(P1, P0, tol);
|
||||
P0 = P1;
|
||||
A0 = A1;
|
||||
} //end while
|
||||
|
||||
// ensure symmetry of X=P0=(P0+P0')/2;
|
||||
blas::gemm("T", "N", half, P1, I, half, P0);
|
||||
X = P0;
|
||||
|
||||
// Check that X is positive definite
|
||||
if (flag_ch == 1) // calc NormCholesky (P0)
|
||||
{
|
||||
lapack_int lpinfo = 0;
|
||||
lapack_int lrows = P0.getRows();
|
||||
lapack_int ldl = P0.getLd();
|
||||
dpotrf("L", &lrows, P0.getData(), &ldl, &lpinfo);
|
||||
if (lpinfo < 0)
|
||||
throw DLPException((int) lpinfo, std::string("DiscLyapFast:Internal error in NormCholesky calculator"));
|
||||
else if (lpinfo > 0)
|
||||
throw DLPException((int) lpinfo, std::string("DiscLyapFast:The matrix is not positive definite in NormCholesky calculator"));
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
#endif //if !defined(DiscLyapFast_INCLUDE)
|
|
@ -0,0 +1,184 @@
|
|||
/*
|
||||
* Copyright (C) 2010 Dynare Team
|
||||
*
|
||||
* This file is part of Dynare.
|
||||
*
|
||||
* Dynare is free software: you can redistribute it and/or modify
|
||||
* it under the terms of the GNU General Public License as published by
|
||||
* the Free Software Foundation, either version 3 of the License, or
|
||||
* (at your option) any later version.
|
||||
*
|
||||
* Dynare is distributed in the hope that it will be useful,
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
* GNU General Public License for more details.
|
||||
*
|
||||
* You should have received a copy of the GNU General Public License
|
||||
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
// Test for InitializeKalmanFilter
|
||||
// Uses fs2000k2e.mod and its ..._dynamic.mexw32
|
||||
|
||||
#include "InitializeKalmanFilter.hh"
|
||||
|
||||
int
|
||||
main(int argc, char **argv)
|
||||
{
|
||||
|
||||
std::string
|
||||
modName("fs2000k2_dynamic.mexw32");
|
||||
const int npar = 7; //(int)mxGetM(mxFldp);
|
||||
const size_t n_endo = 15, n_exo = 2;
|
||||
std::vector<size_t> zeta_fwrd_arg;
|
||||
std::vector<size_t> zeta_back_arg;
|
||||
std::vector<size_t> zeta_mixed_arg;
|
||||
std::vector<size_t> zeta_static_arg;
|
||||
//std::vector<size_t>
|
||||
double qz_criterium = 1.000001; //1.0+1.0e-9;
|
||||
Vector
|
||||
steadyState(n_endo), deepParams(npar);
|
||||
|
||||
double dYSparams [] = {
|
||||
1.0002, 0.9933, 1.0070, 1.0000,
|
||||
2.7186, 1.0073, 18.9822, 0.8608,
|
||||
0.3167, 0.8610, 1.0085, 0.9917,
|
||||
1.3559, 1.0085, 0.9929
|
||||
};
|
||||
double vcov1[] = {
|
||||
0.1960e-3, 0.0,
|
||||
0.0, 0.0250e-3
|
||||
};
|
||||
|
||||
double vcov[] = {
|
||||
0.0013, 0,
|
||||
0, 0.0001
|
||||
};
|
||||
|
||||
Matrix
|
||||
ll_incidence(3, n_endo); // leads and lags indices
|
||||
double inllincidence[] = {
|
||||
1, 5, 0,
|
||||
2, 6, 20,
|
||||
0, 7, 21,
|
||||
0, 8, 0,
|
||||
0, 9, 0,
|
||||
0, 10, 0,
|
||||
3, 11, 0,
|
||||
0, 12, 0,
|
||||
0, 13, 0,
|
||||
0, 14, 0,
|
||||
0, 15, 0,
|
||||
0, 16, 0,
|
||||
4, 17, 0,
|
||||
0, 18, 0,
|
||||
0, 19, 22,
|
||||
};
|
||||
MatrixView
|
||||
llincidence(inllincidence, 3, n_endo, 3); // leads and lags indices
|
||||
ll_incidence = llincidence;
|
||||
|
||||
double dparams[] = {
|
||||
0.3560,
|
||||
0.9930,
|
||||
0.0085,
|
||||
1.0002,
|
||||
0.1290,
|
||||
0.6500,
|
||||
0.0100
|
||||
};
|
||||
|
||||
VectorView
|
||||
modParamsVW(dparams, npar, 1);
|
||||
deepParams = modParamsVW;
|
||||
VectorView
|
||||
steadyStateVW(dYSparams, n_endo, 1);
|
||||
steadyState = steadyStateVW;
|
||||
std::cout << "Vector deepParams: " << std::endl << deepParams << std::endl;
|
||||
std::cout << "MatrixVw llincidence: " << std::endl << llincidence << std::endl;
|
||||
std::cout << "Matrix ll_incidence: " << std::endl << ll_incidence << std::endl;
|
||||
std::cout << "Vector steadyState: " << std::endl << steadyState << std::endl;
|
||||
|
||||
// Set zeta vectors [0:(n-1)] from Matlab indices [1:n] so that:
|
||||
// order_var = [ stat_var(:); pred_var(:); both_var(:); fwrd_var(:)];
|
||||
size_t statc[] = { 4, 5, 6, 8, 9, 10, 11, 12, 14};
|
||||
size_t back[] = {1, 7, 13};
|
||||
size_t both[] = {2};
|
||||
size_t fwd[] = { 3, 15};
|
||||
for (int i = 0; i < 9; ++i)
|
||||
zeta_static_arg.push_back(statc[i]-1);
|
||||
for (int i = 0; i < 3; ++i)
|
||||
zeta_back_arg.push_back(back[i]-1);
|
||||
for (int i = 0; i < 1; ++i)
|
||||
zeta_mixed_arg.push_back(both[i]-1);
|
||||
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<size_t> riv;
|
||||
for (size_t i = 0; i < nriv; ++i)
|
||||
riv.push_back(sriv[i]-1);
|
||||
std::vector<size_t> ric;
|
||||
for (size_t i = 0; i < nric; ++i)
|
||||
ric.push_back(sric[i]-1);
|
||||
|
||||
size_t nobs = 2;
|
||||
size_t svarobs[] = {8, 7};
|
||||
std::vector<size_t> varobs;
|
||||
for (size_t i = 0; i < nobs; ++i)
|
||||
varobs.push_back(svarobs[i]-1);
|
||||
|
||||
Matrix
|
||||
T(riv.size(), riv.size()), R(riv.size(), n_exo), 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;
|
||||
|
||||
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<size_t> 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<size_t> 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
|
||||
dataView(&lyapunov_tol, 1, 1, 1); // dummy
|
||||
MatrixView
|
||||
yView(&lyapunov_tol, 1, 1, 1); // dummy
|
||||
const Vector
|
||||
xparams1(0); // dummy
|
||||
double penalty = 1e8;
|
||||
|
||||
InitializeKalmanFilter
|
||||
initializeKalmanFilter(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);
|
||||
|
||||
std::cout << "Initilise KF with Q: " << std::endl << Q << std::endl;
|
||||
std::cout << "and Z" << std::endl << Z << std::endl;
|
||||
|
||||
initializeKalmanFilter.initialize(steadyState, deepParams, xparams1, R, Z, Q, T, Pstar, Pinf,
|
||||
penalty, dataView, yView, info);
|
||||
|
||||
std::cout << "Matrix T: " << std::endl << T << std::endl;
|
||||
std::cout << "Matrix R: " << std::endl << R << std::endl;
|
||||
std::cout << "Matrix Pstar: " << std::endl << Pstar << std::endl;
|
||||
std::cout << "Matrix Pinf: " << std::endl << Pinf << std::endl;
|
||||
|
||||
}
|
||||
|
Loading…
Reference in New Issue