From 698a44c98a8030e6f42ce4b57a3befeb5d88dd84 Mon Sep 17 00:00:00 2001
From: Johannes Pfeifer <jpfeifer@gmx.de>
Date: Mon, 10 Aug 2015 21:36:48 +0200
Subject: [PATCH] Add option for storing contemporaneous correlation

---
 doc/dynare.texi                  | 43 ++++++++++++++++++++++++++------
 matlab/compute_moments_varendo.m |  5 ++++
 matlab/covariance_mc_analysis.m  | 28 +++++++++++++++++++--
 matlab/disp_moments.m            |  3 +++
 matlab/disp_th_moments.m         |  3 +++
 matlab/global_initialization.m   |  1 +
 matlab/posterior_analysis.m      |  2 +-
 matlab/prior_analysis.m          |  2 +-
 preprocessor/DynareBison.yy      |  5 +++-
 preprocessor/DynareFlex.ll       |  1 +
 tests/TeX/fs2000_corr_ME.mod     |  6 ++---
 11 files changed, 84 insertions(+), 15 deletions(-)

diff --git a/doc/dynare.texi b/doc/dynare.texi
index 6a553ff45e..08cc205d53 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -2189,7 +2189,7 @@ possible types of simulations in stochastic mode:
 @itemize
 
 @item
-@ref{stoch_simul}, if the @code{periods} options is specified
+@ref{stoch_simul}, if the @code{periods} option is specified
 
 @item
 @ref{forecast} as the initial point at which the forecasts are computed
@@ -2419,7 +2419,7 @@ Moreover, as only states enter the recursive policy functions, all values specif
 @itemize
 
 @item
-in @ref{stoch_simul}, if the @code{periods} options is specified. Note that this only affects the starting point for the simulation, but not for the impulse response functions.
+in @ref{stoch_simul}, if the @code{periods} option is specified. Note that this only affects the starting point for the simulation, but not for the impulse response functions.
 
 @item
 in @ref{forecast} as the initial point at which the forecasts are computed
@@ -3925,6 +3925,11 @@ Tolerance for the suppression of small terms in the display of decision rules. R
 smaller than @code{dr_display_tol} are not displayed. 
 Default value: @code{1e-6}.
 
+@item contemporaneous_correlation
+@anchor{contemporaneous_correlation} 
+Saves the contemporaneous correlation between the endogenous variables in @code{oo_.contemporaneous_correlation}.
+Requires the @code{nocorr}-option not to be set.
+
 @end table
 
 @outputhead
@@ -3942,6 +3947,8 @@ If options @code{irf} is different from zero, sets @code{oo_.irfs}
 the global workspace (this latter way of accessing the IRFs is
 deprecated and will disappear in a future version).
 
+If the option @code{contemporaneous_correlation} is different from 0, sets 
+{oo_.contemporaneous_correlation}, which is described below.
 
 @customhead{Example 1}
 
@@ -3971,14 +3978,15 @@ response functions on 60 periods for variables @code{y} and @code{k}.
 @defvr {MATLAB/Octave variable} oo_.mean
 After a run of @code{stoch_simul}, contains the mean of the endogenous
 variables. Contains theoretical mean if the @code{periods} option is
-not present, and empirical mean otherwise. The variables are arranged
+not present, and simulated mean otherwise. The variables are arranged
 in declaration order.
 @end defvr
 
 @defvr {MATLAB/Octave variable} oo_.var
 After a run of @code{stoch_simul}, contains the variance-covariance of
 the endogenous variables. Contains theoretical variance if the
-@code{periods} option is not present (or an approximation thereof for @code{order=2}), and empirical variance
+@code{periods} option is not present (or an approximation thereof for @code{order=2}), 
+and simulated variance
 otherwise. The variables are arranged in declaration order.
 @end defvr
 
@@ -3990,7 +3998,7 @@ number of the matrix in the cell array corresponds to the order of
 autocorrelation. The option @code{ar} specifies the number of
 autocorrelation matrices available. Contains theoretical
 autocorrelations if the @code{periods} option is not present (or an approximation thereof for @code{order=2}), and
-empirical autocorrelations otherwise. The field is only created if stationary variables are present.
+simulated autocorrelations otherwise. The field is only created if stationary variables are present.
 
 The element @code{oo_.autocorr@{i@}(k,l)} is equal to the correlation
 between @math{y^k_t} and @math{y^l_{t-i}}, where @math{y^k}
@@ -4024,13 +4032,18 @@ If a second order approximation has been requested, contains the
 vector of the mean correction terms.
 @end table
 
