Commit fddee8e1 authored by Johannes Pfeifer's avatar Johannes Pfeifer
Browse files

Bugfixes for correlated shocks

Uses preprocessing capabilities introduced in 07137e80

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
parent 07137e80
......@@ -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
......
......@@ -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)
......
......@@ -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
%------------------------------------------------------------------------------
......
......@@ -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');
......
......@@ -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
......@@ -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;
......@@ -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)];
......
//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
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment