Merge branch 'master' into ecb-master

Fixed conflicts in:
	matlab/backward/backward_model_irf.m
	matlab/modules/dseries
	matlab/utilities/general/isint.m
	preprocessor/DynareBison.yy
	preprocessor/DynareFlex.ll
	preprocessor/ModFile.cc
	preprocessor/ParsingDriver.hh

doc/dynare.texi

@@ -363,7 +363,7 @@ both non-profit and for-profit purposes. Most of the source files are covered
 by the GNU General Public Licence (GPL) version 3 or later (there are some
 exceptions to this, see the file @file{license.txt} in Dynare distribution). It
 is available for the Windows, macOS, and Linux platforms and is fully
-documented through a user guide and a reference manual. Part of Dynare is
+documented in this reference manual. Part of Dynare is
 programmed in C++, while the rest is written using the
 @uref{http://www.mathworks.com/products/matlab/, MATLAB} programming language.
 The latter implies that commercially-available MATLAB software is required in
@@ -407,31 +407,42 @@ providing financial support.
 The present document is the reference manual for Dynare. It documents
 all commands and features in a systematic fashion.
 
-New users should rather begin with Dynare User Guide (@cite{Mancini
-(2007)}), distributed with Dynare and also available from the
-@uref{http://www.dynare.org,official Dynare web site}.
-
 Other useful sources of information include the
-@uref{http://www.dynare.org,Dynare wiki} and the
+@uref{http://www.dynare.org/DynareWiki,old Dynare wiki}, the
-@uref{http://www.dynare.org/phpBB3, Dynare forums}.
+@uref{https://github.com/DynareTeam/dynare/wiki,new Dynare wiki}, the documentation
+section of the @uref{http://www.dynare.org/documentation-and-support,Dynare
+website} and the @uref{https://forum.dynare.org/,Dynare forum}.
 
 @node Citing Dynare in your research
 @section Citing Dynare in your research
 
-If you would like to refer to Dynare in a research article, the
+You should cite Dynare if you use it in your research. The recommended way to
-recommended way is to cite the present manual, as follows:
+do this is to cite the present manual as:
 
 @quotation
-Stéphane Adjemian, Houtan Bastani, Michel Juillard, Frédéric Karamé,
+Stéphane Adjemian, Houtan Bastani, Michel Juillard, Frédéric Karamé, Junior
-Ferhat Mihoubi, George Perendia, Johannes Pfeifer, Marco Ratto and
+Maih, Ferhat Mihoubi, George Perendia, Johannes Pfeifer, Marco Ratto and
-Sébastien Villemot (2011), ``Dynare: Reference Manual, Version 4,''
+Sébastien Villemot (2011), ``Dynare: Reference Manual, Version 4,'' @i{Dynare
-@i{Dynare Working Papers}, 1, CEPREMAP
+Working Papers}, 1, CEPREMAP
 @end quotation
 
-Note that citing the Dynare Reference Manual in your research is a
+@noindent For convenience, you can copy and paste the following into your BibTeX file:
-good way to help the Dynare project.
 
-If you want to give a URL, use the address of the Dynare website:
+@verbatim
+@TechReport{Adjemianetal2011,
+  author      = {Adjemian, St\'ephane and Bastani, Houtan and Juillard, Michel and
+                 Karam\'e, Fr\'ederic and Maih, Junior and Mihoubi, Ferhat and
+                 Perendia, George and Pfeifer, Johannes and Ratto, Marco and
+                 Villemot, S\'ebastien},
+  title       = {Dynare: Reference Manual Version 4},
+  year        = {2011},
+  institution = {CEPREMAP},
+  type        = {Dynare Working Papers},
+  number      = {1},
+}
+@end verbatim
+
+@noindent If you want to give a URL, use the address of the Dynare website:
 @uref{http://www.dynare.org}.
 
 @node Installation and configuration
@@ -447,7 +458,7 @@ If you want to give a URL, use the address of the Dynare website:
 @node Software requirements
 @section Software requirements
 
-Packaged versions of Dynare are available for Windows XP/Vista/7/8,
+Packaged versions of Dynare are available for Windows XP/Vista/7/8/10,
 @uref{http://www.debian.org,Debian GNU/Linux},
 @uref{http://www.ubuntu.com/,Ubuntu} and macOS 10.8 or later.  Dynare should
 work on other systems, but some compilation steps are necessary in that case.
@@ -460,7 +471,7 @@ In order to run Dynare, you need one of the following:
 MATLAB version 7.5 (R2007b) or above (MATLAB R2009b 64-bit for macOS);
 
 @item
-GNU Octave version 3.6 or above.
+GNU Octave version 4.2.1 or above.
 @end itemize
 
 Packages of GNU Octave can be downloaded on the
@@ -4382,7 +4393,17 @@ period(s). The periods must be strictly positive. Conditional variances are give
 decomposition provides the decomposition of the effects of shocks upon
 impact. The results are stored in
 @code{oo_.conditional_variance_decomposition}
-(@pxref{oo_.conditional_variance_decomposition}). The variance decomposition is only conducted, if theoretical moments are requested, @i{i.e.} using the @code{periods=0}-option. In case of @code{order=2}, Dynare provides a second-order accurate approximation to the true second moments based on the linear terms of the second-order solution (see @cite{Kim, Kim, Schaumburg and Sims (2008)}). Note that the unconditional variance decomposition (@i{i.e.} at horizon infinity) is automatically conducted if theoretical moments are requested and if @code{nodecomposition} is not set (@pxref{oo_.variance_decomposition})
+(@pxref{oo_.conditional_variance_decomposition}).
+In the presence of measurement error, the @code{oo_.conditional_variance_decomposition} field will contain
+the variance contribution after measurement error has been taken out, i.e. the decomposition will 
+be conducted of the actual as opposed to the measured variables. The variance decomposition of the 
+measured variables will be stored in @code{oo_.conditional_variance_decomposition_ME} (@pxref{oo_.conditional_variance_decomposition_ME}). 
+The variance decomposition is only conducted, if theoretical moments are requested, @i{i.e.} using the @code{periods=0}-option. 
+In case of @code{order=2}, Dynare provides a second-order accurate approximation to the true 
+second moments based on the linear terms of the second-order solution (see @cite{Kim, Kim, 
+Schaumburg and Sims (2008)}). Note that the unconditional variance decomposition (@i{i.e.} at horizon infinity) 
+is automatically conducted if theoretical moments are requested and if @code{nodecomposition} 
+is not set (@pxref{oo_.variance_decomposition})
 
 @item pruning
 Discard higher order terms when iteratively computing simulations of
@@ -4611,9 +4632,25 @@ accurate approximation of the true second moments, see @code{conditional_varianc
 @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{i.e.} at horizon infinity).
-The first dimension corresponds to the endogenous variables (in the order of declaration) and 
+The first dimension corresponds to the endogenous variables (in the order of declaration after the command or in
+@code{M_.endo_names} if not specified) and 
 the second dimension corresponds to exogenous variables (in the order of declaration). 
-Numbers are in percent and sum up to 100 across columns.
+Numbers are in percent and sum up to 100 across columns. In the presence of measurement error, the
+field will contain the variance contribution after measurement error has been taken out, i.e. the decomposition will 
+be conducted of the actual as opposed to the measured variables.
+@end defvr
+
+@anchor{oo_.variance_decomposition_ME}
+@defvr {MATLAB/Octave variable} oo_.variance_decomposition_ME
+Field set after a run of @code{stoch_simul} when requesting theoretical moments (@code{periods=0}) if 
+measurement error is present.  
+It is similar to @ref{oo_.variance_decomposition}, but the decomposition will 
+be conducted of the measured variables. The field contains a matrix with the result 
+of the unconditional variance decomposition (@i{i.e.} at horizon infinity).
+The first dimension corresponds to the observed endogenous variables (in the order of declaration after the command)
+and the second dimension corresponds to exogenous variables (in the order of declaration), with the last column
+corresponding to the contribution of the measurement error. 
+Numbers are in percent and sum up to 100 across columns. 
 @end defvr
 
 @anchor{oo_.conditional_variance_decomposition}
@@ -4622,9 +4659,25 @@ After a run of @code{stoch_simul} with the
 @code{conditional_variance_decomposition} option, contains a
 three-dimensional array with the result of the decomposition. The
 first dimension corresponds to forecast horizons (as declared with the
-option), the second dimension corresponds to endogenous variables (in
+option), the second dimension corresponds to endogenous variables (in the order 
+of declaration after the command or in @code{M_.endo_names} if not specified)), 
+the third dimension corresponds to
+exogenous variables (in the order of declaration). In the presence of measurement error, the
+field will contain the variance contribution after measurement error has been taken out, 
+i.e. the decomposition will be conducted of the actual as opposed to the measured variables.
+@end defvr
+
+@anchor{oo_.conditional_variance_decomposition_ME}
+@defvr {MATLAB/Octave variable} oo_.conditional_variance_decomposition_ME
+Field set after a run of @code{stoch_simul} with the @code{conditional_variance_decomposition} 
+option if measurement error is present. It is similar to @ref{oo_.conditional_variance_decomposition}, but 
+the decomposition will be conducted of the measured variables. 
+It contains a three-dimensional array with the result of the decomposition. The
+first dimension corresponds to forecast horizons (as declared with the
+option), the second dimension corresponds to observed endogenous variables (in
 the order of declaration), the third dimension corresponds to
-exogenous variables (in the order of declaration).
+exogenous variables (in the order of declaration), with the last column
+corresponding to the contribution of the measurement error. 
 @end defvr
 
 @anchor{oo_.contemporaneous_correlation}
@@ -6232,9 +6285,13 @@ positive. Conditional variances are given by @math{var(y_{t+k}|t)}. For
 period 1, the conditional variance decomposition provides the
 decomposition of the effects of shocks upon impact. The results are
 stored in
-@code{oo_.PosteriorTheoreticalMoments.dsge.ConditionalVarianceDecomposition},
+@code{oo_.PosteriorTheoreticalMoments.dsge.ConditionalVarianceDecomposition}. 
-but currently there is no displayed output. Note that this option requires the
+Note that this option requires the
-option @code{moments_varendo} to be specified.
+option @code{moments_varendo} to be specified. In the presence of measurement error, the
+field will contain the variance contribution after measurement error has been taken out, 
+i.e. the decomposition will be conducted of the actual as opposed to the measured variables. 
+The variance decomposition of the measured variables will be stored in 
+@code{oo_.PosteriorTheoreticalMoments.dsge.ConditionalVarianceDecompositionME}. 
 
 @item filtered_vars
 @anchor{filtered_vars} Triggers the computation of the posterior
@@ -7027,14 +7084,29 @@ Auto- and cross-correlation of endogenous variables. Fields are vectors with cor
 
 
 @item VarianceDecomposition
+@anchor{VarianceDecomposition}
 Decomposition of variance (unconditional variance, @i{i.e.} at horizon infinity)@footnote{When the shocks are correlated, it
 is the decomposition of orthogonalized shocks via Cholesky
 decomposition according to the order of declaration of shocks
 (@pxref{Variable declarations})}
 
+@item VarianceDecompositionME
+Same as @ref{VarianceDecomposition}, but contains the decomposition of the
+measured as opposed to the actual variable. The joint contribution of the measurement error
+will be saved in a field named @code{ME}.
+
 @item ConditionalVarianceDecomposition
+@anchor{ConditionalVarianceDecomposition}
 Only if the @code{conditional_variance_decomposition} option has been
-specified
+specified. In the presence of measurement error, the field will contain
+the variance contribution after measurement error has been taken out, i.e. the decomposition will 
+be conducted of the actual as opposed to the measured variables. 
+
+@item ConditionalVarianceDecompositionME
+Only if the @code{conditional_variance_decomposition} option has been
+specified. Same as @ref{ConditionalVarianceDecomposition}, but contains the decomposition of the
+measured as opposed to the actual variable. The joint contribution of the measurement error
+will be saved in a field named @code{ME}.
 
 @end table
 
@@ -15092,10 +15164,6 @@ Lubik, Thomas and Frank Schorfheide (2007): ``Do Central Banks Respond
 to Exchange Rate Movements? A Structural Investigation,'' @i{Journal
 of Monetary Economics}, 54(4), 1069--1087
 
-@item
-Mancini-Griffoli, Tommaso (2007): ``Dynare User Guide: An introduction
-to the solution and estimation of DSGE models''
-
 @item
 Murray, Lawrence M., Emlyn M. Jones and John Parslow (2013): ``On Disturbance State-Space Models and the Particle Marginal
 Metropolis-Hastings Sampler'', @i{SIAM/ASA Journal on Uncertainty Quantification}, 1, 494–521.

matlab/GetOneDraw.m

@@ -1,10 +1,15 @@
-function [xparams, logpost] = GetOneDraw(type)
+function [xparams, logpost] = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_)
-% function [xparams, logpost] = GetOneDraw(type)
+% function [xparams, logpost] = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_)
 % draws one parameter vector and its posterior from MCMC or the prior
 %
 % INPUTS
 %    type:      [string]       'posterior': draw from MCMC draws
 %                              'prior': draw from prior
+%    M_         [structure]     Definition of the model
+%    estim_params_ [structure]  characterizing parameters to be estimated
+%    oo_         [structure]    Storage of results
+%    options_    [structure]    Options
+%    bayestopt_  [structure]    describing the priors
 %
 % OUTPUTS
 %    xparams:   vector of estimated parameters (drawn from posterior or prior distribution)
@@ -36,6 +41,6 @@ switch type
   case 'prior'
     xparams = prior_draw();
     if nargout>1
-        logpost = evaluate_posterior_kernel(xparams');
+        logpost = evaluate_posterior_kernel(xparams',M_,estim_params_,oo_,options_,bayestopt_);
     end
 end

matlab/PosteriorIRF.m

@@ -174,7 +174,7 @@ localVars.type=type;
 if strcmpi(type,'posterior')
     while b<B
         b = b + 1;
-        x(b,:) = GetOneDraw(type);
+        x(b,:) = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_);
     end
 end
 

matlab/PosteriorIRF_core1.m

@@ -141,7 +141,7 @@ while fpar<B
     irun = irun+1;
     irun2 = irun2+1;
     if strcmpi(type,'prior')
-        deep = GetOneDraw(type);
+        deep = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_);
     else
         deep = x(fpar,:);
     end

matlab/backward/backward_model_irf.m

@@ -85,9 +85,6 @@ else
             if initialconditionperiod>=shock.dates(1)
                 error('In experiment n°%s, the shock period must follow %s!', string(initialconditionperiod))
             end
-            if shock.nobs>1
-                error('Shocks over multiple periods not implemented yet!')
-            end
         end
     end
 end
@@ -125,8 +122,14 @@ end
     simul_backward_model_init(initialcondition, periods, options_, M_, oo_, Innovations);
 
 % Get the covariance matrix of the shocks.
-Sigma = M_.Sigma_e + 1e-14*eye(M_.exo_nbr);
+if ~deterministicshockflag
-sigma = transpose(chol(Sigma));
+    if ~nnz(M_.Sigma_e)
+        Sigma = M_.Sigma_e + 1e-14*eye(M_.exo_nbr);
+        sigma = transpose(chol(Sigma));
+    else
+        error('You did not specify the size of the shocks!')
+    end
+end
 
 % Initialization of the returned argument. Each will be a dseries object containing the IRFS for the endogenous variables listed in the third input argument.
 deviations = struct();
@@ -154,10 +157,12 @@ for i=1:length(listofshocks)
     % Add the shock.
     if deterministicshockflag
         shock = listofshocks{i};
-        timid = shock.dates(1)-initialconditionperiod;
+        timid = shock.dates-initialconditionperiod;
         for j=1:shock.vobs
             k = find(strcmp(shock.name{i}, exonames));
-            innovations(timid,:) = innovations(timid,:) + shock.data(1,j);
+            for l=1:length(timid)
+                innovations(timid(l),k) = innovations(timid(l),k) + shock.data(l,j);
+            end
         end
     else
         j = find(strcmp(listofshocks{i}, exonames));

matlab/check_posterior_analysis_data.m

@@ -86,7 +86,7 @@ switch type
   case 'conditional decomposition'
     generic_post_data_file_name = 'PosteriorConditionalVarianceDecomposition';
   otherwise
-    disp('This feature is not yest implemented!')
+    disp('This feature is not yet implemented!')
 end
 pdfinfo = dir([ MetropolisFolder filesep M_.fname '_' generic_post_data_file_name '*']);
 if isempty(pdfinfo)

matlab/compute_moments_varendo.m

@@ -34,7 +34,7 @@ function oo_ = compute_moments_varendo(type,options_,M_,oo_,var_list_)
 % 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)... ');
+fprintf('Estimation::compute_moments_varendo: I''m computing endogenous moments (this may take a while)... \n');
 
 if strcmpi(type,'posterior')
     posterior = 1;
@@ -114,6 +114,7 @@ if M_.exo_nbr > 1
                 end
             end
             title='Posterior mean variance decomposition (in percent)';
+            save_name_string='dsge_post_mean_var_decomp_uncond';
         else
             for i=1:NumberOfEndogenousVariables
                 for j=1:NumberOfExogenousVariables
@@ -122,6 +123,7 @@ if M_.exo_nbr > 1
                 end
             end
             title='Prior mean variance decomposition (in percent)';
+            save_name_string='dsge_prior_mean_var_decomp_uncond';
         end
         title=add_filter_subtitle(title,options_);
         headers = M_.exo_names;
@@ -135,11 +137,57 @@ if M_.exo_nbr > 1
             headers = char(' ',headers);
             labels = deblank(var_list_tex);
             lh = size(labels,2)+2;
-            dyn_latex_table(M_,options_,title,'dsge_post_mean_var_decomp_uncond',headers,labels,100*temp,lh,8,2);
+            dyn_latex_table(M_,options_,title,save_name_string,headers,labels,100*temp,lh,8,2);
         end
         skipline();
     end
-    % CONDITIONAL VARIANCE DECOMPOSITION.
+    skipline();
+    if ~all(M_.H==0)
+        [observable_name_requested_vars,varlist_pos]=intersect(var_list_,options_.varobs,'stable');
+        if ~isempty(observable_name_requested_vars)
+            NumberOfObservedEndogenousVariables=length(observable_name_requested_vars);
+            temp=NaN(NumberOfObservedEndogenousVariables,NumberOfExogenousVariables+1);
+            if posterior
+                for i=1:NumberOfObservedEndogenousVariables
+                    for j=1:NumberOfExogenousVariables
+                        temp(i,j,:)=oo_.PosteriorTheoreticalMoments.dsge.VarianceDecompositionME.Mean.(deblank(observable_name_requested_vars{i,1})).(deblank(M_.exo_names(j,:)));
+                    end
+                    endo_index_varlist=strmatch(deblank(observable_name_requested_vars{i,1}),var_list_,'exact');
+                    oo_ = posterior_analysis('decomposition',var_list_(endo_index_varlist,:),'ME',[],options_,M_,oo_);
+                    temp(i,j+1,:)=oo_.PosteriorTheoreticalMoments.dsge.VarianceDecompositionME.Mean.(deblank(observable_name_requested_vars{i,1})).('ME');
+                end
+                title='Posterior mean variance decomposition (in percent) with measurement error';
+                save_name_string='dsge_post_mean_var_decomp_uncond_ME';
+            else
+                for i=1:NumberOfObservedEndogenousVariables
+                    for j=1:NumberOfExogenousVariables
+                        temp(i,j,:)=oo_.PriorTheoreticalMoments.dsge.VarianceDecompositionME.Mean.(deblank(observable_name_requested_vars(i,:))).(deblank(M_.exo_names(j,:)));
+                    end
+                    endo_index_varlist=strmatch(deblank(observable_name_requested_vars{i,1}),var_list_,'exact');
+                    oo_ = prior_analysis('decomposition',var_list_(endo_index_varlist,:),'ME',[],options_,M_,oo_);
+                    temp(i,j+1,:)=oo_.PriorTheoreticalMoments.dsge.VarianceDecompositionME.Mean.(deblank(observable_name_requested_vars{i,1})).('ME');
+                end
+                title='Prior mean variance decomposition (in percent) with measurement error';
+                save_name_string='dsge_prior_mean_var_decomp_uncond_ME';
+            end
+            title=add_filter_subtitle(title,options_);
+            headers = M_.exo_names;
+            headers(M_.exo_names_orig_ord,:) = headers;
+            headers = char(' ',headers,'ME');
+            lh = size(deblank(var_list_),2)+2;
+            dyntable(options_,title,headers,deblank(char(observable_name_requested_vars)),100* ...
+                temp,lh,8,2);
+            if options_.TeX
+                headers=M_.exo_names_tex;
+                headers = char(' ',headers,'ME');
+                labels = deblank(var_list_tex(varlist_pos,:));
+                lh = size(labels,2)+2;
+                dyn_latex_table(M_,options_,title,save_name_string,headers,labels,100*temp,lh,8,2);
+            end
+            skipline();
+        end
+    end
+% CONDITIONAL VARIANCE DECOMPOSITION.
     if Steps
         temp=NaN(NumberOfEndogenousVariables,NumberOfExogenousVariables,length(Steps));
         if posterior
@@ -150,6 +198,7 @@ if M_.exo_nbr > 1
                 end
             end
             title='Posterior mean conditional variance decomposition (in percent)';
+            save_name_string='dsge_post_mean_var_decomp_cond_h';
         else
             for i=1:NumberOfEndogenousVariables
                 for j=1:NumberOfExogenousVariables
@@ -158,6 +207,7 @@ if M_.exo_nbr > 1
                 end
             end
             title='Prior mean conditional variance decomposition (in percent)';
+            save_name_string='dsge_prior_mean_var_decomp_cond_h';
         end
         for step_iter=1:length(Steps)
             title_print=[title, ' Period ' int2str(Steps(step_iter))];
@@ -172,10 +222,56 @@ if M_.exo_nbr > 1
                 headers = char(' ',headers);
                 labels = deblank(var_list_tex);
                 lh = size(labels,2)+2;
-                dyn_latex_table(M_,options_,title_print,['dsge_post_mean_var_decomp_cond_h',int2str(Steps(step_iter))],headers,labels,100*temp(:,:,step_iter),lh,8,2);
+                dyn_latex_table(M_,options_,title_print,[save_name_string,int2str(Steps(step_iter))],headers,labels,100*temp(:,:,step_iter),lh,8,2);
             end
         end
         skipline();
+        if ~all(M_.H==0)
+            if ~isempty(observable_name_requested_vars)
+                NumberOfObservedEndogenousVariables=length(observable_name_requested_vars);
+                temp=NaN(NumberOfObservedEndogenousVariables,NumberOfExogenousVariables+1,length(Steps));
+                if posterior
+                    for i=1:NumberOfObservedEndogenousVariables
+                        for j=1:NumberOfExogenousVariables
+                            temp(i,j,:)=oo_.PosteriorTheoreticalMoments.dsge.ConditionalVarianceDecompositionME.Mean.(deblank(observable_name_requested_vars{i,1})).(deblank(M_.exo_names(j,:)));
+                        end
+                        endo_index_varlist=strmatch(deblank(observable_name_requested_vars{i,1}),var_list_,'exact');
+                        oo_ = posterior_analysis('conditional decomposition',endo_index_varlist,'ME',Steps,options_,M_,oo_);
+                        temp(i,j+1,:)=oo_.PosteriorTheoreticalMoments.dsge.ConditionalVarianceDecompositionME.Mean.(deblank(observable_name_requested_vars{i,1})).('ME');
+                    end                    
+                    title='Posterior mean conditional variance decomposition (in percent) with measurement error';
+                    save_name_string='dsge_post_mean_var_decomp_ME_cond_h';
+                else
+                   for i=1:NumberOfObservedEndogenousVariables
+                        for j=1:NumberOfExogenousVariables
+                            temp(i,j,:)=oo_.PriorTheoreticalMoments.dsge.ConditionalVarianceDecompositionME.Mean.(deblank(observable_name_requested_vars(i,:))).(deblank(M_.exo_names(j,:)));
+                        end
+                        endo_index_varlist=strmatch(deblank(observable_name_requested_vars{i,1}),var_list_,'exact');
+                        oo_ = prior_analysis('conditional decomposition',endo_index_varlist,'ME',Steps,options_,M_,oo_);
+                        temp(i,j+1,:)=oo_.PriorTheoreticalMoments.dsge.ConditionalVarianceDecompositionME.Mean.(deblank(observable_name_requested_vars{i,1})).('ME');
+                    end                    
+                    title='Prior mean conditional variance decomposition (in percent) with measurement error';
+                    save_name_string='dsge_prior_mean_var_decomp_ME_cond_h';
+                end
+                for step_iter=1:length(Steps)
+                    title_print=[title, ' Period ' int2str(Steps(step_iter))];
+                    headers = M_.exo_names;
+                    headers(M_.exo_names_orig_ord,:) = headers;
+                    headers = char(' ',headers,'ME');
+                    lh = size(deblank(var_list_),2)+2;
+                    dyntable(options_,title_print,headers,deblank(char(observable_name_requested_vars)),100* ...
+                        temp(:,:,step_iter),lh,8,2);
+                    if options_.TeX
+                        headers=M_.exo_names_tex;
+                        headers = char(' ',headers,'ME');
+                        labels = deblank(var_list_tex(varlist_pos,:));
+                        lh = size(labels,2)+2;
+                        dyn_latex_table(M_,options_,title_print,[save_name_string,int2str(Steps(step_iter))],headers,labels,100*temp(:,:,step_iter),lh,8,2);
+                    end
+                end
+                skipline();
+            end
+        end        
     end
 end
 

matlab/conditional_variance_decomposition.m

@@ -1,20 +1,21 @@
-function ConditionalVarianceDecomposition = conditional_variance_decomposition(StateSpaceModel, Steps, SubsetOfVariables,sigma_e_is_diagonal)
+function [ConditionalVarianceDecomposition, ConditionalVarianceDecomposition_ME]= conditional_variance_decomposition(StateSpaceModel, Steps, SubsetOfVariables,sigma_e_is_diagonal)
 % This function computes the conditional variance decomposition of a given state space model
 % for a subset of endogenous variables.
 %
 % INPUTS
 %   StateSpaceModel     [structure]   Specification of the state space model.
 %   Steps               [integer]     1*h vector of dates.
-%   SubsetOfVariables   [integer]     1*q vector of indices.
+%   SubsetOfVariables   [integer]     1*q vector of indices (declaration order).
 %
 % OUTPUTS
 %   ConditionalVarianceDecomposition  [double] [n h p] array, where
 %                                                    n is equal to length(SubsetOfVariables)
 %                                                    h is the number of Steps
 %                                                    p is the number of state innovations and
-% SPECIAL REQUIREMENTS
+%   ConditionalVarianceDecomposition_ME  [double] [m h p] array, where
-%
+%                                                    m is equal to length(intersect(SubsetOfVariables,varobs))
-% [1] In this version, absence of measurement errors is assumed...
+%                                                    h is the number of Steps
+%                                                    p is the number of state innovations and
 
 % Copyright (C) 2010-2017 Dynare Team
 %
@@ -82,4 +83,22 @@ for i=1:number_of_state_innovations
     for h = 1:length(Steps)
         ConditionalVarianceDecomposition(:,h,i) = squeeze(ConditionalVariance(:,h,i))./SumOfVariances(:,h);
     end
+end
+
+% get intersection of requested variables and observed variables with
+% Measurement error
+
+if ~all(StateSpaceModel.measurement_error==0)
+    [observable_pos,index_subset,index_observables]=intersect(SubsetOfVariables,StateSpaceModel.observable_pos,'stable');
+    ME_Variance=diag(StateSpaceModel.measurement_error);
+    
+    ConditionalVarianceDecomposition_ME = zeros(length(observable_pos),length(Steps),number_of_state_innovations+1);
+    for i=1:number_of_state_innovations
+        for h = 1:length(Steps)
+            ConditionalVarianceDecomposition_ME(:,h,i) = squeeze(ConditionalVariance(index_subset,h,i))./(SumOfVariances(index_subset,h)+ME_Variance(index_observables));
+        end
+    end
+    ConditionalVarianceDecomposition_ME(:,:,number_of_state_innovations+1)=1-sum(ConditionalVarianceDecomposition_ME(:,:,1:number_of_state_innovations),3);
+else
+    ConditionalVarianceDecomposition_ME=[];    
 end

matlab/conditional_variance_decomposition_ME_mc_analysis.m

@@ -0,0 +1,136 @@
+function oo_ = ...
+    conditional_variance_decomposition_ME_mc_analysis(NumberOfSimulations, type, dname, fname, Steps, exonames, exo, var_list, endogenous_variable_index, mh_conf_sig, oo_,options_)
+% This function analyses the (posterior or prior) distribution of the
+% endogenous variables' conditional variance decomposition with measurement error.
+%
+% INPUTS
+%   NumberOfSimulations     [integer]           scalar, number of simulations.
+%   type                    [string]            'prior' or 'posterior'
+%   dname                   [string]            directory name where to save
+%   fname                   [string]            name of the mod-file
+%   Steps                   [integers]          horizons at which to conduct decomposition
+%   exonames                [string]            (n_exo*char_length) character array with names of exogenous variables
+%   exo                     [string]            name of current exogenous
+%                                               variable
+%   var_list                [string]            (n_endo*char_length) character array with name
+%                                               of endogenous variables
+%   endogenous_variable_index [integer]         index of the current
+%                                               endogenous variable
+%   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.
+%
+% OUTPUTS
+%   oo_          [structure]        Dynare structure where the results are saved.
+
+% Copyright (C) 2017 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/>.
+
+if strcmpi(type,'posterior')
+    TYPE = 'Posterior';
+    PATH = [dname '/metropolis/'];
+else
+    TYPE = 'Prior';
+    PATH = [dname '/prior/moments/'];
+end
+
+% $$$ indx = check_name(vartan,var);
+% $$$ if isempty(indx)
+% $$$     disp([ type '_analysis:: ' var ' is not a stationary endogenous variable!'])
+% $$$     return
+% $$$ end
+% $$$ endogenous_variable_index = sum(1:indx);
+exogenous_variable_index = check_name(exonames,exo);
+if isempty(exogenous_variable_index)
+    if isequal(exo,'ME')
+       exogenous_variable_index=size(exonames,1)+1; 
+    else
+        disp([ type '_analysis:: ' exo ' is not a declared exogenous variable!'])
+    return
+    end
+end
+
+[observable_pos_requested_vars,index_subset,index_observables]=intersect(var_list,options_.varobs,'stable');
+matrix_pos=strmatch(var_list(endogenous_variable_index,:),var_list(index_subset,:),'exact');
+name_1 = deblank(var_list(endogenous_variable_index,:));
+name_2 = deblank(exo);
+name = [ name_1 '.' name_2 ];
+
+if isfield(oo_, [ TYPE 'TheoreticalMoments' ])
+    temporary_structure = oo_.([TYPE 'TheoreticalMoments']);
+    if isfield(temporary_structure,'dsge')
+        temporary_structure = oo_.([TYPE 'TheoreticalMoments']).dsge;
+        if isfield(temporary_structure,'ConditionalVarianceDecompositionME')
+            temporary_structure = oo_.([TYPE 'TheoreticalMoments']).dsge.ConditionalVarianceDecompositionME.Mean;
+            if isfield(temporary_structure,name)
+                if sum(Steps-temporary_structure.(name)(1,:)) == 0
+                    % Nothing (new) to do here...
+                    return
+                end
+            end
+        end
+    end
+end
+
+ListOfFiles = dir([ PATH  fname '_' TYPE 'ConditionalVarianceDecompME*.mat']);
+i1 = 1; tmp = zeros(NumberOfSimulations,length(Steps));
+for file = 1:length(ListOfFiles)
+    load([ PATH ListOfFiles(file).name ]);
+    % 4D-array (endovar,time,exovar,simul)
+    i2 = i1 + size(Conditional_decomposition_array_ME,4) - 1;
+    tmp(i1:i2,:) = transpose(dynare_squeeze(Conditional_decomposition_array_ME(matrix_pos,:,exogenous_variable_index,:)));
+    i1 = i2+1;
+end
+
+p_mean = NaN(1,length(Steps));
+p_median = NaN(1,length(Steps));
+p_variance = NaN(1,length(Steps));
+p_deciles = NaN(9,length(Steps));
+if options_.estimation.moments_posterior_density.indicator
+    p_density = NaN(2^9,2,length(Steps));
+end
+p_hpdinf = NaN(1,length(Steps));
+p_hpdsup = NaN(1,length(Steps));
+for i=1:length(Steps)
+    if options_.estimation.moments_posterior_density.indicator
+        [pp_mean, pp_median, pp_var, hpd_interval, pp_deciles, pp_density] = ...
+            posterior_moments(tmp(:,i),1,mh_conf_sig);
+        p_density(:,:,i) = pp_density;
+    else
+        [pp_mean, pp_median, pp_var, hpd_interval, pp_deciles] = ...
+            posterior_moments(tmp(:,i),0,mh_conf_sig);
+    end
+    p_mean(i) = pp_mean;
+    p_median(i) = pp_median;
+    p_variance(i) = pp_var;
+    p_deciles(:,i) = pp_deciles;
+    p_hpdinf(i) = hpd_interval(1);
+    p_hpdsup(i) = hpd_interval(2);
+end
+
+FirstField = sprintf('%sTheoreticalMoments', TYPE);
+
+oo_.(FirstField).dsge.ConditionalVarianceDecompositionME.Steps = Steps;
+oo_.(FirstField).dsge.ConditionalVarianceDecompositionME.Mean.(name_1).(name_2) = p_mean;
+oo_.(FirstField).dsge.ConditionalVarianceDecompositionME.Median.(name_1).(name_2) = p_median;
+oo_.(FirstField).dsge.ConditionalVarianceDecompositionME.Variance.(name_1).(name_2) = p_variance;
+oo_.(FirstField).dsge.ConditionalVarianceDecompositionME.HPDinf.(name_1).(name_2) = p_hpdinf;
+oo_.(FirstField).dsge.ConditionalVarianceDecompositionME.HPDsup.(name_1).(name_2) = p_hpdsup;
+oo_.(FirstField).dsge.ConditionalVarianceDecompositionME.deciles.(name_1).(name_2)  = p_deciles;
+if options_.estimation.moments_posterior_density.indicator
+    oo_.(FirstField).dsge.ConditionalVarianceDecompositionME.density.(name_1).(name_2) = p_density;
+end

matlab/conditional_variance_decomposition_mc_analysis.m

@@ -56,7 +56,9 @@ end
 % $$$ endogenous_variable_index = sum(1:indx);
 exogenous_variable_index = check_name(exonames,exo);
 if isempty(exogenous_variable_index)
-    disp([ type '_analysis:: ' exo ' is not a declared exogenous variable!'])
+    if ~isequal(exo,'ME')
+        disp([ type '_analysis:: ' exo ' is not a declared exogenous variable!'])        
+    end
     return
 end
 

matlab/disp_moments.m

@@ -48,6 +48,25 @@ end
 
 y = y(ivar,options_.drop+1:end)';
 
+ME_present=0;
+if ~all(M_.H==0)
+    [observable_pos_requested_vars,index_subset,index_observables]=intersect(ivar,options_.varobs_id,'stable');
+    if ~isempty(observable_pos_requested_vars)
+        ME_present=1;
+        i_ME = setdiff([1:size(M_.H,1)],find(diag(M_.H) == 0)); % find ME with 0 variance
+        chol_S = chol(M_.H(i_ME,i_ME)); %decompose rest
+        shock_mat=zeros(options_.periods,size(M_.H,1)); %initialize
+        shock_mat(:,i_ME)=randn(length(i_ME),options_.periods)'*chol_S;
+        y_ME = y(:,index_subset)+shock_mat(options_.drop+1:end,index_observables);       
+        y_ME_only = shock_mat(options_.drop+1:end,index_observables);       
+        m_ME = mean(y_ME);
+        y_ME=get_filtered_time_series(y_ME,m_ME,options_);
+        y_ME_only_filtered=get_filtered_time_series(y_ME_only,mean(y_ME_only),options_);
+        s2_ME = mean(y_ME.*y_ME);
+    end
+end
+
+
 m = mean(y);
 
 % filter series
@@ -151,8 +170,13 @@ if ~options_.nodecomposition
             y_sim_one_shock = simult_(y0,oo_.dr,temp_shock_mat,options_.order);
             y_sim_one_shock=y_sim_one_shock(ivar,1+options_.drop+1:end)';
             y_sim_one_shock=get_filtered_time_series(y_sim_one_shock,mean(y_sim_one_shock),options_);
-            oo_.variance_decomposition(:,i_exo_var(shock_iter))=var(y_sim_one_shock)./s2*100;
+            oo_.variance_decomposition(:,i_exo_var(shock_iter))=var(y_sim_one_shock)./s2*100;            
         end
+        if ME_present
+            oo_.variance_decomposition_ME=oo_.variance_decomposition(index_subset,:)...
+                .*repmat((s2(index_subset)./s2_ME)',1,length(i_exo_var));
+            oo_.variance_decomposition_ME(:,end+1)=var(y_ME_only_filtered)./s2_ME*100;
+        end        
         if ~options_.noprint %options_.nomoments == 0
             skipline()
             title='VARIANCE DECOMPOSITION SIMULATING ONE SHOCK AT A TIME (in percent)';
@@ -164,12 +188,21 @@ if ~options_.nodecomposition
             headers = char(' ',headers);
             lh = size(deblank(M_.endo_names(ivar,:)),2)+2;
             dyntable(options_,title,char(headers,'Tot. lin. contr.'),deblank(M_.endo_names(ivar,:)),[oo_.variance_decomposition sum(oo_.variance_decomposition,2)],lh,8,2);
+            if ME_present
+                headers_ME=char(headers,'ME');
+                dyntable(options_,[title,' WITH MEASUREMENT ERROR'],char(headers_ME,'Tot. lin. contr.'),deblank(M_.endo_names(ivar(index_subset), ...
+                    :)),[oo_.variance_decomposition_ME sum(oo_.variance_decomposition_ME,2)],lh,8,2);
+            end
             if options_.TeX
                 headers=M_.exo_names_tex;
                 headers = char(' ',headers);
                 labels = deblank(M_.endo_names_tex(ivar,:));
                 lh = size(labels,2)+2;
                 dyn_latex_table(M_,options_,title,'sim_var_decomp',char(headers,'Tot. lin. contr.'),labels_TeX,[oo_.variance_decomposition sum(oo_.variance_decomposition,2)],lh,8,2);
+                if ME_present
+                    headers_ME=char(headers,'ME');
+                    dyn_latex_table(M_,options_,[title,' WITH MEASUREMENT ERROR'],'sim_var_decomp_ME',char(headers_ME,'Tot. lin. contr.'),labels_TeX(ivar(index_subset),:),[oo_.variance_decomposition_ME sum(oo_.variance_decomposition_ME,2)],lh,8,2);
+                end
             end
 
             if options_.order == 1

matlab/disp_th_moments.m

@@ -52,9 +52,22 @@ z = [ m sd s2 ];
 oo_.mean = m;
 oo_.var = oo_.gamma_y{1};
 
+ME_present=0;
+if ~all(M_.H==0)
+    [observable_pos_requested_vars,index_subset,index_observables]=intersect(ivar,options_.varobs_id,'stable');
+    if ~isempty(observable_pos_requested_vars)
+        ME_present=1;
+    end
+end
+                
 if size(stationary_vars, 1) > 0
     if ~nodecomposition
         oo_.variance_decomposition=100*oo_.gamma_y{options_.ar+2};
+        if ME_present
+            ME_Variance=diag(M_.H);
+            oo_.variance_decomposition_ME=oo_.variance_decomposition(index_subset,:).*repmat(diag(oo_.var(index_subset,index_subset))./(diag(oo_.var(index_subset,index_subset))+ME_Variance(index_observables)),1,M_.exo_nbr);
+            oo_.variance_decomposition_ME(:,end+1)=100-sum(oo_.variance_decomposition_ME,2);
+        end
     end
     if ~options_.noprint %options_.nomoments == 0
         if options_.order == 2
@@ -88,12 +101,22 @@ if size(stationary_vars, 1) > 0
             dyntable(options_,title,headers,deblank(M_.endo_names(ivar(stationary_vars), ...
                                                               :)),100* ...
                      oo_.gamma_y{options_.ar+2}(stationary_vars,:),lh,8,2);
+            if ME_present
+                [stationary_observables,pos_index_subset]=intersect(index_subset,stationary_vars,'stable');
+                headers_ME=char(headers,'ME');
+                dyntable(options_,[title,' WITH MEASUREMENT ERROR'],headers_ME,deblank(M_.endo_names(ivar(stationary_observables), ...
+                                                  :)),oo_.variance_decomposition_ME(pos_index_subset,:),lh,8,2);
+            end
             if options_.TeX
                 headers=M_.exo_names_tex;
                 headers = char(' ',headers);
                 labels = deblank(M_.endo_names_tex(ivar(stationary_vars),:));
                 lh = size(labels,2)+2;
                 dyn_latex_table(M_,options_,title,'th_var_decomp_uncond',headers,labels,100*oo_.gamma_y{options_.ar+2}(stationary_vars,:),lh,8,2);
+                if ME_present
+                    headers_ME=char(headers,'ME');
+                    dyn_latex_table(M_,options_,[title,' WITH MEASUREMENT ERROR'],'th_var_decomp_uncond_ME',headers_ME,labels,oo_.variance_decomposition_ME(pos_index_subset,:),lh,8,2);
+                end
             end
         end
     end
@@ -106,11 +129,17 @@ if size(stationary_vars, 1) > 0
         [StateSpaceModel.transition_matrix,StateSpaceModel.impulse_matrix] = kalman_transition_matrix(dr,(1:M_.endo_nbr)',M_.nstatic+(1:M_.nspred)',M_.exo_nbr);
         StateSpaceModel.state_innovations_covariance_matrix = M_.Sigma_e;
         StateSpaceModel.order_var = dr.order_var;
-        oo_.conditional_variance_decomposition = conditional_variance_decomposition(StateSpaceModel,conditional_variance_steps,ivar);
+        StateSpaceModel.measurement_error=M_.H;
+        StateSpaceModel.observable_pos=options_.varobs_id;
+        [oo_.conditional_variance_decomposition, oo_.conditional_variance_decomposition_ME]= conditional_variance_decomposition(StateSpaceModel,conditional_variance_steps,ivar);
 
         if options_.noprint == 0
             display_conditional_variance_decomposition(oo_.conditional_variance_decomposition,conditional_variance_steps,...
                                                        ivar,M_,options_);
+            if ME_present
+                display_conditional_variance_decomposition(oo_.conditional_variance_decomposition_ME,conditional_variance_steps,...
+                                                       observable_pos_requested_vars,M_,options_);               
+            end                                                       
         end
     end
 end

matlab/display_conditional_variance_decomposition.m

@@ -30,36 +30,56 @@ function display_conditional_variance_decomposition(conditional_decomposition_ar
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
+if size(conditional_decomposition_array,3)==M_.exo_nbr %no ME input
+    shock_number=M_.exo_nbr;
+    headers = M_.exo_names;
+    headers(M_.exo_names_orig_ord,:) = headers;
+    if options_.TeX
+        headers_TeX=char('',deblank(M_.exo_names_tex));
+    end
+    title_addon='';
+elseif size(conditional_decomposition_array,3)==M_.exo_nbr+1 %ME input
+    shock_number=M_.exo_nbr+1;
+    headers = M_.exo_names;
+    headers(M_.exo_names_orig_ord,:) = headers;
+    headers=char(headers,'ME');
+    if options_.TeX
+        headers_TeX=char('',deblank(strvcat(M_.exo_names_tex,'ME')));
+    end
+    title_addon=' - WITH MEASUREMENT ERROR';
+else
+    error('display_conditional_variance_decomposition:: This case should not happen. Please contact the developers')
+end
+
 if options_.order == 2
     skipline()
-    title='APPROXIMATED CONDITIONAL VARIANCE DECOMPOSITION (in percent)';
+    title=['APPROXIMATED CONDITIONAL VARIANCE DECOMPOSITION (in percent)',title_addon];
     disp(title)
 else
     skipline()
-    title='CONDITIONAL VARIANCE DECOMPOSITION (in percent)';
+    title=['CONDITIONAL VARIANCE DECOMPOSITION (in percent)',title_addon];
     disp(title)
 end
 
-vardec_i = zeros(length(SubsetOfVariables),M_.exo_nbr);
+headers = char(' ',headers);
+lh = size(deblank(M_.endo_names(SubsetOfVariables,:)),2)+2;
+if options_.TeX
+    labels_TeX = deblank(M_.endo_names_tex(SubsetOfVariables,:));
+    lh = size(labels_TeX,2)+2;
+end
+ 
+vardec_i = zeros(length(SubsetOfVariables),shock_number);
 
 for i=1:length(Steps)
     disp(['Period ' int2str(Steps(i)) ':'])
-
+    for j=1:shock_number
-    for j=1:M_.exo_nbr
         vardec_i(:,j) = 100*conditional_decomposition_array(:, ...
                                                           i,j);
     end
-    headers = M_.exo_names;
-    headers(M_.exo_names_orig_ord,:) = headers;
-    headers = char(' ',headers);
-    lh = size(deblank(M_.endo_names(SubsetOfVariables,:)),2)+2;
     dyntable(options_,'',headers,...
              deblank(M_.endo_names(SubsetOfVariables,:)),...
              vardec_i,lh,8,2);
     if options_.TeX
-        labels_TeX = deblank(M_.endo_names_tex(SubsetOfVariables,:));
-        headers_TeX=char('',deblank(M_.exo_names_tex));
-        lh = size(labels_TeX,2)+2;
         dyn_latex_table(M_,options_,[title,'; Period ' int2str(Steps(i))],['th_var_decomp_cond_h',int2str(Steps(i))],headers_TeX,labels_TeX,vardec_i,lh,8,2);
     end
 end

matlab/distributions/inverse_gamma_specification.m

@@ -14,8 +14,10 @@ function [s,nu] = inverse_gamma_specification(mu, sigma2, lb, type, use_fzero_fl
 % - s                [double]    scalar, first hyperparameter.
 % - nu               [double]    scalar, second hyperparameter.
 %
-% REMARK
+% REMARKS
-% The call to the matlab's implementation of the secant method is here for testing purpose and should not be used. This routine fails
+% 1. In the Inverse Gamma parameterization with alpha and beta, we have alpha=nu/2 and beta=2/s, where
+%       if X is IG(alpha,beta) then 1/X is Gamma(alpha,1/beta)
+% 2. The call to the matlab's implementation of the secant method is here for testing purpose and should not be used. This routine fails
 % more often in finding an interval for nu containing a signe change because it expands the interval on both sides and eventually
 % violates  the condition nu>2.
 
@@ -60,7 +62,7 @@ if nargin==4 || isempty(use_fzero_flag)
     use_fzero_flag = false;
 else
     if ~isscalar(use_fzero_flag) || ~islogical(use_fzero_flag)
-        error('Fourth input argument must be a scalar logical!')
+        error('Fifth input argument must be a scalar logical!')
     end
 end
 

matlab/dsge_simulated_theoretical_conditional_variance_decomposition.m

@@ -53,8 +53,10 @@ end
 %delete old stale files before creating new ones
 if posterior
     delete_stale_file([M_.dname '/metropolis/' M_.fname '_PosteriorConditionalVarianceDecomposition*'])
+    delete_stale_file([M_.dname '/metropolis/' M_.fname '_PosteriorConditionalVarianceDecompositionME*'])
 else
     delete_stale_file([M_.dname '/prior/moments/' M_.fname '_PriorConditionalVarianceDecomposition*'])
+    delete_stale_file([M_.dname '/prior/moments/' M_.fname '_PriorConditionalVarianceDecompositionME*'])
 end
 
 % Set varlist (vartan)
@@ -80,6 +82,17 @@ NumberOfDrawsFiles = rows(DrawsFiles);
 NumberOfSavedElementsPerSimulation = nvar*M_.exo_nbr*length(Steps);
 MaXNumberOfConditionalDecompLines = ceil(options_.MaxNumberOfBytes/NumberOfSavedElementsPerSimulation/8);
 
+ME_present=0;
+if ~all(M_.H==0)
+    [observable_pos_requested_vars,index_subset,index_observables]=intersect(ivar,options_.varobs_id,'stable');
+    if ~isempty(observable_pos_requested_vars)
+        ME_present=1;
+        nobs_ME=length(observable_pos_requested_vars);
+        NumberOfSavedElementsPerSimulation_ME = nobs_ME*(M_.exo_nbr+1)*length(Steps);
+        MaXNumberOfConditionalDecompLines_ME = ceil(options_.MaxNumberOfBytes/NumberOfSavedElementsPerSimulation_ME/8);
+    end
+end
+
 if SampleSize<=MaXNumberOfConditionalDecompLines
     Conditional_decomposition_array = zeros(nvar,length(Steps),M_.exo_nbr,SampleSize);
     NumberOfConditionalDecompFiles = 1;
@@ -89,15 +102,30 @@ else
     NumberOfConditionalDecompFiles = ceil(SampleSize/MaXNumberOfConditionalDecompLines);
 end
 
+if ME_present
+    if SampleSize<=MaXNumberOfConditionalDecompLines_ME
+        Conditional_decomposition_array_ME = zeros(nobs_ME,length(Steps),M_.exo_nbr+1,SampleSize);
+        NumberOfConditionalDecompFiles_ME = 1;
+    else
+        Conditional_decomposition_array_ME = zeros(nobs_ME,length(Steps),M_.exo_nbr+1,SampleSize);
+        NumberOfLinesInTheLastConditionalDecompFile_ME = mod(SampleSize,MaXNumberOfConditionalDecompLines_ME);
+        NumberOfConditionalDecompFiles_ME = ceil(SampleSize/MaXNumberOfConditionalDecompLines_ME);
+    end    
+    NumberOfConditionalDecompLines_ME = size(Conditional_decomposition_array_ME,4);
+    ConditionalDecompFileNumber_ME = 0;
+end
+
 NumberOfConditionalDecompLines = size(Conditional_decomposition_array,4);
 ConditionalDecompFileNumber = 0;
 
+
 StateSpaceModel.number_of_state_equations = M_.endo_nbr;
 StateSpaceModel.number_of_state_innovations = M_.exo_nbr;
 
 first_call = 1;
 
 linea = 0;
+linea_ME = 0;
 for file = 1:NumberOfDrawsFiles
     if posterior
         load([M_.dname '/metropolis/' DrawsFiles(file).name ]);
@@ -108,6 +136,7 @@ for file = 1:NumberOfDrawsFiles
     NumberOfDraws = rows(pdraws);
     for linee = 1:NumberOfDraws
         linea = linea+1;
+        linea_ME = linea_ME+1;
         if isdrsaved
             M_=set_parameters_locally(M_,pdraws{linee,1});% Needed to update the covariance matrix of the state innovations.
             dr = pdraws{linee,2};
@@ -125,13 +154,19 @@ for file = 1:NumberOfDrawsFiles
             StateSpaceModel.number_of_state_innovations = M_.exo_nbr;
             StateSpaceModel.sigma_e_is_diagonal = M_.sigma_e_is_diagonal;
             StateSpaceModel.order_var = dr.order_var;
+            StateSpaceModel.observable_pos=options_.varobs_id;
             first_call = 0;
             clear('endo_nbr','nstatic','nspred','k');
         end
         [StateSpaceModel.transition_matrix,StateSpaceModel.impulse_matrix] = kalman_transition_matrix(dr,iv,ic,M_.exo_nbr);
         StateSpaceModel.state_innovations_covariance_matrix = M_.Sigma_e;
+        StateSpaceModel.measurement_error=M_.H;
         clear('dr');
-        Conditional_decomposition_array(:,:,:,linea) = conditional_variance_decomposition(StateSpaceModel, Steps, ivar);
+        [ConditionalVarianceDecomposition, ConditionalVarianceDecomposition_ME]=conditional_variance_decomposition(StateSpaceModel, Steps, ivar);
+        Conditional_decomposition_array(:,:,:,linea) =ConditionalVarianceDecomposition;
+        if ME_present
+            Conditional_decomposition_array_ME(:,:,:,linea) =ConditionalVarianceDecomposition_ME;
+        end
         if linea == NumberOfConditionalDecompLines
             ConditionalDecompFileNumber = ConditionalDecompFileNumber + 1;
             linea = 0;
@@ -151,6 +186,28 @@ for file = 1:NumberOfDrawsFiles
                 clear('Conditional_decomposition_array');
             end
         end
+        %with measurement error
+        if ME_present
+            if linea_ME == NumberOfConditionalDecompLines_ME
+                ConditionalDecompFileNumber_ME = ConditionalDecompFileNumber_ME + 1;
+                linea_ME = 0;
+                if posterior
+                    save([M_.dname '/metropolis/' M_.fname '_PosteriorConditionalVarianceDecompME' int2str(ConditionalDecompFileNumber_ME) '.mat' ], ...
+                        'Conditional_decomposition_array_ME');
+                else
+                    save([M_.dname '/prior/moments/' M_.fname '_PriorConditionalVarianceDecompME' int2str(ConditionalDecompFileNumber_ME) '.mat' ], ...
+                        'Conditional_decomposition_array_ME');
+                end
+                if (ConditionalDecompFileNumber_ME==NumberOfConditionalDecompFiles_ME-1)% Prepare last round.
+                    Conditional_decomposition_array_ME = zeros(nobs_ME, length(Steps),M_.exo_nbr+1,NumberOfLinesInTheLastConditionalDecompFile_ME) ;
+                    NumberOfConditionalDecompLines_ME = NumberOfLinesInTheLastConditionalDecompFile_ME;
+                elseif ConditionalDecompFileNumber_ME<NumberOfConditionalDecompFiles_ME-1
+                    Conditional_decomposition_array_ME = zeros(nobs_ME,length(Steps),M_.exo_nbr+1,MaXNumberOfConditionalDecompLines_ME);
+                else
+                    clear('Conditional_decomposition_array_ME');
+                end
+            end
+        end
     end
 end
 

matlab/dsge_simulated_theoretical_variance_decomposition.m

@@ -53,9 +53,10 @@ end
 %delete old stale files before creating new ones
 if posterior
     delete_stale_file([M_.dname '/metropolis/' M_.fname '_PosteriorVarianceDecomposition*']);
+    delete_stale_file([M_.dname '/metropolis/' M_.fname '_PosteriorVarianceDecompME*']);
 else
     delete_stale_file([M_.dname '/prior/moments/' M_.fname '_PriorVarianceDecomposition*']);
-end
+    delete_stale_file([M_.dname '/prior/moments/' M_.fname '_PriorVarianceDecompME*']);end
 
 % Set varlist (vartan)
 if ~posterior
@@ -84,6 +85,17 @@ NumberOfDrawsFiles = rows(DrawsFiles);
 NumberOfSavedElementsPerSimulation = nvar*(nexo+1);
 MaXNumberOfDecompLines = ceil(options_.MaxNumberOfBytes/NumberOfSavedElementsPerSimulation/8);
 
+ME_present=0;
+if ~all(M_.H==0)
+    [observable_pos_requested_vars,index_subset,index_observables]=intersect(ivar,options_.varobs_id,'stable');
+    if ~isempty(observable_pos_requested_vars)
+        ME_present=1;
+        nobs_ME=length(observable_pos_requested_vars);
+        NumberOfSavedElementsPerSimulation_ME = nobs_ME*(nexo+1);
+        MaXNumberOfDecompLines_ME = ceil(options_.MaxNumberOfBytes/NumberOfSavedElementsPerSimulation_ME/8);
+    end
+end
+
 if SampleSize<=MaXNumberOfDecompLines
     Decomposition_array = zeros(SampleSize,nvar*nexo);
     NumberOfDecompFiles = 1;
@@ -96,9 +108,22 @@ end
 NumberOfDecompLines = rows(Decomposition_array);
 DecompFileNumber = 1;
 
+if ME_present
+    if SampleSize<=MaXNumberOfDecompLines_ME
+        Decomposition_array_ME = zeros(SampleSize,nobs_ME*(nexo+1));
+        NumberOfDecompFiles_ME = 1;
+    else
+        Decomposition_array_ME = zeros(MaXNumberOfDecompLines_ME,nobs_ME*(nexo+1));
+        NumberOfLinesInTheLastDecompFile_ME = mod(SampleSize,MaXNumberOfDecompLines_ME);
+        NumberOfDecompFiles_ME = ceil(SampleSize/MaXNumberOfDecompLines_ME);
+    end
+    NumberOfDecompLines_ME = rows(Decomposition_array_ME);
+    DecompFileNumber_ME = 1;
+end
 % Compute total variances (covariances are not saved) and variances
 % implied by each structural shock.
 linea = 0;
+linea_ME = 0;
 only_non_stationary_vars=0;
 for file = 1:NumberOfDrawsFiles
     if posterior
@@ -110,6 +135,7 @@ for file = 1:NumberOfDrawsFiles
     NumberOfDraws = rows(pdraws);
     for linee = 1:NumberOfDraws
         linea = linea+1;
+        linea_ME = linea_ME+1;
         if isdrsaved
             dr = pdraws{linee,2};
         else
@@ -125,11 +151,7 @@ for file = 1:NumberOfDrawsFiles
             end
         end
         if only_non_stationary_vars
-            for i=1:nvar
+             Decomposition_array(linea,:) = NaN;
-                for j=1:nexo
-                    Decomposition_array(linea,(i-1)*nexo+j) = NaN;
-                end
-            end
         else
             tmp = th_autocovariances(dr,ivar,M_,options_,nodecomposition);
             for i=1:nvar
@@ -137,6 +159,17 @@ for file = 1:NumberOfDrawsFiles
                     Decomposition_array(linea,(i-1)*nexo+j) = tmp{2}(i,j);
                 end
             end
+            if ME_present
+                ME_Variance=diag(M_.H);
+                tmp_ME=NaN(nobs_ME,nexo+1);
+                tmp_ME(:,1:end-1)=tmp{2}(index_subset,:).*repmat(diag(tmp{1}(index_subset,index_subset))./(diag(tmp{1}(index_subset,index_subset))+ME_Variance(index_observables)),1,nexo);
+                tmp_ME(:,end)=1-sum(tmp_ME(:,1:end-1),2);
+                for i=1:nobs_ME
+                    for j=1:nexo+1
+                        Decomposition_array_ME(linea,(i-1)*(nexo+1)+j) = tmp_ME(i,j);
+                    end
+                end
+            end
         end
         if linea == NumberOfDecompLines
             if posterior
@@ -156,6 +189,26 @@ for file = 1:NumberOfDrawsFiles
                 clear('Decomposition_array');
             end
         end
+        if ME_present
+            if linea_ME == NumberOfDecompLines_ME
+                if posterior
+                    save([M_.dname '/metropolis/' M_.fname '_PosteriorVarianceDecompME' int2str(DecompFileNumber_ME) '.mat' ],'Decomposition_array_ME');
+                else
+                    save([M_.dname '/prior/moments/' M_.fname '_PriorVarianceDecompME' int2str(DecompFileNumber_ME) '.mat' ],'Decomposition_array_ME');
+                end
+                DecompFileNumber_ME = DecompFileNumber_ME + 1;
+                linea_ME = 0;
+                test = DecompFileNumber_ME-NumberOfDecompFiles_ME;
+                if ~test% Prepare the last round...
+                    Decomposition_array_ME = zeros(NumberOfLinesInTheLastDecompFile_ME,nobs_ME*(nexo+1));
+                    NumberOfDecompLines_ME = NumberOfLinesInTheLastDecompFile_ME;
+                elseif test<0
+                    Decomposition_array_ME = zeros(MaXNumberOfDecompLines_ME,nobs_ME*(nexo+1));
+                else
+                    clear('Decomposition_array_ME');
+                end
+            end
+        end
     end
 end
 

matlab/dynare.m

@@ -59,7 +59,6 @@ if nargin>1
     id = strfind(varargin, 'nopreprocessoroutput');
     if ~all(cellfun(@isempty, id))
         preprocessoroutput = false;
-        varargin(cellfun(@isempty, id) == 0) = [];
     end
 end
 

matlab/dynare_config.m

@@ -125,6 +125,11 @@ if ~exist('struct2array')
     p{end+1} = '/missing/struct2array';
 end
 
+% isfile is missing in Octave and Matlab<R2017b
+if isoctave || matlab_ver_less_than('9.3')
+    p{end+1} = '/missing/isfile';
+end
+
 P = cellfun(@(c)[dynareroot(1:end-1) c], p, 'uni',false);
 
 % Get mex files folder(s)

matlab/dynare_estimation_1.m

@@ -235,7 +235,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
                 [junk1, junk2,junk3, junk4, hh] = feval(objective_function,xparam1, ...
                                                         dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
                 options_.analytic_derivation = ana_deriv_old;
-            elseif ~isnumeric(options_.mode_compute) || ~(isequal(options_.mode_compute,5) && newratflag~=1)
+            elseif ~isnumeric(options_.mode_compute) || ~(isequal(options_.mode_compute,5) && newratflag~=1 && strcmp(func2str(objective_function),'dsge_likelihood'))
                 % with flag==0, we force to use the hessian from outer product gradient of optimizer 5
                 if options_.hessian.use_penalized_objective
                     penalized_objective_function = str2func('penalty_objective_function');
@@ -519,7 +519,7 @@ if (any(bayestopt_.pshape  >0 ) && options_.mh_replic) || ...
             end
             prior_posterior_statistics('posterior',dataset_,dataset_info);
         end
-        xparam1 = get_posterior_parameters('mean');
+        xparam1 = get_posterior_parameters('mean',M_,estim_params_,oo_,options_);
         M_ = set_all_parameters(xparam1,estim_params_,M_);
     end
 end

matlab/dynare_estimation_init.m

@@ -53,15 +53,22 @@ hh = [];
 xparam1 = [];
 
 if isempty(gsa_flag)
-    gsa_flag = 0;
+    gsa_flag = false;
 else
     % Decide if a DSGE or DSGE-VAR has to be estimated.
     if ~isempty(strmatch('dsge_prior_weight',M_.param_names))
         options_.dsge_var = 1;
     end
-    % Get the list of the endogenous variables for which posterior statistics wil be computed
+    if isempty(var_list_)
-    var_list_ = check_list_of_variables(options_, M_, var_list_);
+        var_list_ = check_list_of_variables(options_, M_, var_list_);
-    options_.varlist = var_list_;
+        options_.varlist = var_list_;
+    end
+    if gsa_flag
+        % Get the list of the endogenous variables for which posterior statistics wil be computed.
+        options_.varlist = var_list_;
+    else
+        % This was done in dynare_estimation_1
+    end
 end
 
 if options_.dsge_var && options_.presample~=0
@@ -549,7 +556,7 @@ else
     steadystate_check_flag = 1;
 end
 
-% If the steady state of the observed variables is non zero, set noconstant equal 0 ()
+%check steady state at initial parameters
 M = M_;
 nvx = estim_params_.nvx;
 ncx = estim_params_.ncx;
@@ -565,6 +572,7 @@ if info(1)
     print_info(info, 0, options_);
 end
 
+% If the steady state of the observed variables is non zero, set noconstant equal 0 ()
 if (~options_.loglinear && all(abs(oo_.steady_state(bayestopt_.mfys))<1e-9)) || (options_.loglinear && all(abs(log(oo_.steady_state(bayestopt_.mfys)))<1e-9))
     options_.noconstant = 1;
 else

matlab/dynare_identification.m

@@ -271,15 +271,15 @@ if iload <=0
               case 'posterior_mode'
                 parameters_TeX = 'Posterior mode';
                 disp('Testing posterior mode')
-                params(1,:) = get_posterior_parameters('mode');
+                params(1,:) = get_posterior_parameters('mode',M_,estim_params_,oo_,options_);
               case 'posterior_mean'
                 parameters_TeX = 'Posterior mean';
                 disp('Testing posterior mean')
-                params(1,:) = get_posterior_parameters('mean');
+                params(1,:) = get_posterior_parameters('mean',M_,estim_params_,oo_,options_);
               case 'posterior_median'
                 parameters_TeX = 'Posterior median';
                 disp('Testing posterior median')
-                params(1,:) = get_posterior_parameters('median');
+                params(1,:) = get_posterior_parameters('median',M_,estim_params_,oo_,options_);
               case 'prior_mode'
                 parameters_TeX = 'Prior mode';
                 disp('Testing prior mode')

matlab/evaluate_likelihood.m

@@ -1,10 +1,14 @@
-function [llik,parameters] = evaluate_likelihood(parameters)
+function [llik,parameters] = evaluate_likelihood(parameters,M_,estim_params_,oo_,options_,bayestopt_)
 % Evaluate the logged likelihood at parameters.
 %
 % INPUTS
 %    o parameters  a string ('posterior mode','posterior mean','posterior median','prior mode','prior mean') or a vector of values for
 %                  the (estimated) parameters of the model.
-%
+%    o M_          [structure]  Definition of the model
+%    o estim_params_ [structure] characterizing parameters to be estimated
+%    o oo_         [structure]  Storage of results
+%    o options_    [structure]  Options
+%    o bayestopt_  [structure]  describing the priors
 %
 % OUTPUTS
 %    o ldens       [double]  value of the sample logged density at parameters.
@@ -35,8 +39,6 @@ function [llik,parameters] = evaluate_likelihood(parameters)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-global options_ M_ bayestopt_ oo_ estim_params_
-
 persistent dataset dataset_info
 
 if nargin==0
@@ -46,11 +48,11 @@ end
 if ischar(parameters)
     switch parameters
       case 'posterior mode'
-        parameters = get_posterior_parameters('mode');
+        parameters = get_posterior_parameters('mode',M_,estim_params_,oo_,options_);
       case 'posterior mean'
-        parameters = get_posterior_parameters('mean');
+        parameters = get_posterior_parameters('mean',M_,estim_params_,oo_,options_);
       case 'posterior median'
-        parameters = get_posterior_parameters('median');
+        parameters = get_posterior_parameters('median',M_,estim_params_,oo_,options_);
       case 'prior mode'
         parameters = bayestopt_.p5(:);
       case 'prior mean'
@@ -72,5 +74,5 @@ end
 options_=select_qz_criterium_value(options_);
 
 llik = -dsge_likelihood(parameters,dataset,dataset_info,options_,M_,estim_params_,bayestopt_,prior_bounds(bayestopt_,options_.prior_trunc),oo_);
-ldens = evaluate_prior(parameters);
+ldens = evaluate_prior(parameters,M_,estim_params_,oo_,options_,bayestopt_);
 llik = llik - ldens;

matlab/evaluate_posterior_kernel.m

@@ -1,10 +1,16 @@
-function lpkern = evaluate_posterior_kernel(parameters,llik)
+function lpkern = evaluate_posterior_kernel(parameters,M_,estim_params_,oo_,options_,bayestopt_,llik)
-% Evaluate the prior density at parameters.
+% Evaluate the evaluate_posterior_kernel at parameters.
 %
 % INPUTS
 %    o parameters  a string ('posterior mode','posterior mean','posterior median','prior mode','prior mean') or a vector of values for
 %                  the (estimated) parameters of the model.
-%
+%    o M_         [structure]     Definition of the model
+%    o estim_params_ [structure]  characterizing parameters to be estimated
+%    o oo_         [structure]    Storage of results
+%    o options_    [structure]    Options
+%    o bayestopt_  [structure]    describing the priors
+%    o llik        [double]       value of the logged likelihood if it
+%                                   should not be computed
 %
 % OUTPUTS
 %    o lpkern      [double]  value of the logged posterior kernel.
@@ -34,8 +40,8 @@ function lpkern = evaluate_posterior_kernel(parameters,llik)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-[ldens,parameters] = evaluate_prior(parameters);
+[ldens,parameters] = evaluate_prior(parameters,M_,estim_params_,oo_,options_,bayestopt_);
-if nargin==1
+if nargin==6 %llik provided as an input
-    llik = evaluate_likelihood(parameters);
+    llik = evaluate_likelihood(parameters,M_,estim_params_,oo_,options_,bayestopt_);
 end
 lpkern = ldens+llik;

matlab/evaluate_prior.m

@@ -1,9 +1,14 @@
-function [ldens,parameters] = evaluate_prior(parameters)
+function [ldens,parameters] = evaluate_prior(parameters,M_,estim_params_,oo_,options_,bayestopt_)
 % Evaluate the prior density at parameters.
 %
 % INPUTS
 %    o parameters  a string ('posterior mode','posterior mean','posterior median','prior mode','prior mean') or a vector of values for
 %                  the (estimated) parameters of the model.
+%    o M_          [structure]  Definition of the model
+%    o oo_         [structure]  Storage of results
+%    o options_    [structure]  Options
+%    o bayestopt_  [structure]  describing the priors
+%    o estim_params_ [structure] characterizing parameters to be estimated
 %
 %
 % OUTPUTS
@@ -33,8 +38,6 @@ function [ldens,parameters] = evaluate_prior(parameters)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-global bayestopt_
-
 if nargin==0
     parameters = 'posterior mode';
 end
@@ -42,11 +45,11 @@ end
 if ischar(parameters)
     switch parameters
       case 'posterior mode'
-        parameters = get_posterior_parameters('mode');
+        parameters = get_posterior_parameters('mode',M_,estim_params_,oo_,options_);
       case 'posterior mean'
-        parameters = get_posterior_parameters('mean');
+        parameters = get_posterior_parameters('mean',M_,estim_params_,oo_,options_);
       case 'posterior median'
-        parameters = get_posterior_parameters('median');
+        parameters = get_posterior_parameters('median',M_,estim_params_,oo_,options_);
       case 'prior mode'
         parameters = bayestopt_.p5(:);
       case 'prior mean'

matlab/evaluate_smoother.m

@@ -73,13 +73,13 @@ end
 if ischar(parameters)
     switch parameters
       case 'posterior_mode'
-        parameters = get_posterior_parameters('mode');
+        parameters = get_posterior_parameters('mode',M_,estim_params_,oo_,options_);
       case 'posterior_mean'
-        parameters = get_posterior_parameters('mean');
+        parameters = get_posterior_parameters('mean',M_,estim_params_,oo_,options_);
       case 'posterior_median'
-        parameters = get_posterior_parameters('median');
+        parameters = get_posterior_parameters('median',M_,estim_params_,oo_,options_);
       case 'mle_mode'
-        parameters = get_posterior_parameters('mode','mle_');
+        parameters = get_posterior_parameters('mode',M_,estim_params_,oo_,options_,'mle_');
       case 'prior_mode'
         parameters = bayestopt_.p5(:);
       case 'prior_mean'

matlab/execute_prior_posterior_function.m

@@ -62,17 +62,24 @@ if strcmpi(type,'posterior')
     n_draws=options_.sub_draws;
     prior = false;
 elseif strcmpi(type,'prior')
+    if isempty(bayestopt_)
+        if ~isempty(estim_params_) && ~(isfield(estim_params_,'nvx') && (size(estim_params_.var_exo,1)+size(estim_params_.var_endo,1)+size(estim_params_.corrx,1)+size(estim_params_.corrn,1)+size(estim_params_.param_vals,1))==0)
+            [xparam1,estim_params_,bayestopt_,lb,ub,M_] = set_prior(estim_params_,M_,options_);
+        else
+            error('The prior distributions are not properly set up.')
+        end
+    end
     prior_draw(bayestopt_, options_.prior_trunc);
 else
     error('EXECUTE_POSTERIOR_FUNCTION: Unknown type!')
 end
 
 %get draws for later use
-first_draw=GetOneDraw(type);
+first_draw=GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_);
 parameter_mat=NaN(n_draws,length(first_draw));
 parameter_mat(1,:)=first_draw;
 for draw_iter=2:n_draws
-    parameter_mat(draw_iter,:) = GetOneDraw(type);
+    parameter_mat(draw_iter,:) = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_);
 end
 
 % get output size

matlab/fjaco.m

@@ -1,14 +1,14 @@
 function fjac = fjaco(f,x,varargin)
 
-% FDJAC Computes two-sided finite difference Jacobian
+% FJACO Computes two-sided finite difference Jacobian
 % USAGE
-%   fjac = fdjac(f,x,P1,P2,...)
+%   fjac = fjaco(f,x,P1,P2,...)
 % INPUTS
 %   f         : name of function of form fval = f(x)
 %   x         : evaluation point
 %   P1,P2,... : additional arguments for f (optional)
 % OUTPUT
-%   fjac      : finite differnce Jacobian
+%   fjac      : finite difference Jacobian
 %
 % Copyright (C) 2010-2017 Dynare Team
 %

matlab/get_all_parameters.m

@@ -6,7 +6,7 @@ function xparam1=get_all_parameters(estim_params_,M_)
 % parameter values
 %
 % INPUTS
-%        estim_params_:  Dynare structure describing the estimated parameters.
+%    estim_params_:  Dynare structure describing the estimated parameters.
 %    M_:             Dynare structure describing the model.
 %
 % OUTPUTS

matlab/get_posterior_parameters.m

@@ -1,11 +1,13 @@
-function xparam = get_posterior_parameters(type,field1)
+function xparam = get_posterior_parameters(type,M_,estim_params_,oo_,options_,field1)
 
-% function xparam = get_posterior_parameters(type)
+% function xparam = get_posterior_parameters(type,M_,estim_params_,oo_,options_,field1)
 % Selects (estimated) parameters (posterior mode or posterior mean).
 %
 % INPUTS
-%   o type       [char]     = 'mode' or 'mean'.
+%   o type              [char]     = 'mode' or 'mean'.
-%   o field_1    [char]     optional field like 'mle_'.
+%   o M_:               [structure] Dynare structure describing the model.
+%   o estim_params_:    [structure] Dynare structure describing the estimated parameters.
+%   o field_1           [char]     optional field like 'mle_'.
 %
 % OUTPUTS
 %   o xparam     vector of estimated parameters
@@ -30,9 +32,7 @@ function xparam = get_posterior_parameters(type,field1)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-global estim_params_ oo_ options_ M_
+if nargin<6
-
-if nargin<2
     field1='posterior_';
 end
 nvx = estim_params_.nvx;
@@ -48,7 +48,6 @@ for i=1:nvx
     k1 = estim_params_.var_exo(i,1);
     name1 = deblank(M_.exo_names(k1,:));
     xparam(m) = oo_.([field1 type]).shocks_std.(name1);
-    M_.Sigma_e(k1,k1) = xparam(m)^2;
     m = m+1;
 end
 
@@ -65,8 +64,6 @@ for i=1:ncx
     name1 = deblank(M_.exo_names(k1,:));
     name2 = deblank(M_.exo_names(k2,:));
     xparam(m) = oo_.([field1 type]).shocks_corr.([name1 '_' name2]);
-    M_.Sigma_e(k1,k2) = xparam(m);
-    M_.Sigma_e(k2,k1) = xparam(m);
     m = m+1;
 end
 
@@ -84,10 +81,5 @@ FirstDeep = m;
 for i=1:np
     name1 = deblank(M_.param_names(estim_params_.param_vals(i,1),:));
     xparam(m) = oo_.([field1 type]).parameters.(name1);
-    %assignin('base',name1,xparam(m));% Useless with version 4 (except maybe for users)
     m = m+1;
-end
-
-if np
-    M_.params(estim_params_.param_vals(:,1)) = xparam(FirstDeep:end);
 end

matlab/global_initialization.m

@@ -174,6 +174,11 @@ options_.bandpass.indicator = 0;
 options_.bandpass.passband = [6; 32];
 options_.bandpass.K=12;
 
+options_.irf_opt.diagonal_only = 0;
+options_.irf_opt.stderr_multiples = 0;
+options_.irf_opt.irf_shock_graphtitles = {};
+options_.irf_opt.irf_shocks = [];
+
 % Extended path options
 %
 % Set debug flag
@@ -825,6 +830,8 @@ options_.mcppath.mu0 = [];
 %Figure options
 options_.figures.textwidth=0.8;
 
+options_.varobs_id=[]; %initialize field
+
 % initialize persistent variables in priordens()
 priordens([],[],[],[],[],[],1);
 % initialize persistent variables in dyn_first_order_solver()

matlab/imcforecast.m

@@ -80,16 +80,16 @@ if estimated_model
     if ischar(options_cond_fcst.parameter_set)
         switch options_cond_fcst.parameter_set
           case 'posterior_mode'
-            xparam = get_posterior_parameters('mode');
+            xparam = get_posterior_parameters('mode',M_,estim_params_,oo_,options_);
             graph_title='Posterior Mode';
           case 'posterior_mean'
-            xparam = get_posterior_parameters('mean');
+            xparam = get_posterior_parameters('mean',M_,estim_params_,oo_,options_);
             graph_title='Posterior Mean';
           case 'posterior_median'
-            xparam = get_posterior_parameters('median');
+            xparam = get_posterior_parameters('median',M_,estim_params_,oo_,options_);
             graph_title='Posterior Median';
           case 'mle_mode'
-            xparam = get_posterior_parameters('mode','mle_');
+            xparam = get_posterior_parameters('mode',M_,estim_params_,oo_,options_,'mle_');
             graph_title='ML Mode';
           case 'prior_mode'
             xparam = bayestopt_.p5(:);

matlab/lmmcp/get_complementarity_conditions.m

@@ -44,7 +44,6 @@ if ramsey_policy
               otherwise
                 error('Wrong operator in get_complementarity_conditions')
             end
-            eq_index(i) = 1;
         end
     end
 end
@@ -52,6 +51,10 @@ end
 etags = M.equations_tags;
 for i=1:size(etags,1)
     if strcmp(etags{i,2},'mcp')
+        eq_nbr = etags{i,1};
+        if ramsey_policy
+            eq_nbr = eq_nbr + M.ramsey_eq_nbr;
+        end
         str = etags{i,3};
         kop = strfind(etags{i,3},'<');
         if ~isempty(kop)
@@ -61,8 +64,8 @@ for i=1:size(etags,1)
                                'not recognized'],etags{i,3},strtrim(str(1:kop-1))))
             end
             ub(k) = str2num(str(kop+1:end));
-            eq_index(etags{i,1}) = k;
+            eq_index(eq_nbr) = k;
-            eq_index(k) = etags{i,1};
+            eq_index(k) = eq_nbr;
         else
             kop = strfind(etags{i,3},'>');
             if ~isempty(kop)
@@ -72,8 +75,8 @@ for i=1:size(etags,1)
                                    'not recognized'],etags{i,3},strtrim(str(1:kop-1))))
                 end
                 lb(k) = str2num(str(kop+1:end));
-                eq_index(etags{i,1}) = k;
+                eq_index(eq_nbr) = k;
-                eq_index(k) = etags{i,1};
+                eq_index(k) = eq_nbr;
             else
                 error(sprintf(['Complementarity condition %s can''t be ' ...
                                'parsed'],etags{i,3}))

matlab/lyapunov_solver.m

@@ -41,7 +41,7 @@ function P=lyapunov_solver(T,R,Q,DynareOptions) % --*-- Unitary tests --*--
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
 if DynareOptions.lyapunov_fp == 1
-    P = lyapunov_symm(T,R*Q'*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, 3, DynareOptions.debug);
+    P = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, 3, DynareOptions.debug);
 elseif DynareOptions.lyapunov_db == 1
     [P, errorflag] = disclyap_fast(T,R*Q*R',DynareOptions.lyapunov_doubling_tol);
     if errorflag %use Schur-based method
@@ -183,4 +183,4 @@ end
 %$ end
 %$
 %$ T = all(t);
-%@eof:1
+%@eof:1

matlab/missing/isfile/isfile.m

@@ -27,7 +27,7 @@ function a = isfile(b)
 %! @end deftypefn
 %@eod:
 
-% Copyright (C) 2012 Dynare Team
+% Copyright (C) 2012-2017 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -44,18 +44,45 @@ function a = isfile(b)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-% Original author: stephane DOT adjemian AT univ DASH lemans DOT fr
+stringarrayflag = false;
+cellofstringflag = false;
+n = 1;
+a = false;
 
-[base,ext] = strtok(b,'.');
+if isstring(b) && length(b)>1 && isvector(b)
+    n = length(b);
+    stringarrayflag = true;
+    a = false(size(b));
+end
 
-if isempty(ext)
+if iscell(b) && length(b)>1 && isvector(b)
-    % File has no extension.
+    if all(cellfun(@ischar, b))
-    [status, c] = fileattrib(b);
+        n = length(b);
-    if status
+        cellofstringflag = true;
-        a = ~c.directory;
+        a = false(size(b));
     else
-        a = 0;
+        error('Wrong input argument type!')
+    end
+end
+
+for i=1:n
+    if stringarrayflag
+        d = b(i);
+    elseif cellofstringflag
+        d = b{i};
+    elseif ischar(b) && size(b, 1)==1 
+        d = b;
+    else
+        error('Wrong input argument type!')
+    end
+    [base, ext] = strtok(d, '.');
+    if isempty(ext)
+        % File has no extension.
+        [status, c] = fileattrib(d);
+        if status
+            a(i) = ~c.directory;
+        end
+    else
+        a(i) = isequal(exist(d, 'file'), 2);
     end
-else
-    a = isequal(exist(b,'file'),2);
 end

matlab/mode_check.m

@@ -111,8 +111,15 @@ for plt = 1:nbplt
             end
         end
         xx = x;
-        l1 = max(BoundsInfo.lb(kk),(1-sign(x(kk))*ll)*x(kk)); m1 = 0; %lower bound
+        if x(kk)~=0
-        l2 = min(BoundsInfo.ub(kk),(1+sign(x(kk))*ll)*x(kk)); %upper bound
+            l1 = max(BoundsInfo.lb(kk),(1-sign(x(kk))*ll)*x(kk)); m1 = 0; %lower bound
+            l2 = min(BoundsInfo.ub(kk),(1+sign(x(kk))*ll)*x(kk)); %upper bound
+        else
+            %size info for 0 parameter is missing, use prior standard
+            %deviation
+            l1 = max(BoundsInfo.lb(kk),-BayesInfo.p2(kk)); m1 = 0; %lower bound
+            l2 = min(BoundsInfo.ub(kk),BayesInfo.p2(kk)); %upper bound
+        end            
         binding_lower_bound=0;
         binding_upper_bound=0;
         if isequal(x(kk),BoundsInfo.lb(kk))

matlab/occbin/evaluate_model.m

@@ -0,0 +1,19 @@
+function [r,g1,g2,g3] = evaluate_model(z,x,M,ss)
+
+ll = M.lead_lag_incidence';
+y = z(find(ll(:)));
+
+switch nargout
+  case 1
+    r = feval([M.fname '_dynamic'],y,x, ...
+              M.params, ss, 1);
+  case 2
+    [r,g1] = feval([M.fname '_dynamic'],y,x, ...
+                    M.params, ss, 1);
+  case 3
+    [r,g1,g2] = feval([M.fname '_dynamic'],y,x, ...
+                    M.params, ss, 1);
+  case 4
+    [r,g1,g2,g3] = feval([M.fname '_dynamic'],y,x, ...
+                         M.params, ss, 1);
+end

matlab/occbin/get_coef.m

@@ -0,0 +1,26 @@
+function [coef_y,coef_u] = get_coef(jacobian,M)
+
+ll = M.lead_lag_incidence;
+endo_nbr = M.endo_nbr;
+coef_y = zeros(endo_nbr,3*endo_nbr);
+coef_u = zeros(endo_nbr,M.exo_nbr);
+
+if M.maximum_lag > 0
+    [junk,c1,c2] = find(ll(1,:));
+    coef_y(:,c1) = jacobian(:,c2);
+    [junk,c1,c2] = find(ll(2,:));
+    coef_y(:,c1+endo_nbr) = jacobian(:,c2);
+    if M.maximum_lead > 0
+            [junk,c1,c2] = find(ll(3,:));
+            coef_y(:,c1+2*endo_nbr) = jacobian(:,c2);
+    end
+else
+    [junk,c1,c2] = find(ll(1,:));
+    coef_y(:,c1+endo_nbr) = jacobian(:,c2);
+    if M.maximum_lead > 0
+            [junk,c1,c2] = find(ll(2,:));
+            coef_y(:,c1+2*endo_nbr) = jacobian(:,c2);
+    end
+end
+
+coef_u = jacobian(:,max(ll(end,:))+1:end);

matlab/occbin/get_complementarity_conditions.m

@@ -0,0 +1,72 @@
+function [lb,ub,eq_index] = get_complementarity_conditions(M,ramsey_policy)
+
+% Copyright (C) 2014 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/>.
+
+ub = inf(M.endo_nbr,1);
+lb = -ub;
+eq_index = (1:M.endo_nbr)';
+if ramsey_policy
+    if isfield(M,'ramsey_model_constraints')
+        rc = M.ramsey_model_constraints;
+        for i = 1:length(rc)
+            switch rc{i}{2}
+              case {'>','>='}
+                lb(rc{i}{1}) = eval(rc{i}{3});
+              case {'<','<='}
+                ub(rc{i}{1}) = eval(rc{i}{3});
+              otherwise
+                error('Wrong operator in get_complementarity_conditions')
+            end
+            eq_index(i) = 1;
+        end
+    end
+end
+
+etags = M.equations_tags;
+for i=1:size(etags,1)
+    if strcmp(etags{i,2},'mcp')
+        str = etags{i,3};
+        kop = strfind(etags{i,3},'<');
+        if ~isempty(kop)
+                k = find(strcmp(strtrim(str(1:kop-1)),cellstr(M.endo_names)));
+                if isempty(k)
+                    error(sprintf(['Complementarity condition %s: variable %s is ' ...
+                                   'not recognized',etags{i,3},b{1}]))
+                end
+                ub(k) = str2num(str(kop+1:end));
+                eq_index(etags{i,1}) = k;
+                eq_index(k) = etags{i,1};
+        else
+            kop = strfind(etags{i,3},'>');
+            if ~isempty(kop)
+                k = find(strcmp(strtrim(str(1:kop-1)),cellstr(M.endo_names)));
+                if isempty(k)
+                    error(sprintf(['Complementarity condition %s: variable %s is ' ...
+                                   'not recognized',etags{i},b{1}]))
+                end
+                lb(k) = str2num(str(kop+1:end));
+                eq_index(etags{i,1}) = k;
+                eq_index(k) = etags{i,1};
+            else
+                error(sprintf(['Complementarity condition %s can''t be ' ...
+                               'parsed'],etags{i,3}))
+            end
+        end
+    end
+end
+

matlab/occbin/get_occbin_complementarity_conditions.m

@@ -0,0 +1,75 @@
+function [ivar,ieq,lb,ub] = get_occbin_complementarity_conditions(M,ramsey_policy)
+
+% Copyright (C) 2015 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/>.
+
+nrow = 1;
+if ramsey_policy
+    if isfield(M,'ramsey_model_constraints')
+        rc = M.ramsey_model_constraints;
+        for i = 1:length(rc)
+            switch rc{i}{2}
+              case {'>','>='}
+                ivar(nrow) = rc{i}{1};
+                ieq(nrow) = rc{i}{1};
+                lb(nrow) = eval(rc{i}{3});
+              case {'<','<='}
+                ivar(nrow) = rc{i}{1};
+                ieq(nrow) = rc{i}{1};
+                ub(nrow) = eval(rc{i}{3});
+              otherwise
+                error('Wrong operator in get_complementarity_conditions')
+            end
+            nrow = nrow + 1;
+        end
+    end
+end
+
+etags = M.equations_tags;
+for i=1:size(etags,1)
+    if strcmp(etags{i,2},'mcp')
+        str = etags{i,3};
+        kop = strfind(etags{i,3},'<');
+        if ~isempty(kop)
+                k = find(strcmp(strtrim(str(1:kop-1)),cellstr(M.endo_names)));
+                if isempty(k)
+                    error(sprintf(['Complementarity condition %s: variable %s is ' ...
+                                   'not recognized',etags{i,3},b{1}]))
+                end
+                ivar(nrow) = k;
+                ieq(nrow) = etags{i,1};
+                ub(nrow) = eval(str(kop+1:end));
+        else
+            kop = strfind(etags{i,3},'>');
+            if ~isempty(kop)
+                k = find(strcmp(strtrim(str(1:kop-1)),cellstr(M.endo_names)));
+                if isempty(k)
+                    error(sprintf(['Complementarity condition %s: variable %s is ' ...
+                                   'not recognized',etags{i},b{1}]))
+                end
+                ivar(nrow) = k;
+                ieq(nrow) = etags{i,1};
+                lb(k) = eval(str(kop+1:end));
+            else
+                error(sprintf(['Complementarity condition %s can''t be ' ...
+                               'parsed'],etags{i,3}))
+            end
+        end
+        nrow = nrow + 1;
+    end
+end
+

matlab/occbin/get_occbin_constraints.m

@@ -0,0 +1,88 @@
+function [i_base,i_alt,c_base,c_alt] = get_occbin_constraints(M,steady_state,ramsey_policy)
+
+% Copyright (C) 2015 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/>.
+
+nrow = 1;
+if ramsey_policy
+    if isfield(M,'ramsey_model_constraints')
+        rc = M.ramsey_model_constraints;
+        for i = 1:length(rc)
+            switch rc{i}{2}
+              case {'>','>='}
+                ivar(nrow) = rc{i}{1};
+                ieq(nrow) = rc{i}{1};
+                lb(nrow) = eval(rc{i}{3});
+              case {'<','<='}
+                ivar(nrow) = rc{i}{1};
+                ieq(nrow) = rc{i}{1};
+                ub(nrow) = eval(rc{i}{3});
+              otherwise
+                error('Wrong operator in get_complementarity_conditions')
+            end
+            nrow = nrow + 1;
+        end
+    end
+end
+
+i_base = {};
+i_alt = {};
+etags = M.equations_tags;
+m = 1;
+base = true;
+for i=1:size(etags,1)
+    [iv,boundary,operator] = parse_constraint(etags{i,3},M.endo_names,M.params,M.param_names);
+    if strcmp(etags{i,2},'OCCBIN')
+        if base
+            i_alt{m} = 1:M.eq_nbr;
+            i_alt{m}(etags{i,1}) = [];
+            c_base{m,1} = etags{i,1};
+            c_base{m,2} = iv;
+            c_base{m,3} = boundary - steady_state(iv);
+            c_base{m,4} = operator;
+            base = false;
+        else
+            i_base{m} = 1:M.eq_nbr;
+            i_base{m}(etags{i,1}) = [];
+            c_alt{m,1} = etags{i,1};
+            c_alt{m,2} = iv;
+            c_alt{m,3} = boundary - steady_state(iv);
+            c_alt{m,4} = operator;
+            base = true;
+            m = m + 1;
+        end
+    end
+end
+if ~base
+    error('OCCBIN: constraints must come by pair')
+end
+
+function [iv,boundary,operator] = parse_constraint(str,endo_names,params,param_names)
+delim = {'<=','>=','<','>'};
+[c,operator] = strsplit(str,delim);
+operator = operator{1};
+iv = strmatch(strtrim(c{1}),endo_names);
+% try for a number
+boundary = str2num(strtrim(c{2}));
+% if not a number try for a parameter name
+if isempty(boundary)
+    k = strmatch(strtrim(c{2}),param_names);
+    if isempty(k)
+        error(['OCCBIN: illegal constraint ' str]);
+    end
+    boundary = params(k);
+end

matlab/occbin/get_residuals.m

@@ -0,0 +1,10 @@
+function r = get_residuals(ivar,lb,ub,M,oo)
+
+ss = oo.steady_state;
+for i = 1:length(ivar)
+    % only one is different from zero
+    ss(ivar(i)) = lb(i) + ub(i);
+end
+oo.steady_state = ss;
+
+r = evaluate_model(M,oo);

matlab/occbin/occbin.m

@@ -0,0 +1,19 @@
+function [endo_simul,endo_simul_no_constraint,status] = occbin(M,oo,options)
+% function oo=occbin(M,oo,options) solves linear models with occasionally
+% binding constraints using OCCBIN by L. Guerrieri
+
+status = 1;
+constraint_nbr = sum(strcmp(upper(M.equations_tags(:,2)),'OCCBIN'))/2;
+
+switch(constraint_nbr)
+  case 1
+    [zdatalinear_ zdatapiecewise_ zdatass_ oobase_ ] = ...
+    solve_one_constraint(M,oo,options);
+  case 2
+    [zdatalinear_ zdatapiecewise_ zdatass_ oobase_ ] = ...
+    solve_two_constraints(M,oo,options);
+  otherwise
+    error('OCCBIN can only handle two constraints in a model')
+end
+endo_simul = zdatapiecewise_';
+endo_simul_no_constraint = zdatalinear_';

matlab/occbin/test_constraint.m

@@ -0,0 +1,16 @@
+function b=test_constraint(x,constr)
+b = zeros(size(x,1),size(constr,1));
+for i=1:size(constr,1)
+    switch constr{i,4}
+      case '<'
+        b(:,i) = ~(x(:,constr{i,2}) < constr{i,3});
+      case '<='
+        b(:,i) = ~(x(:,constr{i,2}) <= constr{i,3});
+      case '>'
+        b(:,i) = ~(x(:,constr{i,2}) > constr{i,3});
+      case '>='
+        b(:,i) = ~(x(:,constr{i,2}) >= constr{i,3});
+      otherwise
+        error('OCCBIN: wrong inequality sign')
+    end
+end

matlab/optimization/newrat.m

@@ -154,7 +154,7 @@ while norm(gg)>gtol && check==0 && jit<nit
     if length(find(ig))<nx
         ggx=ggx*0;
         ggx(find(ig))=gg(find(ig));
-        if analytic_derivation
+        if analytic_derivation || ~outer_product_gradient
             hhx=hh;
         else
             hhx = reshape(dum,nx,nx);

matlab/posterior_analysis.m

@@ -63,6 +63,13 @@ switch type
     end
     oo_ = variance_decomposition_mc_analysis(SampleSize,'posterior',M_.dname,M_.fname,...
                                              M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_);
+    if ~all(M_.H==0)
+        if strmatch(arg1,options_.varobs,'exact')
+            [observable_name_requested_vars,index_subset,index_observables]=intersect(vartan,options_.varobs,'stable');            
+            oo_ = variance_decomposition_ME_mc_analysis(SampleSize,'posterior',M_.dname,M_.fname,...
+                                             M_.exo_names,arg2,observable_name_requested_vars,arg1,options_.mh_conf_sig,oo_,options_);
+        end
+    end        
   case 'correlation'
     if nargin==narg1
         [nvar,vartan,NumberOfFiles] = ...
@@ -77,6 +84,12 @@ switch type
     end
     oo_ = conditional_variance_decomposition_mc_analysis(SampleSize,'posterior',M_.dname,M_.fname,...
                                                       arg3,M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_);
+    if ~all(M_.H==0)
+        if strmatch(vartan(arg1,:),options_.varobs,'exact')
+            oo_ = conditional_variance_decomposition_ME_mc_analysis(SampleSize,'posterior',M_.dname,M_.fname,...
+                arg3,M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_);
+        end
+    end        
   otherwise
     disp('Not yet implemented')
 end

matlab/prior_analysis.m

@@ -78,6 +78,12 @@ switch type
     end
     oo_ = conditional_variance_decomposition_mc_analysis(SampleSize,'prior',M_.dname,M_.fname,...
                                                       arg3,M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_);
+    if ~all(M_.H==0)
+        if strmatch(vartan(arg1,:),options_.varobs,'exact')
+            oo_ = conditional_variance_decomposition_ME_mc_analysis(SampleSize,'prior',M_.dname,M_.fname,...
+                arg3,M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_);
+        end
+    end        
   otherwise
     disp('Not yet implemented')
 end

matlab/prior_posterior_statistics.m

@@ -226,7 +226,7 @@ if strcmpi(type,'posterior')
     else
         logpost=NaN(B,1);
         for b=1:B
-            [x(b,:), logpost(b)] = GetOneDraw(type);
+            [x(b,:), logpost(b)] = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_);
         end
     end
     localVars.logpost=logpost;

matlab/prior_posterior_statistics_core.m

@@ -190,14 +190,14 @@ end
 for b=fpar:B
     if strcmpi(type,'prior')
 
-        [deep, logpo] = GetOneDraw(type);
+        [deep, logpo] = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_);
 
     else
         deep = x(b,:);
         if strcmpi(type,'posterior')
             logpo = logpost(b);
         else
-            logpo = evaluate_posterior_kernel(deep');
+            logpo = evaluate_posterior_kernel(deep',M_,estim_params_,oo_,options_,bayestopt_);
         end
     end
     M_ = set_all_parameters(deep,estim_params_,M_);

matlab/set_parameters_locally.m

@@ -56,7 +56,7 @@ end
 % and update offset
 offset = offset + nvx + nvn;
 
-% correlations amonx shocks (ncx)
+% correlations among shocks (ncx)
 if ncx
     corrx = estim_params_.corrx;
     for i=1:ncx

matlab/simult.m

@@ -1,4 +1,4 @@
-function [y_,DynareResults] =simult(y0, dr,DynareModel,DynareOptions,DynareResults)
+function [y_out,DynareResults] =simult(y0, dr,DynareModel,DynareOptions,DynareResults)
 % Simulate a DSGE model (perturbation approach).
 
 %@info:
@@ -25,7 +25,7 @@ function [y_,DynareResults] =simult(y0, dr,DynareModel,DynareOptions,DynareResul
 %! @strong{Outputs}
 %! @sp 1
 %! @table @ @var
-%! @item y_
+%! @item y_out
 %! Matrix of doubles, simulated time series for all the endogenous variables (one per row).
 %! @item DynareResults
 %! Matlab's structure gathering the results (see @ref{oo_}).
@@ -88,6 +88,9 @@ for i=1:replic
     if replic > 1
         fwrite(fh,y_,'float64');
     end
+    if i==1
+        y_out=y_;
+    end
 end
 
 if replic > 1

matlab/stoch_simul.m

@@ -36,6 +36,11 @@ if options_.order == 1
     options_.replic = 1;
 end
 
+if M_.hessian_eq_zero && options_.order~=1
+    options_.order = 1;
+    warning('stoch_simul: using order = 1 because Hessian is equal to zero');
+end
+
 if isempty(options_.qz_criterium)
     options_.qz_criterium = 1+1e-6;
 end
@@ -116,6 +121,19 @@ if ~options_.noprint
         lh = size(labels,2)+2;
         dyn_latex_table(M_,options_,my_title,'covar_ex_shocks',headers,labels,M_.Sigma_e,lh,10,6);
     end
+    if ~all(M_.H==0)
+        my_title='MATRIX OF COVARIANCE OF MEASUREMENT ERRORS';
+        labels = [repmat('SE_',length(options_.varobs),1),char(options_.varobs')];
+        headers = char('Variables',labels);
+        lh = size(labels,2)+2;
+        dyntable(options_,my_title,headers,labels,M_.H,lh,10,6);
+        if options_.TeX
+            labels = deblank(M_.exo_names_tex);
+            headers = char('Variables',labels);
+            lh = size(labels,2)+2;
+            dyn_latex_table(M_,options_,my_title,'covar_ME',headers,labels,M_.H,lh,10,6);
+        end
+    end
     if options_.partial_information
         skipline()
         disp('SOLUTION UNDER PARTIAL INFORMATION')

matlab/utilities/dataset/makedataset.m

@@ -6,8 +6,9 @@ function [DynareDataset, DatasetInfo, newdatainterface] = makedataset(DynareOpti
 % INPUTS
 % ======
 %
-%     DynareOptions [struct] Structure of options built by Dynare's preprocessor.
+%     DynareOptions         [struct]    Structure of options built by Dynare's preprocessor.
-%
+%     initialconditions     [double]    number of lags for VAR and DSGE_VAR
+%     gsa_flag              [integer]   1: GSA, 0: other
 %
 % OUTPUTS
 % =======

matlab/variance_decomposition_ME_mc_analysis.m

@@ -0,0 +1,112 @@
+function oo_ = variance_decomposition_ME_mc_analysis(NumberOfSimulations,type,dname,fname,exonames,exo,vartan,var,mh_conf_sig,oo_,options_)
+% function oo_ = variance_decomposition_ME_mc_analysis(NumberOfSimulations,type,dname,fname,exonames,exo,vartan,var,mh_conf_sig,oo_)
+% This function analyses the (posterior or prior) distribution of the
+% endogenous variables' variance decomposition.
+%
+% INPUTS
+%   NumberOfSimulations     [integer]           scalar, number of simulations.
+%   type                    [string]            'prior' or 'posterior'
+%   dname                   [string]            directory name where to save
+%   fname                   [string]            name of the mod-file
+%   exonames                [string]            (n_exo*char_length) character array with names of exogenous variables
+%   exo                     [string]            name of current exogenous
+%                                               variable
+%   vartan                  [string]            (n_endo*char_length) character array with name
+%                                               of endogenous variables
+%   var                     [integer]         index of the current
+%                                               endogenous variable
+%   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.
+
+
+
+% Copyright (C) 2017 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/>.
+
+if strcmpi(type,'posterior')
+    TYPE = 'Posterior';
+    PATH = [dname '/metropolis/'];
+else
+    TYPE = 'Prior';
+    PATH = [dname '/prior/moments/'];
+end
+
+indx = check_name(vartan,var);
+if isempty(indx)
+    disp([ type '_analysis:: ' var ' is not a stationary endogenous variable!'])
+    return
+end
+jndx = check_name(exonames,exo);
+if isempty(jndx)
+    if isequal(exo,'ME')
+       jndx=size(exonames,1)+1; 
+    else
+        disp([ type '_analysis:: ' exo ' is not a declared exogenous variable!'])
+    return
+    end
+end
+
+var=deblank(var);
+exo=deblank(exo);
+
+name = [ var '.' exo ];
+if isfield(oo_, [ TYPE 'TheoreticalMoments'])
+    temporary_structure = oo_.([TYPE, 'TheoreticalMoments']);
+    if isfield(temporary_structure,'dsge')
+        temporary_structure = oo_.([TYPE, 'TheoreticalMoments']).dsge;
+        if isfield(temporary_structure,'VarianceDecompositionME')
+            temporary_structure = oo_.([TYPE, 'TheoreticalMoments']).dsge.VarianceDecompositionME.Mean;
+            if isfield(temporary_structure,name)
+                % Nothing to do.
+                return
+            end
+        end
+    end
+end
+
+ListOfFiles = dir([ PATH  fname '_' TYPE 'VarianceDecompME*.mat']);
+i1 = 1; tmp = zeros(NumberOfSimulations,1);
+indice = (indx-1)*rows(exonames)+jndx;
+for file = 1:length(ListOfFiles)
+    load([ PATH ListOfFiles(file).name ]);
+    i2 = i1 + rows(Decomposition_array_ME) - 1;
+    tmp(i1:i2) = Decomposition_array_ME(:,indice);
+    i1 = i2+1;
+end
+
+if options_.estimation.moments_posterior_density.indicator
+    [p_mean, p_median, p_var, hpd_interval, p_deciles, density] = ...
+        posterior_moments(tmp,1,mh_conf_sig);
+else
+    [p_mean, p_median, p_var, hpd_interval, p_deciles] = ...
+        posterior_moments(tmp,0,mh_conf_sig);
+end
+
+oo_.([TYPE, 'TheoreticalMoments']).dsge.VarianceDecompositionME.Mean.(var).(exo) = p_mean;
+oo_.([TYPE, 'TheoreticalMoments']).dsge.VarianceDecompositionME.Median.(var).(exo) = p_median;
+oo_.([TYPE, 'TheoreticalMoments']).dsge.VarianceDecompositionME.Variance.(var).(exo) = p_var;
+oo_.([TYPE, 'TheoreticalMoments']).dsge.VarianceDecompositionME.HPDinf.(var).(exo) = hpd_interval(1);
+oo_.([TYPE, 'TheoreticalMoments']).dsge.VarianceDecompositionME.HPDsup.(var).(exo) = hpd_interval(2);
+oo_.([TYPE, 'TheoreticalMoments']).dsge.VarianceDecompositionME.deciles.(var).(exo) = p_deciles;
+if options_.estimation.moments_posterior_density.indicator
+    oo_.([TYPE, 'TheoreticalMoments']).dsge.VarianceDecompositionME.density.(var).(exo) = density;
+end

matlab/variance_decomposition_mc_analysis.m

@@ -57,8 +57,10 @@ if isempty(indx)
 end
 jndx = check_name(exonames,exo);
 if isempty(jndx)
-    disp([ type '_analysis:: ' exo ' is not a declared exogenous variable!'])
+        if ~isequal(exo,'ME')
-    return
+            disp([ type '_analysis:: ' exo ' is not a declared exogenous variable!'])
+        end
+        return
 end
 
 var=deblank(var);

preprocessor/ComputingTasks.cc

@@ -2053,7 +2053,7 @@ PlannerObjectiveStatement::getPlannerObjective() const
 void
 PlannerObjectiveStatement::computingPass()
 {
-  model_tree->computingPass(eval_context_t(), false, true, true, none, false, false);
+  model_tree->computingPass(eval_context_t(), false, true, true, none, false, false, false);
   computing_pass_called = true;
 }
 
@@ -4686,3 +4686,141 @@ Smoother2histvalStatement::writeJsonOutput(ostream &output) const
     }
   output << "}";
 }
+
+GMMEstimationStatement::GMMEstimationStatement(const SymbolList &symbol_list_arg,
+                                               const OptionsList &options_list_arg) :
+  symbol_list(symbol_list_arg),
+  options_list(options_list_arg)
+{
+}
+
+void
+GMMEstimationStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
+{
+  symbol_list.writeOutput("var_list_", output);
+  options_list.writeOutput(output);
+  output << "[M_, oo_, estim_params_, bayestopt_, dataset_, dataset_info] = "
+         << "GMM_SMM_estimation_core(var_list_, M_, options_, oo_, estim_params_, bayestopt_, dataset_, dataset_info, 'GMM');" << endl;
+}
+
+void
+GMMEstimationStatement::writeJsonOutput(ostream &output) const
+{
+  output << "{\"statementName\": \"gmm_estimation\"";
+  if (options_list.getNumberOfOptions())
+    {
+      output << ", ";
+      options_list.writeJsonOutput(output);
+    }
+  if (!symbol_list.empty())
+    {
+      output << ", ";
+      symbol_list.writeJsonOutput(output);
+    }
+  output << "}";
+}
+
+SMMEstimationStatement::SMMEstimationStatement(const SymbolList &symbol_list_arg,
+                                               const OptionsList &options_list_arg) :
+  symbol_list(symbol_list_arg),
+  options_list(options_list_arg)
+{
+}
+
+void
+SMMEstimationStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
+{
+  symbol_list.writeOutput("var_list_", output);
+  options_list.writeOutput(output);
+  output << "[M_, oo_, estim_params_, bayestopt_, dataset_, dataset_info] = "
+         << "GMM_SMM_estimation_core(var_list_, M_, options_, oo_, estim_params_, bayestopt_, dataset_, dataset_info, 'SMM');" << endl;
+}
+
+void
+SMMEstimationStatement::writeJsonOutput(ostream &output) const
+{
+  output << "{\"statementName\": \"smm_estimation\"";
+  if (options_list.getNumberOfOptions())
+    {
+      output << ", ";
+      options_list.writeJsonOutput(output);
+    }
+  if (!symbol_list.empty())
+    {
+      output << ", ";
+      symbol_list.writeJsonOutput(output);
+    }
+  output << "}";
+}
+
+GenerateIRFsStatement::GenerateIRFsStatement(const OptionsList &options_list_arg,
+                                             const vector<string> &generate_irf_names_arg,
+                                             const vector<map<string, double> > &generate_irf_elements_arg) :
+  options_list(options_list_arg),
+  generate_irf_names(generate_irf_names_arg),
+  generate_irf_elements(generate_irf_elements_arg)
+{
+}
+
+void
+GenerateIRFsStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
+{
+  options_list.writeOutput(output);
+
+  if (generate_irf_names.empty())
+    return;
+
+  output << "options_.irf_opt.irf_shock_graphtitles = { ";
+  for (vector<string>::const_iterator it = generate_irf_names.begin();
+       it != generate_irf_names.end(); it++)
+    output << "'" << *it << "'; ";
+  output << "};" << endl;
+
+  output << "options_.irf_opt.irf_shocks = zeros(M_.exo_nbr, "
+         << generate_irf_names.size() << ");" << endl;
+
+  for (size_t i = 0; i < generate_irf_names.size(); i++)
+    {
+      map<string, double> m = generate_irf_elements[i];
+      for (map<string, double>::const_iterator it = m.begin();
+         it != m.end(); it++)
+        output << "options_.irf_opt.irf_shocks(M_.exo_names == '"
+               << it->first << "', " << i + 1 << ") = "
+               << it->second << ";" << endl;
+    }
+}
+
+void
+GenerateIRFsStatement::writeJsonOutput(ostream &output) const
+{
+  output << "{\"statementName\": \"generate_irfs\"";
+  if (options_list.getNumberOfOptions())
+    {
+      output << ", ";
+      options_list.writeJsonOutput(output);
+    }
+
+  if (!generate_irf_names.empty())
+    {
+      output << ", \"irf_elements\": [";
+      for (size_t i = 0; i < generate_irf_names.size(); i++)
+        {
+          output << "{\"name\": \"" << generate_irf_names[i] << "\", \"shocks\": [";
+          map<string, double> m = generate_irf_elements[i];
+          size_t idx = 0;
+          for (map<string, double>::const_iterator it = m.begin();
+               it != m.end(); it++, idx++)
+            {
+              output << "{\"exogenous_variable\": \"" << it->first << "\", "
+                     << "\"exogenous_variable_value\": \"" << it->second << "\"}";
+              if (idx + 1 < m.size())
+                output << ", ";
+            }
+          output << "]}";
+          if (i + 1 < generate_irf_names.size())
+            output << ", ";
+        }
+      output << "]";
+    }
+  output << "}";
+}

preprocessor/ComputingTasks.hh

@@ -1138,4 +1138,41 @@ public:
   virtual void writeJsonOutput(ostream &output) const;
 };
 
+class GMMEstimationStatement : public Statement
+{
+private:
+  const SymbolList symbol_list;
+  const OptionsList options_list;
+public:
+  GMMEstimationStatement(const SymbolList &symbol_list_arg, const OptionsList &options_list_arg);
+  virtual void writeOutput(ostream &output, const string &basename, bool minimal_workspace) const;
+  virtual void writeJsonOutput(ostream &output) const;
+};
+
+class SMMEstimationStatement : public Statement
+{
+private:
+  const SymbolList symbol_list;
+  const OptionsList options_list;
+public:
+  SMMEstimationStatement(const SymbolList &symbol_list_arg, const OptionsList &options_list_arg);
+  virtual void writeOutput(ostream &output, const string &basename, bool minimal_workspace) const;
+  virtual void writeJsonOutput(ostream &output) const;
+};
+
+class GenerateIRFsStatement : public Statement
+{
+public:
+private:
+  const OptionsList options_list;
+  const vector<string> generate_irf_names;
+  const vector<map<string, double> > generate_irf_elements;
+public:
+  GenerateIRFsStatement(const OptionsList &options_list_arg,
+                        const vector<string> &generate_irf_names_arg,
+                        const vector<map<string, double> > &generate_irf_elements_arg);
+  virtual void writeOutput(ostream &output, const string &basename, bool minimal_workspace) const;
+  virtual void writeJsonOutput(ostream &output) const;
+};
+
 #endif

preprocessor/DynamicModel.cc

@@ -2539,6 +2539,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia
 
       DynamicOutput << "#=" << endl << comments.str() << "=#" << endl
                     << "  @assert size(g2) == (" << nrows << ", " << hessianColsNbr << ")" << endl
+                    << "  fill!(g2, 0.0)" << endl
                     << "  dynamic!(y, x, params, steady_state, it_, residual, g1)" << endl;
       if (second_derivatives.size())
         DynamicOutput << model_local_vars_output.str()
@@ -2562,6 +2563,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia
 
       DynamicOutput << "#=" << endl << comments.str() << "=#" << endl
                     << "  @assert size(g3) == (" << nrows << ", " << ncols << ")" << endl
+                    << "  fill!(g3, 0.0)" << endl
                     << "  dynamic!(y, x, params, steady_state, it_, residual, g1, g2)" << endl;
       if (third_derivatives.size())
         DynamicOutput << model_local_vars_output.str()
@@ -3270,7 +3272,7 @@ DynamicModel::addEquationsForVar(map<string, pair<SymbolList, int> > var_model_i
 void
 DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivatives, int paramsDerivsOrder,
                             const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll,
-                            bool bytecode)
+                            bool bytecode, const bool nopreprocessoroutput)
 {
   assert(jacobianExo || !(hessian || thirdDerivatives || paramsDerivsOrder));
 
@@ -3293,19 +3295,22 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
     }
 
   // Launch computations
-  cout << "Computing dynamic model derivatives:" << endl
+  if (!nopreprocessoroutput)
-       << " - order 1" << endl;
+    cout << "Computing dynamic model derivatives:" << endl
+         << " - order 1" << endl;
   computeJacobian(vars);
 
   if (hessian)
     {
-      cout << " - order 2" << endl;
+      if (!nopreprocessoroutput)
+        cout << " - order 2" << endl;
       computeHessian(vars);
     }
 
   if (paramsDerivsOrder > 0)
     {
-      cout << " - derivatives of Jacobian/Hessian w.r. to parameters" << endl;
+      if (!nopreprocessoroutput)
+        cout << " - derivatives of Jacobian/Hessian w.r. to parameters" << endl;
       computeParamsDerivatives(paramsDerivsOrder);
 
       if (!no_tmp_terms)
@@ -3314,7 +3319,8 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
 
   if (thirdDerivatives)
     {
-      cout << " - order 3" << endl;
+      if (!nopreprocessoroutput)
+        cout << " - order 3" << endl;
       computeThirdDerivatives(vars);
     }
 
@@ -3336,7 +3342,8 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
 
       equation_type_and_normalized_equation = equationTypeDetermination(first_order_endo_derivatives, variable_reordered, equation_reordered, mfs);
 
-      cout << "Finding the optimal block decomposition of the model ...\n";
+      if (!nopreprocessoroutput)
+        cout << "Finding the optimal block decomposition of the model ...\n";
 
       lag_lead_vector_t equation_lag_lead, variable_lag_lead;
 
@@ -3858,7 +3865,7 @@ DynamicModel::replaceMyEquations(DynamicModel &dynamic_model) const
 }
 
 void
-DynamicModel::computeRamseyPolicyFOCs(const StaticModel &static_model)
+DynamicModel::computeRamseyPolicyFOCs(const StaticModel &static_model, const bool nopreprocessoroutput)
 {
   // Add aux LM to constraints in equations
   // equation[i]->lhs = rhs becomes equation[i]->MULT_(i+1)*(lhs-rhs) = 0
@@ -3869,8 +3876,8 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel &static_model)
       assert(substeq != NULL);
       equations[i] = substeq;
     }
-
+  if (!nopreprocessoroutput)
-  cout << "Ramsey Problem: added " << i << " Multipliers." << endl;
+    cout << "Ramsey Problem: added " << i << " Multipliers." << endl;
 
   // Add Planner Objective to equations to include in computeDerivIDs
   assert(static_model.equations.size() == 1);
@@ -5632,84 +5639,43 @@ DynamicModel::writeJsonOutput(ostream &output) const
   writeJsonXrefs(output);
 }
 
+void
+DynamicModel::writeJsonXrefsHelper(ostream &output, const map<pair<int, int>, set<int> > &xrefs) const
+{
+  for (map<pair<int, int>, set<int> >::const_iterator it = xrefs.begin();
+       it != xrefs.end(); it++)
+    {
+      if (it != xrefs.begin())
+        output << ", ";
+      output << "{\"name\": \"" << symbol_table.getName(it->first.first) << "\""
+             << ", \"shift\": " << it->first.second
+             << ", \"equations\": [";
+      for (set<int>::const_iterator it1 = it->second.begin();
+           it1 != it->second.end(); it1++)
+        {
+          if (it1 != it->second.begin())
+            output << ", ";
+          output << *it1 + 1;
+        }
+      output << "]}";
+    }
+}
+
 void
 DynamicModel::writeJsonXrefs(ostream &output) const
 {
   output << "\"xrefs\": {"
          << "\"parameters\": [";
-  for (map<pair<int, int>, set<int> >::const_iterator it = xref_param.begin();
+  writeJsonXrefsHelper(output, xref_param);
-       it != xref_param.end(); it++)
-    {
-      if (it != xref_param.begin())
-        output << ", ";
-      output << "{\"parameter\": \"" << symbol_table.getName(it->first.first) << "\""
-             << ", \"equations\": [";
-      for (set<int>::const_iterator it1 = it->second.begin();
-           it1 != it->second.end(); it1++)
-        {
-          if (it1 != it->second.begin())
-            output << ", ";
-          output << *it1 + 1;
-        }
-      output << "]}";
-    }
   output << "]"
          << ", \"endogenous\": [";
-  for (map<pair<int, int>, set<int> >::const_iterator it = xref_endo.begin();
+  writeJsonXrefsHelper(output, xref_endo);
-       it != xref_endo.end(); it++)
-    {
-      if (it != xref_endo.begin())
-        output << ", ";
-      output << "{\"endogenous\": \"" << symbol_table.getName(it->first.first) << "\""
-             << ", \"shift\": " << it->first.second
-             << ", \"equations\": [";
-      for (set<int>::const_iterator it1 = it->second.begin();
-           it1 != it->second.end(); it1++)
-        {
-          if (it1 != it->second.begin())
-            output << ", ";
-          output << *it1 + 1;
-        }
-      output << "]}";
-    }
   output << "]"
          << ", \"exogenous\": [";
-  for (map<pair<int, int>, set<int> >::const_iterator it = xref_exo.begin();
+    writeJsonXrefsHelper(output, xref_exo);
-       it != xref_exo.end(); it++)
-    {
-      if (it != xref_exo.begin())
-        output << ", ";
-      output << "{\"exogenous\": \"" << symbol_table.getName(it->first.first) << "\""
-             << ", \"shift\": " << it->first.second
-             << ", \"equations\": [";
-      for (set<int>::const_iterator it1 = it->second.begin();
-           it1 != it->second.end(); it1++)
-        {
-          if (it1 != it->second.begin())
-            output << ", ";
-          output << *it1 + 1;
-        }
-      output << "]}";
-    }
   output << "]"
          << ", \"exogenous_deterministic\": [";
-  for (map<pair<int, int>, set<int> >::const_iterator it = xref_exo_det.begin();
+  writeJsonXrefsHelper(output, xref_exo_det);
-       it != xref_exo_det.end(); it++)
-    {
-      if (it != xref_exo_det.begin())
-        output << ", ";
-      output << "{\"exogenous_det\": \"" << symbol_table.getName(it->first.first) << "\""
-             << ", \"shift\": " << it->first.second
-             << ", \"equations\": [";
-      for (set<int>::const_iterator it1 = it->second.begin();
-           it1 != it->second.end(); it1++)
-        {
-          if (it1 != it->second.begin())
-            output << ", ";
-          output << *it1 + 1;
-        }
-      output << "]}";
-    }
   output << "]}" << endl;
 }
 

preprocessor/DynamicModel.hh

@@ -248,7 +248,7 @@ public:
     \param no_tmp_terms if true, no temporary terms will be computed in the dynamic files
   */
   void computingPass(bool jacobianExo, bool hessian, bool thirdDerivatives, int paramsDerivsOrder,
-                     const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll, bool bytecode);
+                     const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll, bool bytecode, const bool nopreprocessoroutput);
   //! Writes model initialization and lead/lag incidence matrix to output
   void writeOutput(ostream &output, const string &basename, bool block, bool byte_code, bool use_dll, int order, bool estimation_present, bool compute_xrefs, bool julia) const;
 
@@ -269,6 +269,7 @@ public:
 
   //! Write cross reference output if the xref maps have been filed
   void writeJsonXrefs(ostream &output) const;
+  void writeJsonXrefsHelper(ostream &output, const map<pair<int, int>, set<int> > &xrefs) const;
 
   //! Return true if the hessian is equal to zero
   inline bool checkHessianZero() const;
@@ -306,7 +307,7 @@ public:
   void cloneDynamic(DynamicModel &dynamic_model) const;
 
   //! Replaces model equations with derivatives of Lagrangian w.r.t. endogenous
-  void computeRamseyPolicyFOCs(const StaticModel &static_model);
+  void computeRamseyPolicyFOCs(const StaticModel &static_model, const bool nopreprocessoroutput);
   //! Replaces the model equations in dynamic_model with those in this model
   void replaceMyEquations(DynamicModel &dynamic_model) const;
 

preprocessor/DynareBison.yy

@@ -113,9 +113,9 @@ class ParsingDriver;
 %token CPF_WEIGHTS AMISANOTRISTANI MURRAYJONESPARSLOW WRITE_EQUATION_TAGS METHOD
 %token NONLINEAR_FILTER_INITIALIZATION FILTER_ALGORITHM PROPOSAL_APPROXIMATION CUBATURE UNSCENTED MONTECARLO DISTRIBUTION_APPROXIMATION
 %token <string_val> NAME
-%token USE_PENALIZED_OBJECTIVE_FOR_HESSIAN INIT_STATE RESCALE_PREDICTION_ERROR_COVARIANCE
+%token USE_PENALIZED_OBJECTIVE_FOR_HESSIAN INIT_STATE RESCALE_PREDICTION_ERROR_COVARIANCE GENERATE_IRFS
 %token NAN_CONSTANT NO_STATIC NOBS NOCONSTANT NODISPLAY NOCORR NODIAGNOSTIC NOFUNCTIONS NO_HOMOTOPY
-%token NOGRAPH POSTERIOR_NOGRAPH POSTERIOR_GRAPH NOMOMENTS NOPRINT NORMAL_PDF SAVE_DRAWS MODEL_NAME
+%token NOGRAPH POSTERIOR_NOGRAPH POSTERIOR_GRAPH NOMOMENTS NOPRINT NORMAL_PDF SAVE_DRAWS MODEL_NAME STDERR_MULTIPLES DIAGONAL_ONLY
 %token OBSERVATION_TRENDS OPTIM OPTIM_WEIGHTS ORDER OSR OSR_PARAMS MAX_DIM_COVA_GROUP ADVANCED OUTFILE OUTVARS OVERWRITE
 %token PARALLEL_LOCAL_FILES PARAMETERS PARAMETER_SET PARTIAL_INFORMATION PERIODS PERIOD PLANNER_OBJECTIVE PLOT_CONDITIONAL_FORECAST PLOT_PRIORS PREFILTER PRESAMPLE
 %token PERFECT_FORESIGHT_SETUP PERFECT_FORESIGHT_SOLVER NO_POSTERIOR_KERNEL_DENSITY FUNCTION
@@ -151,7 +151,7 @@ class ParsingDriver;
 %token VLISTLOG VLISTPER SPECTRAL_DENSITY
 %token RESTRICTION RESTRICTION_FNAME CROSS_RESTRICTIONS NLAGS CONTEMP_REDUCED_FORM REAL_PSEUDO_FORECAST
 %token DUMMY_OBS NSTATES INDXSCALESSTATES NO_BAYESIAN_PRIOR SPECIFICATION SIMS_ZHA
-%token <string_val> ALPHA BETA ABAND NINV CMS NCMS CNUM GAMMA INV_GAMMA INV_GAMMA1 INV_GAMMA2 NORMAL UNIFORM EPS PDF FIG DR NONE PRIOR PRIOR_VARIANCE HESSIAN IDENTITY_MATRIX DIRICHLET
+%token <string_val> ALPHA BETA ABAND NINV CMS NCMS CNUM GAMMA INV_GAMMA INV_GAMMA1 INV_GAMMA2 NORMAL UNIFORM EPS PDF FIG DR NONE PRIOR PRIOR_VARIANCE HESSIAN IDENTITY_MATRIX DIRICHLET DIAGONAL OPTIMAL
 %token GSIG2_LMDM Q_DIAG FLAT_PRIOR NCSK NSTD WEIBULL WEIBULL_PDF
 %token INDXPARR INDXOVR INDXAP APBAND INDXIMF IMFBAND INDXFORE FOREBAND INDXGFOREHAT INDXGIMFHAT
 %token INDXESTIMA INDXGDLS EQ_MS FILTER_COVARIANCE FILTER_DECOMPOSITION SMOOTHED_STATE_UNCERTAINTY
@@ -168,13 +168,17 @@ class ParsingDriver;
 %token SHOCK_DRAWS FREE_PARAMETERS MEDIAN DATA_OBS_NBR NEIGHBORHOOD_WIDTH PVALUE_KS PVALUE_CORR
 %token FILTERED_PROBABILITIES REAL_TIME_SMOOTHED PRIOR_FUNCTION POSTERIOR_FUNCTION SAMPLING_DRAWS
 %token PROPOSAL_TYPE PROPOSAL_UPPER_BOUND PROPOSAL_LOWER_BOUND PROPOSAL_DRAWS USE_MEAN_CENTER
-%token ADAPTIVE_MH_DRAWS THINNING_FACTOR COEFFICIENTS_PRIOR_HYPERPARAMETERS
+%token ADAPTIVE_MH_DRAWS THINNING_FACTOR COEFFICIENTS_PRIOR_HYPERPARAMETERS SMM_ESTIMATION GMM_ESTIMATION
 %token CONVERGENCE_STARTING_VALUE CONVERGENCE_ENDING_VALUE CONVERGENCE_INCREMENT_VALUE
 %token MAX_ITERATIONS_STARTING_VALUE MAX_ITERATIONS_INCREMENT_VALUE MAX_BLOCK_ITERATIONS
 %token MAX_REPEATED_OPTIMIZATION_RUNS FUNCTION_CONVERGENCE_CRITERION SAVE_REALTIME
 %token PARAMETER_CONVERGENCE_CRITERION NUMBER_OF_LARGE_PERTURBATIONS NUMBER_OF_SMALL_PERTURBATIONS
 %token NUMBER_OF_POSTERIOR_DRAWS_AFTER_PERTURBATION MAX_NUMBER_OF_STAGES
 %token RANDOM_FUNCTION_CONVERGENCE_CRITERION RANDOM_PARAMETER_CONVERGENCE_CRITERION
+%token CENTERED_MOMENTS AUTOLAG RECURSIVE_ORDER_ESTIMATION BARTLETT_KERNEL_LAG WEIGHTING_MATRIX PENALIZED_ESTIMATOR VERBOSE
+%token SIMULATION_MULTIPLE SEED BOUNDED_SHOCK_SUPPORT
+%token ANALYTICAL_GIRF IRF_IN_PERCENT EMAS_GIRF EMAS_DROP EMAS_TOLF EMAS_MAX_ITER
+
 %token <vector_string_val> SYMBOL_VEC
 
 %type <node_val> expression expression_or_empty
@@ -280,6 +284,7 @@ statement : parameters
           | external_function
           | steady_state_model
           | trend_var
+          | generate_irfs
           | log_trend_var
           | ms_estimation
           | ms_simulation
@@ -299,6 +304,8 @@ statement : parameters
           | perfect_foresight_solver
           | prior_function
           | posterior_function
+          | gmm_estimation
+          | smm_estimation
           | shock_groups
           ;
 
@@ -1215,6 +1222,100 @@ perfect_foresight_solver_options : o_stack_solve_algo
                                  | o_pf_tolx
                                  ;
 
+gmm_smm_common_option : o_datafile
+                      | o_nobs
+                      | o_first_obs
+                      | o_optim
+                      | o_mode_file
+                      | o_mode_compute
+                      | o_prior_trunc
+                      | o_loglinear
+                      | o_logdata
+                      | o_relative_irf
+                      | o_irf
+                      | o_tex
+                      | o_xls_sheet
+                      | o_xls_range
+                      | o_solve_algo
+                      | o_plot_priors
+                      | o_aim_solver
+                      | o_selected_variables_only
+                      | o_irf_shocks
+                      | o_sylvester
+                      | o_sylvester_fixed_point_tol
+                      | o_lyapunov
+                      | o_lyapunov_fixed_point_tol
+                      | o_lyapunov_doubling_tol
+                      | o_dr
+                      | o_dr_cycle_reduction_tol
+                      | o_dr_logarithmic_reduction_tol
+                      | o_dr_logarithmic_reduction_maxiter
+                      | o_qz_zero_threshold
+                      | o_irf_plot_threshold
+                      | o_consider_all_endogenous
+                      | o_consider_only_observed
+                      | o_dirname
+                      | o_huge_number
+                      | o_silent_optimizer
+                      | o_nograph
+                      | o_nodisplay
+                      | o_graph_format
+                      | o_analytical_girf
+                      | o_irf_in_percent
+                      | o_emas_girf
+                      | o_emas_drop
+                      | o_emas_tolf
+                      | o_emas_max_iter
+                      | o_stderr_multiples
+                      | o_diagonal_only
+                      ;
+
+gmm_estimation : GMM_ESTIMATION '(' gmm_estimation_options_list ')' ';'
+                { driver.gmm_estimation(); }
+               | GMM_ESTIMATION '(' gmm_estimation_options_list ')' symbol_list ';'
+                { driver.gmm_estimation(); }
+               ;
+
+gmm_estimation_options_list : gmm_estimation_option COMMA gmm_estimation_options_list
+                            | gmm_estimation_option
+                            ;
+
+gmm_estimation_option : gmm_smm_common_option
+                       | o_gmm_order
+                       | o_gmm_centered_moments
+                       | o_gmm_autolag
+                       | o_gmm_recursive_order_estimation
+                       | o_gmm_bartlett_kernel_lag
+                       | o_gmm_weighting_matrix
+                       | o_gmm_penalized_estimator
+                       | o_gmm_verbose
+                       ;
+
+smm_estimation : SMM_ESTIMATION '(' smm_estimation_options_list ')' ';'
+                { driver.smm_estimation(); }
+               | SMM_ESTIMATION '(' smm_estimation_options_list ')' symbol_list ';'
+                { driver.smm_estimation(); }
+               ;
+
+smm_estimation_options_list : smm_estimation_option COMMA smm_estimation_options_list
+                            | smm_estimation_option
+                            ;
+
+smm_estimation_option : gmm_smm_common_option
+                       | o_smm_order
+                       | o_smm_centered_moments
+                       | o_smm_autolag
+                       | o_smm_recursive_order_estimation
+                       | o_smm_bartlett_kernel_lag
+                       | o_smm_weighting_matrix
+                       | o_smm_penalized_estimator
+                       | o_smm_verbose
+                       | o_smm_simulation_multiple
+                       | o_smm_drop
+                       | o_smm_seed
+                       | o_smm_bounded_shock_support
+                       ;
+
 prior_function : PRIOR_FUNCTION '(' prior_posterior_function_options_list ')' ';'
                 { driver.prior_posterior_function(true); }
                ;
@@ -1291,6 +1392,14 @@ stoch_simul_primary_options : o_dr_algo
                             | o_irf
                             | o_irf_shocks
                             | o_relative_irf
+                            | o_analytical_girf
+                            | o_irf_in_percent
+                            | o_emas_girf
+                            | o_emas_drop
+                            | o_emas_tolf
+                            | o_emas_max_iter
+                            | o_stderr_multiples
+                            | o_diagonal_only
                             | o_hp_filter
                             | o_hp_ngrid
                             | o_periods
@@ -1975,6 +2084,14 @@ estimation_options : o_datafile
                    | o_keep_kalman_algo_if_singularity_is_detected
                    | o_use_penalized_objective_for_hessian
                    | o_rescale_prediction_error_covariance
+                   | o_analytical_girf
+                   | o_irf_in_percent
+                   | o_emas_girf
+                   | o_emas_drop
+                   | o_emas_tolf
+                   | o_emas_max_iter
+                   | o_stderr_multiples
+                   | o_diagonal_only
                    ;
 
 list_optim_option : QUOTED_STRING COMMA QUOTED_STRING
@@ -2843,6 +2960,38 @@ calib_smoother_option : o_filtered_vars
                       | o_parameter_set
                       ;
 
+generate_irfs : GENERATE_IRFS ';' END ';'
+                { driver.end_generate_irfs(); }
+              | GENERATE_IRFS ';' generate_irfs_element_list END ';'
+                { driver.end_generate_irfs(); }
+              | GENERATE_IRFS '(' generate_irfs_options_list ')' ';' END ';'
+                { driver.end_generate_irfs(); }
+              | GENERATE_IRFS '(' generate_irfs_options_list ')' ';' generate_irfs_element_list END ';'
+                { driver.end_generate_irfs(); }
+              ;
+
+generate_irfs_options_list : generate_irfs_option COMMA generate_irfs_options_list
+                           | generate_irfs_option
+                           ;
+
+generate_irfs_option : o_stderr_multiples
+                     | o_diagonal_only
+                     ;
+
+generate_irfs_element_list : generate_irfs_element_list generate_irfs_element
+                           | generate_irfs_element
+                           ;
+
+generate_irfs_element : NAME COMMA generate_irfs_exog_element_list ';'
+                        { driver.add_generate_irfs_element($1); }
+                      ;
+
+generate_irfs_exog_element_list : generate_irfs_exog_element_list COMMA symbol EQUAL signed_number
+                                  { driver.add_generate_irfs_exog_element($3, $5); }
+                                | symbol EQUAL signed_number
+                                  { driver.add_generate_irfs_exog_element($1, $3); }
+                                ;
+
 extended_path : EXTENDED_PATH ';'
                 { driver.extended_path(); }
               | EXTENDED_PATH '(' extended_path_options_list ')' ';'            
@@ -3238,7 +3387,8 @@ o_bvar_prior_omega : BVAR_PRIOR_OMEGA EQUAL INT_NUMBER { driver.option_num("bvar
 o_bvar_prior_flat : BVAR_PRIOR_FLAT { driver.option_num("bvar_prior_flat", "1"); };
 o_bvar_prior_train : BVAR_PRIOR_TRAIN EQUAL INT_NUMBER { driver.option_num("bvar_prior_train", $3); };
 o_bvar_replic : BVAR_REPLIC EQUAL INT_NUMBER { driver.option_num("bvar_replic", $3); };
-
+o_stderr_multiples : STDERR_MULTIPLES { driver.option_num("irf_opt.stderr_multiples", "1"); };
+o_diagonal_only : DIAGONAL_ONLY { driver.option_num("irf_opt.diagonal_only", "1"); };
 o_number_of_particles : NUMBER_OF_PARTICLES EQUAL INT_NUMBER { driver.option_num("particle.number_of_particles", $3); };
 o_resampling : RESAMPLING EQUAL SYSTEMATIC
               | RESAMPLING EQUAL NONE {driver.option_num("particle.resampling.status.systematic", "0"); driver.option_num("particle.resampling.status.none", "1"); }
@@ -3542,6 +3692,59 @@ o_use_shock_groups : USE_SHOCK_GROUPS { driver.option_str("plot_shock_decomp.use
                    ;
 o_colormap : COLORMAP EQUAL symbol { driver.option_num("plot_shock_decomp.colormap",$3); };
 
+o_gmm_order : ORDER EQUAL INT_NUMBER { driver.option_num("gmm.order", $3); };
+o_smm_order : ORDER EQUAL INT_NUMBER { driver.option_num("smm.order", $3); };
+o_gmm_centered_moments : CENTERED_MOMENTS { driver.option_num("gmm.centered_moments", "1"); };
+o_smm_centered_moments : CENTERED_MOMENTS { driver.option_num("smm.centered_moments", "1"); };
+o_gmm_autolag : AUTOLAG EQUAL vec_int 
+                { driver.option_vec_int("gmm.autolag", $3); }
+              | AUTOLAG EQUAL vec_int_number
+                { driver.option_vec_int("gmm.autolag", $3); }
+              ;
+o_smm_autolag : AUTOLAG EQUAL vec_int 
+                { driver.option_vec_int("smm.autolag", $3); }
+              | AUTOLAG EQUAL vec_int_number
+                { driver.option_vec_int("smm.autolag", $3); }
+              ;
+o_gmm_recursive_order_estimation : RECURSIVE_ORDER_ESTIMATION { driver.option_num("gmm.recursive_estimation", "1"); };
+o_smm_recursive_order_estimation : RECURSIVE_ORDER_ESTIMATION { driver.option_num("smm.recursive_estimation", "1"); };
+o_gmm_bartlett_kernel_lag : BARTLETT_KERNEL_LAG EQUAL INT_NUMBER { driver.option_num("gmm.qLag", $3); };
+o_smm_bartlett_kernel_lag : BARTLETT_KERNEL_LAG EQUAL INT_NUMBER { driver.option_num("smm.qLag", $3); };
+o_gmm_weighting_matrix : WEIGHTING_MATRIX EQUAL OPTIMAL
+                     { driver.option_str("gmm.weighting_matrix", $3); }                                  
+                   | WEIGHTING_MATRIX EQUAL IDENTITY_MATRIX
+                     { driver.option_str("gmm.weighting_matrix", $3); }
+                   | WEIGHTING_MATRIX EQUAL DIAGONAL
+                     { driver.option_str("gmm.weighting_matrix", $3); }
+                   | WEIGHTING_MATRIX EQUAL filename
+                     { driver.option_str("gmm.weighting_matrix", $3); }
+                   ;
+o_smm_weighting_matrix : WEIGHTING_MATRIX EQUAL OPTIMAL
+                     { driver.option_str("smm.weighting_matrix", $3); }                                  
+                   | WEIGHTING_MATRIX EQUAL IDENTITY_MATRIX
+                     { driver.option_str("smm.weighting_matrix", $3); }
+                   | WEIGHTING_MATRIX EQUAL DIAGONAL
+                     { driver.option_str("smm.weighting_matrix", $3); }
+                   | WEIGHTING_MATRIX EQUAL filename
+                     { driver.option_str("smm.weighting_matrix", $3); }
+                   ;
+o_gmm_penalized_estimator : PENALIZED_ESTIMATOR { driver.option_num("gmm.penalized_estimator", "1"); };
+o_smm_penalized_estimator : PENALIZED_ESTIMATOR { driver.option_num("smm.penalized_estimator", "1"); };
+o_gmm_verbose : VERBOSE { driver.option_num("gmm.verbose", "1"); };
+o_smm_verbose : VERBOSE { driver.option_num("smm.verbose", "1"); };
+
+o_smm_simulation_multiple : SIMULATION_MULTIPLE EQUAL INT_NUMBER { driver.option_num("smm.simulation_multiple", $3); };
+o_smm_drop : DROP EQUAL INT_NUMBER { driver.option_num("smm.drop", $3); };
+o_smm_seed : SEED EQUAL INT_NUMBER { driver.option_num("smm.seed", $3); };
+o_smm_bounded_shock_support : BOUNDED_SHOCK_SUPPORT { driver.option_num("smm.bounded_support", "1"); };
+
+o_analytical_girf : ANALYTICAL_GIRF { driver.option_num("irf_opt.analytical_GIRF", "1"); };
+o_irf_in_percent : IRF_IN_PERCENT { driver.option_num("irf_opt.percent", "1"); };
+o_emas_girf : EMAS_GIRF { driver.option_num("irf_opt.ergodic_mean_irf", "1"); };
+o_emas_drop : EMAS_DROP EQUAL INT_NUMBER { driver.option_num("irf_opt.EM.drop", $3); };
+o_emas_tolf : EMAS_TOLF EQUAL non_negative_number { driver.option_num("irf_opt.EM.tolf", $3); };
+o_emas_max_iter : EMAS_MAX_ITER EQUAL INT_NUMBER { driver.option_num("irf_opt.EM.iter", $3); };
+
 range : symbol ':' symbol
         {
           $1->append(":");

preprocessor/DynareFlex.ll

@@ -168,6 +168,8 @@ DATE -?[0-9]+([YyAa]|[Mm]([1-9]|1[0-2])|[Qq][1-4]|[Ww]([1-9]{1}|[1-4][0-9]|5[0-2
 <INITIAL>ms_variance_decomposition {BEGIN DYNARE_STATEMENT; return token::MS_VARIANCE_DECOMPOSITION;}
 <INITIAL>conditional_forecast {BEGIN DYNARE_STATEMENT; return token::CONDITIONAL_FORECAST;}
 <INITIAL>plot_conditional_forecast {BEGIN DYNARE_STATEMENT; return token::PLOT_CONDITIONAL_FORECAST;}
+<INITIAL>gmm_estimation {BEGIN DYNARE_STATEMENT; return token::GMM_ESTIMATION;}
+<INITIAL>smm_estimation {BEGIN DYNARE_STATEMENT; return token::SMM_ESTIMATION;}
 
 <INITIAL>markov_switching {BEGIN DYNARE_STATEMENT; return token::MARKOV_SWITCHING;}
 <INITIAL>svar {BEGIN DYNARE_STATEMENT; return token::SVAR;}
@@ -212,6 +214,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
 <INITIAL>irf_calibration {BEGIN DYNARE_BLOCK; return token::IRF_CALIBRATION;}
 <INITIAL>ramsey_constraints {BEGIN DYNARE_BLOCK; return token::RAMSEY_CONSTRAINTS;}
 <INITIAL>restrictions {BEGIN DYNARE_BLOCK; return token::RESTRICTIONS;}
+<INITIAL>generate_irfs {BEGIN DYNARE_BLOCK; return token::GENERATE_IRFS;}
 
  /* For the semicolon after an "end" keyword */
 <INITIAL>; {return Dynare::parser::token_type (yytext[0]);}
@@ -640,6 +643,30 @@ 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>silent_optimizer {return token::SILENT_OPTIMIZER;}
 <DYNARE_STATEMENT>lmmcp {return token::LMMCP;}
 <DYNARE_STATEMENT>occbin {return token::OCCBIN;}
+<DYNARE_STATEMENT>centered_moments {return token::CENTERED_MOMENTS; }
+<DYNARE_STATEMENT>autolag {return token::AUTOLAG; }
+<DYNARE_STATEMENT>recursive_order_estimation {return token::RECURSIVE_ORDER_ESTIMATION; }
+<DYNARE_STATEMENT>bartlett_kernel_lag {return token::BARTLETT_KERNEL_LAG; }
+<DYNARE_STATEMENT>optimal {
+  yylval->string_val = new string(yytext);
+  return token::OPTIMAL;
+}
+<DYNARE_STATEMENT>diagonal  {
+  yylval->string_val = new string(yytext);
+  return token::DIAGONAL;
+}
+<DYNARE_STATEMENT>weighting_matrix {return token::WEIGHTING_MATRIX; }
+<DYNARE_STATEMENT>penalized_estimator {return token::PENALIZED_ESTIMATOR; }
+<DYNARE_STATEMENT>verbose {return token::VERBOSE; }
+<DYNARE_STATEMENT>simulation_multiple {return token::SIMULATION_MULTIPLE; }
+<DYNARE_STATEMENT>seed {return token::SEED; }
+<DYNARE_STATEMENT>bounded_shock_support {return token::BOUNDED_SHOCK_SUPPORT; }
+<DYNARE_STATEMENT>analytical_girf {return token::ANALYTICAL_GIRF; }
+<DYNARE_STATEMENT>irf_in_percent {return token::IRF_IN_PERCENT; }
+<DYNARE_STATEMENT>emas_girf {return token::EMAS_GIRF; }
+<DYNARE_STATEMENT>emas_drop {return token::EMAS_DROP; }
+<DYNARE_STATEMENT>emas_tolf {return token::EMAS_TOLF; }
+<DYNARE_STATEMENT>emas_max_iter {return token::EMAS_MAX_ITER; }
 
 <DYNARE_STATEMENT>[\$][^$]*[\$] {
   strtok(yytext+1, "$");
@@ -714,6 +741,8 @@ 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>irf_plot_threshold {return token::IRF_PLOT_THRESHOLD;}
 <DYNARE_STATEMENT>no_homotopy {return token::NO_HOMOTOPY;}
 
+<DYNARE_BLOCK>stderr_multiples {return token::STDERR_MULTIPLES;}
+<DYNARE_BLOCK>diagonal_only {return token::DIAGONAL_ONLY;}
 <DYNARE_BLOCK>equation {return token::EQUATION;}
 <DYNARE_BLOCK>exclusion {return token::EXCLUSION;}
 <DYNARE_BLOCK>lag {return token::LAG;}

preprocessor/DynareMain.cc

@@ -46,6 +46,7 @@ void main2(stringstream &in, string &basename, bool debug, bool clear_all, bool
            , bool cygwin, bool msvc, bool mingw
 #endif
            , JsonOutputPointType json, JsonFileOutputType json_output_mode, bool onlyjson, bool jsonderivsimple
+           , bool nopreprocessoroutput
            );
 
 void main1(string &modfile, string &basename, string &modfiletxt, bool debug, bool save_macro, string &save_macro_file,
@@ -120,6 +121,7 @@ main(int argc, char **argv)
   bool onlyjson = false;
   bool jsonderivsimple = false;
   LanguageOutputType language = matlab;
+  bool nopreprocessoroutput = false;
 
   // Parse options
   for (int arg = 2; arg < argc; arg++)
@@ -303,6 +305,8 @@ main(int argc, char **argv)
         json_output_mode = standardout;
       else if (!strcmp(argv[arg], "onlyjson"))
         onlyjson = true;
+      else if (!strcmp(argv[arg], "nopreprocessoroutput"))
+        nopreprocessoroutput = true;
       else if (!strcmp(argv[arg], "jsonderivsimple"))
         jsonderivsimple = true;
       else if (strlen(argv[arg]) >= 4 && !strncmp(argv[arg], "json", 4))
@@ -333,8 +337,9 @@ main(int argc, char **argv)
         }
     }
 
-  cout << "Starting Dynare (version " << PACKAGE_VERSION << ")." << endl
+  if (!nopreprocessoroutput)
-       << "Starting preprocessing of the model file ..." << endl;
+    cout << "Starting Dynare (version " << PACKAGE_VERSION << ")." << endl
+         << "Starting preprocessing of the model file ..." << endl;
 
   // Construct basename (i.e. remove file extension if there is one)
   string basename = argv[1];
@@ -397,7 +402,7 @@ main(int argc, char **argv)
 #if defined(_WIN32) || defined(__CYGWIN32__) || defined(__MINGW32__)
         , cygwin, msvc, mingw
 #endif
-        , json, json_output_mode, onlyjson, jsonderivsimple
+        , json, json_output_mode, onlyjson, jsonderivsimple, nopreprocessoroutput
         );
 
   return EXIT_SUCCESS;

preprocessor/DynareMain2.cc

@@ -35,6 +35,7 @@ main2(stringstream &in, string &basename, bool debug, bool clear_all, bool clear
       , bool cygwin, bool msvc, bool mingw
 #endif
       , JsonOutputPointType json, JsonFileOutputType json_output_mode, bool onlyjson, bool jsonderivsimple
+      , bool nopreprocessoroutput
       )
 {
   ParsingDriver p(warnings, nostrict);
@@ -42,38 +43,40 @@ main2(stringstream &in, string &basename, bool debug, bool clear_all, bool clear
   // Do parsing and construct internal representation of mod file
   ModFile *mod_file = p.parse(in, debug);
   if (json == parsing)
-    mod_file->writeJsonOutput(basename, json, json_output_mode, onlyjson);
+    mod_file->writeJsonOutput(basename, json, json_output_mode, onlyjson, nopreprocessoroutput);
 
   // Run checking pass
   mod_file->checkPass(nostrict, stochastic);
   if (json == checkpass)
-    mod_file->writeJsonOutput(basename, json, json_output_mode, onlyjson);
+    mod_file->writeJsonOutput(basename, json, json_output_mode, onlyjson, nopreprocessoroutput);
 
   // Perform transformations on the model (creation of auxiliary vars and equations)
-  mod_file->transformPass(nostrict, stochastic, compute_xrefs || json == transformpass);
+  mod_file->transformPass(nostrict, stochastic, compute_xrefs || json == transformpass, nopreprocessoroutput);
   if (json == transformpass)
-    mod_file->writeJsonOutput(basename, json, json_output_mode, onlyjson);
+    mod_file->writeJsonOutput(basename, json, json_output_mode, onlyjson, nopreprocessoroutput);
 
   // Evaluate parameters initialization, initval, endval and pounds
-  mod_file->evalAllExpressions(warn_uninit);
+  mod_file->evalAllExpressions(warn_uninit, nopreprocessoroutput);
 
   // Do computations
-  mod_file->computingPass(no_tmp_terms, output_mode, params_derivs_order);
+  mod_file->computingPass(no_tmp_terms, output_mode, params_derivs_order, nopreprocessoroutput);
   if (json == computingpass)
-    mod_file->writeJsonOutput(basename, json, json_output_mode, onlyjson, jsonderivsimple);
+    mod_file->writeJsonOutput(basename, json, json_output_mode, onlyjson, nopreprocessoroutput, jsonderivsimple);
 
   // Write outputs
   if (output_mode != none)
-    mod_file->writeExternalFiles(basename, output_mode, language);
+    mod_file->writeExternalFiles(basename, output_mode, language, nopreprocessoroutput);
   else
     mod_file->writeOutputFiles(basename, clear_all, clear_global, no_log, no_warn, console, nograph,
                                nointeractive, config_file, check_model_changes, minimal_workspace, compute_xrefs
 #if defined(_WIN32) || defined(__CYGWIN32__) || defined(__MINGW32__)
                                , cygwin, msvc, mingw
 #endif
+                               , nopreprocessoroutput
                                );
 
   delete mod_file;
 
-  cout << "Preprocessing completed." << endl;
+  if (!nopreprocessoroutput)
+    cout << "Preprocessing completed." << endl;
 }

preprocessor/ModFile.cc

@@ -53,9 +53,10 @@ ModFile::~ModFile()
 }
 
 void
-ModFile::evalAllExpressions(bool warn_uninit)
+ModFile::evalAllExpressions(bool warn_uninit, const bool nopreprocessoroutput)
 {
-  cout << "Evaluating expressions...";
+  if (!nopreprocessoroutput)
+    cout << "Evaluating expressions...";
 
   // Loop over all statements, and fill global eval context if relevant
   for (vector<Statement *>::const_iterator it = statements.begin(); it != statements.end(); it++)
@@ -76,7 +77,8 @@ ModFile::evalAllExpressions(bool warn_uninit)
   // Evaluate model local variables
   dynamic_model.fillEvalContext(global_eval_context);
 
-  cout << "done" << endl;
+  if (!nopreprocessoroutput)
+    cout << "done" << endl;
 
   // Check if some symbols are not initialized, and give them a zero value then
   for (int id = 0; id <= symbol_table.maxID(); id++)
@@ -340,7 +342,7 @@ ModFile::checkPass(bool nostrict, bool stochastic)
 }
 
 void
-ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs)
+ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const bool nopreprocessoroutput)
 {
   // Save the original model (must be done before any model transformations by preprocessor)
   // - except adl and diff which we always want expanded
@@ -408,7 +410,7 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs)
       if (linear)
         dynamic_model.cloneDynamic(orig_ramsey_dynamic_model);
       dynamic_model.cloneDynamic(ramsey_FOC_equations_dynamic_model);
-      ramsey_FOC_equations_dynamic_model.computeRamseyPolicyFOCs(*planner_objective);
+      ramsey_FOC_equations_dynamic_model.computeRamseyPolicyFOCs(*planner_objective, nopreprocessoroutput);
       ramsey_FOC_equations_dynamic_model.replaceMyEquations(dynamic_model);
       mod_file_struct.ramsey_eq_nbr = dynamic_model.equation_number() - mod_file_struct.orig_eq_nbr;
     }
@@ -517,13 +519,14 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs)
       exit(EXIT_FAILURE);
     }
 
-  if (!mod_file_struct.ramsey_model_present)
+  if (!nopreprocessoroutput)
-    cout << "Found " << dynamic_model.equation_number() << " equation(s)." << endl;
+    if (!mod_file_struct.ramsey_model_present)
-  else
+      cout << "Found " << dynamic_model.equation_number() << " equation(s)." << endl;
-    {
+    else
-      cout << "Found " << mod_file_struct.orig_eq_nbr  << " equation(s)." << endl;
+      {
-      cout << "Found " << dynamic_model.equation_number() << " FOC equation(s) for Ramsey Problem." << endl;
+        cout << "Found " << mod_file_struct.orig_eq_nbr  << " equation(s)." << endl;
-    }
+        cout << "Found " << dynamic_model.equation_number() << " FOC equation(s) for Ramsey Problem." << endl;
+      }
 
   if (symbol_table.exists("dsge_prior_weight"))
     if (mod_file_struct.bayesian_irf_present)
@@ -545,7 +548,7 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs)
 }
 
 void
-ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_derivs_order)
+ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_derivs_order, const bool nopreprocessoroutput)
 {
   // Mod file may have no equation (for example in a standalone BVAR estimation)
   if (dynamic_model.equation_number() > 0)
@@ -569,7 +572,7 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_deri
           if (mod_file_struct.identification_present || mod_file_struct.estimation_analytic_derivation)
             paramsDerivsOrder = params_derivs_order;
           static_model.computingPass(global_eval_context, no_tmp_terms, static_hessian,
-                                     false, paramsDerivsOrder, block, byte_code);
+                                     false, paramsDerivsOrder, block, byte_code, nopreprocessoroutput);
         }
       // Set things to compute for dynamic model
       if (mod_file_struct.perfect_foresight_solver_present || mod_file_struct.check_present
@@ -579,7 +582,7 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_deri
           || mod_file_struct.calib_smoother_present)
         {
           if (mod_file_struct.perfect_foresight_solver_present)
-            dynamic_model.computingPass(true, false, false, none, global_eval_context, no_tmp_terms, block, use_dll, byte_code);
+            dynamic_model.computingPass(true, false, false, none, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput);
           else
             {
               if (mod_file_struct.stoch_simul_present
@@ -604,13 +607,13 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_deri
               int paramsDerivsOrder = 0;
               if (mod_file_struct.identification_present || mod_file_struct.estimation_analytic_derivation)
                 paramsDerivsOrder = params_derivs_order;
-              dynamic_model.computingPass(true, hessian, thirdDerivatives, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code);
+              dynamic_model.computingPass(true, hessian, thirdDerivatives, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput);
               if (linear && mod_file_struct.ramsey_model_present)
-                orig_ramsey_dynamic_model.computingPass(true, true, false, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code);
+                orig_ramsey_dynamic_model.computingPass(true, true, false, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput);
             }
         }
       else // No computing task requested, compute derivatives up to 2nd order by default
-        dynamic_model.computingPass(true, true, false, none, global_eval_context, no_tmp_terms, block, use_dll, byte_code);
+        dynamic_model.computingPass(true, true, false, none, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput);
 
       if ((linear && !mod_file_struct.ramsey_model_present && !dynamic_model.checkHessianZero())
           || (linear && mod_file_struct.ramsey_model_present && !orig_ramsey_dynamic_model.checkHessianZero()))
@@ -646,6 +649,7 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
 #if defined(_WIN32) || defined(__CYGWIN32__)
                           , bool cygwin, bool msvc, bool mingw
 #endif
+                          , const bool nopreprocessoroutput
                           ) const
 {
   ofstream mOutputFile;
@@ -718,7 +722,8 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
   if (param_used_with_lead_lag)
     mOutputFile << "M_.parameter_used_with_lead_lag = true;" << endl;
 
-  cout << "Processing outputs ..." << endl;
+  if (!nopreprocessoroutput)
+    cout << "Processing outputs ..." << endl;
 
   symbol_table.writeOutput(mOutputFile);
 
@@ -960,11 +965,12 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
       steady_state_model.writeSteadyStateFile(basename, mod_file_struct.ramsey_model_present, false);
     }
 
