diff --git a/doc/dynare.texi b/doc/dynare.texi index 4fa3d0abe..368c641b6 100644 --- a/doc/dynare.texi +++ b/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/phpBB3, Dynare forums}. +@uref{http://www.dynare.org/DynareWiki,old Dynare wiki}, the +@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 -recommended way is to cite the present manual, as follows: +You should cite Dynare if you use it in your research. The recommended way to +do this is to cite the present manual as: @quotation -Stéphane Adjemian, Houtan Bastani, Michel Juillard, Frédéric Karamé, -Ferhat Mihoubi, George Perendia, Johannes Pfeifer, Marco Ratto and -Sébastien Villemot (2011), ``Dynare: Reference Manual, Version 4,'' -@i{Dynare Working Papers}, 1, CEPREMAP +Stéphane Adjemian, Houtan Bastani, Michel Juillard, Frédéric Karamé, Junior +Maih, Ferhat Mihoubi, George Perendia, Johannes Pfeifer, Marco Ratto and +Sébastien Villemot (2011), ``Dynare: Reference Manual, Version 4,'' @i{Dynare +Working Papers}, 1, CEPREMAP @end quotation -Note that citing the Dynare Reference Manual in your research is a -good way to help the Dynare project. +@noindent For convenience, you can copy and paste the following into your BibTeX file: -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}, -but currently there is no displayed output. Note that this option requires the -option @code{moments_varendo} to be specified. +@code{oo_.PosteriorTheoreticalMoments.dsge.ConditionalVarianceDecomposition}. +Note that this option requires the +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. diff --git a/matlab/GetOneDraw.m b/matlab/GetOneDraw.m index f8e2a0d63..e4fba3a43 100644 --- a/matlab/GetOneDraw.m +++ b/matlab/GetOneDraw.m @@ -1,10 +1,15 @@ -function [xparams, logpost] = GetOneDraw(type) -% function [xparams, logpost] = GetOneDraw(type) +function [xparams, logpost] = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_) +% 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 \ No newline at end of file diff --git a/matlab/PosteriorIRF.m b/matlab/PosteriorIRF.m index de0f625d2..476b721b6 100644 --- a/matlab/PosteriorIRF.m +++ b/matlab/PosteriorIRF.m @@ -174,7 +174,7 @@ localVars.type=type; if strcmpi(type,'posterior') while b=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); -sigma = transpose(chol(Sigma)); +if ~deterministicshockflag + 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)); diff --git a/matlab/check_posterior_analysis_data.m b/matlab/check_posterior_analysis_data.m index 75aa78a7f..38d10aaae 100644 --- a/matlab/check_posterior_analysis_data.m +++ b/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) diff --git a/matlab/compute_moments_varendo.m b/matlab/compute_moments_varendo.m index c072f9c10..fd2203e3f 100644 --- a/matlab/compute_moments_varendo.m +++ b/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 . -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 diff --git a/matlab/conditional_variance_decomposition.m b/matlab/conditional_variance_decomposition.m index 177684ce6..7b85af755 100644 --- a/matlab/conditional_variance_decomposition.m +++ b/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 -% -% [1] In this version, absence of measurement errors is assumed... +% ConditionalVarianceDecomposition_ME [double] [m h p] array, where +% m is equal to length(intersect(SubsetOfVariables,varobs)) +% 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 \ No newline at end of file diff --git a/matlab/conditional_variance_decomposition_ME_mc_analysis.m b/matlab/conditional_variance_decomposition_ME_mc_analysis.m new file mode 100644 index 000000000..0742e4547 --- /dev/null +++ b/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 . + +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 \ No newline at end of file diff --git a/matlab/conditional_variance_decomposition_mc_analysis.m b/matlab/conditional_variance_decomposition_mc_analysis.m index c320554eb..d66bd0623 100644 --- a/matlab/conditional_variance_decomposition_mc_analysis.m +++ b/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 diff --git a/matlab/disp_moments.m b/matlab/disp_moments.m index 6e55bb571..be56d2c70 100644 --- a/matlab/disp_moments.m +++ b/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 diff --git a/matlab/disp_th_moments.m b/matlab/disp_th_moments.m index f3227f4f4..5ef074780 100644 --- a/matlab/disp_th_moments.m +++ b/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 diff --git a/matlab/display_conditional_variance_decomposition.m b/matlab/display_conditional_variance_decomposition.m index f0c115ad4..dd029bde6 100644 --- a/matlab/display_conditional_variance_decomposition.m +++ b/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 . +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:M_.exo_nbr + for j=1:shock_number 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 \ No newline at end of file diff --git a/matlab/distributions/inverse_gamma_specification.m b/matlab/distributions/inverse_gamma_specification.m index 410d4e52b..a5b0a0a2c 100644 --- a/matlab/distributions/inverse_gamma_specification.m +++ b/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 -% The call to the matlab's implementation of the secant method is here for testing purpose and should not be used. This routine fails +% REMARKS +% 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 diff --git a/matlab/dsge_simulated_theoretical_conditional_variance_decomposition.m b/matlab/dsge_simulated_theoretical_conditional_variance_decomposition.m index 4a2604d5b..c3b1902fa 100644 --- a/matlab/dsge_simulated_theoretical_conditional_variance_decomposition.m +++ b/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_ME1 id = strfind(varargin, 'nopreprocessoroutput'); if ~all(cellfun(@isempty, id)) preprocessoroutput = false; - varargin(cellfun(@isempty, id) == 0) = []; end end diff --git a/matlab/dynare_config.m b/matlab/dynare_config.m index 8ebe8da50..77237955d 100644 --- a/matlab/dynare_config.m +++ b/matlab/dynare_config.m @@ -125,6 +125,11 @@ if ~exist('struct2array') p{end+1} = '/missing/struct2array'; end +% isfile is missing in Octave and Matlab0 ) && 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 diff --git a/matlab/dynare_estimation_init.m b/matlab/dynare_estimation_init.m index 3cbbf6d60..b87b76873 100644 --- a/matlab/dynare_estimation_init.m +++ b/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 - var_list_ = check_list_of_variables(options_, M_, var_list_); - options_.varlist = var_list_; + if isempty(var_list_) + var_list_ = check_list_of_variables(options_, M_, 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 diff --git a/matlab/dynare_identification.m b/matlab/dynare_identification.m index 7b5c1e556..bc16725fa 100644 --- a/matlab/dynare_identification.m +++ b/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') diff --git a/matlab/evaluate_likelihood.m b/matlab/evaluate_likelihood.m index a76534adf..20da8eb4c 100644 --- a/matlab/evaluate_likelihood.m +++ b/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 . -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; \ No newline at end of file diff --git a/matlab/evaluate_posterior_kernel.m b/matlab/evaluate_posterior_kernel.m index 997962127..c4132fb38 100644 --- a/matlab/evaluate_posterior_kernel.m +++ b/matlab/evaluate_posterior_kernel.m @@ -1,10 +1,16 @@ -function lpkern = evaluate_posterior_kernel(parameters,llik) -% Evaluate the prior density at parameters. +function lpkern = evaluate_posterior_kernel(parameters,M_,estim_params_,oo_,options_,bayestopt_,llik) +% 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 . -[ldens,parameters] = evaluate_prior(parameters); -if nargin==1 - llik = evaluate_likelihood(parameters); +[ldens,parameters] = evaluate_prior(parameters,M_,estim_params_,oo_,options_,bayestopt_); +if nargin==6 %llik provided as an input + llik = evaluate_likelihood(parameters,M_,estim_params_,oo_,options_,bayestopt_); end lpkern = ldens+llik; \ No newline at end of file diff --git a/matlab/evaluate_prior.m b/matlab/evaluate_prior.m index 01cede921..e340667c3 100644 --- a/matlab/evaluate_prior.m +++ b/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 . -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' diff --git a/matlab/evaluate_smoother.m b/matlab/evaluate_smoother.m index c6d3d7113..c22cc0635 100644 --- a/matlab/evaluate_smoother.m +++ b/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' diff --git a/matlab/execute_prior_posterior_function.m b/matlab/execute_prior_posterior_function.m index 7e77a6ae8..cb7884e83 100644 --- a/matlab/execute_prior_posterior_function.m +++ b/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 diff --git a/matlab/fjaco.m b/matlab/fjaco.m index 5f5de0783..9b6f68abc 100644 --- a/matlab/fjaco.m +++ b/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 % diff --git a/matlab/get_all_parameters.m b/matlab/get_all_parameters.m index 289b6a679..d2aa30c61 100644 --- a/matlab/get_all_parameters.m +++ b/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 diff --git a/matlab/get_posterior_parameters.m b/matlab/get_posterior_parameters.m index e10616d20..0417956ed 100644 --- a/matlab/get_posterior_parameters.m +++ b/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 field_1 [char] optional field like 'mle_'. +% o type [char] = 'mode' or 'mean'. +% 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 . -global estim_params_ oo_ options_ M_ - -if nargin<2 +if nargin<6 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 \ No newline at end of file diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index 762a23281..077ad3795 100644 --- a/matlab/global_initialization.m +++ b/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() diff --git a/matlab/imcforecast.m b/matlab/imcforecast.m index 66841319a..3f77a3f2f 100644 --- a/matlab/imcforecast.m +++ b/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(:); diff --git a/matlab/lmmcp/get_complementarity_conditions.m b/matlab/lmmcp/get_complementarity_conditions.m index fb48ece9a..8459c0514 100644 --- a/matlab/lmmcp/get_complementarity_conditions.m +++ b/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(k) = etags{i,1}; + eq_index(eq_nbr) = k; + 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(k) = etags{i,1}; + eq_index(eq_nbr) = k; + eq_index(k) = eq_nbr; else error(sprintf(['Complementarity condition %s can''t be ' ... 'parsed'],etags{i,3})) diff --git a/matlab/lyapunov_solver.m b/matlab/lyapunov_solver.m index 98d782cdf..2bafb065a 100644 --- a/matlab/lyapunov_solver.m +++ b/matlab/lyapunov_solver.m @@ -41,7 +41,7 @@ function P=lyapunov_solver(T,R,Q,DynareOptions) % --*-- Unitary tests --*-- % along with Dynare. If not, see . 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 \ No newline at end of file +%@eof:1 diff --git a/matlab/utilities/general/isfile.m b/matlab/missing/isfile/isfile.m similarity index 57% rename from matlab/utilities/general/isfile.m rename to matlab/missing/isfile/isfile.m index 59f118710..ece12f43f 100644 --- a/matlab/utilities/general/isfile.m +++ b/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 . -% 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) - % File has no extension. - [status, c] = fileattrib(b); - if status - a = ~c.directory; +if iscell(b) && length(b)>1 && isvector(b) + if all(cellfun(@ischar, b)) + n = length(b); + cellofstringflag = true; + 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 \ No newline at end of file diff --git a/matlab/mode_check.m b/matlab/mode_check.m index 2df26ff0c..3131f10cc 100644 --- a/matlab/mode_check.m +++ b/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 - l2 = min(BoundsInfo.ub(kk),(1+sign(x(kk))*ll)*x(kk)); %upper bound + if x(kk)~=0 + 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)) diff --git a/matlab/occbin/evaluate_model.m b/matlab/occbin/evaluate_model.m new file mode 100644 index 000000000..f01d58a08 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/matlab/occbin/get_coef.m b/matlab/occbin/get_coef.m new file mode 100644 index 000000000..f2dbde5a6 --- /dev/null +++ b/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); \ No newline at end of file diff --git a/matlab/occbin/get_complementarity_conditions.m b/matlab/occbin/get_complementarity_conditions.m new file mode 100644 index 000000000..cf5a40866 --- /dev/null +++ b/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 . + +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 + diff --git a/matlab/occbin/get_occbin_complementarity_conditions.m b/matlab/occbin/get_occbin_complementarity_conditions.m new file mode 100644 index 000000000..29c865587 --- /dev/null +++ b/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 . + +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 + diff --git a/matlab/occbin/get_occbin_constraints.m b/matlab/occbin/get_occbin_constraints.m new file mode 100644 index 000000000..6723f4542 --- /dev/null +++ b/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 . + +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 diff --git a/matlab/occbin/get_residuals.m b/matlab/occbin/get_residuals.m new file mode 100644 index 000000000..f86d05cfb --- /dev/null +++ b/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); \ No newline at end of file diff --git a/matlab/occbin/occbin.m b/matlab/occbin/occbin.m new file mode 100644 index 000000000..ea0780d21 --- /dev/null +++ b/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_'; \ No newline at end of file diff --git a/matlab/occbin/test_constraint.m b/matlab/occbin/test_constraint.m new file mode 100644 index 000000000..96e41839a --- /dev/null +++ b/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 \ No newline at end of file diff --git a/matlab/optimization/newrat.m b/matlab/optimization/newrat.m index c554a22dc..9f14d9f1c 100644 --- a/matlab/optimization/newrat.m +++ b/matlab/optimization/newrat.m @@ -154,7 +154,7 @@ while norm(gg)>gtol && check==0 && jit 1 fwrite(fh,y_,'float64'); end + if i==1 + y_out=y_; + end end if replic > 1 diff --git a/matlab/stoch_simul.m b/matlab/stoch_simul.m index 81f8b886a..67dcf2a85 100644 --- a/matlab/stoch_simul.m +++ b/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') diff --git a/matlab/utilities/dataset/makedataset.m b/matlab/utilities/dataset/makedataset.m index 85d86d841..c5c244b54 100644 --- a/matlab/utilities/dataset/makedataset.m +++ b/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 % ======= diff --git a/matlab/variance_decomposition_ME_mc_analysis.m b/matlab/variance_decomposition_ME_mc_analysis.m new file mode 100644 index 000000000..d5c873df1 --- /dev/null +++ b/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 . + +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 \ No newline at end of file diff --git a/matlab/variance_decomposition_mc_analysis.m b/matlab/variance_decomposition_mc_analysis.m index 4621f25d3..c3acd14a6 100644 --- a/matlab/variance_decomposition_mc_analysis.m +++ b/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!']) - return + if ~isequal(exo,'ME') + disp([ type '_analysis:: ' exo ' is not a declared exogenous variable!']) + end + return end var=deblank(var); diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc index 1af68165f..2daf075dc 100644 --- a/preprocessor/ComputingTasks.cc +++ b/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 &generate_irf_names_arg, + const vector > &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::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 m = generate_irf_elements[i]; + for (map::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 m = generate_irf_elements[i]; + size_t idx = 0; + for (map::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 << "}"; +} diff --git a/preprocessor/ComputingTasks.hh b/preprocessor/ComputingTasks.hh index 3f09c1a7c..8432be9fa 100644 --- a/preprocessor/ComputingTasks.hh +++ b/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 generate_irf_names; + const vector > generate_irf_elements; +public: + GenerateIRFsStatement(const OptionsList &options_list_arg, + const vector &generate_irf_names_arg, + const vector > &generate_irf_elements_arg); + virtual void writeOutput(ostream &output, const string &basename, bool minimal_workspace) const; + virtual void writeJsonOutput(ostream &output) const; +}; + #endif diff --git a/preprocessor/DynamicModel.cc b/preprocessor/DynamicModel.cc index aba9a8740..19072ce71 100644 --- a/preprocessor/DynamicModel.cc +++ b/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 > 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 - << " - order 1" << endl; + if (!nopreprocessoroutput) + 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; } - - cout << "Ramsey Problem: added " << i << " Multipliers." << endl; + if (!nopreprocessoroutput) + 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, set > &xrefs) const +{ + for (map, set >::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::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, set >::const_iterator it = xref_param.begin(); - it != xref_param.end(); it++) - { - if (it != xref_param.begin()) - output << ", "; - output << "{\"parameter\": \"" << symbol_table.getName(it->first.first) << "\"" - << ", \"equations\": ["; - for (set::const_iterator it1 = it->second.begin(); - it1 != it->second.end(); it1++) - { - if (it1 != it->second.begin()) - output << ", "; - output << *it1 + 1; - } - output << "]}"; - } + writeJsonXrefsHelper(output, xref_param); output << "]" << ", \"endogenous\": ["; - for (map, set >::const_iterator it = xref_endo.begin(); - 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::const_iterator it1 = it->second.begin(); - it1 != it->second.end(); it1++) - { - if (it1 != it->second.begin()) - output << ", "; - output << *it1 + 1; - } - output << "]}"; - } + writeJsonXrefsHelper(output, xref_endo); output << "]" << ", \"exogenous\": ["; - for (map, set >::const_iterator it = xref_exo.begin(); - 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::const_iterator it1 = it->second.begin(); - it1 != it->second.end(); it1++) - { - if (it1 != it->second.begin()) - output << ", "; - output << *it1 + 1; - } - output << "]}"; - } + writeJsonXrefsHelper(output, xref_exo); output << "]" << ", \"exogenous_deterministic\": ["; - for (map, set >::const_iterator it = xref_exo_det.begin(); - 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::const_iterator it1 = it->second.begin(); - it1 != it->second.end(); it1++) - { - if (it1 != it->second.begin()) - output << ", "; - output << *it1 + 1; - } - output << "]}"; - } + writeJsonXrefsHelper(output, xref_exo_det); output << "]}" << endl; } diff --git a/preprocessor/DynamicModel.hh b/preprocessor/DynamicModel.hh index 9bc5c2959..c9eaec7ff 100644 --- a/preprocessor/DynamicModel.hh +++ b/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, set > &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; diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy index a6bfaf79c..821e97f2b 100644 --- a/preprocessor/DynareBison.yy +++ b/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 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 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 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 SYMBOL_VEC %type 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(":"); diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll index 034bda4ae..668db31e4 100644 --- a/preprocessor/DynareFlex.ll +++ b/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 ms_variance_decomposition {BEGIN DYNARE_STATEMENT; return token::MS_VARIANCE_DECOMPOSITION;} conditional_forecast {BEGIN DYNARE_STATEMENT; return token::CONDITIONAL_FORECAST;} plot_conditional_forecast {BEGIN DYNARE_STATEMENT; return token::PLOT_CONDITIONAL_FORECAST;} +gmm_estimation {BEGIN DYNARE_STATEMENT; return token::GMM_ESTIMATION;} +smm_estimation {BEGIN DYNARE_STATEMENT; return token::SMM_ESTIMATION;} markov_switching {BEGIN DYNARE_STATEMENT; return token::MARKOV_SWITCHING;} 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 irf_calibration {BEGIN DYNARE_BLOCK; return token::IRF_CALIBRATION;} ramsey_constraints {BEGIN DYNARE_BLOCK; return token::RAMSEY_CONSTRAINTS;} restrictions {BEGIN DYNARE_BLOCK; return token::RESTRICTIONS;} +generate_irfs {BEGIN DYNARE_BLOCK; return token::GENERATE_IRFS;} /* For the semicolon after an "end" keyword */ ; {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 silent_optimizer {return token::SILENT_OPTIMIZER;} lmmcp {return token::LMMCP;} occbin {return token::OCCBIN;} +centered_moments {return token::CENTERED_MOMENTS; } +autolag {return token::AUTOLAG; } +recursive_order_estimation {return token::RECURSIVE_ORDER_ESTIMATION; } +bartlett_kernel_lag {return token::BARTLETT_KERNEL_LAG; } +optimal { + yylval->string_val = new string(yytext); + return token::OPTIMAL; +} +diagonal { + yylval->string_val = new string(yytext); + return token::DIAGONAL; +} +weighting_matrix {return token::WEIGHTING_MATRIX; } +penalized_estimator {return token::PENALIZED_ESTIMATOR; } +verbose {return token::VERBOSE; } +simulation_multiple {return token::SIMULATION_MULTIPLE; } +seed {return token::SEED; } +bounded_shock_support {return token::BOUNDED_SHOCK_SUPPORT; } +analytical_girf {return token::ANALYTICAL_GIRF; } +irf_in_percent {return token::IRF_IN_PERCENT; } +emas_girf {return token::EMAS_GIRF; } +emas_drop {return token::EMAS_DROP; } +emas_tolf {return token::EMAS_TOLF; } +emas_max_iter {return token::EMAS_MAX_ITER; } [\$][^$]*[\$] { 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 irf_plot_threshold {return token::IRF_PLOT_THRESHOLD;} no_homotopy {return token::NO_HOMOTOPY;} +stderr_multiples {return token::STDERR_MULTIPLES;} +diagonal_only {return token::DIAGONAL_ONLY;} equation {return token::EQUATION;} exclusion {return token::EXCLUSION;} lag {return token::LAG;} diff --git a/preprocessor/DynareMain.cc b/preprocessor/DynareMain.cc index d37b2314d..eb2891a6e 100644 --- a/preprocessor/DynareMain.cc +++ b/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 - << "Starting preprocessing of the model file ..." << endl; + if (!nopreprocessoroutput) + 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; diff --git a/preprocessor/DynareMain2.cc b/preprocessor/DynareMain2.cc index 6f0e11391..d61df45b5 100644 --- a/preprocessor/DynareMain2.cc +++ b/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; } diff --git a/preprocessor/ModFile.cc b/preprocessor/ModFile.cc index 2ab5a4b78..416de5043 100644 --- a/preprocessor/ModFile.cc +++ b/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::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) - 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; - } + if (!nopreprocessoroutput) + if (!mod_file_struct.ramsey_model_present) + 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; + } 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) - { - case parsing: - cout << "JSON written after Parsing step." << endl; - break; - case checkpass: - cout << "JSON written after Check step." << endl; - break; - case transformpass: - cout << "JSON written after Transform step." << endl; - break; - case computingpass: - cout << "JSON written after Computing step." << endl; - break; - case nojson: - cerr << "ModFile::writeJsonOutput: should not arrive here." << endl; - exit(EXIT_FAILURE); - } + if (!nopreprocessoroutput) + switch (json) + { + case parsing: + cout << "JSON written after Parsing step." << endl; + break; + case checkpass: + cout << "JSON written after Check step." << endl; + break; + case transformpass: + cout << "JSON written after Transform step." << endl; + break; + case computingpass: + cout << "JSON written after Computing step." << endl; + break; + case nojson: + cerr << "ModFile::writeJsonOutput: should not arrive here." << endl; + exit(EXIT_FAILURE); + } if (onlyjson) exit(EXIT_SUCCESS); diff --git a/preprocessor/ModFile.hh b/preprocessor/ModFile.hh index 216542cbf..ac678b5b5 100644 --- a/preprocessor/ModFile.hh +++ b/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 diff --git a/preprocessor/ParsingDriver.cc b/preprocessor/ParsingDriver.cc index fbade4c56..5d2639f00 100644 --- a/preprocessor/ParsingDriver.cc +++ b/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::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) { diff --git a/preprocessor/ParsingDriver.hh b/preprocessor/ParsingDriver.hh index 252d79ee4..e7f54b0a9 100644 --- a/preprocessor/ParsingDriver.hh +++ b/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 generate_irf_names; + vector > generate_irf_elements; + map generate_irf_exos; //! Temporary storage for argument list of external function stack > 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 diff --git a/preprocessor/StaticModel.cc b/preprocessor/StaticModel.cc index 37133f8a3..ea6ed05e1 100644 --- a/preprocessor/StaticModel.cc +++ b/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 - << " - order 1" << endl; + if (!nopreprocessoroutput) + 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() diff --git a/preprocessor/StaticModel.hh b/preprocessor/StaticModel.hh index 09abf2067..68b6db3d8 100644 --- a/preprocessor/StaticModel.hh +++ b/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, diff --git a/tests/Makefile.am b/tests/Makefile.am index 891eb42d2..65ef2be3b 100644 --- a/tests/Makefile.am +++ b/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 \ diff --git a/tests/TeX/fs2000_corr_ME.mod b/tests/TeX/fs2000_corr_ME.mod index 33123e54b..b91acf92c 100644 --- a/tests/TeX/fs2000_corr_ME.mod +++ b/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); diff --git a/tests/conditional_variance_decomposition/example1.mod b/tests/conditional_variance_decomposition/example1.mod index 84429c5b9..455489c83 100644 --- a/tests/conditional_variance_decomposition/example1.mod +++ b/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 \ No newline at end of file diff --git a/tests/example2.mod b/tests/example2.mod index 6022f3e8a..f3f7fd940 100644 --- a/tests/example2.mod +++ b/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 \ No newline at end of file diff --git a/tests/stochastic-backward-models/backward_linear.mod b/tests/stochastic-backward-models/backward_linear.mod deleted file mode 100644 index b3d514647..000000000 --- a/tests/stochastic-backward-models/backward_linear.mod +++ /dev/null @@ -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; diff --git a/tests/stochastic-backward-models/solow_cd.mod b/tests/stochastic-backward-models/solow_cd.mod index e1e54aa5e..7302199ee 100644 --- a/tests/stochastic-backward-models/solow_cd.mod +++ b/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); diff --git a/tests/stochastic-backward-models/solow_cd_with_steadystate.mod b/tests/stochastic-backward-models/solow_cd_with_steadystate.mod index 500cc137c..906bd8d46 100644 --- a/tests/stochastic-backward-models/solow_cd_with_steadystate.mod +++ b/tests/stochastic-backward-models/solow_cd_with_steadystate.mod @@ -33,35 +33,31 @@ model; PhysicalCapitalStock = (1-delta)*PhysicalCapitalStock(-1) + s*Output; end; -histval; - Efficiency(0) = .5; - EfficiencyGrowth(0) = 1.00; - Population(0) = 1; - PopulationGrowth(0) = 1.04; - PhysicalCapitalStock(0) = 15; -end; +d = dseries([.5 1 1 1.04 15], 2000Q1, {'Efficiency'; 'EfficiencyGrowth'; 'Population'; 'PopulationGrowth'; 'PhysicalCapitalStock'}); -LongRunEfficiency = M_.endo_histval(1)*M_.endo_histval(2)^(rho_x/(1-rho_x)); -LongRunPopulation = M_.endo_histval(3)*M_.endo_histval(4)^(rho_n/(1-rho_n)); +LongRunEfficiency = d.Efficiency*d.EfficiencyGrowth^(rho_x/(1-rho_x)); +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) - error('Wrong long run level!') +if abs(IntensiveCapitalStock.data(end)-LongRunIntensiveCapitalStock)>1e-3 + error('Wrong long run level!') end diff --git a/tests/stochastic-backward-models/solow_ces.mod b/tests/stochastic-backward-models/solow_ces.mod index 235ff37b8..71d90717f 100644 --- a/tests/stochastic-backward-models/solow_ces.mod +++ b/tests/stochastic-backward-models/solow_ces.mod @@ -35,17 +35,11 @@ model; PhysicalCapitalStock = (1-delta)*PhysicalCapitalStock(-1) + s*Output; end; -histval; - Efficiency(0) = 1; - EfficiencyGrowth(0) = 1.02; - Population(0) = 1; - PopulationGrowth(0) = 1.02; - PhysicalCapitalStock(0) = 1; -end; +d = dseries([1 1.02 1 1.02 1], 2000Q1, {'Efficiency'; 'EfficiencyGrowth'; 'Population'; 'PopulationGrowth'; 'PhysicalCapitalStock'}); 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); \ No newline at end of file