Merge pull request #511 from JohannesPfeifer/Correlated_errors_preprocessor

Bugfixes for correlated shocks

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
 
 

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)

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
 
 %------------------------------------------------------------------------------

matlab/print_info.m

@@ -114,6 +114,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');

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

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;

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)];

tests/estimation/fs2000_calibrated_covariance.mod

+//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