-  cout << "done" << endl;
+  if (!nopreprocessoroutput)
+    cout << "done" << endl;
 }
 
 void
-ModFile::writeExternalFiles(const string &basename, FileOutputType output, LanguageOutputType language) const
+ModFile::writeExternalFiles(const string &basename, FileOutputType output, LanguageOutputType language, const bool nopreprocessoroutput) const
 {
   switch (language)
     {
@@ -975,7 +981,7 @@ ModFile::writeExternalFiles(const string &basename, FileOutputType output, Langu
       writeExternalFilesCC(basename, output);
       break;
     case julia:
-      writeExternalFilesJulia(basename, output);
+      writeExternalFilesJulia(basename, output, nopreprocessoroutput);
       break;
     default:
       cerr << "This case shouldn't happen. Contact the authors of Dynare" << endl;
@@ -1193,7 +1199,7 @@ ModFile::writeModelCC(const string &basename) const
 }
 
 void
-ModFile::writeExternalFilesJulia(const string &basename, FileOutputType output) const
+ModFile::writeExternalFilesJulia(const string &basename, FileOutputType output, const bool nopreprocessoroutput) const
 {
   ofstream jlOutputFile;
   if (basename.size())
@@ -1270,7 +1276,8 @@ ModFile::writeExternalFilesJulia(const string &basename, FileOutputType output)
     jlOutputFile << "model_.h = zeros(Float64, 1, 1)" << endl
                  << "model_.correlation_matrix_me = ones(Float64, 1, 1)" << endl;
 
-  cout << "Processing outputs ..." << endl;
+  if (!nopreprocessoroutput)
+    cout << "Processing outputs ..." << endl;
   symbol_table.writeJuliaOutput(jlOutputFile);
 
   if (dynamic_model.equation_number() > 0)
@@ -1314,11 +1321,12 @@ ModFile::writeExternalFilesJulia(const string &basename, FileOutputType output)
                << "end" << endl
                << "end" << endl;
   jlOutputFile.close();
-  cout << "done" << endl;
+  if (!nopreprocessoroutput)
+    cout << "done" << endl;
 }
 
 void
-ModFile::writeJsonOutput(const string &basename, JsonOutputPointType json, JsonFileOutputType json_output_mode, bool onlyjson, bool jsonderivsimple)
+ModFile::writeJsonOutput(const string &basename, JsonOutputPointType json, JsonFileOutputType json_output_mode, bool onlyjson, const bool nopreprocessoroutput, bool jsonderivsimple)
 {
   if (json == nojson)
     return;
@@ -1342,24 +1350,25 @@ ModFile::writeJsonOutput(const string &basename, JsonOutputPointType json, JsonF
     cout << "}" << endl
          << "//-- END JSON --// " << endl;
 
-  switch (json)
+  if (!nopreprocessoroutput)
-    {
+    switch (json)
-    case parsing:
+      {
-      cout << "JSON written after Parsing step." << endl;
+      case parsing:
-      break;
+        cout << "JSON written after Parsing step." << endl;
-    case checkpass:
+        break;
-      cout << "JSON written after Check step." << endl;
+      case checkpass:
-      break;
+        cout << "JSON written after Check step." << endl;
-    case transformpass:
+        break;
-      cout << "JSON written after Transform step." << endl;
+      case transformpass:
-      break;
+        cout << "JSON written after Transform step." << endl;
-    case computingpass:
+        break;
-      cout << "JSON written after Computing step." << endl;
+      case computingpass:
-      break;
+        cout << "JSON written after Computing step." << endl;
-    case nojson:
+        break;
-      cerr << "ModFile::writeJsonOutput: should not arrive here." << endl;
+      case nojson:
-      exit(EXIT_FAILURE);
+        cerr << "ModFile::writeJsonOutput: should not arrive here." << endl;
-    }
+        exit(EXIT_FAILURE);
+      }
 
   if (onlyjson)
     exit(EXIT_SUCCESS);

preprocessor/ModFile.hh

@@ -128,17 +128,17 @@ public:
   void addStatementAtFront(Statement *st);
   //! Evaluate all the statements
   /*! \param warn_uninit Should a warning be displayed for uninitialized endogenous/exogenous/parameters ? */
-  void evalAllExpressions(bool warn_uninit);
+  void evalAllExpressions(bool warn_uninit, const bool nopreprocessoroutput);
   //! Do some checking and fills mod_file_struct
   /*! \todo add check for number of equations and endogenous if ramsey_policy is present */
   void checkPass(bool nostrict, bool stochastic);
   //! Perform some transformations on the model (creation of auxiliary vars and equations)
   /*! \param compute_xrefs if true, equation cross references will be computed */
-  void transformPass(bool nostrict, bool stochastic, bool compute_xrefs);
+  void transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const bool nopreprocessoroutput);
   //! Execute computations
   /*! \param no_tmp_terms if true, no temporary terms will be computed in the static and dynamic files */
   /*! \param params_derivs_order compute this order of derivs wrt parameters */
-  void computingPass(bool no_tmp_terms, FileOutputType output, int params_derivs_order);
+  void computingPass(bool no_tmp_terms, FileOutputType output, int params_derivs_order, const bool nopreprocessoroutput);
   //! Writes Matlab/Octave output files
   /*!
     \param basename The base name used for writing output files. Should be the name of the mod file without its extension
@@ -157,11 +157,12 @@ public:
 #if defined(_WIN32) || defined(__CYGWIN32__) || defined(__MINGW32__)
                         , bool cygwin, bool msvc, bool mingw
 #endif
+                        , const bool nopreprocessoroutput
                         ) const;
-  void writeExternalFiles(const string &basename, FileOutputType output, LanguageOutputType language) const;
+  void writeExternalFiles(const string &basename, FileOutputType output, LanguageOutputType language, const bool nopreprocessoroutput) const;
   void writeExternalFilesC(const string &basename, FileOutputType output) const;
   void writeExternalFilesCC(const string &basename, FileOutputType output) const;
-  void writeExternalFilesJulia(const string &basename, FileOutputType output) const;
+  void writeExternalFilesJulia(const string &basename, FileOutputType output, const bool nopreprocessoroutput) const;
   //! Writes C output files only => No further Matlab processing
   void writeCOutputFiles(const string &basename) const;
   void writeModelC(const string &basename) const;
@@ -174,7 +175,7 @@ public:
   //! Initially created to enable Julia to work with .mod files
   //! Potentially outputs ModFile after the various parts of processing (parsing, checkPass, transformPass, computingPass)
   //! Allows user of other host language platforms (python, fortran, etc) to provide support for dynare .mod files
-  void writeJsonOutput(const string &basename, JsonOutputPointType json, JsonFileOutputType json_output_mode, bool onlyjson, bool jsonderivsimple = false);
+  void writeJsonOutput(const string &basename, JsonOutputPointType json, JsonFileOutputType json_output_mode, bool onlyjson, const bool nopreprocessoroutput, bool jsonderivsimple = false);
 };
 
 #endif // ! MOD_FILE_HH

preprocessor/ParsingDriver.cc

@@ -746,6 +746,46 @@ ParsingDriver::homotopy_val(string *name, expr_t val1, expr_t val2)
   delete name;
 }
 
+void
+ParsingDriver::end_generate_irfs()
+{
+  mod_file->addStatement(new GenerateIRFsStatement(options_list, generate_irf_names, generate_irf_elements));
+
+  generate_irf_elements.clear();
+  generate_irf_names.clear();
+  options_list.clear();
+}
+
+void
+ParsingDriver::add_generate_irfs_element(string *name)
+{
+  for (vector<string>::const_iterator it = generate_irf_names.begin();
+       it != generate_irf_names.end(); it++)
+    if (*it == *name)
+      error("Names in the generate_irfs block must be unique but you entered '"
+            + *name + "' more than once.");
+
+  generate_irf_names.push_back(*name);
+  generate_irf_elements.push_back(generate_irf_exos);
+
+  generate_irf_exos.clear();
+
+  delete name;
+}
+
+void
+ParsingDriver::add_generate_irfs_exog_element(string *exo, string *value)
+{
+  check_symbol_is_exogenous(exo);
+  if (generate_irf_exos.find(*exo) != generate_irf_exos.end())
+    error("You have set the exogenous variable " + *exo + " twice.");
+
+  generate_irf_exos[*exo] = atof(value->c_str());
+
+  delete exo;
+  delete value;
+}
+
 void
 ParsingDriver::forecast()
 {
@@ -1809,6 +1849,21 @@ ParsingDriver::check_symbol_is_endogenous_or_exogenous(string *name)
     }
 }
 
+void
+ParsingDriver::check_symbol_is_exogenous(string *name)
+{
+  check_symbol_existence(*name);
+  int symb_id = mod_file->symbol_table.getID(*name);
+  switch (mod_file->symbol_table.getType(symb_id))
+    {
+    case eExogenous:
+    case eExogenousDet:
+      break;
+    default:
+      error(*name + " is not exogenous.");
+    }
+}
+
 void
 ParsingDriver::set_std_prior(string *name, string *subsample_name)
 {
@@ -3221,6 +3276,22 @@ ParsingDriver::perfect_foresight_solver()
   options_list.clear();
 }
 
+void
+ParsingDriver::gmm_estimation()
+{
+  mod_file->addStatement(new GMMEstimationStatement(symbol_list, options_list));
+  symbol_list.clear();
+  options_list.clear();
+}
+
+void
+ParsingDriver::smm_estimation()
+{
+  mod_file->addStatement(new SMMEstimationStatement(symbol_list, options_list));
+  symbol_list.clear();
+  options_list.clear();
+}
+
 void
 ParsingDriver::prior_posterior_function(bool prior_func)
 {

preprocessor/ParsingDriver.hh

@@ -93,6 +93,9 @@ private:
   //! Checks that a given symbol exists and is a endogenous or exogenous, and stops with an error message if it isn't
   void check_symbol_is_endogenous_or_exogenous(string *name);
 
+  //! Checks that a given symbol exists and is a exogenous, and stops with an error message if it isn't
+  void check_symbol_is_exogenous(string *name);
+
   //! Checks for symbol existence in model block. If it doesn't exist, an error message is stored to be printed at
   //! the end of the model block
   void check_symbol_existence_in_model_block(const string &name);
@@ -196,7 +199,10 @@ private:
       Ri_TYPE
     };
   SvarRestrictionType svar_restriction_type;
-
+  //! Temporary storage for generate_irfs
+  vector<string> generate_irf_names;
+  vector<map<string, double> > generate_irf_elements;
+  map<string, double> generate_irf_exos;
   //! Temporary storage for argument list of external function
   stack<vector<expr_t> >  stack_external_function_args;
   //! Temporary storage for parameters in joint prior statement
@@ -553,6 +559,10 @@ public:
   void add_lower_cholesky();
   //! Svar_Global_Identification_Check Statement
   void add_svar_global_identification_check();
+  //! generate_irfs Block
+  void end_generate_irfs();
+  void add_generate_irfs_element(string *name);
+  void add_generate_irfs_exog_element(string *exo, string *value);
   //! Forecast Statement
   void forecast();
   void set_trends();
@@ -810,6 +820,10 @@ public:
   void add_VAR_covariance_pair_restriction(string *name11, string *name12, string *name21, string *name22);
   //! Runs VAR estimation process
   void run_var_estimation();
+  //! GMM Estimation statement
+  void gmm_estimation();
+  //! SMM Estimation statement
+  void smm_estimation();
 };
 
 #endif // ! PARSING_DRIVER_HH

preprocessor/StaticModel.cc

@@ -1047,7 +1047,7 @@ StaticModel::collect_first_order_derivatives_endogenous()
 }
 
 void
-StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms, bool hessian, bool thirdDerivatives, int paramsDerivsOrder, bool block, bool bytecode)
+StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms, bool hessian, bool thirdDerivatives, int paramsDerivsOrder, bool block, bool bytecode, const bool nopreprocessoroutput)
 {
   initializeVariablesAndEquations();
 
@@ -1077,27 +1077,31 @@ StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms
     }
 
   // Launch computations
-  cout << "Computing static model derivatives:" << endl
+  if (!nopreprocessoroutput)
-       << " - order 1" << endl;
+    cout << "Computing static model derivatives:" << endl
+         << " - order 1" << endl;
   first_derivatives.clear();
 
   computeJacobian(vars);
 
   if (hessian)
     {
-      cout << " - order 2" << endl;
+      if (!nopreprocessoroutput)
+        cout << " - order 2" << endl;
       computeHessian(vars);
     }
 
   if (thirdDerivatives)
     {
-      cout << " - order 3" << endl;
+      if (!nopreprocessoroutput)
+        cout << " - order 3" << endl;
       computeThirdDerivatives(vars);
     }
 
   if (paramsDerivsOrder > 0)
     {
-      cout << " - derivatives of Jacobian/Hessian w.r. to parameters" << endl;
+      if (!nopreprocessoroutput)
+        cout << " - derivatives of Jacobian/Hessian w.r. to parameters" << endl;
       computeParamsDerivatives(paramsDerivsOrder);
 
       if (!no_tmp_terms)
@@ -1122,7 +1126,8 @@ StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms
 
       equation_type_and_normalized_equation = equationTypeDetermination(first_order_endo_derivatives, variable_reordered, equation_reordered, mfs);
 
-      cout << "Finding the optimal block decomposition of the model ...\n";
+      if (!nopreprocessoroutput)
+        cout << "Finding the optimal block decomposition of the model ...\n";
 
       lag_lead_vector_t equation_lag_lead, variable_lag_lead;
 
@@ -1550,6 +1555,7 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) c
 
       StaticOutput << "#=" << endl << comments.str() << "=#" << endl
                    << "  @assert size(g2) == (" << equations.size() << ", " << g2ncols << ")" << endl
+                   << "  fill!(g2, 0.0)" << endl
                    << "  static!(y, x, params, residual, g1)" << endl;
       if (second_derivatives.size())
         StaticOutput << model_local_vars_output.str()
@@ -1573,6 +1579,7 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) c
 
       StaticOutput << "#=" << endl << comments.str() << "=#" << endl
                    << "  @assert size(g3) == (" << nrows << ", " << ncols << ")" << endl
+                   << "  fill!(g3, 0.0)" << endl
                    << "  static!(y, x, params, residual, g1, g2)" << endl;
       if (third_derivatives.size())
         StaticOutput << model_local_vars_output.str()

preprocessor/StaticModel.hh

@@ -161,7 +161,7 @@ public:
     \param hessian whether 2nd derivatives w.r. to exo, exo_det and endo should be computed
     \param paramsDerivsOrder order of derivatives w.r. to a pair (endo/exo/exo_det, parameter) to be computed
   */
-  void computingPass(const eval_context_t &eval_context, bool no_tmp_terms, bool hessian, bool thirdDerivatices, int paramsDerivsOrder, bool block, bool bytecode);
+  void computingPass(const eval_context_t &eval_context, bool no_tmp_terms, bool hessian, bool thirdDerivatices, int paramsDerivsOrder, bool block, bool bytecode, const bool nopreprocessoroutput);
 
   //! Adds informations for simulation in a binary file for a block decomposed model
   void Write_Inf_To_Bin_File_Block(const string &static_basename, const string &bin_basename, const int &num,

tests/Makefile.am

@@ -251,7 +251,6 @@ MODFILES = \
 	stochastic-backward-models/solow_cd.mod \
 	stochastic-backward-models/solow_ces.mod \
 	stochastic-backward-models/solow_cd_with_steadystate.mod \
-	stochastic-backward-models/backward_linear.mod \
 	deterministic_simulations/purely_forward/ar1.mod \
 	deterministic_simulations/purely_forward/nk.mod \
 	deterministic_simulations/purely_backward/ar1.mod \

tests/TeX/fs2000_corr_ME.mod

@@ -182,6 +182,8 @@ trace_plot(options_,M_,estim_params_,'StructuralShock',1,'e_a')
 
 shock_decomposition y W R;
 
+stoch_simul(order=1,irf=20,graph_format=eps,periods=0,contemporaneous_correlation,conditional_variance_decomposition=[1,3]);
+
 collect_latex_files;
 
 //identification(advanced=1,max_dim_cova_group=3,prior_mc=250);

tests/conditional_variance_decomposition/example1.mod

@@ -41,6 +41,86 @@ var u; stderr 0.009;
 //var e, u = phi*0.009*0.009;
 end;
 
-stoch_simul(conditional_variance_decomposition = 100,irf=0);
+varobs a y;
+
+
+@#for order in [1,2]
+options_.order=@{order};
+stoch_simul(conditional_variance_decomposition = 100,irf=0) a y k;
+if max(abs(sum(oo_.variance_decomposition,2)-100))>1e-6
+    error(['Variance decomposition at order ',num2str(options_.order),' does not work'])
+end
 
 stoch_simul(conditional_variance_decomposition = [1 2 3 5 10 100],irf=0) a y k;
+if max(max(abs(sum(oo_.conditional_variance_decomposition,3)-1)))>1e-6
+    error(['Conditional variance decomposition at order ',num2str(options_.order),' does not work'])
+end
+
+shocks;
+var y; stderr 0.01;
+var a; stderr 0.009;
+end;
+
+stoch_simul(conditional_variance_decomposition = [1 2 3 5 10 200],irf=0) a y k;
+if max(max(abs(sum(oo_.conditional_variance_decomposition,3)-1)))>1e-6
+    error(['Conditional variance decomposition at order ',num2str(options_.order),' does not work'])
+end
+if max(max(abs(sum(oo_.conditional_variance_decomposition_ME,3)-1)))>1e-6
+    error(['Conditional variance decomposition at order ',num2str(options_.order),' does not work'])
+end
+
+nvar = size(var_list_,1);
+SubsetOfVariables=zeros(nvar,1);
+for i=1:nvar
+    i_tmp = strmatch(var_list_(i,:),M_.endo_names,'exact');
+    SubsetOfVariables(i) = i_tmp;
+end
+
+[observable_pos,index_observables,index_subset]=intersect(SubsetOfVariables,options_.varobs_id,'stable');
+y_pos=strmatch('y',var_list_,'exact');
+y_pos_varobs=strmatch('y',options_.varobs,'exact');
+a_pos_varobs=strmatch('a',options_.varobs,'exact');
+
+if (oo_.conditional_variance_decomposition_ME(index_observables(y_pos_varobs),end,end)+oo_.var(y_pos,y_pos)/(oo_.var(y_pos,y_pos)+M_.H(y_pos_varobs,y_pos_varobs)))-1>1e-5...
+        || abs(oo_.conditional_variance_decomposition_ME(index_observables(a_pos_varobs),1,end)-0.5)>1e-5    
+    error(['Conditional variance decomposition at order ',num2str(options_.order),' with ME does not work'])
+end
+
+if (oo_.variance_decomposition_ME(index_observables(y_pos_varobs),end)/100+oo_.var(y_pos,y_pos)/(oo_.var(y_pos,y_pos)+M_.H(y_pos_varobs,y_pos_varobs)))-1>1e-5   
+    error(['Unconditional variance decomposition at order ',num2str(options_.order),' with ME does not work'])
+end
+
+shocks;
+var y; stderr 0;
+end;
+
+@#endfor
+
+%% do simulated moments
+@#for order in [1,2]
+options_.order=@{order};
+
+shocks;
+var y; stderr 0.01;
+var a; stderr 0.009;
+end;
+
+stoch_simul(irf=0,periods=1000000) a y k;
+if max(abs(sum(oo_.variance_decomposition,2)-100))>2
+    error(['Variance decomposition at order ',num2str(options_.order),' does not work'])
+end
+
+[observable_pos,index_observables,index_subset]=intersect(SubsetOfVariables,options_.varobs_id,'stable');
+y_pos=strmatch('y',var_list_,'exact');
+y_pos_varobs=strmatch('y',options_.varobs,'exact');
+a_pos_varobs=strmatch('a',options_.varobs,'exact');
+
+if (oo_.variance_decomposition_ME(index_observables(y_pos_varobs),end)/100+oo_.var(y_pos,y_pos)/(oo_.var(y_pos,y_pos)+M_.H(y_pos_varobs,y_pos_varobs)))-1>5e-4
+    error(['Unconditional variance decomposition at order ',num2str(options_.order),' with ME does not work'])
+end
+
+shocks;
+var y; stderr 0;
+end;
+
+@#endfor

tests/example2.mod

@@ -41,3 +41,11 @@ var u = 0.009^2;
 end;
 
 stoch_simul(periods=2000, drop=200);
+
+%% test that simul_replic does not affect simulated moments
+moments_temp=oo_.var;
+set_dynare_seed('default');
+stoch_simul(periods=2000, drop=200,simul_replic=2);
+if ~isequal(moments_temp,oo_.var)
+    error('Simul_replic affects simulated moments')
+end

tests/stochastic-backward-models/backward_linear.mod

@@ -1,41 +0,0 @@
-var y gdp gdp_pot g pi P rs;
-varexo e_y e_g e_p;
-
-parameters alpha_1 alpha_2 beta_1 beta_2 theta gamma_1 gamma_2 pi_tar g_ss rr_bar;
-alpha_1 = 0.9;
-alpha_2 = -0.1;
-beta_1 = 0.8;
-beta_2 = 0.1;
-theta = 0.9;
-gamma_1 = 1.5;
-gamma_2 = 0.5;
-pi_tar = 2;
-g_ss = 0.005;
-rr_bar = 1;
-
-model(linear);
-gdp = gdp_pot + y;
-y = alpha_1*y(-1) + alpha_2*(rs - pi - rr_bar) + e_y;
-gdp_pot = gdp_pot(-1) + g;
-g = (1 - theta)*g_ss + theta*g(-1) + e_g;
-pi = (1 - beta_1)*pi_tar + beta_1*pi(-1) + beta_2*y + e_p;
-P = pi/400 + P(-1);
-rs = rr_bar + pi_tar + gamma_1*(pi(-1) - pi_tar) + gamma_2*y(-1);
-end;
-
-histval;
-gdp_pot(0) = 1;
-P(0) = 1;
-g(0) = g_ss;
-pi(0) = pi_tar;
-end;
-
-oo_.steadystate = NaN(7,1);
-
-oo_ = simul_backward_model(M_.endo_histval, 10, options_, M_, oo_, zeros(11,3));
-
-err1 = norm(abs(oo_.endo_simul([1 4 5 7],2:end) - repmat([0 0.005 2 3]',1,10)));
-err2 = norm(abs(oo_.endo_simul([2 3 6],2:end) - repmat(linspace(1.005,1.05,10),3,1))); 
-if err1 > 1e-14 || err2 > 1e-14;
-   error('Error in backward_linear.mod');   
-end;

tests/stochastic-backward-models/solow_cd.mod

@@ -54,4 +54,6 @@ shocks;
     var e_n = 0.001;
 end;
 
-oo_ = simul_backward_nonlinear_model([], 5000, options_, M_, oo_);
+initialconditions = dseries([1 1.02 1 1.02 1], 2000Q1, {'Efficiency'; 'EfficiencyGrowth'; 'Population'; 'PopulationGrowth'; 'PhysicalCapitalStock'});
+
+simulations = simul_backward_model(initialconditions, 5000);

tests/stochastic-backward-models/solow_cd_with_steadystate.mod

@@ -33,35 +33,31 @@ model;
     PhysicalCapitalStock = (1-delta)*PhysicalCapitalStock(-1) + s*Output;
 end;
 
-histval;
+d = dseries([.5 1 1 1.04 15], 2000Q1, {'Efficiency'; 'EfficiencyGrowth'; 'Population'; 'PopulationGrowth'; 'PhysicalCapitalStock'});
-    Efficiency(0) = .5;
-    EfficiencyGrowth(0) = 1.00;
-    Population(0) = 1;
-    PopulationGrowth(0) = 1.04;
-    PhysicalCapitalStock(0) = 15;
-end;
 
-LongRunEfficiency = M_.endo_histval(1)*M_.endo_histval(2)^(rho_x/(1-rho_x));
+LongRunEfficiency = d.Efficiency*d.EfficiencyGrowth^(rho_x/(1-rho_x));
-LongRunPopulation = M_.endo_histval(3)*M_.endo_histval(4)^(rho_n/(1-rho_n));
+LongRunPopulation = d.Population*d.PopulationGrowth^(rho_n/(1-rho_n));
 LongRunEfficiencyGrowth = EfficiencyGrowth_ss; 
 LongRunPopulationGrowth = PopulationGrowth_ss;
-LongRunIntensiveCapitalStock = LongRunEfficiencyGrowth*LongRunPopulationGrowth*(s/(LongRunEfficiencyGrowth*LongRunPopulationGrowth-1+delta))^(1/(1-alpha));
+LongRunIntensiveCapitalStock = (s/(LongRunEfficiencyGrowth*LongRunPopulationGrowth-1+delta))^(1/(1-alpha)); //LongRunEfficiencyGrowth*LongRunPopulationGrowth*
 
 precision = 1e-6;
 T = 5*floor(log(precision)/log(max(rho_x, rho_n)));
 
-oo_ = simul_backward_model(M_.endo_histval, T, options_, M_, oo_, zeros(T+1,2));
+e = dseries(zeros(T, 2), 2000Q2, {'e_x'; 'e_n'});
 
-if abs(oo_.endo_simul(1,end)-LongRunEfficiency)>1e-10
+simulations = simul_backward_model(d, T, e);
+
+if abs(simulations.Efficiency.data(end)-LongRunEfficiency.data)>1e-10
    error('Wrong long run level!')
 end
 
-if abs(oo_.endo_simul(3,end)-LongRunPopulation)>1e-10
+if abs(simulations.Population.data(end)-LongRunPopulation.data)>1e-10
    error('Wrong long run level!')
 end
 
-IntensiveCapitalStock = oo_.endo_simul(6,1:end)./(oo_.endo_simul(1,1:end).*oo_.endo_simul(3,1:end));
+IntensiveCapitalStock = simulations.PhysicalCapitalStock/(simulations.Efficiency*simulations.Population);
 
-if abs(IntensiveCapitalStock(end)-LongRunIntensiveCapitalStock>1e-10)
+if abs(IntensiveCapitalStock.data(end)-LongRunIntensiveCapitalStock)>1e-3
-   error('Wrong long run level!')
+    error('Wrong long run level!')
 end

tests/stochastic-backward-models/solow_ces.mod

@@ -35,17 +35,11 @@ model;
     PhysicalCapitalStock = (1-delta)*PhysicalCapitalStock(-1) + s*Output;
 end;
 
-histval;
+d = dseries([1 1.02 1 1.02 1], 2000Q1, {'Efficiency'; 'EfficiencyGrowth'; 'Population'; 'PopulationGrowth'; 'PhysicalCapitalStock'});
-    Efficiency(0) = 1;
-    EfficiencyGrowth(0) = 1.02;
-    Population(0) = 1;
-    PopulationGrowth(0) = 1.02;
-    PhysicalCapitalStock(0) = 1;
-end;
 
 shocks;
     var e_x = 0.005;
     var e_n = 0.001;
 end;
 
-oo_ = simul_backward_nonlinear_model([], 5000, options_, M_, oo_);
+simulations = simul_backward_model(d, 5000);