Merge remote branch 'george/master'
Sébastien Villemot
commit
c38f4893ef
|
|
@ -104,6 +104,10 @@ extern "C" {
|
|||
void dpotrf(LACHAR uplo, CONST_LAINT n, LADOU a, CONST_LAINT lda,
|
||||
LAINT info);
|
||||
|
||||
#define dppsv FORTRAN_WRAPPER(dppsv)
|
||||
void dppsv(LACHAR uplo, CONST_LAINT n, CONST_LAINT m, LADOU a, LADOU b, CONST_LAINT ldb,
|
||||
LAINT info);
|
||||
|
||||
#define dpotri FORTRAN_WRAPPER(dpotri)
|
||||
void dpotri(LACHAR uplo, CONST_LAINT n, LADOU a, CONST_LAINT lda,
|
||||
LAINT info);
|
||||
|
|
|
|||
|
|
@ -24,6 +24,7 @@
|
|||
///////////////////////////////////////////////////////////
|
||||
|
||||
#include "KalmanFilter.hh"
|
||||
#include "LapackBindings.hh"
|
||||
|
||||
KalmanFilter::~KalmanFilter()
|
||||
{
|
||||
|
|
@ -37,13 +38,14 @@ KalmanFilter::KalmanFilter(const std::string &dynamicDllFile, size_t n_endo, siz
|
|||
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()),
|
||||
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),
|
||||
oldKFinv(zeta_varobs_back_mixed.size(), 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(dynamicDllFile, n_endo, n_exo, zeta_fwrd_arg, zeta_back_arg, zeta_mixed_arg,
|
||||
zeta_static_arg, zeta_varobs_back_mixed, qz_criterium_arg, lyapunov_tol_arg, info)
|
||||
zeta_static_arg, zeta_varobs_back_mixed, qz_criterium_arg, lyapunov_tol_arg, info),
|
||||
FUTP(varobs_arg.size()*(varobs_arg.size()+1)/2)
|
||||
{
|
||||
a_init.setAll(0.0);
|
||||
Z.setAll(0.0);
|
||||
|
|
@ -65,7 +67,7 @@ KalmanFilter::compute_zeta_varobs_back_mixed(const std::vector<size_t> &zeta_bac
|
|||
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(),
|
||||
varobs_sorted.begin(), varobs_sorted.end(),
|
||||
back_inserter(zeta_varobs_back_mixed));
|
||||
return zeta_varobs_back_mixed;
|
||||
}
|
||||
|
|
@ -83,7 +85,12 @@ KalmanFilter::compute(const MatrixConstView &dataView, VectorView &steadyState,
|
|||
initKalmanFilter.initialize(steadyState, deepParams, R, Q, RQRt, T,
|
||||
penalty, dataView, detrendedDataView, info);
|
||||
|
||||
return filter(detrendedDataView, H, vll, start, info);
|
||||
double lik= filter(detrendedDataView, H, vll, start, info);
|
||||
|
||||
if (info != 0)
|
||||
return penalty;
|
||||
else
|
||||
return lik;
|
||||
|
||||
};
|
||||
|
||||
|
|
@ -95,7 +102,6 @@ KalmanFilter::filter(const MatrixView &detrendedDataView, const Matrix &H, Vect
|
|||
{
|
||||
double loglik=0.0, ll, logFdet, Fdet;
|
||||
size_t p = Finv.getRows();
|
||||
|
||||
bool nonstationary = true;
|
||||
for (size_t t = 0; t < detrendedDataView.getCols(); ++t)
|
||||
{
|
||||
|
|
@ -111,13 +117,57 @@ KalmanFilter::filter(const MatrixView &detrendedDataView, const Matrix &H, Vect
|
|||
|
||||
// Finv=inv(F)
|
||||
mat::set_identity(Finv);
|
||||
Finverter.invMult("N", F, Finv); // F now contains its LU decomposition!
|
||||
// Pack F upper trinagle as vector
|
||||
for (size_t i=1;i<=p;++i)
|
||||
for(size_t j=i;j<=p;++j)
|
||||
FUTP(i + (j-1)*j/2 -1) = F(i-1,j-1);
|
||||
|
||||
info=lapack::choleskySolver(FUTP, Finv, "U"); // F now contains its Chol decomposition!
|
||||
if (info<0)
|
||||
{
|
||||
std::cout << "Info:" << info << std::endl;
|
||||
std::cout << "t:" << t << std::endl;
|
||||
return 0;
|
||||
}
|
||||
if (info>0)
|
||||
{
|
||||
//enforce Pstar symmetry with P=(P+P')/2=0.5P+0.5P'
|
||||
mat::transpose(Ptmp, Pstar);
|
||||
mat::add(Pstar,Ptmp);
|
||||
for (size_t i=0;i<Pstar.getCols();++i)
|
||||
for(size_t j=0;j<Pstar.getCols();++j)
|
||||
Pstar(i,j)*=0.5;
|
||||
|
||||
// K=PZ'
|
||||
blas::gemm("N", "T", 1.0, Pstar, Z, 0.0, K);
|
||||
|
||||
//F=ZPZ' +H = ZK+H
|
||||
F = H;
|
||||
blas::gemm("N", "N", 1.0, Z, K, 1.0, F);
|
||||
|
||||
// Finv=inv(F)
|
||||
mat::set_identity(Finv);
|
||||
// Pack F upper trinagle as vector
|
||||
for (size_t i=1;i<=p;++i)
|
||||
for(size_t j=i;j<=p;++j)
|
||||
FUTP(i + (j-1)*j/2 -1) = F(i-1,j-1);
|
||||
|
||||
info=lapack::choleskySolver(FUTP, Finv, "U"); // F now contains its Chol decomposition!
|
||||
if (info!=0)
|
||||
{
|
||||
return 0;
|
||||
}
|
||||
|
||||
}
|
||||
// KFinv gain matrix
|
||||
blas::gemm("N", "N", 1.0, K, Finv, 0.0, KFinv);
|
||||
// deteminant of F:
|
||||
Fdet = 1;
|
||||
for (size_t d = 0; d < p; ++d)
|
||||
Fdet *= (-F(d, d));
|
||||
for (size_t d = 1; d <= p; ++d)
|
||||
Fdet *= FUTP(d + (d-1)*d/2 -1);
|
||||
Fdet *=Fdet;//*pow(-1.0,p);
|
||||
// for (size_t d = 0; d < p; ++d)
|
||||
// Fdet *= (-F(d, d));
|
||||
|
||||
logFdet=log(fabs(Fdet));
|
||||
|
||||
|
|
@ -131,6 +181,11 @@ KalmanFilter::filter(const MatrixView &detrendedDataView, const Matrix &H, Vect
|
|||
// 3) Pt+1= Ptmp*T' +RQR'
|
||||
Pstar = RQRt;
|
||||
blas::gemm("N", "T", 1.0, Ptmp, T, 1.0, Pstar);
|
||||
//enforce Pstar symmetry with P=(P+P')/2=0.5P+0.5P'
|
||||
//blas::gemm("N", "T", 0.5, Ptmp, T, 0.5, Pstar);
|
||||
//mat::transpose(Ptmp, Pstar);
|
||||
//mat::add(Pstar,Ptmp);
|
||||
|
||||
if (t>0)
|
||||
nonstationary = mat::isDiff(KFinv, oldKFinv, riccati_tol);
|
||||
oldKFinv=KFinv;
|
||||
|
|
@ -139,6 +194,7 @@ KalmanFilter::filter(const MatrixView &detrendedDataView, const Matrix &H, Vect
|
|||
// err= Yt - Za
|
||||
MatrixConstView yt(detrendedDataView, 0, t, detrendedDataView.getRows(), 1); // current observation vector
|
||||
vt = yt;
|
||||
|
||||
blas::gemm("N", "N", -1.0, Z, a_init, 1.0, vt);
|
||||
// at+1= T(at+ KFinv *err)
|
||||
blas::gemm("N", "N", 1.0, KFinv, vt, 1.0, a_init);
|
||||
|
|
|
|||
|
|
@ -70,12 +70,12 @@ private:
|
|||
Matrix RQRt, Ptmp; //mm*mm variance-covariance matrix of variable disturbances
|
||||
Matrix F, Finv; // nob*nob F=ZPZt +H an inv(F)
|
||||
Matrix K, KFinv, oldKFinv; // mm*nobs K=PZt and K*Finv gain matrices
|
||||
LUSolver Finverter; // matrix inversion algorithm
|
||||
Matrix a_init, a_new; // state vector
|
||||
Matrix vt; // current observation error vectors
|
||||
Matrix vtFinv, vtFinvVt; // intermeiate observation error *Finv vector
|
||||
double riccati_tol;
|
||||
InitializeKalmanFilter initKalmanFilter; //Initialise KF matrices
|
||||
Vector FUTP; // F upper triangle packed as vector FUTP(i + (j-1)*j/2) = F(i,j) for 1<=i<=j;
|
||||
|
||||
// Method
|
||||
double filter(const MatrixView &detrendedDataView, const Matrix &H, VectorView &vll, size_t start, int &info);
|
||||
|
|
|
|||
|
|
@ -30,7 +30,7 @@ namespace lapack
|
|||
// calc Cholesky Decomposition (Mat A, char "U"pper/"L"ower)
|
||||
template<class Mat>
|
||||
inline int
|
||||
choleskyDecomp(Mat A, const char *UL)
|
||||
choleskyDecomp(Mat &A, const char *UL)
|
||||
{
|
||||
assert(A.getCols() == A.getRows());
|
||||
lapack_int lpinfo = 0;
|
||||
|
|
@ -42,12 +42,24 @@ namespace lapack
|
|||
}
|
||||
|
||||
|
||||
/**
|
||||
* "DSTEQR computes all eigenvalues and, optionally, eigenvectors of a sym-
|
||||
* metric tridiagonal matrix using the implicit QL or QR method. The eigen-
|
||||
* vectors of a full or band symmetric matrix can also be found if DSYTRD or
|
||||
* DSPTRD or DSBTRD has been used to reduce this matrix to tridiagonal form."
|
||||
*/
|
||||
// calc Cholesky Decomposition based solution X to A*X=B
|
||||
// for A pos. def. and symmetric supplied as uppper/lower triangle
|
||||
// packed in a vector if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
|
||||
// solutino X is returned in B and if B_in=I then B_out=X=inv(A)
|
||||
// A_out contains uupper or lower Cholesky decomposition
|
||||
template<class VecA, class MatB>
|
||||
inline int
|
||||
choleskySolver(VecA &A, MatB &B, const char *UL)
|
||||
{
|
||||
//assert(A.getCols() == A.getRows());
|
||||
lapack_int lpinfo = 0;
|
||||
lapack_int size = B.getRows();
|
||||
lapack_int bcols = B.getCols();
|
||||
lapack_int ldl = B.getLd();
|
||||
dppsv(UL, &size, &bcols, A.getData(), B.getData(), &ldl, &lpinfo);
|
||||
int info = (int) lpinfo;
|
||||
return info;
|
||||
}
|
||||
|
||||
|
||||
} // End of namespace
|
||||
|
|
|
|||
|
|
@ -132,6 +132,7 @@ logposterior(const VectorConstView &estParams, const MatrixConstView &data,
|
|||
double qz_criterium = *mxGetPr(mxGetField(options_, 0, "qz_criterium"));
|
||||
double lyapunov_tol = *mxGetPr(mxGetField(options_, 0, "lyapunov_complex_threshold"));
|
||||
double riccati_tol = *mxGetPr(mxGetField(options_, 0, "riccati_tol"));
|
||||
size_t presample = (size_t) *mxGetPr(mxGetField(options_, 0, "presample"));
|
||||
|
||||
std::vector<size_t> varobs;
|
||||
const mxArray *varobs_mx = mxGetField(options_, 0, "varobs_id");
|
||||
|
|
@ -175,16 +176,17 @@ logposterior(const VectorConstView &estParams, const MatrixConstView &data,
|
|||
Vector deepParams(n_param);
|
||||
deepParams = VectorConstView(mxGetPr(mxGetField(M_, 0, "params")), n_param, 1);
|
||||
Matrix Q(n_exo);
|
||||
Q = MatrixConstView(mxGetPr(mxGetField(M_, 0, "Sigma_e")), n_exo, n_exo, 1);
|
||||
Q = MatrixConstView(mxGetPr(mxGetField(M_, 0, "Sigma_e")), n_exo, n_exo, n_exo);
|
||||
|
||||
Matrix H(n_varobs);
|
||||
const mxArray *H_mx = mxGetField(M_, 0, "H");
|
||||
if (mxGetM(H_mx) == 1 && mxGetN(H_mx) == 1 && *mxGetPr(H_mx) == 0)
|
||||
H.setAll(0.0);
|
||||
else
|
||||
H = MatrixConstView(mxGetPr(mxGetField(M_, 0, "H")), n_varobs, n_varobs, 1);
|
||||
H = MatrixConstView(mxGetPr(mxGetField(M_, 0, "H")), n_varobs, n_varobs, n_varobs);
|
||||
|
||||
// Compute the posterior
|
||||
double logPD = lpd.compute(steadyState, estParams2, deepParams, data, Q, H, 0, info);
|
||||
double logPD =lpd.compute(steadyState, estParams2, deepParams, data, Q, H, presample, info);
|
||||
|
||||
// Cleanups
|
||||
for (std::vector<EstimatedParameter>::iterator it = estParamsInfo.begin();
|
||||
|
|
|
|||
|
|
@ -0,0 +1,345 @@
|
|||
function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations)
|
||||
% function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations)
|
||||
% Evaluates the posterior kernel of a dsge model.
|
||||
%
|
||||
% INPUTS
|
||||
% xparam1 [double] vector of model parameters.
|
||||
% gend [integer] scalar specifying the number of observations.
|
||||
% data [double] matrix of data
|
||||
% data_index [cell] cell of column vectors
|
||||
% number_of_observations [integer]
|
||||
% no_more_missing_observations [integer]
|
||||
% OUTPUTS
|
||||
% fval : value of the posterior kernel at xparam1.
|
||||
% cost_flag : zero if the function returns a penalty, one otherwise.
|
||||
% ys : steady state of original endogenous variables
|
||||
% trend_coeff :
|
||||
% info : vector of informations about the penalty:
|
||||
% 41: one (many) parameter(s) do(es) not satisfied the lower bound
|
||||
% 42: one (many) parameter(s) do(es) not satisfied the upper bound
|
||||
%
|
||||
% SPECIAL REQUIREMENTS
|
||||
%
|
||||
|
||||
% Copyright (C) 2004-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/>.
|
||||
|
||||
global bayestopt_ estim_params_ options_ trend_coeff_ M_ oo_
|
||||
fval = [];
|
||||
ys = [];
|
||||
trend_coeff = [];
|
||||
cost_flag = 1;
|
||||
nobs = size(options_.varobs,1);
|
||||
%------------------------------------------------------------------------------
|
||||
% 1. Get the structural parameters & define penalties
|
||||
%------------------------------------------------------------------------------
|
||||
if options_.mode_compute ~= 1 & any(xparam1 < bayestopt_.lb)
|
||||
k = find(xparam1 < bayestopt_.lb);
|
||||
fval = bayestopt_.penalty+sum((bayestopt_.lb(k)-xparam1(k)).^2);
|
||||
cost_flag = 0;
|
||||
info = 41;
|
||||
return;
|
||||
end
|
||||
if options_.mode_compute ~= 1 & any(xparam1 > bayestopt_.ub)
|
||||
k = find(xparam1 > bayestopt_.ub);
|
||||
fval = bayestopt_.penalty+sum((xparam1(k)-bayestopt_.ub(k)).^2);
|
||||
cost_flag = 0;
|
||||
info = 42;
|
||||
return;
|
||||
end
|
||||
Q = M_.Sigma_e;
|
||||
H = M_.H;
|
||||
for i=1:estim_params_.nvx
|
||||
k =estim_params_.var_exo(i,1);
|
||||
Q(k,k) = xparam1(i)*xparam1(i);
|
||||
end
|
||||
offset = estim_params_.nvx;
|
||||
if estim_params_.nvn
|
||||
for i=1:estim_params_.nvn
|
||||
k = estim_params_.var_endo(i,1);
|
||||
H(k,k) = xparam1(i+offset)*xparam1(i+offset);
|
||||
end
|
||||
offset = offset+estim_params_.nvn;
|
||||
end
|
||||
if estim_params_.ncx
|
||||
for i=1:estim_params_.ncx
|
||||
k1 =estim_params_.corrx(i,1);
|
||||
k2 =estim_params_.corrx(i,2);
|
||||
Q(k1,k2) = xparam1(i+offset)*sqrt(Q(k1,k1)*Q(k2,k2));
|
||||
Q(k2,k1) = Q(k1,k2);
|
||||
end
|
||||
[CholQ,testQ] = chol(Q);
|
||||
if testQ %% The variance-covariance matrix of the structural innovations is not definite positive.
|
||||
%% We have to compute the eigenvalues of this matrix in order to build the penalty.
|
||||
a = diag(eig(Q));
|
||||
k = find(a < 0);
|
||||
if k > 0
|
||||
fval = bayestopt_.penalty+sum(-a(k));
|
||||
cost_flag = 0;
|
||||
info = 43;
|
||||
return
|
||||
end
|
||||
end
|
||||
offset = offset+estim_params_.ncx;
|
||||
end
|
||||
if estim_params_.ncn
|
||||
for i=1:estim_params_.ncn
|
||||
k1 = options_.lgyidx2varobs(estim_params_.corrn(i,1));
|
||||
k2 = options_.lgyidx2varobs(estim_params_.corrn(i,2));
|
||||
H(k1,k2) = xparam1(i+offset)*sqrt(H(k1,k1)*H(k2,k2));
|
||||
H(k2,k1) = H(k1,k2);
|
||||
end
|
||||
[CholH,testH] = chol(H);
|
||||
if testH
|
||||
a = diag(eig(H));
|
||||
k = find(a < 0);
|
||||
if k > 0
|
||||
fval = bayestopt_.penalty+sum(-a(k));
|
||||
cost_flag = 0;
|
||||
info = 44;
|
||||
return
|
||||
end
|
||||
end
|
||||
offset = offset+estim_params_.ncn;
|
||||
end
|
||||
if estim_params_.np > 0
|
||||
M_.params(estim_params_.param_vals(:,1)) = xparam1(offset+1:end);
|
||||
end
|
||||
M_.Sigma_e = Q;
|
||||
M_.H = H;
|
||||
%------------------------------------------------------------------------------
|
||||
% 2. call model setup & reduction program
|
||||
%------------------------------------------------------------------------------
|
||||
[T,R,SteadyState,info] = dynare_resolve(bayestopt_.restrict_var_list,...
|
||||
bayestopt_.restrict_columns,...
|
||||
bayestopt_.restrict_aux);
|
||||
if info(1) == 1 || info(1) == 2 || info(1) == 5
|
||||
fval = bayestopt_.penalty+1;
|
||||
cost_flag = 0;
|
||||
return
|
||||
elseif info(1) == 3 || info(1) == 4 || info(1)==6 ||info(1) == 19 || info(1) == 20 || info(1) == 21
|
||||
fval = bayestopt_.penalty+info(2);
|
||||
cost_flag = 0;
|
||||
return
|
||||
end
|
||||
bayestopt_.mf = bayestopt_.mf1;
|
||||
if options_.noconstant
|
||||
constant = zeros(nobs,1);
|
||||
else
|
||||
if options_.loglinear
|
||||
constant = log(SteadyState(bayestopt_.mfys));
|
||||
else
|
||||
constant = SteadyState(bayestopt_.mfys);
|
||||
end
|
||||
end
|
||||
if bayestopt_.with_trend
|
||||
trend_coeff = zeros(nobs,1);
|
||||
t = options_.trend_coeffs;
|
||||
for i=1:length(t)
|
||||
if ~isempty(t{i})
|
||||
trend_coeff(i) = evalin('base',t{i});
|
||||
end
|
||||
end
|
||||
trend = repmat(constant,1,gend)+trend_coeff*[1:gend];
|
||||
else
|
||||
trend = repmat(constant,1,gend);
|
||||
end
|
||||
start = options_.presample+1;
|
||||
np = size(T,1);
|
||||
mf = bayestopt_.mf;
|
||||
no_missing_data_flag = (number_of_observations==gend*nobs);
|
||||
%------------------------------------------------------------------------------
|
||||
% 3. Initial condition of the Kalman filter
|
||||
%------------------------------------------------------------------------------
|
||||
T
|
||||
R
|
||||
Q
|
||||
R*Q*R'
|
||||
pause
|
||||
options_.lik_init = 1;
|
||||
kalman_algo = options_.kalman_algo;
|
||||
if options_.lik_init == 1 % Kalman filter
|
||||
if kalman_algo ~= 2
|
||||
kalman_algo = 1;
|
||||
end
|
||||
Pstar = lyapunov_symm(T,R*Q*R',options_.qz_criterium,options_.lyapunov_complex_threshold);
|
||||
Pinf = [];
|
||||
elseif options_.lik_init == 2 % Old Diffuse Kalman filter
|
||||
if kalman_algo ~= 2
|
||||
kalman_algo = 1;
|
||||
end
|
||||
Pstar = options_.Harvey_scale_factor*eye(np);
|
||||
Pinf = [];
|
||||
elseif options_.lik_init == 3 % Diffuse Kalman filter
|
||||
if kalman_algo ~= 4
|
||||
kalman_algo = 3;
|
||||
end
|
||||
[QT,ST] = schur(T);
|
||||
e1 = abs(ordeig(ST)) > 2-options_.qz_criterium;
|
||||
[QT,ST] = ordschur(QT,ST,e1);
|
||||
k = find(abs(ordeig(ST)) > 2-options_.qz_criterium);
|
||||
nk = length(k);
|
||||
nk1 = nk+1;
|
||||
Pinf = zeros(np,np);
|
||||
Pinf(1:nk,1:nk) = eye(nk);
|
||||
Pstar = zeros(np,np);
|
||||
B = QT'*R*Q*R'*QT;
|
||||
for i=np:-1:nk+2
|
||||
if ST(i,i-1) == 0
|
||||
if i == np
|
||||
c = zeros(np-nk,1);
|
||||
else
|
||||
c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+...
|
||||
ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i);
|
||||
end
|
||||
q = eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i);
|
||||
Pstar(nk1:i,i) = q\(B(nk1:i,i)+c);
|
||||
Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)';
|
||||
else
|
||||
if i == np
|
||||
c = zeros(np-nk,1);
|
||||
c1 = zeros(np-nk,1);
|
||||
else
|
||||
c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+...
|
||||
ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i)+...
|
||||
ST(i,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1);
|
||||
c1 = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i-1,i+1:end)')+...
|
||||
ST(i-1,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1)+...
|
||||
ST(i-1,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i);
|
||||
end
|
||||
q = [eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i) -ST(nk1:i,nk1:i)*ST(i,i-1);...
|
||||
-ST(nk1:i,nk1:i)*ST(i-1,i) eye(i-nk)-ST(nk1:i,nk1:i)*ST(i-1,i-1)];
|
||||
z = q\[B(nk1:i,i)+c;B(nk1:i,i-1)+c1];
|
||||
Pstar(nk1:i,i) = z(1:(i-nk));
|
||||
Pstar(nk1:i,i-1) = z(i-nk+1:end);
|
||||
Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)';
|
||||
Pstar(i-1,nk1:i-2) = Pstar(nk1:i-2,i-1)';
|
||||
i = i - 1;
|
||||
end
|
||||
end
|
||||
if i == nk+2
|
||||
c = ST(nk+1,:)*(Pstar(:,nk+2:end)*ST(nk1,nk+2:end)')+ST(nk1,nk1)*ST(nk1,nk+2:end)*Pstar(nk+2:end,nk1);
|
||||
Pstar(nk1,nk1)=(B(nk1,nk1)+c)/(1-ST(nk1,nk1)*ST(nk1,nk1));
|
||||
end
|
||||
Z = QT(mf,:);
|
||||
R1 = QT'*R;
|
||||
[QQ,RR,EE] = qr(Z*ST(:,1:nk),0);
|
||||
k = find(abs(diag([RR; zeros(nk-size(Z,1),size(RR,2))])) < 1e-8);
|
||||
if length(k) > 0
|
||||
k1 = EE(:,k);
|
||||
dd =ones(nk,1);
|
||||
dd(k1) = zeros(length(k1),1);
|
||||
Pinf(1:nk,1:nk) = diag(dd);
|
||||
end
|
||||
end
|
||||
if kalman_algo == 2
|
||||
end
|
||||
kalman_tol = options_.kalman_tol;
|
||||
riccati_tol = options_.riccati_tol;
|
||||
mf = bayestopt_.mf1;
|
||||
Y = data-trend;
|
||||
Pstar
|
||||
pause
|
||||
%------------------------------------------------------------------------------
|
||||
% 4. Likelihood evaluation
|
||||
%------------------------------------------------------------------------------
|
||||
if (kalman_algo==1)% Multivariate Kalman Filter
|
||||
if no_missing_data_flag
|
||||
LIK = kalman_filter(T,R,Q,H,Pstar,Y,start,mf,kalman_tol,riccati_tol);
|
||||
else
|
||||
LIK = ...
|
||||
missing_observations_kalman_filter(T,R,Q,H,Pstar,Y,start,mf,kalman_tol,riccati_tol, ...
|
||||
data_index,number_of_observations,no_more_missing_observations);
|
||||
end
|
||||
if isinf(LIK)
|
||||
kalman_algo = 2;
|
||||
end
|
||||
end
|
||||
if (kalman_algo==2)% Univariate Kalman Filter
|
||||
no_correlation_flag = 1;
|
||||
if length(H)==1 & H == 0
|
||||
H = zeros(nobs,1);
|
||||
else
|
||||
if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
|
||||
H = diag(H);
|
||||
else
|
||||
no_correlation_flag = 0;
|
||||
end
|
||||
end
|
||||
if no_correlation_flag
|
||||
LIK = univariate_kalman_filter(T,R,Q,H,Pstar,Y,start,mf,kalman_tol,riccati_tol,data_index,number_of_observations,no_more_missing_observations);
|
||||
else
|
||||
LIK = univariate_kalman_filter_corr(T,R,Q,H,Pstar,Y,start,mf,kalman_tol,riccati_tol,data_index,number_of_observations,no_more_missing_observations);
|
||||
end
|
||||
end
|
||||
if (kalman_algo==3)% Multivariate Diffuse Kalman Filter
|
||||
if no_missing_data_flag
|
||||
LIK = diffuse_kalman_filter(ST,R1,Q,H,Pinf,Pstar,Y,start,Z,kalman_tol, ...
|
||||
riccati_tol);
|
||||
else
|
||||
LIK = missing_observations_diffuse_kalman_filter(ST,R1,Q,H,Pinf, ...
|
||||
Pstar,Y,start,Z,kalman_tol,riccati_tol,...
|
||||
data_index,number_of_observations,...
|
||||
no_more_missing_observations);
|
||||
end
|
||||
if isinf(LIK)
|
||||
kalman_algo = 4;
|
||||
end
|
||||
end
|
||||
if (kalman_algo==4)% Univariate Diffuse Kalman Filter
|
||||
no_correlation_flag = 1;
|
||||
if length(H)==1 & H == 0
|
||||
H = zeros(nobs,1);
|
||||
else
|
||||
if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
|
||||
H = diag(H);
|
||||
else
|
||||
no_correlation_flag = 0;
|
||||
end
|
||||
end
|
||||
if no_correlation_flag
|
||||
LIK = univariate_diffuse_kalman_filter(ST,R1,Q,H,Pinf,Pstar,Y, ...
|
||||
start,Z,kalman_tol,riccati_tol,data_index,...
|
||||
number_of_observations,no_more_missing_observations);
|
||||
else
|
||||
LIK = univariate_diffuse_kalman_filter_corr(ST,R1,Q,H,Pinf,Pstar, ...
|
||||
Y,start,Z,kalman_tol,riccati_tol,...
|
||||
data_index,number_of_observations,...
|
||||
no_more_missing_observations);
|
||||
end
|
||||
end
|
||||
if isnan(LIK)
|
||||
cost_flag = 0;
|
||||
return
|
||||
end
|
||||
if imag(LIK)~=0
|
||||
likelihood = bayestopt_.penalty;
|
||||
else
|
||||
likelihood = LIK;
|
||||
end
|
||||
% ------------------------------------------------------------------------------
|
||||
% Adds prior if necessary
|
||||
% ------------------------------------------------------------------------------
|
||||
lnprior = priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
|
||||
fval = (likelihood-lnprior);
|
||||
likelihood
|
||||
lnprior
|
||||
fval
|
||||
pause
|
||||
LIKDLL=logposterior(xparam1,Y,mexext)
|
||||
pause
|
||||
options_.kalman_algo = kalman_algo;
|
||||
|
|
@ -0,0 +1,967 @@
|
|||
C =[
|
||||
-7.4073
|
||||
-6.1860
|
||||
-6.5983
|
||||
-5.6088
|
||||
-5.0547
|
||||
-4.4774
|
||||
-3.8081
|
||||
-3.8425
|
||||
-2.4178
|
||||
-1.9835
|
||||
-1.0395
|
||||
-0.1583
|
||||
-0.0397
|
||||
0.3505
|
||||
-0.1879
|
||||
-0.0067
|
||||
0.0478
|
||||
-1.2247
|
||||
-1.4349
|
||||
-0.7973
|
||||
-0.0461
|
||||
0.5844
|
||||
1.1372
|
||||
1.3801
|
||||
1.8023
|
||||
2.2972
|
||||
2.0469
|
||||
2.5435
|
||||
2.8169
|
||||
3.2007
|
||||
2.6705
|
||||
3.0518
|
||||
3.2445
|
||||
3.8443
|
||||
3.8525
|
||||
4.9494
|
||||
4.2770
|
||||
4.9532
|
||||
5.1441
|
||||
3.7124
|
||||
3.9880
|
||||
3.6926
|
||||
2.6005
|
||||
1.8679
|
||||
1.9085
|
||||
1.5563
|
||||
1.2308
|
||||
0.3264
|
||||
-0.2208
|
||||
-0.2483
|
||||
-0.4082
|
||||
-1.0315
|
||||
-1.6030
|
||||
-1.5499
|
||||
-1.3777
|
||||
-2.1675
|
||||
-2.5138
|
||||
-2.8820
|
||||
-2.6958
|
||||
-2.4719
|
||||
-1.9854
|
||||
-1.7954
|
||||
-2.2362
|
||||
-1.0595
|
||||
-0.8808
|
||||
-0.8548
|
||||
-1.2839
|
||||
-0.1363
|
||||
0.2104
|
||||
0.8810
|
||||
0.3555
|
||||
0.4766
|
||||
1.3269
|
||||
1.4506
|
||||
1.4308
|
||||
1.6263
|
||||
1.9842
|
||||
2.3948
|
||||
2.8710
|
||||
3.0177
|
||||
2.9305
|
||||
3.1739
|
||||
3.7380
|
||||
3.8285
|
||||
3.3342
|
||||
3.7447
|
||||
3.7830
|
||||
3.1039
|
||||
2.8413
|
||||
3.0338
|
||||
0.3669
|
||||
0.0847
|
||||
0.0104
|
||||
0.2115
|
||||
-0.6649
|
||||
-0.9625
|
||||
-0.7330
|
||||
-0.8664
|
||||
-1.4441
|
||||
-1.0179
|
||||
-1.2729
|
||||
-1.9539
|
||||
-1.4427
|
||||
-2.0371
|
||||
-1.9764
|
||||
-2.5654
|
||||
-2.8570
|
||||
-2.5842
|
||||
-3.0427
|
||||
-2.8312
|
||||
-2.3320
|
||||
-2.2768
|
||||
-2.1816
|
||||
-2.1043
|
||||
-1.8969
|
||||
-2.2388
|
||||
-2.1679
|
||||
-2.1172
|
||||
];
|
||||
|
||||
E =[
|
||||
0.6263
|
||||
0.7368
|
||||
0.7477
|
||||
1.0150
|
||||
0.6934
|
||||
0.4135
|
||||
0.3845
|
||||
0.2380
|
||||
0.2853
|
||||
0.5999
|
||||
0.8622
|
||||
1.2116
|
||||
1.4921
|
||||
1.5816
|
||||
1.7259
|
||||
1.6276
|
||||
1.2422
|
||||
0.8084
|
||||
0.4710
|
||||
-0.3704
|
||||
-0.6427
|
||||
-0.5323
|
||||
-0.5562
|
||||
-0.3651
|
||||
-0.4356
|
||||
-0.7164
|
||||
-0.5816
|
||||
-0.4635
|
||||
-0.8456
|
||||
-0.9708
|
||||
-0.7138
|
||||
-0.7499
|
||||
-0.6941
|
||||
-0.6656
|
||||
-0.2912
|
||||
-0.1650
|
||||
0.0774
|
||||
0.2307
|
||||
0.4484
|
||||
0.4942
|
||||
0.4653
|
||||
0.2196
|
||||
0.1736
|
||||
-0.1595
|
||||
-0.3918
|
||||
-0.4611
|
||||
-0.8493
|
||||
-0.7384
|
||||
-1.0604
|
||||
-1.2166
|
||||
-1.7187
|
||||
-1.6932
|
||||
-1.7830
|
||||
-1.7035
|
||||
-2.2079
|
||||
-2.3769
|
||||
-2.2511
|
||||
-2.1093
|
||||
-2.4638
|
||||
-2.4027
|
||||
-2.1313
|
||||
-1.9199
|
||||
-1.7941
|
||||
-1.4661
|
||||
-1.2269
|
||||
-1.0392
|
||||
-1.0725
|
||||
-0.7156
|
||||
-0.4778
|
||||
-0.4233
|
||||
-0.0409
|
||||
0.1620
|
||||
0.4280
|
||||
0.5873
|
||||
1.0323
|
||||
1.3420
|
||||
1.6902
|
||||
2.0680
|
||||
2.8219
|
||||
3.2511
|
||||
3.2930
|
||||
3.5633
|
||||
3.8992
|
||||
3.6874
|
||||
3.2849
|
||||
3.1614
|
||||
2.6221
|
||||
2.5067
|
||||
1.9223
|
||||
1.1777
|
||||
0.4483
|
||||
-0.0661
|
||||
-0.4424
|
||||
-0.9000
|
||||
-1.1478
|
||||
-1.2047
|
||||
-1.1412
|
||||
-1.2383
|
||||
-1.1048
|
||||
-0.9716
|
||||
-0.9287
|
||||
-1.0057
|
||||
-1.0827
|
||||
-1.0200
|
||||
-1.0072
|
||||
-1.1740
|
||||
-1.2809
|
||||
-1.1086
|
||||
-0.9866
|
||||
-0.8947
|
||||
-0.5875
|
||||
-0.2329
|
||||
0.1493
|
||||
0.4906
|
||||
0.8400
|
||||
1.0720
|
||||
1.2648
|
||||
1.5431
|
||||
];
|
||||
|
||||
I =[
|
||||
2.6617
|
||||
2.4325
|
||||
1.9592
|
||||
3.2530
|
||||
2.9949
|
||||
3.7918
|
||||
4.7444
|
||||
4.8289
|
||||
5.5983
|
||||
7.8923
|
||||
9.4297
|
||||
9.5010
|
||||
10.0150
|
||||
10.0413
|
||||
9.6046
|
||||
6.4766
|
||||
5.9647
|
||||
3.0114
|
||||
0.5683
|
||||
-2.1226
|
||||
-2.1855
|
||||
-0.8329
|
||||
-1.5207
|
||||
-1.3419
|
||||
-1.7897
|
||||
-0.1476
|
||||
0.4675
|
||||
-1.6516
|
||||
-1.5419
|
||||
-1.3050
|
||||
-1.2451
|
||||
-0.7815
|
||||
-0.7796
|
||||
-0.3612
|
||||
-2.4072
|
||||
1.1162
|
||||
1.1383
|
||||
3.4132
|
||||
5.0356
|
||||
2.8016
|
||||
2.1734
|
||||
0.9366
|
||||
-0.7050
|
||||
-1.5021
|
||||
-2.9868
|
||||
-6.0237
|
||||
-6.2589
|
||||
-6.9138
|
||||
-8.2340
|
||||
-9.2589
|
||||
-9.2465
|
||||
-9.6988
|
||||
-9.7782
|
||||
-10.5645
|
||||
-10.7544
|
||||
-13.1583
|
||||
-12.2718
|
||||
-12.0131
|
||||
-13.5983
|
||||
-12.3579
|
||||
-10.9146
|
||||
-11.1572
|
||||
-12.4935
|
||||
-9.4393
|
||||
-8.5535
|
||||
-7.3723
|
||||
-10.0169
|
||||
-6.6088
|
||||
-5.2045
|
||||
-4.1024
|
||||
-2.8472
|
||||
-1.3139
|
||||
0.0477
|
||||
1.5629
|
||||
3.6947
|
||||
4.0327
|
||||
4.1320
|
||||
7.1400
|
||||
9.1036
|
||||
8.5609
|
||||
7.6576
|
||||
8.8022
|
||||
8.9611
|
||||
10.0871
|
||||
9.4797
|
||||
9.3964
|
||||
10.0363
|
||||
8.6340
|
||||
6.6522
|
||||
4.4471
|
||||
0.2854
|
||||
-2.1879
|
||||
-2.9879
|
||||
-4.1021
|
||||
-2.7713
|
||||
-2.2281
|
||||
-1.2908
|
||||
-0.3250
|
||||
0.6534
|
||||
0.3942
|
||||
0.3534
|
||||
-0.1532
|
||||
-1.7936
|
||||
0.4909
|
||||
0.3634
|
||||
0.4290
|
||||
-0.9709
|
||||
0.1942
|
||||
0.6103
|
||||
1.4426
|
||||
2.7225
|
||||
1.7525
|
||||
3.2780
|
||||
3.5985
|
||||
4.9011
|
||||
5.3312
|
||||
6.4402
|
||||
6.6529
|
||||
];
|
||||
|
||||
L =[
|
||||
0.6263
|
||||
0.7368
|
||||
0.7477
|
||||
1.0150
|
||||
0.6934
|
||||
0.4135
|
||||
0.3845
|
||||
0.2380
|
||||
0.2853
|
||||
0.5999
|
||||
0.8622
|
||||
1.2116
|
||||
1.4921
|
||||
1.5816
|
||||
1.7259
|
||||
1.6276
|
||||
1.2422
|
||||
0.8084
|
||||
0.4710
|
||||
-0.3704
|
||||
-0.6427
|
||||
-0.5323
|
||||
-0.5562
|
||||
-0.3651
|
||||
-0.4356
|
||||
-0.7164
|
||||
-0.5816
|
||||
-0.4635
|
||||
-0.8456
|
||||
-0.9708
|
||||
-0.7138
|
||||
-0.7499
|
||||
-0.6941
|
||||
-0.6656
|
||||
-0.2912
|
||||
-0.1650
|
||||
0.0774
|
||||
0.2307
|
||||
0.4484
|
||||
0.4942
|
||||
0.4653
|
||||
0.2196
|
||||
0.1736
|
||||
-0.1595
|
||||
-0.3918
|
||||
-0.4611
|
||||
-0.8493
|
||||
-0.7384
|
||||
-1.0604
|
||||
-1.2166
|
||||
-1.7187
|
||||
-1.6932
|
||||
-1.7830
|
||||
-1.7035
|
||||
-2.2079
|
||||
-2.3769
|
||||
-2.2511
|
||||
-2.1093
|
||||
-2.4638
|
||||
-2.4027
|
||||
-2.1313
|
||||
-1.9199
|
||||
-1.7941
|
||||
-1.4661
|
||||
-1.2269
|
||||
-1.0392
|
||||
-1.0725
|
||||
-0.7156
|
||||
-0.4778
|
||||
-0.4233
|
||||
-0.0409
|
||||
0.1620
|
||||
0.4280
|
||||
0.5873
|
||||
1.0323
|
||||
1.3420
|
||||
1.6902
|
||||
2.0680
|
||||
2.8219
|
||||
3.2511
|
||||
3.2930
|
||||
3.5633
|
||||
3.8992
|
||||
3.6874
|
||||
3.2849
|
||||
3.1614
|
||||
2.6221
|
||||
2.5067
|
||||
1.9223
|
||||
1.1777
|
||||
0.4483
|
||||
-0.0661
|
||||
-0.4424
|
||||
-0.9000
|
||||
-1.1478
|
||||
-1.2047
|
||||
-1.1412
|
||||
-1.2383
|
||||
-1.1048
|
||||
-0.9716
|
||||
-0.9287
|
||||
-1.0057
|
||||
-1.0827
|
||||
-1.0200
|
||||
-1.0072
|
||||
-1.1740
|
||||
-1.2809
|
||||
-1.1086
|
||||
-0.9866
|
||||
-0.8947
|
||||
-0.5875
|
||||
-0.2329
|
||||
0.1493
|
||||
0.4906
|
||||
0.8400
|
||||
1.0720
|
||||
1.2648
|
||||
1.5431
|
||||
];
|
||||
|
||||
PIE =[
|
||||
-1.0113
|
||||
-0.8305
|
||||
0.2332
|
||||
-0.8746
|
||||
-0.7978
|
||||
-0.9220
|
||||
-0.2487
|
||||
-0.7749
|
||||
-0.5460
|
||||
-0.5347
|
||||
0.5050
|
||||
-0.0334
|
||||
0.6756
|
||||
0.8791
|
||||
0.7267
|
||||
1.0997
|
||||
1.1750
|
||||
1.1927
|
||||
0.4420
|
||||
0.5357
|
||||
0.0345
|
||||
0.0196
|
||||
0.3371
|
||||
0.9379
|
||||
1.2160
|
||||
0.3393
|
||||
0.5813
|
||||
0.7410
|
||||
0.3374
|
||||
0.2616
|
||||
0.4025
|
||||
0.4799
|
||||
0.5981
|
||||
-0.1523
|
||||
0.4458
|
||||
0.2182
|
||||
0.9793
|
||||
0.7562
|
||||
1.0064
|
||||
0.8203
|
||||
0.6966
|
||||
0.3352
|
||||
0.6581
|
||||
0.6111
|
||||
0.9833
|
||||
1.1991
|
||||
0.9562
|
||||
0.3868
|
||||
0.2939
|
||||
0.2471
|
||||
0.8331
|
||||
0.0715
|
||||
0.3910
|
||||
0.3301
|
||||
0.2547
|
||||
-0.2702
|
||||
-0.2998
|
||||
-0.1953
|
||||
-0.2293
|
||||
-0.3284
|
||||
0.0480
|
||||
-0.0374
|
||||
0.3253
|
||||
-0.3434
|
||||
-0.3892
|
||||
-0.7178
|
||||
-0.4758
|
||||
-0.6794
|
||||
-0.8505
|
||||
-0.3512
|
||||
-0.4436
|
||||
-0.5101
|
||||
-0.4574
|
||||
-0.2696
|
||||
-0.1047
|
||||
-0.5745
|
||||
-0.2989
|
||||
-0.0063
|
||||
0.0088
|
||||
-0.1184
|
||||
-0.1506
|
||||
-0.4073
|
||||
0.2674
|
||||
0.2896
|
||||
0.0669
|
||||
0.1166
|
||||
-0.1699
|
||||
-0.2518
|
||||
-0.0562
|
||||
-0.3269
|
||||
-0.0703
|
||||
-0.1046
|
||||
-0.4888
|
||||
-0.3524
|
||||
-0.2485
|
||||
-0.5870
|
||||
-0.4546
|
||||
-0.3970
|
||||
-0.2353
|
||||
-0.0352
|
||||
-0.2171
|
||||
-0.3754
|
||||
-0.4322
|
||||
-0.4572
|
||||
-0.4903
|
||||
-0.4518
|
||||
-0.6435
|
||||
-0.6304
|
||||
-0.4148
|
||||
-0.2892
|
||||
-0.4318
|
||||
-0.6010
|
||||
-0.4148
|
||||
-0.4315
|
||||
-0.3531
|
||||
-0.8053
|
||||
-0.4680
|
||||
-0.4263
|
||||
];
|
||||
|
||||
R =[
|
||||
-1.0750
|
||||
-1.1540
|
||||
-1.3682
|
||||
-1.4569
|
||||
-1.3490
|
||||
-1.4011
|
||||
-1.6486
|
||||
-1.6968
|
||||
-1.6976
|
||||
-1.2567
|
||||
-1.1392
|
||||
-0.7783
|
||||
-0.3021
|
||||
-0.0435
|
||||
0.0066
|
||||
-0.0043
|
||||
0.1029
|
||||
-0.0628
|
||||
-0.5358
|
||||
-0.9627
|
||||
-1.1079
|
||||
-1.0918
|
||||
-0.9966
|
||||
-0.6223
|
||||
-0.3616
|
||||
-0.2711
|
||||
-0.0997
|
||||
-0.2810
|
||||
-0.3710
|
||||
-0.3167
|
||||
-0.5301
|
||||
-0.5826
|
||||
-0.3194
|
||||
-0.2713
|
||||
-0.5287
|
||||
-0.2432
|
||||
0.1098
|
||||
0.5349
|
||||
0.7094
|
||||
0.8415
|
||||
0.6226
|
||||
0.7376
|
||||
0.9316
|
||||
1.4370
|
||||
1.5853
|
||||
1.4267
|
||||
1.1783
|
||||
1.2046
|
||||
0.9689
|
||||
0.7918
|
||||
0.6315
|
||||
0.5950
|
||||
0.6853
|
||||
0.7171
|
||||
0.5887
|
||||
0.4873
|
||||
0.4027
|
||||
0.3489
|
||||
0.2934
|
||||
0.3060
|
||||
0.1741
|
||||
0.0348
|
||||
0.0771
|
||||
-0.1005
|
||||
-0.1518
|
||||
-0.1104
|
||||
-0.0681
|
||||
-0.0059
|
||||
0.0256
|
||||
0.0404
|
||||
-0.1721
|
||||
-0.2002
|
||||
0.0015
|
||||
0.1249
|
||||
0.3738
|
||||
0.4320
|
||||
0.5579
|
||||
0.8186
|
||||
0.8727
|
||||
0.7356
|
||||
0.7243
|
||||
0.8635
|
||||
0.9058
|
||||
0.7656
|
||||
0.7936
|
||||
0.8631
|
||||
0.9074
|
||||
0.9547
|
||||
1.2045
|
||||
1.0850
|
||||
0.9178
|
||||
0.5242
|
||||
0.3178
|
||||
0.1472
|
||||
0.0227
|
||||
-0.0799
|
||||
-0.0611
|
||||
-0.0140
|
||||
0.1132
|
||||
0.1774
|
||||
0.0782
|
||||
0.0436
|
||||
-0.1596
|
||||
-0.2691
|
||||
-0.2895
|
||||
-0.3791
|
||||
-0.4020
|
||||
-0.4166
|
||||
-0.4037
|
||||
-0.3636
|
||||
-0.4075
|
||||
-0.4311
|
||||
-0.4470
|
||||
-0.5111
|
||||
-0.6274
|
||||
-0.7261
|
||||
-0.6974
|
||||
-0.5012
|
||||
];
|
||||
|
||||
W =[
|
||||
-14.8791
|
||||
-13.2300
|
||||
-13.5037
|
||||
-13.0249
|
||||
-11.2546
|
||||
-10.0148
|
||||
-8.8586
|
||||
-8.5739
|
||||
-7.7851
|
||||
-6.7136
|
||||
-5.5878
|
||||
-4.6881
|
||||
-3.8039
|
||||
-3.0366
|
||||
-2.7342
|
||||
-1.3135
|
||||
-0.7387
|
||||
-0.1131
|
||||
-0.2769
|
||||
0.8696
|
||||
1.8855
|
||||
2.3667
|
||||
2.4942
|
||||
3.2049
|
||||
3.9682
|
||||
5.1500
|
||||
4.7047
|
||||
4.7827
|
||||
5.3377
|
||||
5.6614
|
||||
5.2813
|
||||
5.2967
|
||||
5.5175
|
||||
6.1526
|
||||
5.6627
|
||||
6.0694
|
||||
6.5824
|
||||
6.9032
|
||||
6.7849
|
||||
6.6896
|
||||
6.6201
|
||||
6.9933
|
||||
5.8959
|
||||
6.7419
|
||||
6.9999
|
||||
6.4009
|
||||
5.5083
|
||||
5.1054
|
||||
5.2813
|
||||
4.5790
|
||||
3.9589
|
||||
3.8599
|
||||
3.8978
|
||||
2.7957
|
||||
3.2480
|
||||
1.4634
|
||||
1.9219
|
||||
1.8398
|
||||
1.9279
|
||||
1.8316
|
||||
1.6092
|
||||
1.2741
|
||||
0.2031
|
||||
-0.0236
|
||||
-0.1004
|
||||
-0.3034
|
||||
-1.0273
|
||||
-0.2205
|
||||
0.0458
|
||||
0.2386
|
||||
-0.0977
|
||||
-0.3145
|
||||
-0.1416
|
||||
-0.7009
|
||||
-0.9082
|
||||
-0.8802
|
||||
-0.5644
|
||||
-0.5852
|
||||
-0.5346
|
||||
0.0652
|
||||
0.1301
|
||||
0.3444
|
||||
-0.3592
|
||||
0.8096
|
||||
0.9644
|
||||
1.0289
|
||||
1.2781
|
||||
1.2298
|
||||
2.2134
|
||||
2.0808
|
||||
0.4925
|
||||
0.6506
|
||||
0.5531
|
||||
0.2456
|
||||
-0.5351
|
||||
-0.8183
|
||||
-0.8967
|
||||
-0.7268
|
||||
-1.0738
|
||||
-1.2844
|
||||
-1.4338
|
||||
-1.6995
|
||||
-1.7085
|
||||
-2.2889
|
||||
-2.1018
|
||||
-2.4273
|
||||
-2.4609
|
||||
-2.1407
|
||||
-2.3847
|
||||
-3.1689
|
||||
-4.5581
|
||||
-4.1027
|
||||
-4.2436
|
||||
-4.8836
|
||||
-5.9660
|
||||
-4.9971
|
||||
-5.2386
|
||||
-5.6618
|
||||
];
|
||||
|
||||
Y =[
|
||||
-4.9347
|
||||
-4.6205
|
||||
-5.2198
|
||||
-4.5937
|
||||
-3.8015
|
||||
-3.6643
|
||||
-2.7239
|
||||
-2.7524
|
||||
-2.0634
|
||||
-1.0112
|
||||
0.0530
|
||||
0.7623
|
||||
1.7927
|
||||
2.1486
|
||||
2.4866
|
||||
2.1456
|
||||
2.1671
|
||||
-0.0254
|
||||
-1.6716
|
||||
-1.9673
|
||||
-1.6109
|
||||
-1.0292
|
||||
-0.1222
|
||||
0.7329
|
||||
1.1234
|
||||
2.0603
|
||||
1.7998
|
||||
1.4820
|
||||
1.1732
|
||||
1.6424
|
||||
1.5382
|
||||
2.1399
|
||||
2.0127
|
||||
2.7210
|
||||
2.4966
|
||||
3.5249
|
||||
3.6237
|
||||
4.2011
|
||||
4.5634
|
||||
3.3442
|
||||
2.7761
|
||||
1.9812
|
||||
1.3779
|
||||
1.4616
|
||||
1.3029
|
||||
0.7594
|
||||
0.3695
|
||||
0.0832
|
||||
-0.8118
|
||||
-1.4557
|
||||
-1.4850
|
||||
-1.2346
|
||||
-1.5696
|
||||
-1.3785
|
||||
-0.7682
|
||||
-2.0308
|
||||
-1.7778
|
||||
-1.7801
|
||||
-2.1711
|
||||
-1.7469
|
||||
-1.3413
|
||||
-1.3352
|
||||
-2.4390
|
||||
-1.2125
|
||||
-1.1695
|
||||
-1.0891
|
||||
-2.4753
|
||||
-1.3503
|
||||
-0.9412
|
||||
-0.1470
|
||||
0.0026
|
||||
0.1108
|
||||
0.6890
|
||||
1.3520
|
||||
1.6018
|
||||
2.0667
|
||||
1.7625
|
||||
2.6658
|
||||
3.4048
|
||||
3.2507
|
||||
3.4251
|
||||
3.2174
|
||||
3.1903
|
||||
3.3396
|
||||
3.1358
|
||||
2.8625
|
||||
3.3546
|
||||
2.4609
|
||||
1.9534
|
||||
0.9962
|
||||
-0.7904
|
||||
-1.1672
|
||||
-1.2586
|
||||
-1.3593
|
||||
-1.3443
|
||||
-0.9413
|
||||
-0.6023
|
||||
-0.4516
|
||||
-0.5129
|
||||
-0.8741
|
||||
-1.0784
|
||||
-1.4091
|
||||
-1.3627
|
||||
-1.5731
|
||||
-1.6037
|
||||
-1.8814
|
||||
-2.1482
|
||||
-1.3597
|
||||
-1.1855
|
||||
-1.1122
|
||||
-0.8424
|
||||
-0.9747
|
||||
-1.1385
|
||||
-1.4548
|
||||
-1.4284
|
||||
-1.4633
|
||||
-1.0621
|
||||
-0.7871
|
||||
];
|
||||
|
|
@ -0,0 +1,182 @@
|
|||
//options_.usePartInfo=1;
|
||||
|
||||
var MC E EF R_KF QF CF IF YF LF PIEF WF RF R_K Q C I Y L PIE W R EE_A PIE_BAR EE_B EE_G EE_L EE_I KF K one BIGTHETA;
|
||||
|
||||
varexo E_A E_B E_G E_L E_I ETA_R E_PIE_BAR ETA_Q ETA_P ETA_W ;
|
||||
|
||||
parameters
|
||||
xi_e lambda_w alpha czcap beta phi_i tau sig_c hab ccs cinvs phi_y gamma_w xi_w gamma_p xi_p sig_l r_dpi
|
||||
r_pie r_dy r_y rho rho_a rho_pb rho_b rho_g rho_l rho_i LMP ;
|
||||
|
||||
|
||||
|
||||
alpha=.30;
|
||||
beta=.99;
|
||||
tau=0.025;
|
||||
ccs=0.6;
|
||||
cinvs=.22; //% alpha*(tau+ctrend)/R_K R_K=ctrend/beta-1+tau
|
||||
lambda_w = 0.5;
|
||||
phi_i= 6.771;
|
||||
sig_c= 1.353;
|
||||
hab= 0.573;
|
||||
xi_w= 0.737;
|
||||
sig_l= 2.400;
|
||||
xi_p= 0.908;
|
||||
xi_e= 0.599;
|
||||
gamma_w= 0.763;
|
||||
gamma_p= 0.469;
|
||||
czcap= 0.169;
|
||||
phi_y= 1.408;
|
||||
r_pie= 1.684;
|
||||
r_dpi= 0.14;
|
||||
rho= 0.961;
|
||||
r_y= 0.099;
|
||||
r_dy= 0.159;
|
||||
rho_a= 0.823;
|
||||
rho_b= 0.855;
|
||||
rho_g= 0.949;
|
||||
rho_l= 0.889;
|
||||
rho_i= 0.927;
|
||||
rho_pb= 0.924;
|
||||
LMP = 0.0 ; //NEW.
|
||||
|
||||
model(linear, use_dll);
|
||||
CF = (1/(1+hab))*(CF(1)+hab*CF(-1))-((1-hab)/((1+hab)*sig_c))*(RF-PIEF(1)-EE_B) ;
|
||||
0 = alpha*R_KF+(1-alpha)*WF -EE_A ;
|
||||
PIEF = 0*one;
|
||||
IF = (1/(1+beta))* (( IF(-1) + beta*(IF(1)))+(1/phi_i)*QF)+0*ETA_Q+EE_I ;
|
||||
QF = -(RF-PIEF(1))+(1-beta*(1-tau))*((0+czcap)/czcap)*R_KF(1)+beta*(1-tau)*QF(1) +0*EE_I ;
|
||||
KF = (1-tau)*KF(-1)+tau*IF(-1) ;
|
||||
YF = (ccs*CF+cinvs*IF)+EE_G ;
|
||||
|
||||
YF = 1*phi_y*( alpha*KF+alpha*(1/czcap)*R_KF+(1-alpha)*LF+EE_A ) ;
|
||||
WF = (sig_c/(1-hab))*(CF-hab*CF(-1)) + sig_l*LF - EE_L ;
|
||||
LF = R_KF*((1+czcap)/czcap)-WF+KF ;
|
||||
EF = EF(-1)+EF(1)-EF+(LF-EF)*((1-xi_e)*(1-xi_e*beta)/(xi_e));
|
||||
|
||||
C = (hab/(1+hab))*C(-1)+(1/(1+hab))*C(1)-((1-hab)/((1+hab)*sig_c))*(R-PIE(1)-EE_B) ;
|
||||
I = (1/(1+beta))* (( I(-1) + beta*(I(1)))+(1/phi_i)*Q )+1*ETA_Q+1*EE_I ;
|
||||
Q = -(R-PIE(1))+(1-beta*(1-tau))*((0+czcap)/czcap)*R_K(1)+beta*(1-tau)*Q(1) +EE_I*0+0*ETA_Q ;
|
||||
K = (1-tau)*K(-1)+tau*I(-1) ;
|
||||
Y = (ccs*C+cinvs*I)+ EE_G ;
|
||||
Y = phi_y*( alpha*K+alpha*(1/czcap)*R_K+(1-alpha)*L ) +phi_y*EE_A ;
|
||||
PIE = (1/(1+beta*gamma_p))*
|
||||
(
|
||||
(beta)*(PIE(1)) +(gamma_p)*(PIE(-1))
|
||||
+((1-xi_p)*(1-beta*xi_p)/(xi_p))*(MC)
|
||||
) + ETA_P ;
|
||||
|
||||
MC = alpha*R_K+(1-alpha)*W -EE_A;
|
||||
W = (1/(1+beta))*(beta*W(+1)+W(-1))
|
||||
+(beta/(1+beta))*(PIE(+1))
|
||||
-((1+beta*gamma_w)/(1+beta))*(PIE)
|
||||
+(gamma_w/(1+beta))*(PIE(-1))
|
||||
-(1/(1+beta))*(((1-beta*xi_w)*(1-xi_w))/(((1+(((1+lambda_w)*sig_l)/(lambda_w))))*xi_w))*(W-sig_l*L-(sig_c/(1-hab))*(C-hab*C(-1))+EE_L)
|
||||
+ETA_W;
|
||||
L = R_K*((1+czcap)/czcap)-W+K ;
|
||||
|
||||
// R = r_dpi*(PIE-PIE(-1))
|
||||
// +(1-rho)*(r_pie*(PIE(-1)-PIE_BAR)+r_y*(Y-YF))
|
||||
// +r_dy*(Y-YF-(Y(-1)-YF(-1)))
|
||||
// +rho*(R(-1)-PIE_BAR)
|
||||
// +PIE_BAR
|
||||
// +ETA_R;
|
||||
|
||||
|
||||
R =
|
||||
|
||||
r_dpi*(PIE-PIE(-1))
|
||||
|
||||
+(1-rho)*(r_pie*(BIGTHETA)+r_y*(Y-YF))
|
||||
+r_dy*(Y-YF-(Y(-1)-YF(-1)))
|
||||
+rho*(R(-1)-PIE_BAR)
|
||||
+PIE_BAR
|
||||
+ETA_R;
|
||||
|
||||
|
||||
E = E(-1)+E(1)-E+(L-E)*((1-xi_e)*(1-xi_e*beta)/(xi_e));
|
||||
|
||||
|
||||
EE_A = (rho_a)*EE_A(-1) + E_A;
|
||||
PIE_BAR = rho_pb*PIE_BAR(-1)+ E_PIE_BAR ;
|
||||
EE_B = rho_b*EE_B(-1) + E_B ;
|
||||
EE_G = rho_g*EE_G(-1) + E_G ;
|
||||
EE_L = rho_l*EE_L(-1) + E_L ;
|
||||
EE_I = rho_i*EE_I(-1) + E_I ;
|
||||
one = 0*one(-1) ;
|
||||
|
||||
LMP*BIGTHETA(1) = BIGTHETA - (PIE(-1) - PIE_BAR) ;
|
||||
|
||||
end;
|
||||
|
||||
|
||||
shocks;
|
||||
var E_A; stderr 0.598;
|
||||
var E_B; stderr 0.336;
|
||||
var E_G; stderr 0.325;
|
||||
var E_I; stderr 0.085;
|
||||
var E_L; stderr 3.520;
|
||||
var ETA_P; stderr 0.160;
|
||||
var ETA_W; stderr 0.289;
|
||||
var ETA_R; stderr 0.081;
|
||||
var ETA_Q; stderr 0.604;
|
||||
var E_PIE_BAR; stderr 0.017;
|
||||
end;
|
||||
|
||||
//stoch_simul(irf=20) Y C PIE R W R_K L Q I K ;
|
||||
|
||||
// stoch_simul generates what kind of standard errors for the shocks ?
|
||||
|
||||
//steady;
|
||||
//check;
|
||||
//stoch_simul(periods=200,irf=20,simul_seed=3) Y C PIE MC R W R_K E L I ;
|
||||
|
||||
//datatomfile('ddd',[]);
|
||||
|
||||
// new syntax
|
||||
|
||||
estimated_params;
|
||||
// PARAM NAME, INITVAL, LB, UB, PRIOR_SHAPE, PRIOR_P1, PRIOR_P2, PRIOR_P3, PRIOR_P4, JSCALE
|
||||
// PRIOR_SHAPE: BETA_PDF, GAMMA_PDF, NORMAL_PDF, INV_GAMMA_PDF
|
||||
stderr E_A,0.543,0.01,4,INV_GAMMA_PDF,0.4,2;
|
||||
stderr E_PIE_BAR,0.072,0.001,4,INV_GAMMA_PDF,0.02,10;
|
||||
stderr E_B,0.2694,0.01,4,INV_GAMMA_PDF,0.2,2;
|
||||
stderr E_G,0.3052,0.01,4,INV_GAMMA_PDF,0.3,2;
|
||||
stderr E_L,1.4575,0.1,6,INV_GAMMA_PDF,1,2;
|
||||
stderr E_I,0.1318,0.01,4,INV_GAMMA_PDF,0.1,2;
|
||||
stderr ETA_R,0.1363,0.01,4,INV_GAMMA_PDF,0.1,2;
|
||||
stderr ETA_Q,0.4842,0.01,4,INV_GAMMA_PDF,0.4,2;
|
||||
stderr ETA_P,0.1731,0.01,4,INV_GAMMA_PDF,0.15,2;
|
||||
stderr ETA_W,0.2462,0.1,4,INV_GAMMA_PDF,0.25,2;
|
||||
rho_a,.9722,.1,.9999,BETA_PDF,0.85,0.1;
|
||||
rho_pb,.85,.1,.999,BETA_PDF,0.85,0.1;
|
||||
rho_b,.7647,.1,.99,BETA_PDF,0.85,0.1;
|
||||
rho_g,.9502,.1,.9999,BETA_PDF,0.85,0.1;
|
||||
rho_l,.9542,.1,.9999,BETA_PDF,0.85,0.1;
|
||||
rho_i,.6705,.1,.99,BETA_PDF,0.85,0.1;
|
||||
phi_i,5.2083,1,15,NORMAL_PDF,4,1.5;
|
||||
sig_c,0.9817,0.25,3,NORMAL_PDF,1,0.375;
|
||||
hab,0.5612,0.3,0.95,BETA_PDF,0.7,0.1;
|
||||
xi_w,0.7661,0.3,0.9,BETA_PDF,0.75,0.05;
|
||||
sig_l,1.7526,0.5,5,NORMAL_PDF,2,0.75;
|
||||
xi_p,0.8684,0.3,0.95,BETA_PDF,0.75,0.05;
|
||||
xi_e,0.5724,0.1,0.95,BETA_PDF,0.5,0.15;
|
||||
gamma_w,0.6202,0.1,0.99,BETA_PDF,0.75,0.15;
|
||||
gamma_p,0.6638,0.1,0.99,BETA_PDF,0.75,0.15;
|
||||
czcap,0.2516,0.01,2,NORMAL_PDF,0.2,0.075;
|
||||
phi_y,1.3011,1.001,2,NORMAL_PDF,1.45,0.125;
|
||||
r_pie,1.4616,1.2,2,NORMAL_PDF,1.7,0.1;
|
||||
r_dpi,0.1144,0.01,0.5,NORMAL_PDF,0.3,0.1;
|
||||
rho,0.8865,0.5,0.99,BETA_PDF,0.8,0.10;
|
||||
r_y,0.0571,0.01,0.2,NORMAL_PDF,0.125,0.05;
|
||||
r_dy,0.2228,0.05,0.5,NORMAL_PDF,0.0625,0.05;
|
||||
end;
|
||||
|
||||
varobs Y C I E PIE W R;
|
||||
|
||||
//estimation(datafile=rawdata_euromodel_1,presample=40, first_obs=1, nobs=118, lik_init=2, mode_compute=1,mh_replic=0);
|
||||
estimation(datafile=rawdata_euromodel_1,presample=40, first_obs=1, nobs=118,mh_jscale=0.2,mh_replic=150000, mode_check); //
|
||||
|
||||
|
||||
//stoch_simul(periods=200,irf=20,simul_seed=3) Y C PIE R W R_K L Q I K ;
|
||||
|
||||
Loading…
Reference in New Issue