Merge branch 'master' into remove-dynDate-class

README.md

@@ -257,7 +257,7 @@ Configure and make:
     - Next to `Command Line Tools`, click on `Install`
 - Download [MacOSX10.6.sdk.zip](http://www.jamesgeorge.org/uploads/MacOSX10.6.sdk.zip) and unzip it in `/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs`. Change the owner to be `root` and the group to be `wheel`
 
-- Install [Homebrew](http://mxcl.github.io/homebrew/)
+- Install [Homebrew](http://mxcl.github.io/homebrew/) and [Homebrew Science](https://github.com/Homebrew/homebrew-science)
 - Install the following brews:
 ```
     brew install automake
@@ -269,7 +269,12 @@ Configure and make:
     brew install slicot --with-default-integer-8
 ```
 - **(Optional)** To compile Dynare mex files for use on Octave, first install Octave following the [Simple Installation Instructions](http://wiki.octave.org/Octave_for_MacOS_X#Simple_Installation_Instructions_3). Then, you will probably also want to install graphicsmagick via Homebrew with `brew install graphicsmagick`.
-- **(Optional)** To compile Dynare's documentation, first install the latest version of [MacTeX](http://www.tug.org/mactex/). Then install `doxygen` and `latex2html` via Homebrew with the following commands:
+- **(Optional)** To compile Dynare's documentation, first install the latest version of [MacTeX](http://www.tug.org/mactex/). Then install `doxygen`, `latex2html` and `texi2html` via Homebrew with the following commands:
+```
+    brew install doxygen
+    brew install latex2html
+    brew install texi2html
+```
 - **(On OS X 10.7 Only)** Copy [FlexLexer.h](http://www.dynare.org/DynareWiki/BuildingDynareFromSource?action=AttachFile&do=view&target=FlexLexer.h) into the `preprocessor` directory (there was an error in the `FlexLexer.h` file distributed with 10.7)
 - Finally, switch to the root dynare directory. Ensure your path contains `/usr/bin:/bin:/usr/sbin:/sbin:/usr/local/bin:/opt/X11/bin:/usr/texbin:/usr/local/sbin`. Run:
     - `autoconf -si`

doc/dynare.texi

@@ -752,7 +752,7 @@ Activate console mode. In addition to the behavior of
 computations.
 
 @item nograph
-Activate the @code{nograph} option (@xref{nograph}), so that Dynare will not produce any
+Activate the @code{nograph} option (@pxref{nograph}), so that Dynare will not produce any
 graph
 
 @item nointeractive
@@ -4400,6 +4400,7 @@ the total number of Metropolis draws available. Default:
 The fraction of initially generated parameter vectors to be dropped as a burnin before using posterior simulations. Default: @code{0.5}
 
 @item mh_jscale = @var{DOUBLE}
+@anchor{mh_jscale}
 The scale to be used for the jumping distribution in
 Metropolis-Hastings algorithm. The default value is rarely
 satisfactory. This option must be tuned to obtain, ideally, an
@@ -4496,6 +4497,26 @@ value of that function as the posterior mode.
 @noindent
 Default value is @code{4}.
 
+@item mcmc_jumping_covariance = hessian|prior_variance|identity_matrix|@var{FILENAME}
+Tells Dynare which covariance to use for the proposal density of the MCMC sampler. @code{mcmc_jumping_covariance} can be one of the following:
+
+@table @code
+@item hessian 
+Uses the Hessian matrix computed at the mode.
+
+@item prior_variance
+Uses the prior variances. No infinite prior variances are allowed in this case.
+
+@item identity_matrix
+Uses an identity matrix.
+
+@item @var{FILENAME}
+Loads an arbitrary user-specified covariance matrix from @code{@var{FILENAME}.mat}. The covariance matrix must be saved in a variable named @code{jumping_covariance}, must be square, positive definite, and have the same dimension as the number of estimated parameters.
+
+@end table
+@noindent
+Note that the covariance matrices are still scaled with @ref{mh_jscale}. Default value is @code{hessian}.
+
 @item mode_check
 Tells Dynare to plot the posterior density for values around the
 computed mode for each estimated parameter in turn. This is helpful to

matlab/@dseries/merge.m

@@ -51,7 +51,11 @@ tex = [B.tex; C.tex];
 A.tex = tex(IBC); 
 A.vobs=length(IBC);
 
-if B.init >= C.init
+if B.nobs == 0
+    A = C;
+elseif C.nobs == 0
+    A = B;
+elseif B.init >= C.init
     diff = B.init - C.init;
     A.nobs = max(B.nobs + diff, C.nobs);
     A.data = NaN(A.nobs, A.vobs);

matlab/check_consistency_covariances.m

+function correct_flag=check_consistency_covariances(Covariance_matrix)
+% function check_consistency_covariances(Covariance_matrix)
+% checks consistency of covariance matrices by checking whether the
+% covariances imply correlations bigger than 1.
+%
+% Outputs: correct_flag         [scalar] 0 if not consistent, 1 otherwise
+% Inputs: Covariance_matrix     [matrix] covariance matrix to be checked
+
+% Copyright (C) 2013 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/>.
+
+%compute theoretical bound by assuming correlation of 1
+bound=diag(sqrt(diag(Covariance_matrix)))*ones(size(Covariance_matrix))*diag(sqrt(diag(Covariance_matrix)));
+correct_flag=1;
+if any(any(tril(Covariance_matrix,-1)>bound)) || any(any(tril(Covariance_matrix,-1)<-bound))
+    correct_flag=0;
+end
\ No newline at end of file

matlab/check_for_calibrated_covariances.m

+function estim_params=check_for_calibrated_covariances(xparam1,estim_params,M)
+% function check_for_calibrated_covariances(xparam1,estim_params,M)
+% find calibrated covariances to consider during estimation
+% Inputs
+%   -xparam1        [vector] parameters to be estimated
+%   -estim_params   [structure] describing parameters to be estimated
+%   -M              [structure] describing the model
+%
+% Outputs
+%   -estim_params   [structure] describing parameters to be estimated
+%
+% Notes: M is local to this function and not updated when calling
+% set_all_parameters
+
+% Copyright (C) 2013 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/>.
+Sigma_e_calibrated=M.Sigma_e;
+H_calibrated=M.H;
+%check covariance for structural errors
+covariance_pos=find(tril(Sigma_e_calibrated,-1)); %off-diagonal elements set by covariances before updating correlation matrix to reflect estimated covariances
+covariance_pos_ME=find(tril(H_calibrated,-1)); %off-diagonal elements set by covariances before updating correlation matrix to reflect estimated covariances
+
+%locally updated M
+M = set_all_parameters(xparam1,estim_params,M);
+
+correlation_pos=find(tril(M.Correlation_matrix,-1)); %off-diagonal elements set by correlations after accounting for estimation
+calibrated_covariance_pos=covariance_pos(~ismember(covariance_pos,correlation_pos));
+if any(calibrated_covariance_pos)
+    [rows, columns]=ind2sub(size(M.Sigma_e),calibrated_covariance_pos); %find linear indices of lower triangular covariance entries
+    estim_params.calibrated_covariances.position=[calibrated_covariance_pos;sub2ind(size(M.Sigma_e),columns,rows)]; %get linear entries of upper triangular parts
+    estim_params.calibrated_covariances.cov_value=Sigma_e_calibrated(estim_params.calibrated_covariances.position);  
+end
+
+correlation_pos_ME=find(tril(M.Correlation_matrix_ME,-1)); %off-diagonal elements set by correlations after accounting for estimation
+calibrated_covariance_pos_ME=covariance_pos_ME(~ismember(covariance_pos_ME,correlation_pos_ME));
+if any(calibrated_covariance_pos_ME)
+    [rows, columns]=ind2sub(size(M.H),calibrated_covariance_pos_ME); %find linear indices of lower triangular covariance entries
+    estim_params.calibrated_covariances_ME.position=[calibrated_covariance_pos_ME;sub2ind(size(M.H),columns,rows)]; %get linear entries of upper triangular parts
+    estim_params.calibrated_covariances_ME.cov_value=H_calibrated(estim_params.calibrated_covariances_ME.position);  
+end
+

matlab/check_model.m

@@ -30,4 +30,13 @@ end
 if (DynareModel.exo_det_nbr > 0) && (DynareModel.maximum_lag > 1 || DynareModel.maximum_lead > 1)
     error(['Exogenous deterministic variables are currently only allowed in' ...
            ' models with leads and lags on only one period'])
+end
+
+if ~check_consistency_covariances(DynareModel.Sigma_e)
+    error('The specified covariances for the structural errors are not consistent with the variances as they imply a correlation larger than +-1')
+end
+if ~isequal(DynareModel.H,0)
+    if ~check_consistency_covariances(DynareModel.H)
+        error('The specified covariances for the measurement errors are not consistent with the variances as they imply a correlation larger than +-1')
+    end
 end
\ No newline at end of file

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.m

@@ -53,7 +53,8 @@ end
 
 M_.dname = dname;
 
-if options_.mode_compute && options_.analytic_derivation,
+if (isnumeric(options_.mode_compute) && options_.mode_compute && options_.analytic_derivation) ... %no user supplied function
+        || (~isnumeric(options_.mode_compute) && options_.analytic_derivation) % user supplied function
     analytic_derivation0=options_.analytic_derivation;
     options_.analytic_derivation=1;
 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)
@@ -602,8 +580,6 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
 
         [xparam1, fval, nacc, nfcnev, nobds, ier, t, vm] = sa(objective_function,xparam1,maxy,rt_,epsilon,ns,nt ...
                                                               ,neps,maxevl,LB,UB,c,idisp ,t,vm,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
-      case 'prior'
-        hh = diag(bayestopt_.p2.^2);
       otherwise
         if ischar(options_.mode_compute)
             [xparam1, fval, retcode ] = feval(options_.mode_compute,objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
@@ -611,8 +587,8 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
             error(['dynare_estimation:: mode_compute = ' int2str(options_.mode_compute) ' option is unknown!'])
         end
     end
-    if ~isequal(options_.mode_compute,6) && ~isequal(options_.mode_compute,'prior')
-        if options_.cova_compute == 1
+    if ~isequal(options_.mode_compute,6) %always already computes covariance matrix
+        if options_.cova_compute == 1 %user did not request covariance not to be computed
             if options_.analytic_derivation && strcmp(func2str(objective_function),'dsge_likelihood'),
                 ana_deriv = options_.analytic_derivation;
                 options_.analytic_derivation = 2;
@@ -638,6 +614,35 @@ if options_.cova_compute == 0
     hh = [];%NaN(length(xparam1),length(xparam1));
 end
 
+switch options_.MCMC_jumping_covariance
+    case 'hessian' %Baseline
+        %do nothing and use hessian from mode_compute
+    case 'prior_variance' %Use prior variance
+        if any(isinf(bayestopt_.p2))
+            error('Infinite prior variances detected. You cannot use the prior variances as the proposal density, if some variances are Inf.')
+        else            
+            hh = diag(1./(bayestopt_.p2.^2));
+        end
+    case 'identity_matrix' %Use identity
+        hh = eye(nx);   
+    otherwise %user specified matrix in file
+        try 
+            load(options_.MCMC_jumping_covariance,'jumping_covariance')
+            hh=jumping_covariance;
+        catch
+            error(['No matrix named ''jumping_covariance'' could be found in ',options_.MCMC_jumping_covariance,'.mat'])
+        end
+        [nrow, ncol]=size(hh);
+        if ~isequal(nrow,ncol) && ~isequal(nrow,nx) %check if square and right size
+            error(['jumping_covariance matrix must be square and have ',num2str(nx),' rows and columns'])
+        end
+        try %check for positive definiteness
+            chol(hh);
+        catch
+            error(['Specified jumping_covariance is not positive definite'])
+        end
+end
+
 if ~options_.mh_posterior_mode_estimation && options_.cova_compute
     try
         chol(hh);
@@ -646,7 +651,7 @@ if ~options_.mh_posterior_mode_estimation && options_.cova_compute
         disp('POSTERIOR KERNEL OPTIMIZATION PROBLEM!')
         disp(' (minus) the hessian matrix at the "mode" is not positive definite!')
         disp('=> posterior variance of the estimated parameters are not positive.')
-        disp('You should  try  to change the initial values of the parameters using')
+        disp('You should try to change the initial values of the parameters using')
         disp('the estimated_params_init block, or use another optimization routine.')
         params_at_bound=find(xparam1==ub | xparam1==lb);
         if ~isempty(params_at_bound)
@@ -767,11 +772,15 @@ if (any(bayestopt_.pshape  >0 ) && options_.mh_replic) || ...
         CutSample(M_, options_, estim_params_);
         %% Estimation of the marginal density from the Mh draws:
         if options_.mh_replic
-            [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_);
+            [marginal,oo_] = marginal_density(M_, options_, estim_params_, ...
+                                              oo_);
+            % Store posterior statistics by parameter name
             oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bayestopt_, oo_);
             if ~options_.nograph
                 oo_ = PlotPosteriorDistributions(estim_params_, M_, options_, bayestopt_, oo_);
             end
+            % Store posterior mean in a vector and posterior variance in
+            % a matrix
             [oo_.posterior.metropolis.mean,oo_.posterior.metropolis.Variance] ...
                 = GetPosteriorMeanVariance(M_,options_.mh_drop);
         else

matlab/global_initialization.m

@@ -381,6 +381,7 @@ options_.mh_nblck = 2;
 options_.mh_recover = 0;
 options_.mh_replic = 20000;
 options_.recursive_estimation_restart = 0;
+options_.MCMC_jumping_covariance='hessian';
 
 options_.mode_compute = 4;
 options_.mode_file = '';

matlab/model_diagnostics.m

@@ -85,6 +85,7 @@ end
 %
 % singular Jacobian of static model
 %
+singularity_problem = 0;
 if ~isfield(M,'block_structure_stat')
     nb = 1;
 else
@@ -107,6 +108,7 @@ for b=1:nb
     end
     rank_jacob = rank(jacob);
     if rank_jacob < size(jacob,1)
+        singularity_problem = 1;
         disp(['model_diagnostic: the Jacobian of the static model is ' ...
               'singular'])
         disp(['there is ' num2str(endo_nbr-rank_jacob) ...
@@ -143,4 +145,9 @@ for b=1:nb
         end
     end
 end
+if singularity_problem 
+    fprint('The presence of a singularity problem typically indicates that there is one\n')
+    fprint('redundant equation entered in the model block, while another non-redundant equation\n')
+    fprint('is missing. The problem often derives from Walras Law.\n')
+end
 

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

@@ -36,7 +36,8 @@ if ~noprint
       case 2
         error(['The generalized Schur (QZ) decomposition failed. ' ...
                'For more information, see the documentation for Lapack function dgges: info=' ...
-               int2str(info(2)) ', n=' int2str(info(3))])
+               int2str(info(2)) ', n=' int2str(info(3)) ...
+               '. You can also run model_diagnostics to get more information on what may cause this problem.'])
       case 3
         error(['Blanchard Kahn conditions are not satisfied: no stable' ...
                ' equilibrium'])
@@ -113,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