From fddee8e1db602a1dd51a43d1ffa965eeeb9b8528 Mon Sep 17 00:00:00 2001
From: Johannes Pfeifer <jpfeifer@gmx,de>
Date: Wed, 30 Oct 2013 14:26:30 +0100
Subject: [PATCH] Bugfixes for correlated shocks
Uses preprocessing capabilities introduced in 07137e804b34a0e615398cb7cb5876f4a9027274
Fixes #392 and #494. Also fixes a bug in the checking for positive definiteness of covariance matrices in likelihood functions
Allows for calibrated covariances by reading them out and setting them after covariance matrix has been reconstructed from correlation and variances.
Adds unit test
---
matlab/dsge_likelihood.m | 24 +++++-
matlab/dynare_estimation_1.m | 34 ++------
matlab/non_linear_dsge_likelihood.m | 26 +++++-
matlab/print_info.m | 4 +
matlab/set_all_parameters.m | 45 ++++++----
matlab/set_parameters.m | 20 +++--
matlab/set_prior.m | 4 +-
.../fs2000_calibrated_covariance.mod | 86 +++++++++++++++++++
8 files changed, 187 insertions(+), 56 deletions(-)
create mode 100644 tests/estimation/fs2000_calibrated_covariance.mod
diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m
index 62bf68c1fa..8920413aa6 100644
--- a/matlab/dsge_likelihood.m
+++ b/matlab/dsge_likelihood.m
@@ -206,7 +206,7 @@ Q = Model.Sigma_e;
H = Model.H;
% Test if Q is positive definite.
-if ~issquare(Q) && EstimatedParameters.ncx
+if ~issquare(Q) || EstimatedParameters.ncx || isfield(EstimatedParameters,'calibrated_covariances')
[Q_is_positive_definite, penalty] = ispd(Q);
if ~Q_is_positive_definite
fval = objective_function_penalty_base+penalty;
@@ -214,10 +214,20 @@ if ~issquare(Q) && EstimatedParameters.ncx
info = 43;
return
end
+ if isfield(EstimatedParameters,'calibrated_covariances')
+ correct_flag=check_consistency_covariances(Q);
+ if ~correct_flag
+ penalty = sum(Q(EstimatedParameters.calibrated_covariances.position).^2);
+ fval = objective_function_penalty_base+penalty;
+ exit_flag = 0;
+ info = 71;
+ return
+ end
+ end
end
% Test if H is positive definite.
-if ~issquare(H) && EstimatedParameters.ncn
+if ~issquare(H) || EstimatedParameters.ncn || isfield(EstimatedParameters,'calibrated_covariances_ME')
[H_is_positive_definite, penalty] = ispd(H);
if ~H_is_positive_definite
fval = objective_function_penalty_base+penalty;
@@ -225,6 +235,16 @@ if ~issquare(H) && EstimatedParameters.ncn
info = 44;
return
end
+ if isfield(EstimatedParameters,'calibrated_covariances_ME')
+ correct_flag=check_consistency_covariances(H);
+ if ~correct_flag
+ penalty = sum(H(EstimatedParameters.calibrated_covariances_ME.position).^2);
+ fval = objective_function_penalty_base+penalty;
+ exit_flag = 0;
+ info = 72;
+ return
+ end
+ end
end
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index 2c8783bfd2..da865af160 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -80,38 +80,16 @@ end
[dataset_,xparam1, hh, M_, options_, oo_, estim_params_,bayestopt_] = dynare_estimation_init(var_list_, dname, [], M_, options_, oo_, estim_params_, bayestopt_);
+%check for calibrated covariances before updating parameters
+if ~isempty(estim_params_)
+ estim_params_=check_for_calibrated_covariances(xparam1,estim_params_,M_);
+end
% Set sigma_e_is_diagonal flag (needed if the shocks block is not declared in the mod file).
M_.sigma_e_is_diagonal = 1;
-if estim_params_.ncx || any(nnz(tril(M_.Sigma_e,-1)))
+if estim_params_.ncx || any(nnz(tril(M_.Correlation_matrix,-1))) || isfield(estim_params_,'calibrated_covariances')
M_.sigma_e_is_diagonal = 0;
end
-% Set the correlation matrix if necessary.
-if ~isequal(estim_params_.ncx,nnz(tril(M_.Sigma_e,-1)))
- M_.Correlation_matrix = diag(1./sqrt(diag(M_.Sigma_e)))*M_.Sigma_e*diag(1./sqrt(diag(M_.Sigma_e)));
- % Remove NaNs appearing because of variances calibrated to zero.
- if any(isnan(M_.Correlation_matrix))
- zero_variance_idx = find(~diag(M_.Sigma_e));
- for i=1:length(zero_variance_idx)
- M_.Correlation_matrix(zero_variance_idx(i),:) = 0;
- M_.Correlation_matrix(:,zero_variance_idx(i)) = 0;
- end
- end
-end
-
-% Set the correlation matrix of measurement errors if necessary.
-if ~isequal(estim_params_.ncn,nnz(tril(M_.H,-1)))
- M_.Correlation_matrix_ME = diag(1./sqrt(diag(M_.H)))*M_.H*diag(1./sqrt(diag(M_.H)));
- % Remove NaNs appearing because of variances calibrated to zero.
- if any(isnan(M_.Correlation_matrix_ME))
- zero_variance_idx = find(~diag(M_.H));
- for i=1:length(zero_variance_idx)
- M_.Correlation_matrix_ME(zero_variance_idx(i),:) = 0;
- M_.Correlation_matrix_ME(:,zero_variance_idx(i)) = 0;
- end
- end
-end
-
data = dataset_.data;
rawdata = dataset_.rawdata;
data_index = dataset_.missing.aindex;
@@ -151,7 +129,7 @@ if exist(mh_scale_fname)
end
if ~isempty(estim_params_)
- set_parameters(xparam1);
+ M_ = set_all_parameters(xparam1,estim_params_,M_);
end
% compute sample moments if needed (bvar-dsge)
diff --git a/matlab/non_linear_dsge_likelihood.m b/matlab/non_linear_dsge_likelihood.m
index 41ba13fdad..474b8dbb1c 100644
--- a/matlab/non_linear_dsge_likelihood.m
+++ b/matlab/non_linear_dsge_likelihood.m
@@ -163,7 +163,7 @@ Model = set_all_parameters(xparam1,EstimatedParameters,Model);
Q = Model.Sigma_e;
H = Model.H;
-if ~isscalar(Q) && EstimatedParameters.ncx
+if ~issquare(Q) || EstimatedParameters.ncx || isfield(EstimatedParameters,'calibrated_covariances')
[Q_is_positive_definite, penalty] = ispd(Q);
if ~Q_is_positive_definite
fval = objective_function_penalty_base+penalty;
@@ -171,9 +171,20 @@ if ~isscalar(Q) && EstimatedParameters.ncx
info = 43;
return
end
+ if isfield(EstimatedParameters,'calibrated_covariances')
+ correct_flag=check_consistency_covariances(Q);
+ if ~correct_flag
+ penalty = sum(Q(EstimatedParameters.calibrated_covariances.position).^2);
+ fval = objective_function_penalty_base+penalty;
+ exit_flag = 0;
+ info = 71;
+ return
+ end
+ end
+
end
-if ~isscalar(H) && EstimatedParameters.ncn
+if ~issquare(H) || EstimatedParameters.ncn || isfield(EstimatedParameters,'calibrated_covariances_ME')
[H_is_positive_definite, penalty] = ispd(H);
if ~H_is_positive_definite
fval = objective_function_penalty_base+penalty;
@@ -181,6 +192,17 @@ if ~isscalar(H) && EstimatedParameters.ncn
info = 44;
return
end
+ if isfield(EstimatedParameters,'calibrated_covariances_ME')
+ correct_flag=check_consistency_covariances(H);
+ if ~correct_flag
+ penalty = sum(H(EstimatedParameters.calibrated_covariances_ME.position).^2);
+ fval = objective_function_penalty_base+penalty;
+ exit_flag = 0;
+ info = 72;
+ return
+ end
+ end
+
end
%------------------------------------------------------------------------------
diff --git a/matlab/print_info.m b/matlab/print_info.m
index 157275fb4b..691bda9300 100644
--- a/matlab/print_info.m
+++ b/matlab/print_info.m
@@ -113,6 +113,10 @@ if ~noprint
error(['Discretionary policy: some eigenvalues greater than options_.qz_criterium. Model potentially unstable.']);
case 63
error(['Discretionary policy: NaN elements are present in the solution. Procedure failed.']);
+ case 71
+ error(['Calibrated covariance of the structural errors implies correlation larger than +-1.']);
+ case 72
+ error(['Calibrated covariance of the measurement errors implies correlation larger than +-1.']);
% Aim Code Conversions by convertAimCodeToInfo.m
case 102
error('Aim: roots not correctly computed by real_schur');
diff --git a/matlab/set_all_parameters.m b/matlab/set_all_parameters.m
index 656ac34edf..9fb3ef4bd6 100644
--- a/matlab/set_all_parameters.m
+++ b/matlab/set_all_parameters.m
@@ -56,9 +56,11 @@ nvn = estim_params.nvn;
ncn = estim_params.ncn;
np = estim_params.np;
Sigma_e = M.Sigma_e;
+Correlation_matrix = M.Correlation_matrix;
H = M.H;
-
-% setting shocks variance
+Correlation_matrix_ME = M.Correlation_matrix_ME;
+% setting shocks variance on the diagonal of Covariance matrix; used later
+% for updating covariances
if nvx
var_exo = estim_params.var_exo;
for i=1:nvx
@@ -69,7 +71,8 @@ end
% update offset
offset = nvx;
-% setting measument error variance
+% setting measument error variance; on the diagonal of Covariance matrix; used later
+% for updating covariances
if nvn
for i=1:nvn
k = estim_params.nvn_observable_correspondence(i,1);
@@ -81,38 +84,42 @@ end
offset = nvx+nvn;
% setting shocks covariances
-if ~isempty(M.Correlation_matrix)
- Sigma_e = diag(sqrt(diag(Sigma_e)))*M.Correlation_matrix*diag(sqrt(diag(Sigma_e))); % use of old correlation matrix is correct due to the diagonal structure and later only using the hence correctly updated diagonal entries of Sigma_e
-end
if ncx
corrx = estim_params.corrx;
for i=1:ncx
k1 = corrx(i,1);
k2 = corrx(i,2);
- M.Correlation_matrix(k1,k2) = xparam1(i+offset);
- M.Correlation_matrix(k2,k1) = M.Correlation_matrix(k1,k2);
- Sigma_e(k1,k2) = xparam1(i+offset)*sqrt(Sigma_e(k1,k1)*Sigma_e(k2,k2));
- Sigma_e(k2,k1) = Sigma_e(k1,k2);
+ Correlation_matrix(k1,k2) = xparam1(i+offset);
+ Correlation_matrix(k2,k1) = Correlation_matrix(k1,k2);
end
end
-
+%build covariance matrix from correlation matrix and variances already on
+%diagonal
+Sigma_e = diag(sqrt(diag(Sigma_e)))*Correlation_matrix*diag(sqrt(diag(Sigma_e)));
+%if calibrated covariances, set them now to their stored value
+if isfield(estim_params,'calibrated_covariances')
+ Sigma_e(estim_params.calibrated_covariances.position)=estim_params.calibrated_covariances.cov_value;
+end
% update offset
offset = nvx+nvn+ncx;
+
% setting measurement error covariances
-if ~isempty(M.Correlation_matrix_ME)
- H = diag(sqrt(diag(H)))*M.Correlation_matrix_ME*diag(sqrt(diag(H)));
-end
if ncn
corrn_observable_correspondence = estim_params.corrn_observable_correspondence;
for i=1:ncn
k1 = corrn_observable_correspondence(i,1);
k2 = corrn_observable_correspondence(i,2);
- M.Correlation_matrix_ME(k1,k2) = xparam1(i+offset);
- M.Correlation_matrix_ME(k2,k1) = M.Correlation_matrix_ME(k1,k2);
- H(k1,k2) = xparam1(i+offset)*sqrt(H(k1,k1)*H(k2,k2));
- H(k2,k1) = H(k1,k2);
+ Correlation_matrix_ME(k1,k2) = xparam1(i+offset);
+ Correlation_matrix_ME(k2,k1) = Correlation_matrix_ME(k1,k2);
end
end
+%build covariance matrix from correlation matrix and variances already on
+%diagonal
+H = diag(sqrt(diag(H)))*Correlation_matrix_ME*diag(sqrt(diag(H)));
+%if calibrated covariances, set them now to their stored value
+if isfield(estim_params,'calibrated_covariances_ME')
+ H(estim_params.calibrated_covariances_ME.position)=estim_params.calibrated_covariances_ME.cov_value;
+end
% update offset
offset = nvx+ncx+nvn+ncn;
@@ -125,7 +132,9 @@ end
% updating matrices in M
if nvx || ncx
M.Sigma_e = Sigma_e;
+ M.Correlation_matrix=Correlation_matrix;
end
if nvn || ncn
M.H = H;
+ M.Correlation_matrix_ME=Correlation_matrix_ME;
end
\ No newline at end of file
diff --git a/matlab/set_parameters.m b/matlab/set_parameters.m
index d1e81afd81..58e167d554 100644
--- a/matlab/set_parameters.m
+++ b/matlab/set_parameters.m
@@ -14,7 +14,7 @@ function set_parameters(xparam1)
% SPECIAL REQUIREMENTS
% none
-% Copyright (C) 2003-2009 Dynare Team
+% Copyright (C) 2003-2013 Dynare Team
%
% This file is part of Dynare.
%
@@ -39,9 +39,11 @@ nvn = estim_params_.nvn;
ncn = estim_params_.ncn;
np = estim_params_.np;
Sigma_e = M_.Sigma_e;
+Correlation_matrix = M_.Correlation_matrix;
offset = 0;
-% stderrs of the exogenous shocks
+% setting shocks variance on the diagonal of Covariance matrix; used later
+% for updating covariances
if nvx
var_exo = estim_params_.var_exo;
for i=1:nvx
@@ -58,10 +60,17 @@ if ncx
for i=1:ncx
k1 = corrx(i,1);
k2 = corrx(i,2);
- Sigma_e(k1,k2) = xparam1(i+offset)*sqrt(Sigma_e(k1,k1)*Sigma_e(k2,k2));
- Sigma_e(k2,k1) = Sigma_e(k1,k2);
+ Correlation_matrix(k1,k2) = xparam1(i+offset);
+ Correlation_matrix(k2,k1) = Correlation_matrix(k1,k2);
end
end
+%build covariance matrix from correlation matrix and variances already on
+%diagonal
+Sigma_e = diag(sqrt(diag(Sigma_e)))*Correlation_matrix*diag(sqrt(diag(Sigma_e)));
+if isfield(estim_params,'calibrated_covariances')
+ Sigma_e(estim_params.calibrated_covariances.position)=estim_params.calibrated_covariances.cov_value;
+end
+
% and update offset
offset = offset + ncx + ncn;
@@ -70,4 +79,5 @@ if np
M_.params(estim_params_.param_vals(:,1)) = xparam1(offset+1:end);
end
-M_.Sigma_e = Sigma_e;
\ No newline at end of file
+M_.Sigma_e = Sigma_e;
+M_.Correlation_matrix=Correlation_matrix;
diff --git a/matlab/set_prior.m b/matlab/set_prior.m
index f82223ec85..9f2c40285b 100644
--- a/matlab/set_prior.m
+++ b/matlab/set_prior.m
@@ -75,9 +75,10 @@ if nvx
end
if nvn
estim_params_.nvn_observable_correspondence=NaN(nvn,1); % stores number of corresponding observable
- if isequal(M_.H,0)
+ if isequal(M_.H,0) %if no previously set measurement error, initialize H
nvarobs = size(options_.varobs,1);
M_.H = zeros(nvarobs,nvarobs);
+ M_.Correlation_matrix_ME = eye(nvarobs);
end
for i=1:nvn
obsi_ = strmatch(deblank(M_.endo_names(estim_params_.var_endo(i,1),:)),deblank(options_.varobs),'exact');
@@ -116,6 +117,7 @@ if ncn
if isequal(M_.H,0)
nvarobs = size(options_.varobs,1);
M_.H = zeros(nvarobs,nvarobs);
+ M_.Correlation_matrix_ME = eye(nvarobs);
end
xparam1 = [xparam1; estim_params_.corrn(:,3)];
ub = [ub; max(min(estim_params_.corrn(:,5),1),-1)];
diff --git a/tests/estimation/fs2000_calibrated_covariance.mod b/tests/estimation/fs2000_calibrated_covariance.mod
new file mode 100644
index 0000000000..20d34c84ad
--- /dev/null
+++ b/tests/estimation/fs2000_calibrated_covariance.mod
@@ -0,0 +1,86 @@
+//two covariances are calibrated and one covariance of the ME. One of the calibrated covariances is superseded by estimation
+
+var m P c e W R k d n l gy_obs gp_obs y dA;
+varexo e_a e_m e_f;
+
+parameters alp bet gam mst rho psi del;
+
+alp = 0.33;
+bet = 0.99;
+gam = 0.003;
+mst = 1.011;
+rho = 0.7;
+psi = 0.787;
+del = 0.02;
+
+model;
+dA = exp(gam+e_a);
+log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
+-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
+W = l/n;
+-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
+R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
+1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
+c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
+P*c = m;
+m-1+d = l;
+e = exp(e_a);
+y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
+gy_obs = dA*y/y(-1);
+gp_obs = (P/P(-1))*m(-1)/dA;
+end;
+
+initval;
+k = 6;
+m = mst;
+P = 2.25;
+c = 0.45;
+e = 1;
+W = 4;
+R = 1.02;
+d = 0.85;
+n = 0.19;
+l = 0.86;
+y = 0.6;
+gy_obs = exp(gam);
+gp_obs = exp(-gam);
+dA = exp(gam);
+end;
+
+varobs gp_obs gy_obs;
+
+shocks;
+var e_a; stderr 0.01;
+var e_f; stderr 0.01;
+var gy_obs; stderr 0.01;
+var gy_obs, gp_obs= 0.00005;
+var e_a, e_m = 0.00005;
+var e_m, e_f = 0.00001;
+end;
+
+steady;
+check;
+
+estimated_params;
+stderr e_m, 0.008862;
+corr e_a, e_m, 0.5;
+stderr gp_obs, 0.035449;
+end;
+
+estimated_params_init;
+stderr e_m, 0.5;
+corr e_a, e_m, 0.5;
+stderr gp_obs, 0.5;
+end;
+
+estimation(order=1,datafile=fsdat_simul,nobs=192, loglinear, mh_replic=2002, mh_nblocks=1, mh_jscale=0.8);
+
+if ~isequal(M_.Sigma_e(3,2),1e-5) || ~isequal(M_.Sigma_e(2,3),1e-5)
+ error('Problem in setting calibrated covariance of structural shocks')
+end
+if isequal(M_.Sigma_e(2,1),5e-5) || isequal(M_.Sigma_e(1,2),5e-5)
+ error('Problem in overriding calibrated covariance of structural shocks by estimated correlation')
+end
+if ~isequal(M_.H(2,1),5e-5) || ~isequal(M_.H(1,2),5e-5)
+ error('Problem in setting calibrated covariance of measurement errors')
+end
\ No newline at end of file
--
GitLab