-In case of @code{order=2}, the theoretical second moments are a second order accurate approximation of the true second moments, see @code{conditional_variance_decomposition}.
+In case of @code{order=2}, the theoretical second moments are a second order 
+accurate approximation of the true second moments, see @code{conditional_variance_decomposition}.
 
 @end defvr
 
 @anchor{oo_.variance_decomposition}
 @defvr {MATLAB/Octave variable} oo_.variance_decomposition
-After a run of @code{stoch_simul} when requesting theoretical moments (@code{periods=0}), contains a matrix with the result of the unconditional variance decomposition (i.e. at horizon infinity). The first dimension corresponds to the endogenous variables (in the order of declaration) and the second dimension corresponds to exogenous variables (in the order of declaration). Numbers are in percent and sum up to 100 across columns.
+After a run of @code{stoch_simul} when requesting theoretical moments (@code{periods=0}), 
+contains a matrix with the result of the unconditional variance decomposition (i.e. at horizon infinity). 
+The first dimension corresponds to the endogenous variables (in the order of declaration) and 
+the second dimension corresponds to exogenous variables (in the order of declaration). 
+Numbers are in percent and sum up to 100 across columns.
 @end defvr
 
 @anchor{oo_.conditional_variance_decomposition}
@@ -4044,6 +4057,15 @@ the order of declaration), the third dimension corresponds to
 exogenous variables (in the order of declaration).
 @end defvr
 
+@anchor{oo_.contemporaneous_correlation}
+@defvr {MATLAB/Octave variable} oo_.contemporaneous_correlation
+After a run of @code{stoch_simul} with the
+@code{contemporaneous_correlation} option, contains theoretical contemporaneous correlations if the
+@code{periods} option is not present (or an approximation thereof for @code{order=2}), 
+and simulated contemporaneous correlations otherwise. The variables are arranged in declaration order.
+@end defvr
+
+
 @defvr {MATLAB/Octave variable} oo_.irfs
 After a run of @code{stoch_simul} with option @code{irf} different
 from zero, contains the impulse responses, with the following naming
@@ -5380,6 +5402,10 @@ variables. Results are stored in
 @code{oo_.PosteriorTheoreticalMoments} (@pxref{oo_.PosteriorTheoreticalMoments}). The number of lags in the autocorrelation function is
 controlled by the @code{ar} option.
 
+@item contemporaneous_correlation
+@xref{contemporaneous_correlation}. Results are stored in @code{oo_.PosteriorTheoreticalMoments}. 
+Note that the @code{nocorr}-option has no effect.
+
 @item conditional_variance_decomposition = @var{INTEGER}
 See below.
 
@@ -5924,6 +5950,9 @@ where @var{THEORETICAL_MOMENT} is one of the following:
 @item covariance
 Variance-covariance of endogenous variables
 
+@item contemporaneous_correlation
+Contemporaneous correlation of endogenous variables when the @ref{contemporaneous_correlation} option is specified.
+
 @item correlation
 Auto- and cross-correlation of endogenous variables. Fields are vectors with correlations from 1 up to order @code{options_.ar}
 
diff --git a/matlab/compute_moments_varendo.m b/matlab/compute_moments_varendo.m
index cf91019c9a..75952c78ad 100644
--- a/matlab/compute_moments_varendo.m
+++ b/matlab/compute_moments_varendo.m
@@ -33,6 +33,9 @@ function oo_ = compute_moments_varendo(type,options_,M_,oo_,var_list_)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.    
 
+
+fprintf('Estimation::compute_moments_varendo: I''m computing endogenous moments (this may take a while)... ');
+
 if strcmpi(type,'posterior')
     posterior = 1;
     if nargin==4
@@ -126,3 +129,5 @@ if M_.exo_nbr > 1
         end
     end
 end
+
+fprintf(' Done!\n');
diff --git a/matlab/covariance_mc_analysis.m b/matlab/covariance_mc_analysis.m
index d98992f0eb..59cc635f2b 100644
--- a/matlab/covariance_mc_analysis.m
+++ b/matlab/covariance_mc_analysis.m
@@ -1,4 +1,4 @@
-function oo_ = covariance_mc_analysis(NumberOfSimulations,type,dname,fname,vartan,nvar,var1,var2,mh_conf_sig,oo_)
+function oo_ = covariance_mc_analysis(NumberOfSimulations,type,dname,fname,vartan,nvar,var1,var2,mh_conf_sig,oo_,options_)
 % This function analyses the (posterior or prior) distribution of the
 % endogenous variables' covariance matrix.
 % 
@@ -14,6 +14,7 @@ function oo_ = covariance_mc_analysis(NumberOfSimulations,type,dname,fname,varta
 %   mh_conf_sig             [double]            2 by 1 vector with upper
 %                                               and lower bound of HPD intervals
 %   oo_                     [structure]         Dynare structure where the results are saved.
+%   options_                [structure]         Dynare options structure 
 %
 % OUTPUTS
 %   oo_                     [structure]        Dynare structure where the results are saved.
@@ -86,10 +87,21 @@ end
 
 ListOfFiles = dir([ PATH  fname '_' TYPE '2ndOrderMoments*.mat']);
 i1 = 1; tmp = zeros(NumberOfSimulations,1);
+if options_.contemporaneous_correlation
+    tmp_corr_mat = zeros(NumberOfSimulations,1);
+    cov_pos=symmetric_matrix_index(indx1,indx2,nvar);
+    var_pos_1=symmetric_matrix_index(indx1,indx1,nvar);
+    var_pos_2=symmetric_matrix_index(indx2,indx2,nvar);
+end
 for file = 1:length(ListOfFiles)
     load([ PATH ListOfFiles(file).name ]);
     i2 = i1 + rows(Covariance_matrix) - 1;
     tmp(i1:i2) = Covariance_matrix(:,symmetric_matrix_index(indx1,indx2,nvar));
+    if options_.contemporaneous_correlation
+        temp=Covariance_matrix(:,cov_pos)./(sqrt(Covariance_matrix(:,var_pos_1)).*sqrt(Covariance_matrix(:,var_pos_2)));
+        temp(Covariance_matrix(:,cov_pos)==0)=0; %filter out 0 correlations that would result in 0/0
+        tmp_corr_mat(i1:i2)=temp;
+    end    
     i1 = i2+1;
 end
 name = [var1 '.' var2];
@@ -111,4 +123,16 @@ else
     eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.HPDsup.' name ' = NaN;']);
     eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.deciles.' name ' = NaN;']);
     eval(['oo_.' TYPE 'TheoreticalMoments.dsge.covariance.density.' name ' = NaN;']);
-end
\ No newline at end of file
+end
+
+if options_.contemporaneous_correlation
+    [p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ...
+        posterior_moments(tmp_corr_mat,1,mh_conf_sig);
+    eval(['oo_.' TYPE 'TheoreticalMoments.dsge.contemporeaneous_correlation.Mean.' name ' = p_mean;']);
+    eval(['oo_.' TYPE 'TheoreticalMoments.dsge.contemporeaneous_correlation.Median.' name ' = p_median;']);
+    eval(['oo_.' TYPE 'TheoreticalMoments.dsge.contemporeaneous_correlation.Variance.' name ' = p_var;']);
+    eval(['oo_.' TYPE 'TheoreticalMoments.dsge.contemporeaneous_correlation.HPDinf.' name ' = hpd_interval(1);']);
+    eval(['oo_.' TYPE 'TheoreticalMoments.dsge.contemporeaneous_correlation.HPDsup.' name ' = hpd_interval(2);']);
+    eval(['oo_.' TYPE 'TheoreticalMoments.dsge.contemporeaneous_correlation.deciles.' name ' = p_deciles;']);
+    eval(['oo_.' TYPE 'TheoreticalMoments.dsge.contemporeaneous_correlation.density.' name ' = density;']);
+end
diff --git a/matlab/disp_moments.m b/matlab/disp_moments.m
index 478936e965..cd7ebf2b9c 100644
--- a/matlab/disp_moments.m
+++ b/matlab/disp_moments.m
@@ -69,6 +69,9 @@ end
 
 if options_.nocorr == 0
     corr = (y'*y/size(y,1))./(s'*s);
+    if options_.contemporaneous_correlation 
+        oo_.contemporaneous_correlation = corr;
+    end
     if options_.noprint == 0
         title = 'CORRELATION OF SIMULATED VARIABLES';
         if options_.hp_filter
diff --git a/matlab/disp_th_moments.m b/matlab/disp_th_moments.m
index 68a24d03e3..b6edede955 100644
--- a/matlab/disp_th_moments.m
+++ b/matlab/disp_th_moments.m
@@ -117,6 +117,9 @@ end
 
 if options_.nocorr == 0 && size(stationary_vars, 1) > 0
     corr = oo_.gamma_y{1}(i1,i1)./(sd(i1)*sd(i1)');
+    if options_.contemporaneous_correlation 
+        oo_.contemporaneous_correlation = corr;
+    end
     if ~options_.noprint,
         skipline()
         if options_.order == 2
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index c42edf4167..44694a40ab 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -480,6 +480,7 @@ end
 options_.filter_covariance = 0;
 options_.filter_decomposition = 0;
 options_.selected_variables_only = 0;
+options_.contemporaneous_correlation = 0;
 options_.initialize_estimated_parameters_with_the_prior_mode = 0;
 options_.estimation_dll = 0;
 % Misc
diff --git a/matlab/posterior_analysis.m b/matlab/posterior_analysis.m
index 0382f923fa..5209cb3f8c 100644
--- a/matlab/posterior_analysis.m
+++ b/matlab/posterior_analysis.m
@@ -55,7 +55,7 @@ switch type
             dsge_simulated_theoretical_covariance(SampleSize,M_,options_,oo_,'posterior');
     end
     oo_ = covariance_mc_analysis(SampleSize,'posterior',M_.dname,M_.fname,...
-                                 vartan,nvar,arg1,arg2,options_.mh_conf_sig,oo_);          
+                                 vartan,nvar,arg1,arg2,options_.mh_conf_sig,oo_,options_);          
   case 'decomposition'
     if nargin==narg1
         [nvar,vartan,NumberOfFiles] = ...
diff --git a/matlab/prior_analysis.m b/matlab/prior_analysis.m
index 69c1ca7d98..b942a509fa 100644
--- a/matlab/prior_analysis.m
+++ b/matlab/prior_analysis.m
@@ -56,7 +56,7 @@ switch type
             dsge_simulated_theoretical_covariance(SampleSize,M_,options_,oo_,'prior');
     end
     oo_ = covariance_mc_analysis(SampleSize,'prior',M_.dname,M_.fname,...
-                                 vartan,nvar,arg1,arg2,options_.mh_conf_sig,oo_);          
+                                 vartan,nvar,arg1,arg2,options_.mh_conf_sig,oo_,options_);          
   case 'decomposition'
     if nargin==narg1
         [nvar,vartan,NumberOfFiles] = ...
diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy
index 789286082b..392d3d3efe 100644
--- a/preprocessor/DynareBison.yy
+++ b/preprocessor/DynareBison.yy
@@ -106,7 +106,7 @@ class ParsingDriver;
 %token LYAPUNOV_FIXED_POINT_TOL LYAPUNOV_DOUBLING_TOL LYAPUNOV_SQUARE_ROOT_SOLVER_TOL LOG_DEFLATOR LOG_TREND_VAR LOG_GROWTH_FACTOR MARKOWITZ MARGINAL_DENSITY MAX MAXIT
 %token MFS MH_CONF_SIG MH_DROP MH_INIT_SCALE MH_JSCALE MH_MODE MH_NBLOCKS MH_REPLIC MH_RECOVER POSTERIOR_MAX_SUBSAMPLE_DRAWS MIN MINIMAL_SOLVING_PERIODS
 %token MODE_CHECK MODE_CHECK_NEIGHBOURHOOD_SIZE MODE_CHECK_SYMMETRIC_PLOTS MODE_CHECK_NUMBER_OF_POINTS MODE_COMPUTE MODE_FILE MODEL MODEL_COMPARISON MODEL_INFO MSHOCKS ABS SIGN
-%token MODEL_DIAGNOSTICS MODIFIEDHARMONICMEAN MOMENTS_VARENDO DIFFUSE_FILTER SUB_DRAWS TAPER_STEPS GEWEKE_INTERVAL MCMC_JUMPING_COVARIANCE MOMENT_CALIBRATION
+%token MODEL_DIAGNOSTICS MODIFIEDHARMONICMEAN MOMENTS_VARENDO CONTEMPORANEOUS_CORRELATION DIFFUSE_FILTER SUB_DRAWS TAPER_STEPS GEWEKE_INTERVAL MCMC_JUMPING_COVARIANCE MOMENT_CALIBRATION
 %token NUMBER_OF_PARTICLES RESAMPLING SYSTEMATIC GENERIC RESAMPLING_THRESHOLD RESAMPLING_METHOD KITAGAWA STRATIFIED SMOOTH
 %token FILTER_ALGORITHM PROPOSAL_APPROXIMATION CUBATURE UNSCENTED MONTECARLO DISTRIBUTION_APPROXIMATION
 %token <string_val> NAME
@@ -1068,6 +1068,7 @@ stoch_simul_primary_options : o_dr_algo
                             | o_drop
                             | o_ar
                             | o_nocorr
+                            | o_contemporaneous_correlation
                             | o_nofunctions
                             | o_nomoments
                             | o_nograph
@@ -1672,6 +1673,7 @@ estimation_options : o_datafile
                    | o_forecast
                    | o_smoother
                    | o_moments_varendo
+                   | o_contemporaneous_correlation
                    | o_filtered_vars
                    | o_kalman_algo
                    | o_kalman_tol
@@ -2706,6 +2708,7 @@ o_tex : TEX { driver.option_num("TeX", "1"); };
 o_forecast : FORECAST EQUAL INT_NUMBER { driver.option_num("forecast", $3); };
 o_smoother : SMOOTHER { driver.option_num("smoother", "1"); };
 o_moments_varendo : MOMENTS_VARENDO { driver.option_num("moments_varendo", "1"); };
+o_contemporaneous_correlation : CONTEMPORANEOUS_CORRELATION { driver.option_num("contemporaneous_correlation", "1"); };
 o_filtered_vars : FILTERED_VARS { driver.option_num("filtered_vars", "1"); };
 o_relative_irf : RELATIVE_IRF { driver.option_num("relative_irf", "1"); };
 o_kalman_algo : KALMAN_ALGO EQUAL INT_NUMBER { driver.option_num("kalman_algo", $3); };
diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll
index 81bcf01225..ce022df743 100644
--- a/preprocessor/DynareFlex.ll
+++ b/preprocessor/DynareFlex.ll
@@ -298,6 +298,7 @@ DATE -?[0-9]+([YyAa]|[Mm]([1-9]|1[0-2])|[Qq][1-4]|[Ww]([1-9]{1}|[1-4][0-9]|5[0-2
 <DYNARE_STATEMENT>dsge_var 	{return token::DSGE_VAR;}
 <DYNARE_STATEMENT>dsge_varlag 	{return token::DSGE_VARLAG;}
 <DYNARE_STATEMENT>moments_varendo {return token::MOMENTS_VARENDO;}
+<DYNARE_STATEMENT>contemporaneous_correlation {return token::CONTEMPORANEOUS_CORRELATION;}
 <DYNARE_STATEMENT>posterior_max_subsample_draws	{return token::POSTERIOR_MAX_SUBSAMPLE_DRAWS;}
 <DYNARE_STATEMENT>filtered_vars	{return token::FILTERED_VARS;}
 <DYNARE_STATEMENT>filter_step_ahead	{return token::FILTER_STEP_AHEAD;}
diff --git a/tests/TeX/fs2000_corr_ME.mod b/tests/TeX/fs2000_corr_ME.mod
index f8bbe3c4ad..bb19e5fef4 100644
--- a/tests/TeX/fs2000_corr_ME.mod
+++ b/tests/TeX/fs2000_corr_ME.mod
@@ -123,7 +123,7 @@ end;
 
 steady;
 
-stoch_simul(order=1,irf=20,graph_format=eps);
+stoch_simul(order=1,irf=20,graph_format=eps,contemporaneous_correlation);
 
 write_latex_original_model;
 write_latex_static_model;
@@ -143,7 +143,7 @@ stderr gy_obs, 1;
 corr gp_obs, gy_obs,0;
 end;
 
-estimation(order=1,datafile='../fs2000/fsdat_simul',mode_check,smoother,filter_decomposition,forecast = 8,filtered_vars,filter_step_ahead=[1,3],irf=20) m P c e W R k d y gy_obs;
+estimation(order=1,datafile='../fs2000/fsdat_simul',mode_check,smoother,filter_decomposition,forecast = 8,filtered_vars,filter_step_ahead=[1,3],irf=20,contemporaneous_correlation) m P c e W R k d y gy_obs;
 
 
 
@@ -159,7 +159,7 @@ stderr gp_obs, inv_gamma_pdf, 0.001, inf;
 //corr gp_obs, gy_obs,normal_pdf, 0, 0.2;
 end;
 
-estimation(mode_compute=9,order=1,datafile='../fs2000/fsdat_simul',mode_check,smoother,filter_decomposition,mh_replic=2002, mh_nblocks=2, mh_jscale=0.8,forecast = 8,bayesian_irf,filtered_vars,filter_step_ahead=[1,3],irf=20,moments_varendo) m P c e W R k d y;
+estimation(mode_compute=9,order=1,datafile='../fs2000/fsdat_simul',mode_check,smoother,filter_decomposition,mh_replic=2002, mh_nblocks=2, mh_jscale=0.8,forecast = 8,bayesian_irf,filtered_vars,filter_step_ahead=[1,3],irf=20,moments_varendo,contemporaneous_correlation) m P c e W R k d y;
 shock_decomposition y W R;
 
 collect_LaTeX_Files(M_);
-- 
GitLab