From 241fd074246269e94f6c0f8cfb3f04bf1fce0881 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Mon, 16 Sep 2013 17:19:56 +0200 Subject: [PATCH 01/80] Add Geweke 1992 convergence diagnostics --- matlab/McMCDiagnostics.m | 59 ++++++++++++++++-- matlab/dynare_estimation_1.m | 4 +- matlab/geweke_chi2_test.m | 74 ++++++++++++++++++++++ matlab/geweke_moments.m | 109 +++++++++++++++++++++++++++++++++ matlab/global_initialization.m | 4 ++ 5 files changed, 242 insertions(+), 8 deletions(-) create mode 100644 matlab/geweke_chi2_test.m create mode 100644 matlab/geweke_moments.m diff --git a/matlab/McMCDiagnostics.m b/matlab/McMCDiagnostics.m index d4a26cb8c..6491be936 100644 --- a/matlab/McMCDiagnostics.m +++ b/matlab/McMCDiagnostics.m @@ -1,4 +1,4 @@ -function McMCDiagnostics(options_, estim_params_, M_) +function oo_ = McMCDiagnostics(options_, estim_params_, M_, oo_) % function McMCDiagnostics % Computes convergence tests % @@ -8,7 +8,7 @@ function McMCDiagnostics(options_, estim_params_, M_) % M_ [structure] % % OUTPUTS -% none +% oo_ [structure] % % SPECIAL REQUIREMENTS % none @@ -38,10 +38,6 @@ MhDirectoryName = CheckPath('metropolis',M_.dname); TeX = options_.TeX; nblck = options_.mh_nblck; -% Brooks and Gelman tests need more than one block -if nblck == 1 - return; -end npar = estim_params_.nvx; npar = npar + estim_params_.nvn; npar = npar + estim_params_.ncx; @@ -56,6 +52,57 @@ NumberOfMcFilesPerBlock = size(dir([MhDirectoryName ,filesep, M_.fname '_mh*_blc % check if all previous files are there for block 1 check_presence_consecutive_MC_files(MhDirectoryName,M_.fname,1) +if nblck == 1 % Brooks and Gelman tests need more than one block + convergence_diagnostics_geweke=zeros(npar,4+2*length(options_.convergence.geweke.taper_steps)); + if any(options_.convergence.geweke.geweke_interval<0) || any(options_.convergence.geweke.geweke_interval>1) || length(any(options_.convergence.geweke.geweke_interval<0))~=2 ... + || (options_.convergence.geweke.geweke_interval(2)-options_.convergence.geweke.geweke_interval(1)<0) + fprintf('\nCONVERGENCE DIAGNOSTICS: Invalid option for geweke_interval. Using the default of [0.2 0.5].\n') + options_.convergence.geweke.geweke_interval=[0.2 0.5]; + end + first_obs_begin_sample = max(1,ceil(options_.mh_drop*options_.mh_replic)); + last_obs_begin_sample = first_obs_begin_sample+round(options_.convergence.geweke.geweke_interval(1)*options_.mh_replic*options_.mh_drop); + first_obs_end_sample = first_obs_begin_sample+round(options_.convergence.geweke.geweke_interval(2)*options_.mh_replic*options_.mh_drop); + param_name=[]; + for jj=1:npar + param_name = strvcat(param_name,get_the_name(jj,options_.TeX,M_,estim_params_,options_)); + end + fprintf('\nGeweke (1992) Convergence Tests, based on means of draws %d to %d vs %d to %d.\n',first_obs_begin_sample,last_obs_begin_sample,first_obs_end_sample,options_.mh_replic); + fprintf('p-values are for Chi2-test for equality of means.\n'); + Geweke_header={'Parameter', 'Post. Mean', 'Post. Std', 'p-val No Taper'}; + print_string=['%',num2str(size(param_name,2)+3),'s \t %12.3f \t %12.3f \t %12.3f']; + print_string_header=['%',num2str(size(param_name,2)+3),'s \t %12s \t %12s \t %12s']; + for ii=1:length(options_.convergence.geweke.taper_steps) + Geweke_header=[Geweke_header, ['p-val ' num2str(options_.convergence.geweke.taper_steps(ii)),'% Taper']]; + print_string=[print_string,'\t %12.3f']; + print_string_header=[print_string_header,'\t %12s']; + end + print_string=[print_string,'\n']; + print_string_header=[print_string_header,'\n']; + fprintf(print_string_header,Geweke_header{1,:}); + for jj=1:npar + startline=0; + for n = 1:NumberOfMcFilesPerBlock + load([MhDirectoryName '/' M_.fname '_mh',int2str(n),'_blck1.mat'],'x2'); + nx2 = size(x2,1); + param_draws(startline+(1:nx2),1) = x2(:,jj); + startline = startline + nx2; + end + [results_vec, results_struct] = geweke_moments(param_draws,options_); + convergence_diagnostics_geweke(jj,:)=results_vec; + + param_draws1 = param_draws(first_obs_begin_sample:last_obs_begin_sample,:); + param_draws2 = param_draws(first_obs_end_sample:end,:); + [results_vec1] = geweke_moments(param_draws1,options_); + [results_vec2] = geweke_moments(param_draws2,options_); + + results_struct = geweke_chi2_test(results_vec1,results_vec2,results_struct,options_); + eval(['oo_.convergence.geweke.',param_name(jj,:),'=results_struct;']) + fprintf(print_string,param_name(jj,:),results_struct.posteriormean,results_struct.posteriorstd,results_struct.prob_chi2_test) + end + skipline(2); + return; +end + for blck = 2:nblck tmp = size(dir([MhDirectoryName ,filesep, M_.fname '_mh*_blck' int2str(blck) '.mat']),1); if tmp~=NumberOfMcFilesPerBlock diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index 7e5071ded..73e24df17 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -611,8 +611,8 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ... CutSample(M_, options_, estim_params_); return else - if ~options_.nodiagnostic && options_.mh_replic > 2000 && options_.mh_nblck > 1 - McMCDiagnostics(options_, estim_params_, M_); + if ~options_.nodiagnostic && options_.mh_replic > 2000 + oo_= McMCDiagnostics(options_, estim_params_, M_,oo_); end %% Here i discard first half of the draws: CutSample(M_, options_, estim_params_); diff --git a/matlab/geweke_chi2_test.m b/matlab/geweke_chi2_test.m new file mode 100644 index 000000000..7ed383139 --- /dev/null +++ b/matlab/geweke_chi2_test.m @@ -0,0 +1,74 @@ +function results_struct = geweke_chi2_test(results1,results2,results_struct,options) +% results_struct = geweke_chi2_test(results1,results2,results_struct,options) +% PURPOSE: computes Geweke's chi-squared test for two sets of MCMC sample draws +% +% INPUTS +% results1 [1 by (4+n_taper*2) vector] vector with post. mean, +% std, NSE_iid, RNE_iid, and tapered NSE and RNE +% for chain part 1 +% results2 [1 by (4+n_taper*2) vector] vector with post. mean, +% std, NSE_iid, RNE_iid, and tapered NSE and RNE +% for chain part 2 +% results_struct [structure] results structure generated by geweke_moments +% Dynareoptions [structure] +% +% OUTPUTS +% results_struct [structure] containing the following fields: +% pooled_mean Pooled mean of the chain parts, weighted +% with precision +% rpooled_nse Pooled NSE of the chain parts, weighted +% with precision +% prob_chi2_test p-value of Chi2-test for equality of +% means in both chain parts +% ----------------------------------------------------------------- + +% +% SPECIAL REQUIREMENTS +% None. + +% Copyright (C) 2013 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . +% +% ------------------------------------------------ +% REFERENCES: Geweke (1992), `Evaluating the accuracy of sampling-based +% approaches to the calculation of posterior moments', in J.O. Berger, +% J.M. Bernardo, A.P. Dawid, and A.F.M. Smith (eds.) Proceedings of +% the Fourth Valencia International Meeting on Bayesian Statistics, +% pp. 169-194, Oxford University Press +% Geweke (1999): `Using simulation methods for Bayesian econometric models: +% Inference, development and communication', Econometric Reviews, 18(1), +% 1-73 + +% written by: Johannes Pfeifer, +% based on code by James P. LeSage, who in turn +% drew on MATLAB programs written by Siddartha Chib + +for k=1:length(options.convergence.geweke.taper_steps)+1; + NSE=[results1(:,3+(k-1)*2) results2(:,3+(k-1)*2)]; + means=[results1(:,1) results2(:,1)]; + diff_Means=means(:,1)-means(:,2); + sum_of_weights=sum(1./(NSE.^2),2); + pooled_mean=sum(means./(NSE.^2),2)./sum_of_weights; + pooled_NSE=1./sqrt(sum_of_weights); + + test_stat=diff_Means.^2./sum(NSE.^2,2); + p = 1-chi2cdf(test_stat,1); + results_struct.pooled_mean(:,k) = pooled_mean; + results_struct.pooled_nse(:,k) = pooled_NSE; + results_struct.prob_chi2_test(:,k) = p; +end; + diff --git a/matlab/geweke_moments.m b/matlab/geweke_moments.m new file mode 100644 index 000000000..242a9d550 --- /dev/null +++ b/matlab/geweke_moments.m @@ -0,0 +1,109 @@ +function [results_vec, results_struct] = geweke_moments(draws,Dynareoptions) +%[results_vec, results_struct] = geweke_moments(draws,Dynareoptions) +% PURPOSE: computes Gewke's convergence diagnostics NSE and RNE +% (numerical std error and relative numerical efficiencies) + +% INPUTS +% draws [ndraws by 1 vector] +% Dynareoptions [structure] +% +% OUTPUTS +% results_vec +% results_struct [structure] containing the following fields: +% posteriormean= posterior parameter mean +% posteriorstd = posterior standard deviation +% nse_iid = nse assuming no serial correlation for variable i +% rne_iid = rne assuming no serial correlation for variable i +% nse_x = nse using x% autocovariance tapered estimate +% rne_x = rne using x% autocovariance taper +% ----------------------------------------------------------------- + +% +% SPECIAL REQUIREMENTS +% None. + +% Copyright (C) 2013 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +% REFERENCES: Geweke (1992), `Evaluating the accuracy of sampling-based +% approaches to the calculation of posterior moments', in J.O. Berger, +% J.M. Bernardo, A.P. Dawid, and A.F.M. Smith (eds.) Proceedings of +% the Fourth Valencia International Meeting on Bayesian Statistics, +% pp. 169-194, Oxford University Press +% Geweke (1999): `Using simulation methods for Bayesian econometric models: +% Inference, development and communication', Econometric Reviews, 18(1), +% 1-73 +% ----------------------------------------------------------------- + +% written by: Johannes Pfeifer, +% based on code by James P. LeSage, who in turn +% drew on MATLAB programs written by Siddartha Chib + + +ndraw = size(draws,1); +n_groups=100; +taper_steps=Dynareoptions.convergence.geweke.taper_steps; +results_vec=zeros(1,4+2*length(taper_steps)); + +ns = floor(ndraw/n_groups); %step_size +n_draws_used = ns*n_groups; %effective number of draws used after rounding down + +window_means= zeros(n_groups,1); +window_uncentered_variances= zeros(n_groups,1); +for ig=1:n_groups; + window_means(ig,1)=sum(draws((ig-1)*ns+1:ig*ns,1))/ns; + window_uncentered_variances(ig,1)=sum(draws((ig-1)*ns+1:ig*ns,1).^2)/ns; +end; %for ig +total_mean=mean(window_means); +total_variance=mean(window_uncentered_variances)-total_mean^2; + +% save posterior means and std deviations to results_struct structure +results_vec(1,1)=total_mean; +results_vec(1,2)=sqrt(total_variance); +results_struct.posteriormean = total_mean; +results_struct.posteriorstd = results_vec(1,2); + +% numerical standard error assuming no serial correlation +NSE=std(draws(1:n_draws_used,1),1)/sqrt(n_draws_used); +% save to results_struct structure +results_vec(1,3)=NSE; +results_vec(1,4)=total_variance/(n_draws_used*NSE^2); +results_struct.nse_iid = NSE; +results_struct.rne_iid = results_vec(1,4); + +%get autocovariance of grouped means +centered_window_means=window_means-total_mean; +autocov_grouped_means=zeros(n_groups,1); +for lag=0:n_groups-1; + autocov_grouped_means(lag+1)=centered_window_means(lag+1:n_groups,1)'*centered_window_means(1:n_groups-lag,1)/100; +end; + +% numerical standard error with tapered autocovariance functions +for taper_index=1:length(taper_steps) + taper=taper_steps(taper_index); + taper_lags=(1:taper-1)'; + taper_lag_weight=1-taper_lags/taper; + tapered_sum_of_covariances=autocov_grouped_means(1)+sum(2*taper_lag_weight.*autocov_grouped_means(1+taper_lags)); + NSE_taper=sqrt(tapered_sum_of_covariances/n_groups); + % save results_struct in structure + results_vec(1,4+taper_index*2-1)=NSE_taper; + results_vec(1,4+taper_index*2)=total_variance/(n_draws_used*NSE_taper^2); + + eval(['results_struct.nse_taper_',num2str(taper),'= NSE_taper;']); + eval(['results_struct.rne_taper_',num2str(taper),'= total_variance/(n_draws_used*NSE_taper^2);']); +end; % end of for mm loop + diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index d06a054d9..c067ef33f 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -577,6 +577,10 @@ options_.osr.verbose=2; % use GPU options_.gpu = 0; +%Geweke convergence diagnostics +options_.convergence.geweke.taper_steps=[4 8 15]; +options_.convergence.geweke.geweke_interval=[0.2 0.5]; + % initialize persistent variables in priordens() priordens([],[],[],[],[],[],1); % initialize persistent variables in dyn_first_order_solver() From 27f1858f1704d64137cc9a1f1501c3cd47613504 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Mon, 16 Sep 2013 19:16:52 +0200 Subject: [PATCH 02/80] Adds documentation of Geweke convergence diagnostics --- doc/dynare.texi | 95 ++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 87 insertions(+), 8 deletions(-) diff --git a/doc/dynare.texi b/doc/dynare.texi index 91758669e..17495f3fa 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -4210,12 +4210,19 @@ graphs of smoothed shocks, smoothed observation errors, smoothed and historical @algorithmshead -The Monte Carlo Markov Chain (MCMC) univariate diagnostics are generated -by the estimation command if @ref{mh_nblocks} is larger than 1, if -@ref{mh_replic} is larger than 2000, and if option @ref{nodiagnostic} is -not used. As described in section 3 of @cite{Brooks and Gelman (1998)} -the convergence diagnostics are based on comparing pooled and within -MCMC moments (Dynare displays the second and third order moments, and +The Monte Carlo Markov Chain (MCMC) diagnostics are generated +by the estimation command if @ref{mh_replic} is larger than 2000 and if +option @ref{nodiagnostic} is not used. If @ref{mh_nblocks} is equal to one, +the convergence diagnostics of @cite{Geweke (1992,1999)} is computed. It uses a +chi square test to compare the means of the first and last draws specified in +@ref{geweke_interval} (@pxref{geweke_interval}) after discarding the burnin of @ref{mh_drop}. The test is +computed using variance estimates under the assumption of no serial correlation +as well as using tapering windows specified in @ref{taper_steps} (@pxref{taper_steps}). +If @ref{mh_nblocks} is larger than 1, the convergence diagnostics of +@cite{Brooks and Gelman (1998)} are used instead. +As described in section 3 of @cite{Brooks and Gelman (1998)} the univariate +convergence diagnostics are based on comparing pooled and within MCMC moments +(Dynare displays the second and third order moments, and the length of the Highest Probability Density interval covering 80% of the posterior distribution). Due to computational reasons, the multivariate convergence diagnostic does not follow @cite{Brooks and @@ -4356,8 +4363,8 @@ the total number of Metropolis draws available. Default: @code{2} @item mh_drop = @var{DOUBLE} -The fraction of initially generated parameter vectors to be dropped -before using posterior simulations. Default: @code{0.5} +@anchor{mh_drop} +The fraction of initially generated parameter vectors to be dropped as a burnin before using posterior simulations. Default: @code{0.5} @item mh_jscale = @var{DOUBLE} The scale to be used for the jumping distribution in @@ -4756,6 +4763,18 @@ Value used to test if a generalized eigenvalue is 0/0 in the generalized Schur decomposition (in which case the model does not admit a unique solution). Default: @code{1e-6}. +@item taper_steps = [@var{INTEGER1} @var{INTEGER2} @dots{}] +@anchor{taper_steps} +Percent tapering used for the spectral window in the @cite{Geweke (1992,1999)} +convergence diagnostics (requires @ref{mh_nblocks}=1). The tapering is used to +take the serial correlation of the posterior draws into account. Default: @code{[4 8 15]}. + +@item geweke_interval = [@var{double} @var{double}] +@anchor{geweke_interval} +Percentage of MCMC draws at the beginning and end of the MCMC chain taken +to compute the @cite{Geweke (1992,1999)} convergence diagnostics (requires @ref{mh_nblocks}=1) +after discarding the first @ref{mh_drop} percent of draws as a burnin. Default: @code{[0.2 0.5]}. + @end table @customhead{Note} @@ -5046,6 +5065,56 @@ Upper/lower bound of the 90\% HPD interval taking into account both parameter an @end defvr +@defvr {MATLAB/Octave variable} oo_.convergence.geweke +@anchor{convergence.geweke} +Variable set by the convergence diagnostics of the @code{estimation} command when used with @ref{mh_nblocks}=1 option (@pxref{mh_nblocks}). + +Fields are of the form: +@example +@code{oo_.convergence.geweke.@var{VARIABLE_NAME}.@var{DIAGNOSTIC_OBJECT}} +@end example +where @var{DIAGNOSTIC_OBJECT} is one of the following: + +@table @code + +@item posteriormean +Mean of the posterior parameter distribution + +@item posteriorstd +Standard deviation of the posterior parameter distribution + +@item nse_iid +Numerical standard error (NSE) under the assumption of iid draws + +@item rne_iid +Relative numerical efficiency (RNE) under the assumption of iid draws + +@item nse_x +Numerical standard error (NSE) when using an x% taper + +@item rne_x +Relative numerical efficiency (RNE) when using an x% taper + +@item pooled_mean +Mean of the parameter when pooling the beginning and end parts of the chain +specified in @ref{geweke_interval} and weighting them with their relative precision. +It is a vector containing the results under the iid assumption followed by the ones +using the @ref{taper_steps} (@pxref{taper_steps}). + +@item pooled_nse +NSE of the parameter when pooling the beginning and end parts of the chain and weighting them with their relative precision. See @code{pooled_mean} + +@item prob_chi2_test +p-value of a chi squared test for equality of means in the beginning and the end +of the MCMC chain. See @code{pooled_mean}. A value above 0.05 indicates that +the null hypothesis of equal means and thus convergence cannot be rejected +at the 5 percent level. Differing values along the @ref{taper_steps} signal +the presence of significant autocorrelation in draws. In this case, the +estimates using a higher tapering are usually more reliable. + +@end table +@end defvr + @deffn Command model_comparison @var{FILENAME}[(@var{DOUBLE})]@dots{}; @deffnx Command model_comparison (marginal_density = laplace | modifiedharmonicmean) @var{FILENAME}[(@var{DOUBLE})]@dots{}; @@ -8686,6 +8755,16 @@ Fernández-Villaverde, Jesús and Juan Rubio-Ramírez (2005): ``Estimating Dynamic Equilibrium Economies: Linear versus Nonlinear Likelihood,'' @i{Journal of Applied Econometrics}, 20, 891--910 +@item +Geweke, John (1992): ``Evaluating the accuracy of sampling-based approaches +to the calculation of posterior moments'', in J.O. Berger, J.M. Bernardo, +A.P. Dawid, and A.F.M. Smith (eds.) Proceedings of the Fourth Valencia +International Meeting on Bayesian Statistics, pp. 169--194, Oxford University Press + +@item +Geweke, John (1999): ``Using simulation methods for Bayesian econometric models: +Inference, development and communication,'' @i{Econometric Reviews}, 18(1), 1--73 + @item Ireland, Peter (2004): ``A Method for Taking Models to the Data,'' @i{Journal of Economic Dynamics and Control}, 28, 1205--26 From 5058e4d00f0ee7d3fcf3e748e691a1fcad2a9ae7 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Mon, 23 Sep 2013 14:31:43 +0200 Subject: [PATCH 03/80] fix typo --- doc/dynare.texi | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/dynare.texi b/doc/dynare.texi index 38f90a9cc..090b805cc 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -4216,10 +4216,10 @@ The Monte Carlo Markov Chain (MCMC) diagnostics are generated by the estimation command if @ref{mh_replic} is larger than 2000 and if option @ref{nodiagnostic} is not used. If @ref{mh_nblocks} is equal to one, the convergence diagnostics of @cite{Geweke (1992,1999)} is computed. It uses a -chi square test to compare the means of the first and last draws specified in -@ref{geweke_interval} (@pxref{geweke_interval}) after discarding the burnin of @ref{mh_drop}. The test is +chi square test to compare the means of the first and last draws specified by +@ref{geweke_interval} after discarding the burnin of @ref{mh_drop}. The test is computed using variance estimates under the assumption of no serial correlation -as well as using tapering windows specified in @ref{taper_steps} (@pxref{taper_steps}). +as well as using tapering windows specified in @ref{taper_steps}. If @ref{mh_nblocks} is larger than 1, the convergence diagnostics of @cite{Brooks and Gelman (1998)} are used instead. As described in section 3 of @cite{Brooks and Gelman (1998)} the univariate @@ -4771,7 +4771,7 @@ Percent tapering used for the spectral window in the @cite{Geweke (1992,1999)} convergence diagnostics (requires @ref{mh_nblocks}=1). The tapering is used to take the serial correlation of the posterior draws into account. Default: @code{[4 8 15]}. -@item geweke_interval = [@var{double} @var{double}] +@item geweke_interval = [@var{DOUBLE} @var{DOUBLE}] @anchor{geweke_interval} Percentage of MCMC draws at the beginning and end of the MCMC chain taken to compute the @cite{Geweke (1992,1999)} convergence diagnostics (requires @ref{mh_nblocks}=1) From 30cb09304818fb312d3913b0d9cebdb6dcd259cc Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Mon, 23 Sep 2013 14:27:43 +0200 Subject: [PATCH 04/80] front end for Geweke convergence diagnostics for single chains --- preprocessor/DynareBison.yy | 8 ++++++-- preprocessor/DynareFlex.ll | 2 ++ 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy index 5207b3a04..dd2192a79 100644 --- a/preprocessor/DynareBison.yy +++ b/preprocessor/DynareBison.yy @@ -113,7 +113,7 @@ class ParsingDriver; %token LYAPUNOV_FIXED_POINT_TOL LYAPUNOV_DOUBLING_TOL LYAPUNOV_SQUARE_ROOT_SOLVER_TOL LOG_DEFLATOR LOG_TREND_VAR LOG_GROWTH_FACTOR MARKOWITZ MARGINAL_DENSITY MAX MAXIT %token MFS MH_DROP MH_INIT_SCALE MH_JSCALE MH_MODE MH_NBLOCKS MH_REPLIC MH_RECOVER MIN MINIMAL_SOLVING_PERIODS SOLVE_MAXIT %token MODE_CHECK MODE_CHECK_NEIGHBOURHOOD_SIZE MODE_CHECK_SYMMETRIC_PLOTS MODE_CHECK_NUMBER_OF_POINTS MODE_COMPUTE MODE_FILE MODEL MODEL_COMPARISON MODEL_INFO MSHOCKS ABS SIGN -%token MODEL_DIAGNOSTICS MODIFIEDHARMONICMEAN MOMENTS_VARENDO DIFFUSE_FILTER SUB_DRAWS +%token MODEL_DIAGNOSTICS MODIFIEDHARMONICMEAN MOMENTS_VARENDO DIFFUSE_FILTER SUB_DRAWS TAPER_STEPS GEWEKE_INTERVAL %token NAME %token NAN_CONSTANT NO_STATIC NOBS NOCONSTANT NODISPLAY NOCORR NODIAGNOSTIC NOFUNCTIONS %token NOGRAPH NOMOMENTS NOPRINT NORMAL_PDF @@ -1555,7 +1555,9 @@ estimation_options : o_datafile | o_ar | o_endogenous_prior | o_use_univariate_filters_if_singularity_is_detected - | o_qz_zero_threshold + | o_qz_zero_threshold + | o_taper_steps + | o_geweke_interval ; list_optim_option : QUOTED_STRING COMMA QUOTED_STRING @@ -2388,6 +2390,8 @@ o_noprint : NOPRINT { driver.option_num("noprint", "1"); }; o_xls_sheet : XLS_SHEET EQUAL symbol { driver.option_str("xls_sheet", $3); }; o_xls_range : XLS_RANGE EQUAL range { driver.option_str("xls_range", $3); }; o_filter_step_ahead : FILTER_STEP_AHEAD EQUAL vec_int { driver.option_vec_int("filter_step_ahead", $3); }; +o_taper_steps : TAPER_STEPS EQUAL vec_int { driver.option_vec_int("taper_steps", $3); }; +o_geweke_interval : GEWEKE_INTERVAL EQUAL vec_value { driver.option_num("geweke_interval",$3); }; o_constant : CONSTANT { driver.option_num("noconstant", "0"); }; o_noconstant : NOCONSTANT { driver.option_num("noconstant", "1"); }; o_mh_recover : MH_RECOVER { driver.option_num("mh_recover", "1"); }; diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll index e40450354..2024f625c 100644 --- a/preprocessor/DynareFlex.ll +++ b/preprocessor/DynareFlex.ll @@ -222,6 +222,8 @@ string eofbuff; presample {return token::PRESAMPLE;} lik_algo {return token::LIK_ALGO;} lik_init {return token::LIK_INIT;} +taper_steps {return token::TAPER_STEPS;} +geweke_interval {return token::GEWEKE_INTERVAL;} graph {return token::GRAPH;} nograph {return token::NOGRAPH;} nodisplay {return token::NODISPLAY;} From e75e6a12dbb0537b519e171248ebf2c229cb8dc3 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Tue, 24 Sep 2013 16:01:35 +0200 Subject: [PATCH 05/80] estim_params: remove short-circuit ops, #476 --- preprocessor/ComputingTasks.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc index 129bfff9e..01be77d71 100644 --- a/preprocessor/ComputingTasks.cc +++ b/preprocessor/ComputingTasks.cc @@ -625,16 +625,16 @@ EstimatedParamsInitStatement::writeOutput(ostream &output, const string &basenam { if (symb_type == eExogenous) { - output << "tmp1 = find((estim_params_.corrx(:,1)==" << symb_id << " && estim_params_.corrx(:,2)==" << symbol_table.getTypeSpecificID(it->name2)+1 << ") || " - << "(estim_params_.corrx(:,2)==" << symb_id << " && estim_params_.corrx(:,1)==" << symbol_table.getTypeSpecificID(it->name2)+1 << "));" << endl; + output << "tmp1 = find((estim_params_.corrx(:,1)==" << symb_id << " & estim_params_.corrx(:,2)==" << symbol_table.getTypeSpecificID(it->name2)+1 << ") | " + << "(estim_params_.corrx(:,2)==" << symb_id << " & estim_params_.corrx(:,1)==" << symbol_table.getTypeSpecificID(it->name2)+1 << "));" << endl; output << "estim_params_.corrx(tmp1,3) = "; it->init_val->writeOutput(output); output << ";" << endl; } else if (symb_type == eEndogenous) { - output << "tmp1 = find((estim_params_.corrn(:,1)==" << symb_id << " && estim_params_.corrn(:,2)==" << symbol_table.getTypeSpecificID(it->name2)+1 << ") || " - << "(estim_params_.corrn(:,2)==" << symb_id << " && estim_params_.corrn(:,1)==" << symbol_table.getTypeSpecificID(it->name2)+1 << "));" << endl; + output << "tmp1 = find((estim_params_.corrn(:,1)==" << symb_id << " & estim_params_.corrn(:,2)==" << symbol_table.getTypeSpecificID(it->name2)+1 << ") | " + << "(estim_params_.corrn(:,2)==" << symb_id << " & estim_params_.corrn(:,1)==" << symbol_table.getTypeSpecificID(it->name2)+1 << "));" << endl; output << "estim_params_.corrn(tmp1,3) = "; it->init_val->writeOutput(output); output << ";" << endl; From 2ae0812012239e0cb0ecfbb3f2b96e1388fe0280 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Tue, 24 Sep 2013 16:25:43 +0200 Subject: [PATCH 06/80] Added the possibility to overwrite the time and init member of a dynSeries object (not allowed for freq, nobs and vobs). (cherry picked from commit b0d6e2b7b62809343dc685485f75b0f71070c245) --- matlab/@dynSeries/subsasgn.m | 21 +++++++++++++++++++-- 1 file changed, 19 insertions(+), 2 deletions(-) diff --git a/matlab/@dynSeries/subsasgn.m b/matlab/@dynSeries/subsasgn.m index da176ea21..796892a9f 100644 --- a/matlab/@dynSeries/subsasgn.m +++ b/matlab/@dynSeries/subsasgn.m @@ -97,8 +97,25 @@ switch length(S) end end end - case '.' % Single variable selection. - if ~isequal(S(1).subs,B.name) + case '.' + if isequal(S(1).subs,'init') && isa(B,'dynDate') + % Overwrite the init member... + A.init = B; + % ... and update freq and time members. + A.freq = A.init.freq; + A.time = A.init:(A.init+(A.nobs-1)); + return + elseif isequal(S(1).subs,'time') && isa(B,'dynDates') + % Overwrite the time member... + A.time = B; + % ... and update the freq and init members. + A.init = B(1); + A.freq = A.init.freq; + return + elseif ismember(S(1).subs,{'freq','nobs','vobs'}) + error(['dynSeries::subsasgn: You cannot overwrite ' S(1).subs ' member!']) + elseif ~isequal(S(1).subs,B.name) + % Single variable selection. if ~isequal(S(1).subs,B.name{1}) % Rename a variable. id = strmatch(S(1).subs,A.name); From 3dee1674d540e26ac0a2bfb326aa5603cc327fec Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Tue, 24 Sep 2013 16:25:57 +0200 Subject: [PATCH 07/80] Added unitary tests. (cherry picked from commit 667e3085851273339c55bffdeda84ae6b9041209) --- matlab/@dynSeries/subsasgn.m | 66 +++++++++++++++++++++++++++++++++++- 1 file changed, 65 insertions(+), 1 deletion(-) diff --git a/matlab/@dynSeries/subsasgn.m b/matlab/@dynSeries/subsasgn.m index 796892a9f..921ef81ba 100644 --- a/matlab/@dynSeries/subsasgn.m +++ b/matlab/@dynSeries/subsasgn.m @@ -708,4 +708,68 @@ end %$ t(7) = dyn_assert(ts1.data,[[A(1:2,1); ones(5,1); A(8:end,1)], [A(1:2,2); ones(5,1); A(8:end,2)], A(:,3)],1e-15); %$ end %$ T = all(t); -%@eof:18 \ No newline at end of file +%@eof:18 + +%@test:19 +%$ % Define a datasets. +%$ A = rand(40,3); +%$ +%$ % Instantiate two dynSeries object. +%$ ts1 = dynSeries(A,'1950Q1',{'A1';'A2';'A3'},[]); +%$ +%$ % Instantiate a dynDate object. +%$ dd = dynDate('1952Q1'); +%$ +%$ % modify first object. +%$ try +%$ ts1.init = dd; +%$ t(1) = 1; +%$ catch +%$ t(1) = 0; +%$ end +%$ +%$ % Instantiate a time series object. +%$ if t(1) +%$ t(2) = dyn_assert(ts1.vobs,3); +%$ t(3) = dyn_assert(ts1.nobs,40); +%$ t(4) = dyn_assert(ts1.name{2},'A2'); +%$ t(5) = dyn_assert(ts1.name{1},'A1'); +%$ t(6) = dyn_assert(ts1.name{3},'A3'); +%$ t(7) = dyn_assert(ts1.data,A,1e-15); +%$ t(8) = dyn_assert(isequal(ts1.init,dd),1); +%$ t(9) = dyn_assert(isequal(ts1.time(1),dd),1); +%$ end +%$ T = all(t); +%@eof:19 + +%@test:20 +%$ % Define a datasets. +%$ A = rand(40,3); +%$ +%$ % Instantiate two dynSeries object. +%$ ts1 = dynSeries(A,'1950Q1',{'A1';'A2';'A3'},[]); +%$ +%$ % Instantiate a dynDate object. +%$ dd = dynDate('1952Q1'); +%$ +%$ % modify first object. +%$ try +%$ ts1.time = dd:(dd+(ts1.nobs-1)); +%$ t(1) = 1; +%$ catch +%$ t(1) = 0; +%$ end +%$ +%$ % Instantiate a time series object. +%$ if t(1) +%$ t(2) = dyn_assert(ts1.vobs,3); +%$ t(3) = dyn_assert(ts1.nobs,40); +%$ t(4) = dyn_assert(ts1.name{2},'A2'); +%$ t(5) = dyn_assert(ts1.name{1},'A1'); +%$ t(6) = dyn_assert(ts1.name{3},'A3'); +%$ t(7) = dyn_assert(ts1.data,A,1e-15); +%$ t(8) = dyn_assert(isequal(ts1.init,dd),1); +%$ t(9) = dyn_assert(isequal(ts1.time(1),dd),1); +%$ end +%$ T = all(t); +%@eof:20 \ No newline at end of file From db485235b54ff4638891f72ed8ac12f1acd726d1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Tue, 24 Sep 2013 16:34:06 +0200 Subject: [PATCH 08/80] data, name and tex of the dynSeries class are also private members that cannot be overwritten. (cherry picked from commit 609f94881f33e857811cbfe86c9419a71d806b87) --- matlab/@dynSeries/subsasgn.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/@dynSeries/subsasgn.m b/matlab/@dynSeries/subsasgn.m index 921ef81ba..a31f75391 100644 --- a/matlab/@dynSeries/subsasgn.m +++ b/matlab/@dynSeries/subsasgn.m @@ -112,7 +112,7 @@ switch length(S) A.init = B(1); A.freq = A.init.freq; return - elseif ismember(S(1).subs,{'freq','nobs','vobs'}) + elseif ismember(S(1).subs,{'freq','nobs','vobs','data','name','tex'}) error(['dynSeries::subsasgn: You cannot overwrite ' S(1).subs ' member!']) elseif ~isequal(S(1).subs,B.name) % Single variable selection. From 2fbd75d11c852c2a31a80eb7fd70a9acff5427dd Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Wed, 25 Sep 2013 14:11:58 +0200 Subject: [PATCH 09/80] reporting: add series option tableRowColor --- doc/dynare.texi | 4 ++++ matlab/reports/@report/write.m | 3 +++ matlab/reports/@series/series.m | 2 ++ matlab/reports/@series/write.m | 3 +++ 4 files changed, 12 insertions(+) diff --git a/doc/dynare.texi b/doc/dynare.texi index 1ab025113..d7c888d32 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -8446,6 +8446,10 @@ The zero tolerance. Anything smaller than @code{zerotol} and larger than @code{-zerotol} will be set to zero before being graphed. Default: @math{1e-6} +@item tableRowColor, @code{STRING} +The color that you want the row to be. Predefined values include +@code{LightCyan} and @code{Gray}. Default: white. + @end table @end defmethod diff --git a/matlab/reports/@report/write.m b/matlab/reports/@report/write.m index 03b0dfc9a..0b34bdadb 100644 --- a/matlab/reports/@report/write.m +++ b/matlab/reports/@report/write.m @@ -54,6 +54,9 @@ else fprintf(fid, '\\usepackage{pgfplots}\n'); end +fprintf(fid, '\\usepackage{color, colortbl}\n'); +fprintf(fid, '\\definecolor{LightCyan}{rgb}{0.88,1,1}\n'); +fprintf(fid, '\\definecolor{Gray}{gray}{0.9}\n'); if o.showDate fprintf(fid, '\\usepackage{fancyhdr, datetime}\n'); fprintf(fid, '\\newdateformat{reportdate}{\\THEDAY\\ \\shortmonthname\\ \\THEYEAR}\n'); diff --git a/matlab/reports/@series/series.m b/matlab/reports/@series/series.m index 6493dced7..754d82b0d 100644 --- a/matlab/reports/@series/series.m +++ b/matlab/reports/@series/series.m @@ -51,6 +51,8 @@ o.tableMarkerLimit = 1e-4; o.tableAlignRight = false; +o.tableRowColor = ''; + o.zerotol = 1e-6; if nargin == 1 diff --git a/matlab/reports/@series/write.m b/matlab/reports/@series/write.m index e0f0ffcae..f626b47ae 100644 --- a/matlab/reports/@series/write.m +++ b/matlab/reports/@series/write.m @@ -51,6 +51,9 @@ dataString = ['%.' num2str(precision) 'f']; precision = 10^precision; fprintf(fid, '%% Table Row (series)\n'); +if ~isempty(o.tableRowColor) + fprintf(fid, '\\rowcolor{%s}', o.tableRowColor); +end if o.tableAlignRight fprintf(fid, '\\multicolumn{1}{r}{'); end From 80768beb1f110a5b3aca5a1e0f41c53524e8ea00 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Wed, 25 Sep 2013 14:44:49 +0200 Subject: [PATCH 10/80] Fix a bug similar to #476 in estimated_params_bounds --- preprocessor/ComputingTasks.cc | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc index 01be77d71..458a1db26 100644 --- a/preprocessor/ComputingTasks.cc +++ b/preprocessor/ComputingTasks.cc @@ -703,7 +703,8 @@ EstimatedParamsBoundsStatement::writeOutput(ostream &output, const string &basen { if (symb_type == eExogenous) { - output << "tmp1 = find((estim_params_.corrx(:,1)==" << symb_id << ")) & (estim_params_.corrx(:,2)==" << symbol_table.getTypeSpecificID(it->name2)+1 << ");" << endl; + output << "tmp1 = find((estim_params_.corrx(:,1)==" << symb_id << " & estim_params_.corrx(:,2)==" << symbol_table.getTypeSpecificID(it->name2)+1 << ") | " + << "(estim_params_.corrx(:,2)==" << symb_id << " & estim_params_.corrx(:,1)==" << symbol_table.getTypeSpecificID(it->name2)+1 << "));" << endl; output << "estim_params_.corrx(tmp1,4) = "; it->low_bound->writeOutput(output); @@ -715,7 +716,8 @@ EstimatedParamsBoundsStatement::writeOutput(ostream &output, const string &basen } else if (symb_type == eEndogenous) { - output << "tmp1 = find((estim_params_.corrn(:,1)==" << symb_id << ")) & (estim_params_.corrn(:,2)==" << symbol_table.getTypeSpecificID(it->name2)+1 << ";" << endl; + output << "tmp1 = find((estim_params_.corrn(:,1)==" << symb_id << " & estim_params_.corrn(:,2)==" << symbol_table.getTypeSpecificID(it->name2)+1 << ") | " + << "(estim_params_.corrn(:,2)==" << symb_id << " & estim_params_.corrn(:,1)==" << symbol_table.getTypeSpecificID(it->name2)+1 << "));" << endl; output << "estim_params_.corrn(tmp1,4) = "; it->low_bound->writeOutput(output); From 55bcee2f9d1e7167d22c4bfd005a8dfb636d907f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Wed, 25 Sep 2013 15:27:34 +0200 Subject: [PATCH 11/80] Fix in license file related to changes to Sims' codes --- license.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/license.txt b/license.txt index 2dcfdcf07..95867ea33 100644 --- a/license.txt +++ b/license.txt @@ -40,12 +40,12 @@ License: public-domain pages 472-489 Files: matlab/bfgsi1.m matlab/csolve.m matlab/csminit1.m matlab/numgrad2.m - matlab/numgrad2_.m matlab/numgrad3.m matlab/numgrad3_.m matlab/numgrad5.m + matlab/numgrad3.m matlab/numgrad3_.m matlab/numgrad5.m matlab/numgrad5_.m matlab/csminwel1.m matlab/bvar_density.m matlab/bvar_toolbox.m matlab/partial_information/PI_gensys.m matlab/qzswitch.m matlab/qzdiv.m Copyright: 1993-2009 Christopher Sims - 2006-2012 Dynare Team + 2006-2013 Dynare Team License: GPL-3+ Files: matlab/cmaes.m From 4522c011d6b61f23546dfc874e3d66eb00b12b7e Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Wed, 25 Sep 2013 14:32:15 +0200 Subject: [PATCH 12/80] reporting: default value bug fix --- matlab/reports/@series/series.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/reports/@series/series.m b/matlab/reports/@series/series.m index 754d82b0d..30f8e462f 100644 --- a/matlab/reports/@series/series.m +++ b/matlab/reports/@series/series.m @@ -51,7 +51,7 @@ o.tableMarkerLimit = 1e-4; o.tableAlignRight = false; -o.tableRowColor = ''; +o.tableRowColor = 'white'; o.zerotol = 1e-6; From 6a9e9dfdea9248e67ea576228ca7c82f3e8334a7 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Wed, 25 Sep 2013 15:25:34 +0200 Subject: [PATCH 13/80] reporting: new table option vlineAfterEndOfPeriod --- doc/dynare.texi | 6 +++++- matlab/reports/@table/table.m | 3 +++ matlab/reports/@table/write.m | 20 ++++++++++++++------ 3 files changed, 22 insertions(+), 7 deletions(-) diff --git a/doc/dynare.texi b/doc/dynare.texi index c99e583b3..280695f76 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -8421,7 +8421,7 @@ Display a solid black line at @math{y = 0}. Default: @code{false} @end table @end defmethod -@defmethod Report addTable data, showHlines, precision, range, seriesToUse, title, titleSize, vlineAfter, showVlines +@defmethod Report addTable data, showHlines, precision, range, seriesToUse, title, titleSize, vlineAfter, vlineAfterEndOfPeriod, showVlines Adds a @code{Table} to a @code{Section}. @optionshead @table @code @@ -8450,6 +8450,10 @@ Title for the table. Default: @code{none} @item vlineAfter, @code{dynDate} Show a vertical line after the specified date. Default: @code{empty} +@item vlineAfterEndOfPeriod, @code{BOOLEAN} +Show a vertical line after the end of every period (@i{i.e.} after +every year, after the fourth quarter, etc.). Default: @code{false} + @item showVlines, @code{BOOLEAN} Whether or not to show vertical lines separating the columns. Default: @code{false} @end table diff --git a/matlab/reports/@table/table.m b/matlab/reports/@table/table.m index 55cdde174..ce3808eeb 100644 --- a/matlab/reports/@table/table.m +++ b/matlab/reports/@table/table.m @@ -39,6 +39,7 @@ o.titleSize = 'large'; o.showHlines = false; o.showVlines = false; o.vlineAfter = ''; +o.vlineAfterEndOfPeriod = false; o.data = ''; o.seriesToUse = ''; @@ -87,6 +88,8 @@ assert(isempty(o.vlineAfter) || isa(o.vlineAfter, 'dynDate'), ... if o.showVlines o.vlineAfter = ''; end +assert(islogical(o.vlineAfterEndOfPeriod), ... + '@table.table: vlineAfterEndOfPeriod must be true or false'); valid_title_sizes = {'Huge', 'huge', 'LARGE', 'Large', 'large', 'normalsize', ... 'small', 'footnotesize', 'scriptsize', 'tiny'}; assert(any(strcmp(o.titleSize, valid_title_sizes)), ... diff --git a/matlab/reports/@table/write.m b/matlab/reports/@table/write.m index afb47f16b..82ba6233f 100644 --- a/matlab/reports/@table/write.m +++ b/matlab/reports/@table/write.m @@ -55,12 +55,20 @@ fprintf(fid, '\\begin{tabular}{@{}l'); for i=1:ndates if o.showVlines - fprintf(fid, '|'); - end - fprintf(fid, 'r'); - if ~isempty(o.vlineAfter) - if dates(i) == o.vlineAfter - fprintf(fid, '|'); + fprintf(fid, 'r|'); + else + fprintf(fid, 'r'); + if o.vlineAfterEndOfPeriod + if dates(i).time(2) == dates(i).freq + fprintf(fid, '|'); + end + end + if ~isempty(o.vlineAfter) + if dates(i) == o.vlineAfter + if ~(o.vlineAfterEndOfPeriod && dates(i).time(2) == dates(i).freq) + fprintf(fid, '|'); + end + end end end end From 79b33ca741df8f45804ad3c7af84e6536c9a19be Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Wed, 25 Sep 2013 15:26:49 +0200 Subject: [PATCH 14/80] reporting: doc bug fix --- doc/dynare.texi | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/doc/dynare.texi b/doc/dynare.texi index 280695f76..6f32eb051 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -8460,7 +8460,7 @@ Whether or not to show vertical lines separating the columns. Default: @code{fal @end defmethod @anchor{addSeries} -@defmethod Report addSeries data, graphLineColor, graphLineStyle, graphLineWidth, graphMarker, graphMarkerEdgeColor, graphMarkerFaceColor, graphMarkerSize, tableShowMarkers, tableAlignRight, tableNegColor, tablePosColor, zerotol +@defmethod Report addSeries data, graphLineColor, graphLineStyle, graphLineWidth, graphMarker, graphMarkerEdgeColor, graphMarkerFaceColor, graphMarkerSize, tableRowColor, tableShowMarkers, tableAlignRight, tableNegColor, tablePosColor, zerotol Adds a @code{Series} to a @code{Graph} or a @code{Table}. @optionshead @table @code @@ -8489,6 +8489,10 @@ The face color of the graph marker. Default: @code{`auto'} @item graphMarkerSize, @code{DOUBLE} The size of the graph marker. Default: @code{6} +@item tableRowColor, @code{STRING} +The color that you want the row to be. Predefined values include +@code{LightCyan} and @code{Gray}. Default: @code{white}. + @item tableShowMarkers, @code{BOOLEAN} In a Table, if @code{true}, surround each cell with brackets and color it according to @ref{tableNegColor} and @ref{tablePosColor}. No effect @@ -8519,10 +8523,6 @@ The zero tolerance. Anything smaller than @code{zerotol} and larger than @code{-zerotol} will be set to zero before being graphed. Default: @math{1e-6} -@item tableRowColor, @code{STRING} -The color that you want the row to be. Predefined values include -@code{LightCyan} and @code{Gray}. Default: white. - @end table @end defmethod From b3d835bd48c3cb2917aa8c04e2dd39e8ed50fa40 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Wed, 25 Sep 2013 15:59:50 +0200 Subject: [PATCH 15/80] Provisions for MATLAB 8.2 (R2013b) --- m4/ax_matlab_version.m4 | 3 +++ matlab/dynare_config.m | 4 ++-- windows/dynare.nsi | 12 ++++++------ 3 files changed, 11 insertions(+), 8 deletions(-) diff --git a/m4/ax_matlab_version.m4 b/m4/ax_matlab_version.m4 index 12b5693dc..55eb57036 100644 --- a/m4/ax_matlab_version.m4 +++ b/m4/ax_matlab_version.m4 @@ -22,6 +22,9 @@ AC_REQUIRE([AX_MATLAB]) AC_MSG_CHECKING([for MATLAB version]) if test "x$MATLAB_VERSION" != "x"; then case $MATLAB_VERSION in + *2013b | *2013B) + MATLAB_VERSION="8.2" + ;; *2013a | *2013A) MATLAB_VERSION="8.1" ;; diff --git a/matlab/dynare_config.m b/matlab/dynare_config.m index 2cd7b7487..3180a81a6 100644 --- a/matlab/dynare_config.m +++ b/matlab/dynare_config.m @@ -118,7 +118,7 @@ else addpath(mexpath) end else - mexpath = [dynareroot '../mex/matlab/win32-7.5-8.1']; + mexpath = [dynareroot '../mex/matlab/win32-7.5-8.2']; if exist(mexpath, 'dir') addpath(mexpath) end @@ -138,7 +138,7 @@ else addpath(mexpath) end else - mexpath = [dynareroot '../mex/matlab/win64-7.8-8.1']; + mexpath = [dynareroot '../mex/matlab/win64-7.8-8.2']; if exist(mexpath, 'dir') addpath(mexpath) end diff --git a/windows/dynare.nsi b/windows/dynare.nsi index bd17b0bfe..e2fbc886d 100644 --- a/windows/dynare.nsi +++ b/windows/dynare.nsi @@ -88,9 +88,9 @@ Section "MEX files for MATLAB 32-bit, version 7.3 to 7.4 (R2006b to R2007a)" File ..\mex\matlab\win32-7.3-7.4\*.mexw32 SectionEnd -Section "MEX files for MATLAB 32-bit, version 7.5 to 8.1 (R2007b to R2013a)" - SetOutPath $INSTDIR\mex\matlab\win32-7.5-8.1 - File ..\mex\matlab\win32-7.5-8.1\*.mexw32 +Section "MEX files for MATLAB 32-bit, version 7.5 to 8.2 (R2007b to R2013b)" + SetOutPath $INSTDIR\mex\matlab\win32-7.5-8.2 + File ..\mex\matlab\win32-7.5-8.2\*.mexw32 SectionEnd Section "MEX files for MATLAB 64-bit, version 7.3 to 7.4 (R2006b to R2007a)" @@ -103,9 +103,9 @@ Section "MEX files for MATLAB 64-bit, version 7.5 to 7.7 (R2007b to R2008b)" File ..\mex\matlab\win64-7.5-7.7\*.mexw64 SectionEnd -Section "MEX files for MATLAB 64-bit, version 7.8 to 8.1 (R2009a to R2013a)" - SetOutPath $INSTDIR\mex\matlab\win64-7.8-8.1 - File ..\mex\matlab\win64-7.8-8.1\*.mexw64 +Section "MEX files for MATLAB 64-bit, version 7.8 to 8.2 (R2009a to R2013b)" + SetOutPath $INSTDIR\mex\matlab\win64-7.8-8.2 + File ..\mex\matlab\win64-7.8-8.2\*.mexw64 SectionEnd SectionGroupEnd From 90fa92f8a95b0908aeb49ec84bddd047afa83362 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Thu, 26 Sep 2013 10:19:05 +0200 Subject: [PATCH 16/80] Adds chi2cdf function if Statistics toolbox is missing --- matlab/missing/stats/chi2cdf.m | 37 ++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) create mode 100644 matlab/missing/stats/chi2cdf.m diff --git a/matlab/missing/stats/chi2cdf.m b/matlab/missing/stats/chi2cdf.m new file mode 100644 index 000000000..0a5db9882 --- /dev/null +++ b/matlab/missing/stats/chi2cdf.m @@ -0,0 +1,37 @@ +function CDF = chi2cdf(x, n) +% chi2cdf CDF of the Chi2 distribution +% CDF = chi2cdf(x, n) computes, for each element of X, the +% CDF at X of the chi-square distribution with N degrees of freedom. +% Original file: statistics/distributions/chi2inv.m + +% Copyright (C) 2013 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +if (nargin ~= 2) + error ('chi2cdf: you must give two arguments'); +end + +if (~isscalar (n)) + [retval, x, n] = common_size (x, n); + if (retval > 0) + error ('chi2cdf: x and n must be of common size or scalar'); + end +end + +CDF = gamcdf (x, n / 2, 2); + +end From 27b48720dcf41870e466fb6fbb040f2ffa2d03ec Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Thu, 26 Sep 2013 11:29:40 +0200 Subject: [PATCH 17/80] Clarify manual on different/inconsistent ordering of variables used in description of decision rules The previous description used the same variables to denote both declaration and DR order, thus confusing users. --- doc/dynare.texi | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/doc/dynare.texi b/doc/dynare.texi index 6f32eb051..970c8d9c2 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -3735,7 +3735,7 @@ declaration order. Conversely, k-th declared variable is numbered @vindex M_.nsfwrd @vindex M_.ndynamic Finally, the state variables of the model are the purely backward variables -and the mixed variables. They are orderer in DR-order when they appear in +and the mixed variables. They are ordered in DR-order when they appear in decision rules elements. There are @code{M_.nspred = M_.npred + M_.nboth} such variables. Similarly, one has @code{M_.nsfwrd = M_.nfwrd + M_.nboth}, and @code{M_.ndynamic = M_.nfwrd+M_.nboth+M_.npred}. @@ -3743,7 +3743,7 @@ and @code{M_.ndynamic = M_.nfwrd+M_.nboth+M_.npred}. @node First order approximation @subsection First order approximation -The approximation has the form: +The approximation has the stylized form: @math{y_t = y^s + A y^h_{t-1} + B u_t} @@ -3772,6 +3772,12 @@ endogenous in DR-order. The matrix columns correspond to exogenous variables in declaration order. @end itemize +Of course, the shown form of the approximation is only stylized, because it neglects the required different ordering in @math{y^s} and @math{y^h_t}. The precise form of the approximation that shows the way Dynare deals with differences between declaration and DR-order, is + +@math{y_t(oo_.dr.order_var) = y^s(oo_.dr.order_var) + A (y_{t-1}(oo_.dr.order_var(k2))-y^s(oo_.dr.order_var(k2))) + B u_t} + +where @math{k2} selects the state variables, @math{y_t} and @math{y^s} are in declaration order and the coefficient matrices are in DR-order. Effectively, all variables on the right hand side are brought into DR order for computations and then assigned to @math{y_t} in declaration order. + @node Second order approximation @subsection Second order approximation @@ -3785,7 +3791,7 @@ A y^h_{t-1} + B u_t + 0.5 C where @math{y^s} is the steady state value of @math{y}, @math{y^h_t=y_t-y^s}, and @math{\Delta^2} is the shift effect of the -variance of future shocks. +variance of future shocks. For the reordering required due to differences in declaration and DR order, see the first order approximation. The coefficients of the decision rules are stored in the variables described for first order approximation, plus the following variables: From 79f03773d77c9d6fd8d5ef0bfc6b01c589a33a44 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Thu, 26 Sep 2013 14:37:33 +0200 Subject: [PATCH 18/80] reporting: add new series option tableSubSectionHeader --- doc/dynare.texi | 7 ++++++- matlab/reports/@series/series.m | 1 + matlab/reports/@series/write.m | 15 +++++++++++++-- 3 files changed, 20 insertions(+), 3 deletions(-) diff --git a/doc/dynare.texi b/doc/dynare.texi index 970c8d9c2..b62056b63 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -8466,7 +8466,7 @@ Whether or not to show vertical lines separating the columns. Default: @code{fal @end defmethod @anchor{addSeries} -@defmethod Report addSeries data, graphLineColor, graphLineStyle, graphLineWidth, graphMarker, graphMarkerEdgeColor, graphMarkerFaceColor, graphMarkerSize, tableRowColor, tableShowMarkers, tableAlignRight, tableNegColor, tablePosColor, zerotol +@defmethod Report addSeries data, graphLineColor, graphLineStyle, graphLineWidth, graphMarker, graphMarkerEdgeColor, graphMarkerFaceColor, graphMarkerSize, tableRowColor, tableShowMarkers, tableAlignRight, tableNegColor, tablePosColor, tableSubSectionHeader, zerotol Adds a @code{Series} to a @code{Graph} or a @code{Table}. @optionshead @table @code @@ -8524,6 +8524,11 @@ zero. Default: @code{`red'} The color to use when marking Table data that is greater than zero. Default: @code{`blue'} +@item tableSubSectionHeader, @code{STRING} +A header for a subsection of the table. No data will be associated +with it. It is equivalent to adding an empty series with a +name. Default: @code{''} + @item zerotol, @code{DOUBLE} The zero tolerance. Anything smaller than @code{zerotol} and larger than @code{-zerotol} will be set to zero before being diff --git a/matlab/reports/@series/series.m b/matlab/reports/@series/series.m index 30f8e462f..ce30bf417 100644 --- a/matlab/reports/@series/series.m +++ b/matlab/reports/@series/series.m @@ -49,6 +49,7 @@ o.tableNegColor = 'red'; o.tablePosColor = 'blue'; o.tableMarkerLimit = 1e-4; +o.tableSubSectionHeader = ''; o.tableAlignRight = false; o.tableRowColor = 'white'; diff --git a/matlab/reports/@series/write.m b/matlab/reports/@series/write.m index f626b47ae..586c7125d 100644 --- a/matlab/reports/@series/write.m +++ b/matlab/reports/@series/write.m @@ -37,8 +37,11 @@ assert(isa(dates, 'dynDates')); assert(isint(precision)); %% Validate options provided by user -assert(~isempty(o.data) && isa(o.data, 'dynSeries'), ... - '@series.write: must provide data as a dynSeries'); +assert(ischar(o.tableSubSectionHeader), '@series.write: tableSubSectionHeader must be a string'); +if isempty(o.tableSubSectionHeader) + assert(~isempty(o.data) && isa(o.data, 'dynSeries'), ... + '@series.write: must provide data as a dynSeries'); +end assert(ischar(o.tableNegColor), '@series.write: tableNegColor must be a string'); assert(ischar(o.tablePosColor), '@series.write: tablePosColor must be a string'); @@ -54,6 +57,14 @@ fprintf(fid, '%% Table Row (series)\n'); if ~isempty(o.tableRowColor) fprintf(fid, '\\rowcolor{%s}', o.tableRowColor); end +if ~isempty(o.tableSubSectionHeader) + fprintf(fid, '%s', o.tableSubSectionHeader); + for i=1:size(dates) + fprintf(fid, ' & '); + end + fprintf(fid, '\\\\%%\n'); + return; +end if o.tableAlignRight fprintf(fid, '\\multicolumn{1}{r}{'); end From a9eb95ff2aebc69e8d9cca38a19cf6e4f386ff00 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Thu, 26 Sep 2013 14:45:46 +0200 Subject: [PATCH 19/80] reporting: series: fix spacing of table cell separators and eol's --- matlab/reports/@series/write.m | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/matlab/reports/@series/write.m b/matlab/reports/@series/write.m index 586c7125d..2511b9944 100644 --- a/matlab/reports/@series/write.m +++ b/matlab/reports/@series/write.m @@ -60,7 +60,7 @@ end if ~isempty(o.tableSubSectionHeader) fprintf(fid, '%s', o.tableSubSectionHeader); for i=1:size(dates) - fprintf(fid, ' & '); + fprintf(fid, '&'); end fprintf(fid, '\\\\%%\n'); return; @@ -75,7 +75,7 @@ end data = o.data(dates); data = data.data; for i=1:size(data,1) - fprintf(fid, ' &'); + fprintf(fid, '&'); if o.tableShowMarkers if data(i) < -o.tableMarkerLimit fprintf(fid, '\\color{%s}', o.tableNegColor); @@ -91,5 +91,5 @@ for i=1:size(data,1) fprintf(fid, ']'); end end -fprintf(fid, ' \\\\\n\n'); +fprintf(fid, '\\\\%%\n'); end From 34b3b94c5f98e44ca6e3e8f5c1d48cf9791daae7 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Thu, 26 Sep 2013 14:48:59 +0200 Subject: [PATCH 20/80] reporting: series: tableRowColor add missing check --- matlab/reports/@series/write.m | 1 + 1 file changed, 1 insertion(+) diff --git a/matlab/reports/@series/write.m b/matlab/reports/@series/write.m index 2511b9944..29c47ddcb 100644 --- a/matlab/reports/@series/write.m +++ b/matlab/reports/@series/write.m @@ -45,6 +45,7 @@ end assert(ischar(o.tableNegColor), '@series.write: tableNegColor must be a string'); assert(ischar(o.tablePosColor), '@series.write: tablePosColor must be a string'); +assert(ischar(o.tableRowColor), '@series.write: tableRowColor must be a string'); assert(islogical(o.tableShowMarkers), '@series.write: tableShowMarkers must be true or false'); assert(islogical(o.tableAlignRight), '@series.write: tableAlignRight must be true or false'); assert(isfloat(o.tableMarkerLimit), '@series,write: tableMarkerLimit must be a float'); From ff3774a4a2d167093239773417e474bc98145860 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Thu, 26 Sep 2013 15:21:08 +0200 Subject: [PATCH 21/80] reporting: table: remove unused variable --- matlab/reports/@table/write.m | 1 - 1 file changed, 1 deletion(-) diff --git a/matlab/reports/@table/write.m b/matlab/reports/@table/write.m index 82ba6233f..ed4239099 100644 --- a/matlab/reports/@table/write.m +++ b/matlab/reports/@table/write.m @@ -83,7 +83,6 @@ fprintf(fid, '\\toprule%%\n'); datedata = dates.time; years = unique(datedata(:, 1)); thdr = num2cell(years, size(years, 1)); -lind = nlhc; switch dates.freq case 1 for i=1:size(thdr, 1) From 79d62d85d0e926ccc2d78f42c612e915b86844d8 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Thu, 26 Sep 2013 16:08:13 +0200 Subject: [PATCH 22/80] reporting: fix hline for table --- matlab/reports/@table/write.m | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/matlab/reports/@table/write.m b/matlab/reports/@table/write.m index ed4239099..e00569d9f 100644 --- a/matlab/reports/@table/write.m +++ b/matlab/reports/@table/write.m @@ -121,11 +121,8 @@ switch dates.freq otherwise error('@table.write: invalid dynSeries frequency'); end -fprintf(fid, '\\\\[-10pt]%%\n'); -for i=1:ndates - fprintf(fid, ' & \\hrulefill'); -end -fprintf(fid, '\\\\%%\n'); +fprintf(fid, '\\\\[-2pt]%%\n'); +fprintf(fid, '\\hline%%\n'); fprintf(fid, '%%\n'); % Write Table Data From d1c012d289c7958ebe262b9337e5bd2a9fd3f2dc Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Thu, 26 Sep 2013 16:34:36 +0200 Subject: [PATCH 23/80] reporting: annualAverages option for addTable --- doc/dynare.texi | 7 ++++++- matlab/reports/@series/write.m | 27 +++++++++++++++++++++++++-- matlab/reports/@table/table.m | 3 +++ matlab/reports/@table/write.m | 21 +++++++++++++++++---- 4 files changed, 51 insertions(+), 7 deletions(-) diff --git a/doc/dynare.texi b/doc/dynare.texi index b62056b63..6902abec9 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -8427,11 +8427,16 @@ Display a solid black line at @math{y = 0}. Default: @code{false} @end table @end defmethod -@defmethod Report addTable data, showHlines, precision, range, seriesToUse, title, titleSize, vlineAfter, vlineAfterEndOfPeriod, showVlines +@defmethod Report addTable annualAverages, data, showHlines, precision, range, seriesToUse, title, titleSize, vlineAfter, vlineAfterEndOfPeriod, showVlines Adds a @code{Table} to a @code{Section}. @optionshead @table @code +@item annualAverages, @code{bool} +Compute the average over every year in the table and display it in a +column to the right of the data (one column for every year). Only +works for quarterly data. Default: @code{false} + @item data, @code{dynSeries} @xref{data}. diff --git a/matlab/reports/@series/write.m b/matlab/reports/@series/write.m index 29c47ddcb..ddca0654e 100644 --- a/matlab/reports/@series/write.m +++ b/matlab/reports/@series/write.m @@ -1,5 +1,5 @@ -function o = write(o, fid, dates, precision) -%function o = write(o, fid, dates, precision) +function o = write(o, fid, dates, precision, yrsForAvgs) +%function o = write(o, fid, dates, precision, yrsForAvgs) % Write Table Row % % INPUTS @@ -7,6 +7,8 @@ function o = write(o, fid, dates, precision) % fid [int] file id % dates [dynDates] dates for series slice % precision [float] precision with which to print the data +% yrsForAvgs [bool] the years for which to compute averages +% % % OUTPUTS % o [series] series object @@ -92,5 +94,26 @@ for i=1:size(data,1) fprintf(fid, ']'); end end + +% Calculate and display annual averages +for i=1:length(yrsForAvgs) + slice = o.data(dynDate([num2str(yrsForAvgs(i)) 'q1']):dynDate([num2str(yrsForAvgs(i)) 'q4'])); + avg = sum(slice.data)/size(slice); + fprintf(fid, '&'); + if o.tableShowMarkers + if avg < -o.tableMarkerLimit + fprintf(fid, '\\color{%s}', o.tableNegColor); + elseif avg > o.tableMarkerLimit + fprintf(fid, '\\color{%s}', o.tablePosColor); + end + fprintf(fid, '['); + end + + fprintf(fid, dataString, round(avg*precision)/precision); + + if o.tableShowMarkers + fprintf(fid, ']'); + end +end fprintf(fid, '\\\\%%\n'); end diff --git a/matlab/reports/@table/table.m b/matlab/reports/@table/table.m index ce3808eeb..eb742409c 100644 --- a/matlab/reports/@table/table.m +++ b/matlab/reports/@table/table.m @@ -41,6 +41,8 @@ o.showVlines = false; o.vlineAfter = ''; o.vlineAfterEndOfPeriod = false; +o.annualAverages = false; + o.data = ''; o.seriesToUse = ''; o.range = {}; @@ -73,6 +75,7 @@ end % Check options provided by user assert(ischar(o.title), '@table.table: title must be a string'); +assert(islogical(o.annualAverages), '@table.table: annualAverages must be true or false'); assert(islogical(o.showHlines), '@table.table: showHlines must be true or false'); assert(islogical(o.showVlines), '@table.table: showVlines must be true or false'); assert(isint(o.precision), '@table.table: precision must be an int'); diff --git a/matlab/reports/@table/write.m b/matlab/reports/@table/write.m index e00569d9f..9e881778d 100644 --- a/matlab/reports/@table/write.m +++ b/matlab/reports/@table/write.m @@ -72,6 +72,19 @@ for i=1:ndates end end end +datedata = dates.time; +years = unique(datedata(:, 1)); +if o.annualAverages + yrsForAvgs = years; +else + yrsForAvgs = []; +end +for i=1:length(yrsForAvgs) + fprintf(fid, 'r'); + if o.showVlines + fprintf(fid, '|'); + end +end fprintf(fid, '@{}}%%\n'); if ~isempty(o.title) fprintf(fid, '\\multicolumn{%d}{c}{\\%s %s}\\\\\n', ... @@ -80,8 +93,6 @@ end fprintf(fid, '\\toprule%%\n'); % Column Headers -datedata = dates.time; -years = unique(datedata(:, 1)); thdr = num2cell(years, size(years, 1)); switch dates.freq case 1 @@ -101,7 +112,6 @@ switch dates.freq end end end - for i=1:size(thdr, 1) fprintf(fid, ' & \\multicolumn{%d}{c}{%d}', size(thdr{i,2}, 2), thdr{i,1}); end @@ -121,6 +131,9 @@ switch dates.freq otherwise error('@table.write: invalid dynSeries frequency'); end +for i=1:length(yrsForAvgs) + fprintf(fid, ' & %d', years(i)); +end fprintf(fid, '\\\\[-2pt]%%\n'); fprintf(fid, '\\hline%%\n'); fprintf(fid, '%%\n'); @@ -128,7 +141,7 @@ fprintf(fid, '%%\n'); % Write Table Data ne = o.seriesElements.numSeriesElements(); for i=1:ne - o.seriesElements(i).write(fid, dates, o.precision); + o.seriesElements(i).write(fid, dates, o.precision, yrsForAvgs); if o.showHlines fprintf(fid, '\\hline\n'); end From 48c428d748e3029e5e817b7515a77a5ec5102fe5 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Fri, 27 Sep 2013 11:14:00 +0200 Subject: [PATCH 24/80] reporting: make tableSubSectionHeader prettier when vlineAfterEndOfPeriod is true --- matlab/reports/@series/write.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/reports/@series/write.m b/matlab/reports/@series/write.m index ddca0654e..eb2b6f277 100644 --- a/matlab/reports/@series/write.m +++ b/matlab/reports/@series/write.m @@ -62,7 +62,7 @@ if ~isempty(o.tableRowColor) end if ~isempty(o.tableSubSectionHeader) fprintf(fid, '%s', o.tableSubSectionHeader); - for i=1:size(dates) + for i=1:size(dates)+length(yrsForAvgs) fprintf(fid, '&'); end fprintf(fid, '\\\\%%\n'); From 1f8b4d9a861a8f3288c4e9b2c5f71a71242f2d74 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Mon, 30 Sep 2013 17:02:25 +0200 Subject: [PATCH 25/80] Fix bug in display of parameter names violating the bounds condition --- matlab/dynare_estimation_init.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/dynare_estimation_init.m b/matlab/dynare_estimation_init.m index 793b96c8e..148140bac 100644 --- a/matlab/dynare_estimation_init.m +++ b/matlab/dynare_estimation_init.m @@ -245,7 +245,7 @@ if ~isempty(estim_params_) outside_bound_pars=find(xparam1 < bounds(:,1) | xparam1 > bounds(:,2)); if ~isempty(outside_bound_pars) for ii=1:length(outside_bound_pars) - outside_bound_par_names{ii,1}=get_the_name(ii,0,M_,estim_params_,options_); + outside_bound_par_names{ii,1}=get_the_name(outside_bound_pars(ii),0,M_,estim_params_,options_); end disp_string=[outside_bound_par_names{1,:}]; for ii=2:size(outside_bound_par_names,1) From 7e7cadb8786b74deb072079ca54cb3c64d4f6b95 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Sun, 29 Sep 2013 16:13:14 +0200 Subject: [PATCH 26/80] Check initial values for violation of inverse gamma prior The inverse gamma distribution does not allow for the value 0, but the current check at the lower bound set LB=0 and tested for Date: Mon, 30 Sep 2013 17:10:39 +0200 Subject: [PATCH 27/80] Fixes bug in display of parameters at prior bound --- matlab/dynare_estimation_1.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index 73e24df17..c9d489d76 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -502,7 +502,7 @@ if ~options_.mh_posterior_mode_estimation && options_.cova_compute params_at_bound=find(xparam1==ub | xparam1==lb); if ~isempty(params_at_bound) for ii=1:length(params_at_bound) - params_at_bound_name{ii,1}=get_the_name(ii,0,M_,estim_params_,options_); + params_at_bound_name{ii,1}=get_the_name(params_at_bound(ii),0,M_,estim_params_,options_); end disp_string=[params_at_bound_name{1,:}]; for ii=2:size(params_at_bound_name,1) From 45942b244ccfac517d84ad72d3724baf6fc77479 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Tue, 1 Oct 2013 15:13:18 +0200 Subject: [PATCH 28/80] stop dynare.m if onlymacro is passed, closes #487 --- matlab/dynare.m | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/matlab/dynare.m b/matlab/dynare.m index 06d81f699..baa87a788 100644 --- a/matlab/dynare.m +++ b/matlab/dynare.m @@ -108,6 +108,10 @@ end [status, result] = system(command); disp(result) +if ismember('onlymacro', varargin) + disp('Preprocesser stopped after macroprocessing step because of ''onlymacro'' option.'); + return; +end % post-dynare-prerocessor-hook if exist([fname(1:end-4) '_post_dynare_preprocessor_hook.m'],'file') From 4b060c47d2362e0bb0b2aedea26feb7ef7aaff56 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Wed, 2 Oct 2013 09:48:07 +0200 Subject: [PATCH 29/80] Add more information when IRFs are not displayed due to explosive simulations --- matlab/stoch_simul.m | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/matlab/stoch_simul.m b/matlab/stoch_simul.m index 6bb7ada05..5c6a1e06f 100644 --- a/matlab/stoch_simul.m +++ b/matlab/stoch_simul.m @@ -176,6 +176,12 @@ if options_.irf y=irf(oo_.dr,cs(M_.exo_names_orig_ord,i), options_.irf, options_.drop, ... options_.replic, options_.order); end + if ~options_.noprint && any(any(isnan(y))) && ~options_.pruning && ~(options_.order==1) + fprintf('\nSTOCH_SIMUL: The simulations conducted for generating IRFs to %s were explosive.\n',M_.exo_names(i,:)) + fprintf('STOCH_SIMUL: No IRFs will be displayed. Either reduce the shock size, \n') + fprintf('STOCH_SIMUL: use pruning, or set the approximation order to 1.'); + skipline(2); + end if options_.relative_irf y = 100*y/cs(i,i); end From b3b2cdfe7b53f14990a2e8a021008cbd120b3ace Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Wed, 2 Oct 2013 11:12:32 +0200 Subject: [PATCH 30/80] bug fix for date_number: don't permit weeks >= 53 --- preprocessor/DynareFlex.ll | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll index 2024f625c..35616b1df 100644 --- a/preprocessor/DynareFlex.ll +++ b/preprocessor/DynareFlex.ll @@ -691,7 +691,7 @@ string eofbuff; return token::INT_NUMBER; } -([1-2][0-9]{3}[Mm](([1-9])|(1[0-2])))|([1-2][0-9]{3}[Qq][1-4])|([1-2][0-9]{3}[Ww](([1-9]{1})|([1-5][0-9]))) { +([1-2][0-9]{3}[Mm](([1-9])|(1[0-2])))|([1-2][0-9]{3}[Qq][1-4])|([1-2][0-9]{3}[Ww](([1-9]{1})|([1-4][0-9])|([5][0-2]))) { yylval->string_val = new string(yytext); return token::DATE_NUMBER; } From 9428007bcb311cac5fc65a3a71c33ab2352cdda6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Wed, 2 Oct 2013 15:29:56 +0200 Subject: [PATCH 31/80] Added simpsa routines (new optimization algorithm) --- matlab/simpsa.m | 535 +++++++++++++++++++++++++++++++++++++++++++++ matlab/simpsaget.m | 125 +++++++++++ matlab/simpsaset.m | 172 +++++++++++++++ 3 files changed, 832 insertions(+) create mode 100644 matlab/simpsa.m create mode 100644 matlab/simpsaget.m create mode 100644 matlab/simpsaset.m diff --git a/matlab/simpsa.m b/matlab/simpsa.m new file mode 100644 index 000000000..97e5c2fd3 --- /dev/null +++ b/matlab/simpsa.m @@ -0,0 +1,535 @@ +function [X,FVAL,EXITFLAG,OUTPUT] = simpsa(FUN,X0,LB,UB,OPTIONS,varargin) + +% Finds a minimum of a function of several variables using an algorithm +% that is based on the combination of the non-linear smplex and the simulated +% annealing algorithm (the SIMPSA algorithm, Cardoso et al., 1996). +% In this paper, the algorithm is shown to be adequate for the global optimi- +% zation of an example set of unconstrained and constrained NLP functions. +% +% SIMPSA attempts to solve problems of the form: +% min F(X) subject to: LB <= X <= UB +% X +% +% Algorithm partly is based on paper of Cardoso et al, 1996. +% +% X=SIMPSA(FUN,X0) start at X0 and finds a minimum X to the function FUN. +% FUN accepts input X and returns a scalar function value F evaluated at X. +% X0 may be a scalar, vector, or matrix. +% +% X=SIMPSA(FUN,X0,LB,UB) defines a set of lower and upper bounds on the +% design variables, X, so that a solution is found in the range +% LB <= X <= UB. Use empty matrices for LB and UB if no bounds exist. +% Set LB(i) = -Inf if X(i) is unbounded below; set UB(i) = Inf if X(i) is +% unbounded above. +% +% X=SIMPSA(FUN,X0,LB,UB,OPTIONS) minimizes with the default optimization +% parameters replaced by values in the structure OPTIONS, an argument +% created with the SIMPSASET function. See SIMPSASET for details. +% Used options are TEMP_START, TEMP_END, COOL_RATE, INITIAL_ACCEPTANCE_RATIO, +% MIN_COOLING_FACTOR, MAX_ITER_TEMP_FIRST, MAX_ITER_TEMP_LAST, MAX_ITER_TEMP, +% MAX_ITER_TOTAL, MAX_TIME, MAX_FUN_EVALS, TOLX, TOLFUN, DISPLAY and OUTPUT_FCN. +% Use OPTIONS = [] as a place holder if no options are set. +% +% X=SIMPSA(FUN,X0,LB,UB,OPTIONS,varargin) is used to supply a variable +% number of input arguments to the objective function FUN. +% +% [X,FVAL]=SIMPSA(FUN,X0,...) returns the value of the objective +% function FUN at the solution X. +% +% [X,FVAL,EXITFLAG]=SIMPSA(FUN,X0,...) returns an EXITFLAG that describes the +% exit condition of SIMPSA. Possible values of EXITFLAG and the corresponding +% exit conditions are: +% +% 1 Change in the objective function value less than the specified tolerance. +% 2 Change in X less than the specified tolerance. +% 0 Maximum number of function evaluations or iterations reached. +% -1 Maximum time exceeded. +% +% [X,FVAL,EXITFLAG,OUTPUT]=SIMPSA(FUN,X0,...) returns a structure OUTPUT with +% the number of iterations taken in OUTPUT.nITERATIONS, the number of function +% evaluations in OUTPUT.nFUN_EVALS, the temperature profile in OUTPUT.TEMPERATURE, +% the simplexes that were evaluated in OUTPUT.SIMPLEX and the best one in +% OUTPUT.SIMPLEX_BEST, the costs associated with each simplex in OUTPUT.COSTS and +% from the best simplex at that iteration in OUTPUT.COST_BEST, the amount of time +% needed in OUTPUT.TIME and the options used in OUTPUT.OPTIONS. +% +% See also SIMPSASET, SIMPSAGET + + +% Copyright (C) 2005 Henning Schmidt, FCC, henning@fcc.chalmers.se +% Copyright (C) 2006 Brecht Donckels, BIOMATH, brecht.donckels@ugent.be +% Copyright (C) 2013 Dynare Team. +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +% handle variable input arguments + +if nargin < 5, + OPTIONS = []; + if nargin < 4, + UB = 1e5; + if nargin < 3, + LB = -1e5; + end + end +end + +% check input arguments + +if ~ischar(FUN), + error('''FUN'' incorrectly specified in ''SIMPSA'''); +end +if ~isfloat(X0), + error('''X0'' incorrectly specified in ''SIMPSA'''); +end +if ~isfloat(LB), + error('''LB'' incorrectly specified in ''SIMPSA'''); +end +if ~isfloat(UB), + error('''UB'' incorrectly specified in ''SIMPSA'''); +end +if length(X0) ~= length(LB), + error('''LB'' and ''X0'' have incompatible dimensions in ''SIMPSA'''); +end +if length(X0) ~= length(UB), + error('''UB'' and ''X0'' have incompatible dimensions in ''SIMPSA'''); +end + +% declaration of global variables + +global NDIM nFUN_EVALS TEMP YBEST PBEST + +% set EXITFLAG to default value + +EXITFLAG = -2; + +% determine number of variables to be optimized + +NDIM = length(X0); + +% seed the random number generator + +rand('state',sum(100*clock)); + +% set default options + +DEFAULT_OPTIONS = simpsaset('TEMP_START',[],... % starting temperature (if none provided, an optimal one will be estimated) + 'TEMP_END',1,... % end temperature + 'COOL_RATE',10,... % small values (<1) means slow convergence,large values (>1) means fast convergence + 'INITIAL_ACCEPTANCE_RATIO',0.95,... % when initial temperature is estimated, this will be the initial acceptance ratio in the first round + 'MIN_COOLING_FACTOR',0.9,... % minimum cooling factor (<1) + 'MAX_ITER_TEMP_FIRST',50,... % number of iterations in the preliminary temperature loop + 'MAX_ITER_TEMP_LAST',50,... % number of iterations in the last temperature loop (pure simplex) + 'MAX_ITER_TEMP',10,... % number of iterations in the remaining temperature loops + 'MAX_ITER_TOTAL',2500,... % maximum number of iterations tout court + 'MAX_TIME',2500,... % maximum duration of optimization + 'MAX_FUN_EVALS',2500,... % maximum number of function evaluations + 'TOLX',1e-6,... % maximum difference between best and worst function evaluation in simplex + 'TOLFUN',1e-3,... % maximum difference between the coordinates of the vertices + 'DISPLAY','none',... % 'iter' or 'none' indicating whether user wants feedback + 'OUTPUT_FCN',[]); % string with output function name + +% update default options with supplied options + +OPTIONS = simpsaset(DEFAULT_OPTIONS,OPTIONS); + +% store options in OUTPUT + +OUTPUT.OPTIONS = OPTIONS; + +% initialize simplex +% ------------------ + +% create empty simplex matrix p (location of vertex i in row i) +P = zeros(NDIM+1,NDIM); +% create empty cost vector (cost of vertex i in row i) +Y = zeros(NDIM+1,1); +% set best vertex of initial simplex equal to initial parameter guess +PBEST = X0(:)'; +% calculate cost with best vertex of initial simplex +YBEST = CALCULATE_COST(FUN,PBEST,LB,UB,varargin{:}); + +% initialize temperature loop +% --------------------------- + +% set temperature loop number to one +TEMP_LOOP_NUMBER = 1; + +% if no TEMP_START is supplied, the initial temperature is estimated in the first +% loop as described by Cardoso et al., 1996 (recommended) + +% therefore, the temperature is set to YBEST*1e5 in the first loop +if isempty(OPTIONS.TEMP_START), + TEMP = abs(YBEST)*1e5; +else + TEMP = OPTIONS.TEMP_START; +end + +% initialize OUTPUT structure +% --------------------------- + +OUTPUT.TEMPERATURE = zeros(OPTIONS.MAX_ITER_TOTAL,1); +OUTPUT.SIMPLEX = zeros(NDIM+1,NDIM,OPTIONS.MAX_ITER_TOTAL); +OUTPUT.SIMPLEX_BEST = zeros(OPTIONS.MAX_ITER_TOTAL,NDIM); +OUTPUT.COSTS = zeros(OPTIONS.MAX_ITER_TOTAL,NDIM+1); +OUTPUT.COST_BEST = zeros(OPTIONS.MAX_ITER_TOTAL,1); + +% initialize iteration data +% ------------------------- + +% start timer +tic +% set number of function evaluations to one +nFUN_EVALS = 1; +% set number of iterations to zero +nITERATIONS = 0; + +% temperature loop: run SIMPSA till stopping criterion is met +% ----------------------------------------------------------- + +while 1, + + % detect if termination criterium was met + % --------------------------------------- + + % if a termination criterium was met, the value of EXITFLAG should have changed + % from its default value of -2 to -1, 0, 1 or 2 + + if EXITFLAG ~= -2, + break + end + + % set MAXITERTEMP: maximum number of iterations at current temperature + % -------------------------------------------------------------------- + + if TEMP_LOOP_NUMBER == 1, + MAXITERTEMP = OPTIONS.MAX_ITER_TEMP_FIRST*NDIM; + % The initial temperature is estimated (is requested) as described in + % Cardoso et al. (1996). Therefore, we need to store the number of + % successful and unsuccessful moves, as well as the increase in cost + % for the unsuccessful moves. + if isempty(OPTIONS.TEMP_START), + [SUCCESSFUL_MOVES,UNSUCCESSFUL_MOVES,UNSUCCESSFUL_COSTS] = deal(0); + end + elseif TEMP < OPTIONS.TEMP_END, + TEMP = 0; + MAXITERTEMP = OPTIONS.MAX_ITER_TEMP_LAST*NDIM; + else + MAXITERTEMP = OPTIONS.MAX_ITER_TEMP*NDIM; + end + + % construct initial simplex + % ------------------------- + + % 1st vertex of initial simplex + P(1,:) = PBEST; + Y(1) = CALCULATE_COST(FUN,P(1,:),LB,UB,varargin{:}); + + % if output function given then run output function to plot intermediate result + if ~isempty(OPTIONS.OUTPUT_FCN), + feval(OPTIONS.OUTPUT_FCN,P(1,:),Y(1)); + end + + % remaining vertices of simplex + for k = 1:NDIM, + % copy first vertex in new vertex + P(k+1,:) = P(1,:); + % alter new vertex + P(k+1,k) = LB(k)+rand*(UB(k)-LB(k)); + % calculate value of objective function at new vertex + Y(k+1) = CALCULATE_COST(FUN,P(k+1,:),LB,UB,varargin{:}); + end + + % store information on what step the algorithm just did + ALGOSTEP = 'initial simplex'; + + % add NDIM+1 to number of function evaluations + nFUN_EVALS = nFUN_EVALS + NDIM; + + % note: + % dimensions of matrix P: (NDIM+1) x NDIM + % dimensions of vector Y: (NDIM+1) x 1 + + % give user feedback if requested + if strcmp(OPTIONS.DISPLAY,'iter'), + if nITERATIONS == 0, + disp(' Nr Iter Nr Fun Eval Min function Best function TEMP Algorithm Step'); + else + disp(sprintf('%5.0f %5.0f %12.6g %15.6g %12.6g %s',nITERATIONS,nFUN_EVALS,Y(1),YBEST,TEMP,'best point')); + end + end + + % run full metropolis cycle at current temperature + % ------------------------------------------------ + + % initialize vector COSTS, needed to calculate new temperature using cooling + % schedule as described by Cardoso et al. (1996) + COSTS = zeros((NDIM+1)*MAXITERTEMP,1); + + % initialize ITERTEMP to zero + + ITERTEMP = 0; + + % start + + for ITERTEMP = 1:MAXITERTEMP, + + % add one to number of iterations + nITERATIONS = nITERATIONS + 1; + + % Press and Teukolsky (1991) add a positive logarithmic distributed variable, + % proportional to the control temperature T to the function value associated with + % every vertex of the simplex. Likewise,they subtract a similar random variable + % from the function value at every new replacement point. + % Thus, if the replacement point corresponds to a lower cost, this method always + % accepts a true down hill step. If, on the other hand, the replacement point + % corresponds to a higher cost, an uphill move may be accepted, depending on the + % relative COSTS of the perturbed values. + % (taken from Cardoso et al.,1996) + + % add random fluctuations to function values of current vertices + YFLUCT = Y+TEMP*abs(log(rand(NDIM+1,1))); + + % reorder YFLUCT, Y and P so that the first row corresponds to the lowest YFLUCT value + help = sortrows([YFLUCT,Y,P],1); + YFLUCT = help(:,1); + Y = help(:,2); + P = help(:,3:end); + + % store temperature at current iteration + OUTPUT.TEMPERATURE(nITERATIONS) = TEMP; + + % store information about simplex at the current iteration + OUTPUT.SIMPLEX(:,:,nITERATIONS) = P; + OUTPUT.SIMPLEX_BEST(nITERATIONS,:) = PBEST; + + % store cost function value of best vertex in current iteration + OUTPUT.COSTS(nITERATIONS,:) = Y; + OUTPUT.COST_BEST(nITERATIONS) = YBEST; + + if strcmp(OPTIONS.DISPLAY,'iter'), + disp(sprintf('%5.0f %5.0f %12.6g %15.6g %12.6g %s',nITERATIONS,nFUN_EVALS,Y(1),YBEST,TEMP,ALGOSTEP)); + end + + % if output function given then run output function to plot intermediate result + if ~isempty(OPTIONS.OUTPUT_FCN), + feval(OPTIONS.OUTPUT_FCN,P(1,:),Y(1)); + end + + % end the optimization if one of the stopping criteria is met + %% 1. difference between best and worst function evaluation in simplex is smaller than TOLFUN + %% 2. maximum difference between the coordinates of the vertices in simplex is less than TOLX + %% 3. no convergence,but maximum number of iterations has been reached + %% 4. no convergence,but maximum time has been reached + + if (abs(max(Y)-min(Y)) < OPTIONS.TOLFUN) && (TEMP_LOOP_NUMBER ~= 1), + if strcmp(OPTIONS.DISPLAY,'iter'), + disp('Change in the objective function value less than the specified tolerance (TOLFUN).') + end + EXITFLAG = 1; + break; + end + + if (max(max(abs(P(2:NDIM+1,:)-P(1:NDIM,:)))) < OPTIONS.TOLX) && (TEMP_LOOP_NUMBER ~= 1), + if strcmp(OPTIONS.DISPLAY,'iter'), + disp('Change in X less than the specified tolerance (TOLX).') + end + EXITFLAG = 2; + break; + end + + if (nITERATIONS >= OPTIONS.MAX_ITER_TOTAL*NDIM) || (nFUN_EVALS >= OPTIONS.MAX_FUN_EVALS*NDIM*(NDIM+1)), + if strcmp(OPTIONS.DISPLAY,'iter'), + disp('Maximum number of function evaluations or iterations reached.'); + end + EXITFLAG = 0; + break; + end + + if toc/60 > OPTIONS.MAX_TIME, + if strcmp(OPTIONS.DISPLAY,'iter'), + disp('Exceeded maximum time.'); + end + EXITFLAG = -1; + break; + end + + % begin a new iteration + + %% first extrapolate by a factor -1 through the face of the simplex + %% across from the high point,i.e.,reflect the simplex from the high point + [YFTRY,YTRY,PTRY] = AMOTRY(FUN,P,-1,LB,UB,varargin{:}); + + %% check the result + if YFTRY <= YFLUCT(1), + %% gives a result better than the best point,so try an additional + %% extrapolation by a factor 2 + [YFTRYEXP,YTRYEXP,PTRYEXP] = AMOTRY(FUN,P,-2,LB,UB,varargin{:}); + if YFTRYEXP < YFTRY, + P(end,:) = PTRYEXP; + Y(end) = YTRYEXP; + ALGOSTEP = 'reflection and expansion'; + else + P(end,:) = PTRY; + Y(end) = YTRY; + ALGOSTEP = 'reflection'; + end + elseif YFTRY >= YFLUCT(NDIM), + %% the reflected point is worse than the second-highest, so look + %% for an intermediate lower point, i.e., do a one-dimensional + %% contraction + [YFTRYCONTR,YTRYCONTR,PTRYCONTR] = AMOTRY(FUN,P,-0.5,LB,UB,varargin{:}); + if YFTRYCONTR < YFLUCT(end), + P(end,:) = PTRYCONTR; + Y(end) = YTRYCONTR; + ALGOSTEP = 'one dimensional contraction'; + else + %% can't seem to get rid of that high point, so better contract + %% around the lowest (best) point + X = ones(NDIM,NDIM)*diag(P(1,:)); + P(2:end,:) = 0.5*(P(2:end,:)+X); + for k=2:NDIM, + Y(k) = CALCULATE_COST(FUN,P(k,:),LB,UB,varargin{:}); + end + ALGOSTEP = 'multiple contraction'; + end + else + %% if YTRY better than second-highest point, use this point + P(end,:) = PTRY; + Y(end) = YTRY; + ALGOSTEP = 'reflection'; + end + + % the initial temperature is estimated in the first loop from + % the number of successfull and unsuccesfull moves, and the average + % increase in cost associated with the unsuccessful moves + + if TEMP_LOOP_NUMBER == 1 && isempty(OPTIONS.TEMP_START), + if Y(1) > Y(end), + SUCCESSFUL_MOVES = SUCCESSFUL_MOVES+1; + elseif Y(1) <= Y(end), + UNSUCCESSFUL_MOVES = UNSUCCESSFUL_MOVES+1; + UNSUCCESSFUL_COSTS = UNSUCCESSFUL_COSTS+(Y(end)-Y(1)); + end + end + + end + + % stop if previous for loop was broken due to some stop criterion + if ITERTEMP < MAXITERTEMP, + break; + end + + % store cost function values in COSTS vector + COSTS((ITERTEMP-1)*NDIM+1:ITERTEMP*NDIM+1) = Y; + + % calculated initial temperature or recalculate temperature + % using cooling schedule as proposed by Cardoso et al. (1996) + % ----------------------------------------------------------- + + if TEMP_LOOP_NUMBER == 1 && isempty(OPTIONS.TEMP_START), + TEMP = -(UNSUCCESSFUL_COSTS/(SUCCESSFUL_MOVES+UNSUCCESSFUL_MOVES))/log(((SUCCESSFUL_MOVES+UNSUCCESSFUL_MOVES)*OPTIONS.INITIAL_ACCEPTANCE_RATIO-SUCCESSFUL_MOVES)/UNSUCCESSFUL_MOVES); + elseif TEMP_LOOP_NUMBER ~= 0, + STDEV_Y = std(COSTS); + COOLING_FACTOR = 1/(1+TEMP*log(1+OPTIONS.COOL_RATE)/(3*STDEV_Y)); + TEMP = TEMP*min(OPTIONS.MIN_COOLING_FACTOR,COOLING_FACTOR); + end + + % add one to temperature loop number + TEMP_LOOP_NUMBER = TEMP_LOOP_NUMBER+1; + +end + +% return solution +X = PBEST; +FVAL = YBEST; + +% store number of function evaluations +OUTPUT.nFUN_EVALS = nFUN_EVALS; + +% store number of iterations +OUTPUT.nITERATIONS = nITERATIONS; + +% trim OUTPUT data structure +OUTPUT.TEMPERATURE(nITERATIONS+1:end) = []; +OUTPUT.SIMPLEX(:,:,nITERATIONS+1:end) = []; +OUTPUT.SIMPLEX_BEST(nITERATIONS+1:end,:) = []; +OUTPUT.COSTS(nITERATIONS+1:end,:) = []; +OUTPUT.COST_BEST(nITERATIONS+1:end) = []; + +% store the amount of time needed in OUTPUT data structure +OUTPUT.TIME = toc; + +return + +% ============================================================================== + +% AMOTRY FUNCTION +% --------------- + +function [YFTRY,YTRY,PTRY] = AMOTRY(FUN,P,fac,LB,UB,varargin) +% Extrapolates by a factor fac through the face of the simplex across from +% the high point, tries it, and replaces the high point if the new point is +% better. + +global NDIM TEMP + +% calculate coordinates of new vertex +psum = sum(P(1:NDIM,:))/NDIM; +PTRY = psum*(1-fac)+P(end,:)*fac; + +% evaluate the function at the trial point. +YTRY = CALCULATE_COST(FUN,PTRY,LB,UB,varargin{:}); +% substract random fluctuations to function values of current vertices +YFTRY = YTRY-TEMP*abs(log(rand(1))); + +return + +% ============================================================================== + +% COST FUNCTION EVALUATION +% ------------------------ + +function [YTRY] = CALCULATE_COST(FUN,PTRY,LB,UB,varargin) + +global YBEST PBEST NDIM nFUN_EVALS + +for i = 1:NDIM, + % check lower bounds + if PTRY(i) < LB(i), + YTRY = 1e12+(LB(i)-PTRY(i))*1e6; + return + end + % check upper bounds + if PTRY(i) > UB(i), + YTRY = 1e12+(PTRY(i)-UB(i))*1e6; + return + end +end + +% calculate cost associated with PTRY +YTRY = feval(FUN,PTRY,varargin{:}); + +% add one to number of function evaluations +nFUN_EVALS = nFUN_EVALS + 1; + +% save the best point ever +if YTRY < YBEST, + YBEST = YTRY; + PBEST = PTRY; +end + +return diff --git a/matlab/simpsaget.m b/matlab/simpsaget.m new file mode 100644 index 000000000..b5014ac1c --- /dev/null +++ b/matlab/simpsaget.m @@ -0,0 +1,125 @@ +function o = SIMPSAGET(options,name,default,flag) +%SIMPSAGET Get SIMPSA OPTIONS parameters. +% VAL = SIMPSAGET(OPTIONS,'NAME') extracts the value of the named parameter +% from optimization options structure OPTIONS, returning an empty matrix if +% the parameter value is not specified in OPTIONS. It is sufficient to +% type only the leading characters that uniquely identify the +% parameter. Case is ignored for parameter names. [] is a valid OPTIONS +% argument. +% +% VAL = SIMPSAGET(OPTIONS,'NAME',DEFAULT) extracts the named parameter as +% above, but returns DEFAULT if the named parameter is not specified (is []) +% in OPTIONS. For example +% +% val = simpsaget(opts,'TolX',1e-4); +% +% returns val = 1e-4 if the TolX property is not specified in opts. +% +% See also SIMPSASET, SIMPSA +% +% Copyright (C) 2006 Brecht Donckels, BIOMATH, brecht.donckels@ugent.be +% +% This program 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 2 +% of the License, or (at your option) any later version. +% +% This program 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 this program; if not, write to the Free Software +% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, +% USA. + + + +% undocumented usage for fast access with no error checking +if (nargin == 4) && isequal(flag,'fast') + o = getknownfield(options,name,default); + return +end + +if nargin < 2 + error('MATLAB:odeget:NotEnoughInputs','Not enough input arguments.'); +end +if nargin < 3 + default = []; +end + +if ~isempty(options) && ~isa(options,'struct') + error('MATLAB:odeget:Arg1NotODESETstruct',... + 'First argument must be an options structure created with ODESET.'); +end + +if isempty(options) + o = default; + return; +end + +Names = [ + 'TEMP_START ' + 'TEMP_END ' + 'COOL_RATE ' + 'INITIAL_ACCEPTANCE_RATIO ' + 'MIN_COOLING_FACTOR ' + 'MAX_ITER_TEMP_FIRST ' + 'MAX_ITER_TEMP_LAST ' + 'MAX_ITER_TEMP ' + 'MAX_ITER_TOTAL ' + 'MAX_TIME ' + 'MAX_FUN_EVALS ' + 'TOLX ' + 'TOLFUN ' + 'DISPLAY ' + 'OUTPUT_FCN ' + ]; + +names = lower(Names); + +lowName = lower(name); +j = strmatch(lowName,names); +if isempty(j) % if no matches + error('MATLAB:odeget:InvalidPropName',... + ['Unrecognized property name ''%s''. ' ... + 'See ODESET for possibilities.'], name); +elseif length(j) > 1 % if more than one match + % Check for any exact matches (in case any names are subsets of others) + k = strmatch(lowName,names,'exact'); + if length(k) == 1 + j = k; + else + msg = sprintf('Ambiguous property name ''%s'' ', name); + msg = [msg '(' deblank(Names(j(1),:))]; + for k = j(2:length(j))' + msg = [msg ', ' deblank(Names(k,:))]; + end + msg = sprintf('%s).', msg); + error('MATLAB:odeget:AmbiguousPropName', msg); + end +end + +if any(strcmp(fieldnames(options),deblank(Names(j,:)))) + o = options.(deblank(Names(j,:))); + if isempty(o) + o = default; + end +else + o = default; +end + +% -------------------------------------------------------------------------- +function v = getknownfield(s, f, d) +%GETKNOWNFIELD Get field f from struct s, or else yield default d. + +if isfield(s,f) % s could be empty. + v = subsref(s, struct('type','.','subs',f)); + if isempty(v) + v = d; + end +else + v = d; +end + diff --git a/matlab/simpsaset.m b/matlab/simpsaset.m new file mode 100644 index 000000000..ca687891a --- /dev/null +++ b/matlab/simpsaset.m @@ -0,0 +1,172 @@ +function options = simpsaset(varargin) + +%SIMPSASET Create/alter simpsa optimization OPTIONS structure. +% OPTIONS = SIMPSASET('PARAM1',VALUE1,'PARAM2',VALUE2,...) creates an +% optimization options structure OPTIONS in which the named parameters have +% the specified values. Any unspecified parameters are set to [] (parameters +% with value [] indicate to use the default value for that parameter when +% OPTIONS is passed to the optimization function). It is sufficient to type +% only the leading characters that uniquely identify the parameter. Case is +% ignored for parameter names. +% NOTE: For values that are strings, the complete string is required. +% +% OPTIONS = SIMPSASET(OLDOPTS,'PARAM1',VALUE1,...) creates a copy of OLDOPTS +% with the named parameters altered with the specified values. +% +% OPTIONS = SIMPSASET(OLDOPTS,NEWOPTS) combines an existing options structure +% OLDOPTS with a new options structure NEWOPTS. Any parameters in NEWOPTS +% with non-empty values overwrite the corresponding old parameters in +% OLDOPTS. +% +% SIMPSASET with no input arguments and no output arguments displays all +% parameter names and their possible values, with defaults shown in {} +% when the default is the same for all functions that use that option -- use +% SIMPSASET(OPTIMFUNCTION) to see options for a specific function. +% +% OPTIONS = SIMPSASET (with no input arguments) creates an options structure +% OPTIONS where all the fields are set to []. + +% Copyright (C) 2005 Henning Schmidt, FCC, henning@fcc.chalmers.se +% Copyright (C) 2006 Brecht Donckels, BIOMATH, brecht.donckels@ugent.be +% Copyright (C) 2013 Dynare Team. +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + + + +% Print out possible values of properties. +if (nargin == 0) && (nargout == 0) + fprintf(' TEMP_START: [ positive scalar ]\n'); + fprintf(' TEMP_END: [ positive scalar ]\n'); + fprintf(' COOL_RATE: [ positive scalar ]\n'); + fprintf(' INITIAL_ACCEPTANCE_RATIO: [ positive scalar < 1 {0.95} ]\n'); + fprintf(' MIN_COOLING_FACTOR: [ positive scalar < 1 {0.9}]\n'); + fprintf(' MAX_ITER_TEMP_FIRST: [ positive scalar {100} ]\n'); + fprintf(' MAX_ITER_TEMP_LAST: [ positive scalar {20} ]\n'); + fprintf(' MAX_ITER_TOTAL: [ positive scalar {2500} ]\n'); + fprintf(' MAX_TIME: [ positive scalar {2500} ]\n'); + fprintf(' MAX_FUN_EVALS: [ positive scalar {2500} ]\n'); + fprintf(' TOLX: [ positive scalar {1e-6} ]\n'); + fprintf(' TOLFUN: [ positive scalar {1e-6} ]\n'); + fprintf(' DISPLAY: [ ''iter'' or ''none'' {''iter''} ]\n'); + fprintf(' OUTPUT_FCN: [ function_handle ]\n'); + fprintf('\n'); +return; +end + +Names = [ + 'TEMP_START ' + 'TEMP_END ' + 'COOL_RATE ' + 'INITIAL_ACCEPTANCE_RATIO ' + 'MIN_COOLING_FACTOR ' + 'MAX_ITER_TEMP_FIRST ' + 'MAX_ITER_TEMP_LAST ' + 'MAX_ITER_TEMP ' + 'MAX_ITER_TOTAL ' + 'MAX_TIME ' + 'MAX_FUN_EVALS ' + 'TOLX ' + 'TOLFUN ' + 'DISPLAY ' + 'OUTPUT_FCN ' + ]; + +m = size(Names,1); +names = lower(Names); + +% Combine all leading options structures o1, o2, ... in odeset(o1,o2,...). +options = []; +for j = 1:m + options.(deblank(Names(j,:))) = []; +end +i = 1; +while i <= nargin + arg = varargin{i}; + if ischar(arg) % arg is an option name + break; + end + if ~isempty(arg) % [] is a valid options argument + if ~isa(arg,'struct') + error('MATLAB:odeset:NoPropNameOrStruct',... + ['Expected argument %d to be a string property name ' ... + 'or an options structure\ncreated with SIMANSET.'], i); + end + for j = 1:m + if any(strcmp(fieldnames(arg),deblank(Names(j,:)))) + val = arg.(deblank(Names(j,:))); + else + val = []; + end + if ~isempty(val) + options.(deblank(Names(j,:))) = val; + end + end + end + i = i + 1; +end + +% A finite state machine to parse name-value pairs. +if rem(nargin-i+1,2) ~= 0 + error('MATLAB:odeset:ArgNameValueMismatch',... + 'Arguments must occur in name-value pairs.'); +end +expectval = 0; % start expecting a name, not a value +while i <= nargin + arg = varargin{i}; + + if ~expectval + if ~ischar(arg) + error('MATLAB:odeset:NoPropName',... + 'Expected argument %d to be a string property name.', i); + end + + lowArg = lower(arg); + j = strmatch(lowArg,names); + if isempty(j) % if no matches + error('MATLAB:odeset:InvalidPropName',... + 'Unrecognized property name ''%s''.', arg); + elseif length(j) > 1 % if more than one match + % Check for any exact matches (in case any names are subsets of others) + k = strmatch(lowArg,names,'exact'); + if length(k) == 1 + j = k; + else + msg = sprintf('Ambiguous property name ''%s'' ', arg); + msg = [msg '(' deblank(Names(j(1),:))]; + for k = j(2:length(j))' + msg = [msg ', ' deblank(Names(k,:))]; + end + msg = sprintf('%s).', msg); + error('MATLAB:odeset:AmbiguousPropName', msg); + end + end + expectval = 1; % we expect a value next + + else + options.(deblank(Names(j,:))) = arg; + expectval = 0; + + end + i = i + 1; +end + +if expectval + error('MATLAB:odeset:NoValueForProp',... + 'Expected value for property ''%s''.', arg); +end + +end From 3f16129e4958af04c37a75a424ca75a698fa8ca2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Wed, 2 Oct 2013 15:30:55 +0200 Subject: [PATCH 32/80] mode_compute = 10 calls simpsa algorithm. --- matlab/dynare_estimation_1.m | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index c9d489d76..1782a8b7e 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -407,7 +407,9 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation [x, fval, COUNTEVAL, STOPFLAG, OUT, BESTEVER] = cmaes(func2str(objective_function),xparam1,H0,options_.cmaes,dataset_,options_,M_,estim_params_,bayestopt_,oo_); xparam1=BESTEVER.x; disp(sprintf('\n Objective function at mode: %f',fval)) - case 10 + case 10 + [xparam1, fval, exitflag] = simpsa(func2str(objective_function),xparam1,-inf(size(xparam1)),inf(size(xparam1)),[],dataset_,options_,M_,estim_params_,bayestopt_,oo_); + case 11 options_.cova_compute = 0 ; [xparam1,stdh,lb_95,ub_95,med_param] = online_auxiliary_filter(xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_) ; case 101 From fb0ccdd5d26edc41e61581fdd3e1eb572b1ccb00 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Wed, 2 Oct 2013 16:45:16 +0200 Subject: [PATCH 33/80] Provide more sensible lower and upper bounds for simpsa algorithm. --- matlab/dynare_estimation_1.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index 1782a8b7e..e5ad0d0ba 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -408,7 +408,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation xparam1=BESTEVER.x; disp(sprintf('\n Objective function at mode: %f',fval)) case 10 - [xparam1, fval, exitflag] = simpsa(func2str(objective_function),xparam1,-inf(size(xparam1)),inf(size(xparam1)),[],dataset_,options_,M_,estim_params_,bayestopt_,oo_); + [xparam1, fval, exitflag] = simpsa(func2str(objective_function),xparam1,lb,ub,[],dataset_,options_,M_,estim_params_,bayestopt_,oo_); case 11 options_.cova_compute = 0 ; [xparam1,stdh,lb_95,ub_95,med_param] = online_auxiliary_filter(xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_) ; From ab7f0ddedc0bbac3084a4009608e043635412669 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Wed, 2 Oct 2013 17:07:00 +0200 Subject: [PATCH 34/80] Removed the seed set on clock. Changed default values for the options. --- matlab/simpsa.m | 15 +++++---------- 1 file changed, 5 insertions(+), 10 deletions(-) diff --git a/matlab/simpsa.m b/matlab/simpsa.m index 97e5c2fd3..be654568d 100644 --- a/matlab/simpsa.m +++ b/matlab/simpsa.m @@ -120,26 +120,21 @@ EXITFLAG = -2; NDIM = length(X0); -% seed the random number generator - -rand('state',sum(100*clock)); - % set default options - DEFAULT_OPTIONS = simpsaset('TEMP_START',[],... % starting temperature (if none provided, an optimal one will be estimated) - 'TEMP_END',1,... % end temperature + 'TEMP_END',.1,... % end temperature 'COOL_RATE',10,... % small values (<1) means slow convergence,large values (>1) means fast convergence 'INITIAL_ACCEPTANCE_RATIO',0.95,... % when initial temperature is estimated, this will be the initial acceptance ratio in the first round 'MIN_COOLING_FACTOR',0.9,... % minimum cooling factor (<1) 'MAX_ITER_TEMP_FIRST',50,... % number of iterations in the preliminary temperature loop - 'MAX_ITER_TEMP_LAST',50,... % number of iterations in the last temperature loop (pure simplex) + 'MAX_ITER_TEMP_LAST',2000,... % number of iterations in the last temperature loop (pure simplex) 'MAX_ITER_TEMP',10,... % number of iterations in the remaining temperature loops 'MAX_ITER_TOTAL',2500,... % maximum number of iterations tout court 'MAX_TIME',2500,... % maximum duration of optimization - 'MAX_FUN_EVALS',2500,... % maximum number of function evaluations + 'MAX_FUN_EVALS',20000,... % maximum number of function evaluations 'TOLX',1e-6,... % maximum difference between best and worst function evaluation in simplex - 'TOLFUN',1e-3,... % maximum difference between the coordinates of the vertices - 'DISPLAY','none',... % 'iter' or 'none' indicating whether user wants feedback + 'TOLFUN',1e-6,... % maximum difference between the coordinates of the vertices + 'DISPLAY','iter',... % 'iter' or 'none' indicating whether user wants feedback 'OUTPUT_FCN',[]); % string with output function name % update default options with supplied options From 1ff6b2de7b158eec68d0e3e62c913e88f95bb30c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Wed, 2 Oct 2013 17:08:04 +0200 Subject: [PATCH 35/80] Fixed bug caused by the size of the vector of parameters (transposition needed). --- matlab/simpsa.m | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/matlab/simpsa.m b/matlab/simpsa.m index be654568d..260607a24 100644 --- a/matlab/simpsa.m +++ b/matlab/simpsa.m @@ -235,7 +235,7 @@ while 1, % if output function given then run output function to plot intermediate result if ~isempty(OPTIONS.OUTPUT_FCN), - feval(OPTIONS.OUTPUT_FCN,P(1,:),Y(1)); + feval(OPTIONS.OUTPUT_FCN,transpose(P(1,:)),Y(1)); end % remaining vertices of simplex @@ -321,7 +321,7 @@ while 1, % if output function given then run output function to plot intermediate result if ~isempty(OPTIONS.OUTPUT_FCN), - feval(OPTIONS.OUTPUT_FCN,P(1,:),Y(1)); + feval(OPTIONS.OUTPUT_FCN,transpose(P(1,:)),Y(1)); end % end the optimization if one of the stopping criteria is met @@ -449,7 +449,7 @@ while 1, end % return solution -X = PBEST; +X = transpose(PBEST); FVAL = YBEST; % store number of function evaluations @@ -516,7 +516,7 @@ for i = 1:NDIM, end % calculate cost associated with PTRY -YTRY = feval(FUN,PTRY,varargin{:}); +YTRY = feval(FUN,PTRY(:),varargin{:}); % add one to number of function evaluations nFUN_EVALS = nFUN_EVALS + 1; From 581f97badaec545a9370ef4bd1c3a4814f3b9845 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Wed, 2 Oct 2013 17:09:05 +0200 Subject: [PATCH 36/80] Linked simpsa's tolerance options to dynare's defaults. --- matlab/dynare_estimation_1.m | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index e5ad0d0ba..14751c988 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -408,7 +408,8 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation xparam1=BESTEVER.x; disp(sprintf('\n Objective function at mode: %f',fval)) case 10 - [xparam1, fval, exitflag] = simpsa(func2str(objective_function),xparam1,lb,ub,[],dataset_,options_,M_,estim_params_,bayestopt_,oo_); + options = simpsaset('TOLX', options_.dynatol.x,'TOLFUN', options_.dynatol.f) + [xparam1, fval, exitflag] = simpsa(func2str(objective_function),xparam1,lb,ub,options,dataset_,options_,M_,estim_params_,bayestopt_,oo_); case 11 options_.cova_compute = 0 ; [xparam1,stdh,lb_95,ub_95,med_param] = online_auxiliary_filter(xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_) ; From a34afdfdc3382c15c288a115277961a08c55aec7 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Wed, 2 Oct 2013 16:12:31 +0200 Subject: [PATCH 37/80] support negative dates --- preprocessor/DynareFlex.ll | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll index 35616b1df..9e3a9d554 100644 --- a/preprocessor/DynareFlex.ll +++ b/preprocessor/DynareFlex.ll @@ -691,7 +691,7 @@ string eofbuff; return token::INT_NUMBER; } -([1-2][0-9]{3}[Mm](([1-9])|(1[0-2])))|([1-2][0-9]{3}[Qq][1-4])|([1-2][0-9]{3}[Ww](([1-9]{1})|([1-4][0-9])|([5][0-2]))) { +(-?[1-2][0-9]{3}[Mm](([1-9])|(1[0-2])))|(-?[1-2][0-9]{3}[Qq][1-4])|(-?[1-2][0-9]{3}[Ww](([1-9]{1})|([1-4][0-9])|([5][0-2]))) { yylval->string_val = new string(yytext); return token::DATE_NUMBER; } From 845a8736b58bab56a7e54db50dcec0a4c7b710e9 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Wed, 2 Oct 2013 16:16:33 +0200 Subject: [PATCH 38/80] simplify regex --- preprocessor/DynareFlex.ll | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll index 9e3a9d554..2f2b89b4a 100644 --- a/preprocessor/DynareFlex.ll +++ b/preprocessor/DynareFlex.ll @@ -691,7 +691,7 @@ string eofbuff; return token::INT_NUMBER; } -(-?[1-2][0-9]{3}[Mm](([1-9])|(1[0-2])))|(-?[1-2][0-9]{3}[Qq][1-4])|(-?[1-2][0-9]{3}[Ww](([1-9]{1})|([1-4][0-9])|([5][0-2]))) { +(-?[1-2][0-9]{3}[Mm](([1-9])|(1[0-2])))|(-?[1-2][0-9]{3}[Qq][1-4])|(-?[1-2][0-9]{3}[Ww](([1-9]{1})|([1-4][0-9])|(5[0-2]))) { yylval->string_val = new string(yytext); return token::DATE_NUMBER; } From c1305f7f87058a969664cc499c405df36002e68a Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Wed, 2 Oct 2013 16:48:34 +0200 Subject: [PATCH 39/80] support unbounded years --- preprocessor/DynareFlex.ll | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll index 2f2b89b4a..2193ff9d3 100644 --- a/preprocessor/DynareFlex.ll +++ b/preprocessor/DynareFlex.ll @@ -691,7 +691,7 @@ string eofbuff; return token::INT_NUMBER; } -(-?[1-2][0-9]{3}[Mm](([1-9])|(1[0-2])))|(-?[1-2][0-9]{3}[Qq][1-4])|(-?[1-2][0-9]{3}[Ww](([1-9]{1})|([1-4][0-9])|(5[0-2]))) { +(-?[0-9]+[Mm](([1-9])|(1[0-2])))|(-?[0-9]+[Qq][1-4])|(-?[0-9]+[Ww](([1-9]{1})|([1-4][0-9])|(5[0-2]))) { yylval->string_val = new string(yytext); return token::DATE_NUMBER; } From 62c623f6c9ac91987a8d710222c8e8aabd32ed02 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Wed, 2 Oct 2013 16:56:38 +0200 Subject: [PATCH 40/80] remove unnecessary parenthesis --- preprocessor/DynareFlex.ll | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll index 2193ff9d3..6d7f3ab6c 100644 --- a/preprocessor/DynareFlex.ll +++ b/preprocessor/DynareFlex.ll @@ -691,7 +691,7 @@ string eofbuff; return token::INT_NUMBER; } -(-?[0-9]+[Mm](([1-9])|(1[0-2])))|(-?[0-9]+[Qq][1-4])|(-?[0-9]+[Ww](([1-9]{1})|([1-4][0-9])|(5[0-2]))) { +-?[0-9]+[Mm]([1-9]|1[0-2])|-?[0-9]+[Qq][1-4]|-?[0-9]+[Ww]([1-9]{1}|[1-4][0-9]|5[0-2]) { yylval->string_val = new string(yytext); return token::DATE_NUMBER; } From 030fe52affb763f11766801f6deeb797e90365e2 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Wed, 2 Oct 2013 16:41:05 +0200 Subject: [PATCH 41/80] preprocessor: replace dates with dynDates --- preprocessor/ParsingDriver.cc | 3 +- preprocessor/Statement.cc | 69 ++++++++++++++++++++++++++++++++++- preprocessor/Statement.hh | 6 ++- 3 files changed, 74 insertions(+), 4 deletions(-) diff --git a/preprocessor/ParsingDriver.cc b/preprocessor/ParsingDriver.cc index 38bf5a782..37029fdcf 100644 --- a/preprocessor/ParsingDriver.cc +++ b/preprocessor/ParsingDriver.cc @@ -2449,7 +2449,8 @@ ParsingDriver::add_model_var_or_external_function(string *function_name, bool in void ParsingDriver::add_native(const string &s) { - mod_file->addStatement(new NativeStatement(s)); + string ss = string(s); + mod_file->addStatement(new NativeStatement(ss)); } void diff --git a/preprocessor/Statement.cc b/preprocessor/Statement.cc index 59021287b..242c73a8f 100644 --- a/preprocessor/Statement.cc +++ b/preprocessor/Statement.cc @@ -18,6 +18,7 @@ */ #include "Statement.hh" +#include ModFileStructure::ModFileStructure() : check_present(false), @@ -65,11 +66,77 @@ Statement::computingPass() { } -NativeStatement::NativeStatement(const string &native_statement_arg) : +NativeStatement::NativeStatement(string &native_statement_arg) : native_statement(native_statement_arg) { } +void +NativeStatement::computingPass() +{ + using namespace boost::xpressive; + // Return if this is a comment + sregex comment_expr = sregex::compile( "\%.*" ); + match_results results; + if (regex_match(native_statement, results, comment_expr)) + return; + + // Otherwise, look at the line and consider substituting date + size_t idx = -1; + vector apostrophes; + while((idx = native_statement.find("'", idx + 1)) != string::npos) + if (apostrophes.size() < 2) + apostrophes.push_back(idx); + else + if (idx == apostrophes.back() + 1) + apostrophes.pop_back(); + else + apostrophes.push_back(idx); + + if (apostrophes.size() % 2) + { + cerr << native_statement << + " seems to be invalid Matlab syntax (an odd number of apostrophes was encountered)" << endl; + exit(EXIT_FAILURE); + } + + bool skip = false; + string newstr = ""; + int lastidx = 0; + sregex date_expr = sregex::compile( "-?[0-9]+[Mm]([1-9]|1[0-2])|-?[0-9]+[Qq][1-4]|-?[0-9]+[Ww]([1-9]{1}|[1-4][0-9]|5[0-2])" ); + string format( "dynDate('$&')" ); + for (size_t i = 0; i < apostrophes.size(); i++) + if (apostrophes[i] == 0) + skip = true; + else + { + if (skip) + { + skip = false; + newstr.append(native_statement.substr(lastidx, apostrophes[i] - lastidx)); + } + else + { + skip = true; + newstr.append(regex_replace(native_statement.substr(lastidx, apostrophes[i] - lastidx), + date_expr, format)); + } + lastidx = apostrophes[i]; + } + + //Replace last (or only) element + if (apostrophes.empty()) + lastidx = 0; + newstr.append(regex_replace(native_statement.substr(lastidx, native_statement.size() - lastidx), + date_expr, format)); + native_statement = newstr; +} + +void +regexReplace() +{ +} + void NativeStatement::writeOutput(ostream &output, const string &basename) const { diff --git a/preprocessor/Statement.hh b/preprocessor/Statement.hh index cd3fd5cf3..82f64fa79 100644 --- a/preprocessor/Statement.hh +++ b/preprocessor/Statement.hh @@ -121,10 +121,12 @@ public: class NativeStatement : public Statement { private: - const string native_statement; + string native_statement; public: - NativeStatement(const string &native_statement_arg); + NativeStatement(string &native_statement_arg); + virtual void computingPass(); virtual void writeOutput(ostream &output, const string &basename) const; + void regexReplace(); }; class OptionsList From b00bf2662183ae0da9989ac4770750514af7875e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Wed, 2 Oct 2013 17:44:47 +0200 Subject: [PATCH 42/80] Added simpsa optimization algorithm in the manual. --- doc/dynare.texi | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/doc/dynare.texi b/doc/dynare.texi index 6902abec9..d2e3bfcad 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -4457,6 +4457,10 @@ available with @code{mode_compute=7}) @item 9 Uses the CMA-ES (Covariance Matrix Adaptation Evolution Strategy) algorithm, an evolutionary algorithm for difficult non-linear non-convex optimization +@item 10 +Uses the simpsa algorithm, based on the combination of the non-linear simplex and simulated annealing algorithms and proposed by +@cite{Cardoso, Salcedo and Feyo de Azevedo (1996)}. + @item @var{FUNCTION_NAME} It is also possible to give a @var{FUNCTION_NAME} to this option, instead of an @var{INTEGER}. In that case, Dynare takes the return @@ -8737,6 +8741,10 @@ Brooks, Stephen P., and Andrew Gelman (1998): ``General methods for monitoring convergence of iterative simulations,'' @i{Journal of computational and graphical statistics}, 7, pp. 434--455 +@item +Cardoso, Margarida F., R. L. Salcedo and S. Feyo de Azevedo (1996): ``The simplex simulated annealing approach to continuous non-linear optimization'', @i{Computers chem. Engng}, 20(9), 1065-1080 + + @item Collard, Fabrice (2001): ``Stochastic simulations with Dynare: A practical guide'' From 979a55a3341c340b03deee67660cdac58183c3fd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Wed, 2 Oct 2013 18:26:57 +0200 Subject: [PATCH 43/80] Added missing semicolon. --- matlab/dynare_estimation_1.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index 14751c988..5e0789491 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -408,7 +408,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation xparam1=BESTEVER.x; disp(sprintf('\n Objective function at mode: %f',fval)) case 10 - options = simpsaset('TOLX', options_.dynatol.x,'TOLFUN', options_.dynatol.f) + options = simpsaset('TOLX', options_.dynatol.x,'TOLFUN', options_.dynatol.f); [xparam1, fval, exitflag] = simpsa(func2str(objective_function),xparam1,lb,ub,options,dataset_,options_,M_,estim_params_,bayestopt_,oo_); case 11 options_.cova_compute = 0 ; From d9e6985e5b0960742081723bf03513bc06be9aea Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= Date: Wed, 2 Oct 2013 18:42:36 +0200 Subject: [PATCH 44/80] Fix license info for simpsa --- license.txt | 6 ++++++ matlab/simpsaget.m | 25 ++++++++++++------------- 2 files changed, 18 insertions(+), 13 deletions(-) diff --git a/license.txt b/license.txt index 95867ea33..c9be4ebec 100644 --- a/license.txt +++ b/license.txt @@ -58,6 +58,12 @@ Copyright: 2011 Lawrence J. Christiano, Mathias Trabandt and Karl Walentin 2013 Dynare Team License: GPL-3+ +Files: matlab/simpsa.m matlab/simpsaget.m matlab/simpsaset.m +Copyright: 2005 Henning Schmidt, FCC, henning@fcc.chalmers.se + 2006 Brecht Donckels, BIOMATH, brecht.donckels@ugent.be + 2013 Dynare Team +License: GPL-3+ + Files: matlab/missing/stats/normpdf.m matlab/missing/stats/gamcdf.m matlab/missing/stats/common_size.m matlab/missing/stats/chi2inv.m matlab/missing/stats/gaminv.m matlab/missing/stats/gampdf.m diff --git a/matlab/simpsaget.m b/matlab/simpsaget.m index b5014ac1c..c69338918 100644 --- a/matlab/simpsaget.m +++ b/matlab/simpsaget.m @@ -16,24 +16,23 @@ function o = SIMPSAGET(options,name,default,flag) % returns val = 1e-4 if the TolX property is not specified in opts. % % See also SIMPSASET, SIMPSA -% + % Copyright (C) 2006 Brecht Donckels, BIOMATH, brecht.donckels@ugent.be % -% This program 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 2 -% of the License, or (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, +% 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. -% +% GNU General Public License for more details. +% % You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, -% USA. - +% along with Dynare. If not, see . % undocumented usage for fast access with no error checking From 51e4e490d61c148f34cdea808f7a36d2c6f7ba57 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Thu, 3 Oct 2013 09:09:18 +0200 Subject: [PATCH 45/80] remove unnecessary test --- preprocessor/Statement.cc | 4 ---- 1 file changed, 4 deletions(-) diff --git a/preprocessor/Statement.cc b/preprocessor/Statement.cc index 242c73a8f..0a449b854 100644 --- a/preprocessor/Statement.cc +++ b/preprocessor/Statement.cc @@ -123,10 +123,6 @@ NativeStatement::computingPass() } lastidx = apostrophes[i]; } - - //Replace last (or only) element - if (apostrophes.empty()) - lastidx = 0; newstr.append(regex_replace(native_statement.substr(lastidx, native_statement.size() - lastidx), date_expr, format)); native_statement = newstr; From b372973ab375e15820c2ac5749a4cf55d86d01f6 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Thu, 3 Oct 2013 10:48:29 +0200 Subject: [PATCH 46/80] catch leading space before comment --- preprocessor/Statement.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/preprocessor/Statement.cc b/preprocessor/Statement.cc index 0a449b854..8ac208fe4 100644 --- a/preprocessor/Statement.cc +++ b/preprocessor/Statement.cc @@ -76,7 +76,7 @@ NativeStatement::computingPass() { using namespace boost::xpressive; // Return if this is a comment - sregex comment_expr = sregex::compile( "\%.*" ); + sregex comment_expr = sregex::compile( "\\s*\%.*" ); match_results results; if (regex_match(native_statement, results, comment_expr)) return; From 62cad6ff44d7c06e0b54f8f3c09b4326127ce101 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Thu, 3 Oct 2013 11:01:11 +0200 Subject: [PATCH 47/80] Cosmetic changes. --- matlab/dynare_estimation_1.m | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index 5e0789491..cbefcd8a1 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -214,8 +214,7 @@ if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && options_.m return end - -%% Estimation of the posterior mode or likelihood mode +% Estimation of the posterior mode or likelihood mode if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation switch options_.mode_compute case 1 @@ -224,17 +223,16 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation elseif ~user_has_matlab_license('optimization_toolbox') error('Option mode_compute=1 requires the Optimization Toolbox') end - - optim_options = optimset('display','iter','LargeScale','off', ... - 'MaxFunEvals',100000,'TolFun',1e-8,'TolX',1e-6); + % Set default optimization options for fmincon. + optim_options = optimset('display','iter', 'LargeScale','off', 'MaxFunEvals',100000, 'TolFun',1e-8, 'TolX',1e-6); if isfield(options_,'optim_opt') eval(['optim_options = optimset(optim_options,' options_.optim_opt ');']); end if options_.analytic_derivation, optim_options = optimset(optim_options,'GradObj','on','TolX',1e-7); end - [xparam1,fval,exitflag,output,lamdba,grad,hessian_fmincon] = ... - fmincon(objective_function,xparam1,[],[],[],[],lb,ub,[],optim_options,dataset_,options_,M_,estim_params_,bayestopt_,oo_); + [xparam1,fval,exitflag,output,lamdba,grad,hessian_fmincon] = ... + fmincon(objective_function,xparam1,[],[],[],[],lb,ub,[],optim_options,dataset_,options_,M_,estim_params_,bayestopt_,oo_); case 2 error('ESTIMATION: mode_compute=2 option (Lester Ingber''s Adaptive Simulated Annealing) is no longer available') case 3 From 49989504ea9e5a1c4fb8263396bad92a2e6ffbde Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Thu, 3 Oct 2013 11:06:07 +0200 Subject: [PATCH 48/80] Cosmetic changes. --- doc/dynare.texi | 2 +- matlab/dynare_estimation_1.m | 3 +-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/doc/dynare.texi b/doc/dynare.texi index d2e3bfcad..06d7c3904 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -4407,7 +4407,7 @@ Specifies the optimizer for the mode computation: @item 0 The mode isn't computed. When @code{mode_file} option is specified, the -mode is simply read from that file. +mode is simply read from that file. When @code{mode_file} option is not specified, Dynare reports the value of the log posterior (log likelihood) diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index cbefcd8a1..c2cfc8dbf 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -241,7 +241,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation elseif ~exist('OCTAVE_VERSION') && ~user_has_matlab_license('optimization_toolbox') error('Option mode_compute=3 requires the Optimization Toolbox') end - + % Set default optimization options for fminunc. optim_options = optimset('display','iter','MaxFunEvals',100000,'TolFun',1e-8,'TolX',1e-6); if isfield(options_,'optim_opt') eval(['optim_options = optimset(optim_options,' options_.optim_opt ');']); @@ -256,7 +256,6 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation func = @(x) objective_function(x, dataset_,options_,M_,estim_params_,bayestopt_,oo_); [xparam1,fval,exitflag] = fminunc(func,xparam1,optim_options); end - case 4 H0 = 1e-4*eye(nx); crit = 1e-7; From b2db159cdd0b79e3a67992250777050639ee1a5c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Thu, 3 Oct 2013 12:35:06 +0200 Subject: [PATCH 49/80] Added the possibility to pass options for csminwell (mode_compute=4) through the optim option (in the estimation command). --- matlab/dynare_estimation_1.m | 39 +++++++++++++++++++++++++++++++----- 1 file changed, 34 insertions(+), 5 deletions(-) diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index c2cfc8dbf..34bd7bd51 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -257,19 +257,48 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation [xparam1,fval,exitflag] = fminunc(func,xparam1,optim_options); end case 4 + % Set default options. H0 = 1e-4*eye(nx); crit = 1e-7; nit = 1000; verbose = 2; - if options_.analytic_derivation, + numgrad = options_.gradient_method; + epsilon = options_.gradient_epsilon; + % Change some options. + if isfield(options_,'optim_opt') + options_list = strsplit(options_.optim_opt,','); + number_of_options = length(options_list)/2; + o = 1; + while o<=number_of_options + switch strtrim(options_list{2*(o-1)+1}) + case '''MaxIter''' + nit = str2num(options_list{2*(o-1)+2}); + case '''InitialHessian''' + H0 = eval(eval(options_list{2*(o-1)+2})); + case '''TolF''' + crit = str2double(options_list{2*(o-1)+2}); + case '''NumgradAlgorithm''' + numgrad = str2num(options_list{2*(o-1)+2}); + case '''NumgradEpsilon''' + epsilon = str2double(options_list{2*(o-1)+2}); + otherwise + error(['csminwel: Unknown option (' options_list{2*(o-1)+1} ')!']) + end + o = o + 1; + end + end + % Set flag for analytical gradient. + if options_.analytic_derivation analytic_grad=1; else analytic_grad=[]; end - - [fval,xparam1,grad,hessian_csminwel,itct,fcount,retcodehat] = ... - csminwel1(objective_function,xparam1,H0,analytic_grad,crit,nit,options_.gradient_method,options_.gradient_epsilon,dataset_,options_,M_,estim_params_,bayestopt_,oo_); - disp(sprintf('Objective function at mode: %f',fval)) + % Call csminwell. + H0 + [fval,xparam1,grad,hessian_csminwel,itct,fcount,retcodehat] = ... + csminwel1(objective_function, xparam1, H0, analytic_grad, crit, nit, numgrad, epsilon, dataset_, options_, M_, estim_params_, bayestopt_, oo_); + % Disp value at the mode. + disp(sprintf('Objective function at mode: %f',fval)) case 5 if isfield(options_,'hess') flag = options_.hess; From 0cac8a2decb93d6b5d66756760db1db9494f7632 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Thu, 3 Oct 2013 15:43:38 +0200 Subject: [PATCH 50/80] clarify error message --- preprocessor/Statement.cc | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/preprocessor/Statement.cc b/preprocessor/Statement.cc index 8ac208fe4..c60be8432 100644 --- a/preprocessor/Statement.cc +++ b/preprocessor/Statement.cc @@ -95,8 +95,7 @@ NativeStatement::computingPass() if (apostrophes.size() % 2) { - cerr << native_statement << - " seems to be invalid Matlab syntax (an odd number of apostrophes was encountered)" << endl; + cerr << "ERROR: A Matlab Statement has an odd number of apostrophes: " << native_statement << endl; exit(EXIT_FAILURE); } From d3111863fe38cb5b88c004ce5d9abf6320349510 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Thu, 3 Oct 2013 16:19:35 +0200 Subject: [PATCH 51/80] fix regular expression --- preprocessor/Statement.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/preprocessor/Statement.cc b/preprocessor/Statement.cc index c60be8432..5e08c164e 100644 --- a/preprocessor/Statement.cc +++ b/preprocessor/Statement.cc @@ -102,7 +102,7 @@ NativeStatement::computingPass() bool skip = false; string newstr = ""; int lastidx = 0; - sregex date_expr = sregex::compile( "-?[0-9]+[Mm]([1-9]|1[0-2])|-?[0-9]+[Qq][1-4]|-?[0-9]+[Ww]([1-9]{1}|[1-4][0-9]|5[0-2])" ); + sregex date_expr = sregex::compile( "-?[0-9]+[Mm](1[0-2]|[1-9])|-?[0-9]+[Qq][1-4]|-?[0-9]+[Ww]([1-4][0-9]|5[0-2]|[1-9])" ); string format( "dynDate('$&')" ); for (size_t i = 0; i < apostrophes.size(); i++) if (apostrophes[i] == 0) From 06ab26f742a21a85fc3e9399509c3bd104e237b2 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Thu, 3 Oct 2013 15:30:22 +0200 Subject: [PATCH 52/80] handle inline comments --- preprocessor/Statement.cc | 39 +++++++++++++++++++++++++++++++++------ 1 file changed, 33 insertions(+), 6 deletions(-) diff --git a/preprocessor/Statement.cc b/preprocessor/Statement.cc index 5e08c164e..649d99a5b 100644 --- a/preprocessor/Statement.cc +++ b/preprocessor/Statement.cc @@ -101,10 +101,11 @@ NativeStatement::computingPass() bool skip = false; string newstr = ""; - int lastidx = 0; sregex date_expr = sregex::compile( "-?[0-9]+[Mm](1[0-2]|[1-9])|-?[0-9]+[Qq][1-4]|-?[0-9]+[Ww]([1-4][0-9]|5[0-2]|[1-9])" ); string format( "dynDate('$&')" ); - for (size_t i = 0; i < apostrophes.size(); i++) + size_t length, i, lastidx, commentidx; + length = lastidx = 0; + for (i = 0; i < apostrophes.size(); i++) if (apostrophes[i] == 0) skip = true; else @@ -116,14 +117,40 @@ NativeStatement::computingPass() } else { + length = apostrophes[i] - lastidx; + commentidx = native_statement.substr(lastidx, length).find("%", 0); + if (commentidx != string::npos) + length = commentidx; + + newstr.append(regex_replace(native_statement.substr(lastidx, length), date_expr, format)); + if (commentidx != string::npos) + { + lastidx += commentidx; + while (i++ < apostrophes.size()) + { + newstr.append(native_statement.substr(lastidx, apostrophes[i] - lastidx)); + lastidx = apostrophes[i]; + } + native_statement = newstr; + return; + } skip = true; - newstr.append(regex_replace(native_statement.substr(lastidx, apostrophes[i] - lastidx), - date_expr, format)); } lastidx = apostrophes[i]; } - newstr.append(regex_replace(native_statement.substr(lastidx, native_statement.size() - lastidx), - date_expr, format)); + length = native_statement.length() - lastidx; + commentidx = native_statement.substr(lastidx, length).find("%", 0); + if (commentidx != string::npos) + length = commentidx; + + newstr.append(regex_replace(native_statement.substr(lastidx, length), date_expr, format)); + + if (commentidx != string::npos) + { + lastidx += commentidx; + newstr.append(native_statement.substr(lastidx, native_statement.length() - lastidx)); + } + native_statement = newstr; } From 68bdf1aff85ce11a525ba26e226bb401943b8997 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Thu, 3 Oct 2013 16:34:51 +0200 Subject: [PATCH 53/80] remove unused portion of code, code cleanup --- preprocessor/Statement.cc | 51 +++++++++++++-------------------------- 1 file changed, 17 insertions(+), 34 deletions(-) diff --git a/preprocessor/Statement.cc b/preprocessor/Statement.cc index 649d99a5b..4017d9959 100644 --- a/preprocessor/Statement.cc +++ b/preprocessor/Statement.cc @@ -103,43 +103,26 @@ NativeStatement::computingPass() string newstr = ""; sregex date_expr = sregex::compile( "-?[0-9]+[Mm](1[0-2]|[1-9])|-?[0-9]+[Qq][1-4]|-?[0-9]+[Ww]([1-4][0-9]|5[0-2]|[1-9])" ); string format( "dynDate('$&')" ); - size_t length, i, lastidx, commentidx; - length = lastidx = 0; - for (i = 0; i < apostrophes.size(); i++) + size_t lastidx = 0; + for (size_t i = 0; i < apostrophes.size(); i++) if (apostrophes[i] == 0) skip = true; else - { - if (skip) - { - skip = false; - newstr.append(native_statement.substr(lastidx, apostrophes[i] - lastidx)); - } - else - { - length = apostrophes[i] - lastidx; - commentidx = native_statement.substr(lastidx, length).find("%", 0); - if (commentidx != string::npos) - length = commentidx; - - newstr.append(regex_replace(native_statement.substr(lastidx, length), date_expr, format)); - if (commentidx != string::npos) - { - lastidx += commentidx; - while (i++ < apostrophes.size()) - { - newstr.append(native_statement.substr(lastidx, apostrophes[i] - lastidx)); - lastidx = apostrophes[i]; - } - native_statement = newstr; - return; - } - skip = true; - } - lastidx = apostrophes[i]; - } - length = native_statement.length() - lastidx; - commentidx = native_statement.substr(lastidx, length).find("%", 0); + if (skip) + { + newstr.append(native_statement.substr(lastidx, apostrophes[i] - lastidx)); + lastidx = apostrophes[i]; + skip = false; + } + else + { + newstr.append(regex_replace(native_statement.substr(lastidx, apostrophes[i] - lastidx), + date_expr, format)); + lastidx = apostrophes[i]; + skip = true; + } + size_t length = native_statement.length() - lastidx; + size_t commentidx = native_statement.substr(lastidx, length).find("%", 0); if (commentidx != string::npos) length = commentidx; From 05946cd684dd31b8ec7876f3c3aeed39fb07287b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Fri, 4 Oct 2013 11:59:35 +0200 Subject: [PATCH 54/80] Changed names of some options for csminwel (mode_compute=4). --- matlab/dynare_estimation_1.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index 34bd7bd51..0a8194130 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -273,9 +273,9 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation switch strtrim(options_list{2*(o-1)+1}) case '''MaxIter''' nit = str2num(options_list{2*(o-1)+2}); - case '''InitialHessian''' + case '''InitialInverseHessian''' H0 = eval(eval(options_list{2*(o-1)+2})); - case '''TolF''' + case '''TolFun''' crit = str2double(options_list{2*(o-1)+2}); case '''NumgradAlgorithm''' numgrad = str2num(options_list{2*(o-1)+2}); From b685fe6197ef5b07ca111e639658a3aa812c01ef Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Fri, 4 Oct 2013 12:01:36 +0200 Subject: [PATCH 55/80] Updated description of the optim option in the estimation command (mode_compute=1,3,7 and 4). --- doc/dynare.texi | 48 +++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 43 insertions(+), 5 deletions(-) diff --git a/doc/dynare.texi b/doc/dynare.texi index 06d7c3904..bea289596 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -4401,6 +4401,7 @@ hessian (@code{hh}, only if @code{cova_compute=1}) in a file called @file{@var{MODEL_FILENAME}_mode.mat} @item mode_compute = @var{INTEGER} | @var{FUNCTION_NAME} +@anchor{mode_compute} Specifies the optimizer for the mode computation: @table @code @@ -4504,11 +4505,45 @@ parameters. Default: @code{1e-32} Metropolis-Hastings simulations instead of starting from scratch. Shouldn't be used together with @code{mh_recover} -@item optim = (@var{fmincon options}) -Can be used to set options for @code{fmincon}, the optimizing function -of MATLAB Optimization toolbox. Use MATLAB's syntax for these -options. Default: -@code{('display','iter','LargeScale','off','MaxFunEvals',100000,'TolFun',1e-8,'TolX',1e-6)} +@item optim = (@var{NAME}, @var{VALUE}, ...) +A list of @var{NAME} and @var{VALUE} pairs. Can be used to set options for the optimization routines. The set of available options depends on the selected optimization routine (ie on the value of option @ref{mode_compute}): + +@table @code + +@item 1, 3, 7 +Available options are given in the documentation of the MATLAB optimization toolbox or in Octave's documentation. + +@item 4 +Available options are: + +@table @code + +@item 'MaxIter' +Maximum number of iterations. Default: @code{1000} + +@item 'NumgradAlgorithm' +Possible values are @code{2}, @code{3} and @code{5} respectively corresponding to the two, three and five points formula used to compute the gradient of the objective function (see @cite{Abramowitz and Stegun (1964)}). Values @code{13} and @code{15} are more experimental. If perturbations on the right and the left increase the value of the objective function (we minimize this function) then we force the corresponding element of the gradient to be zero. The idea is to temporarly reduce the size of the optimization problem. Default: @code{2}. + +@item 'NumgradEpsilon' +Size of the perturbation used to compute numerically the gradient of the objective function. Default: @code{1e-6} + +@item 'TolFun' +Stopping criteria. Default: @code{1e-7} + +@item 'InitialInverseHessian' +Initial approximation for the inverse of the Hessian matrix of the posterior kernel (or likelihood). Obviously this approximation has to be a square, positive definite and symmetric matrix. Default: @code{'1e-4*eye(nx)'}, where @code{nx} is the number of parameters to be estimated. + +@end table + +@item + +@end table + +@customhead{Example 1} +To change the defaults of csminwel (@code{mode_compute=4}): + +@code{estimation(..., mode_compute=4, optim=('NumgradAlgorithm',3,'TolFun',1e-5), ...);} + @item nodiagnostic @anchor{nodiagnostic} Does not compute the convergence diagnostics for @@ -8719,6 +8754,9 @@ internals --test particle/local_state_iteration @itemize +@item +Milton Abramowitz, Irene A. Stegun (1964): ``Handbook of Mathematical Functions'', Courier Dover Publications + @item Aguiar, Mark and Gopinath, Gita (2004): ``Emerging Market Business Cycles: The Cycle is the Trend,'' @i{NBER Working Paper}, 10734 From 36e3fb496c2d926c1307a34aee5728ab1ccc737f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Fri, 4 Oct 2013 12:17:30 +0200 Subject: [PATCH 56/80] Do not crash if an unknown optimization option is declared (replaced an error by a warning). --- matlab/dynare_estimation_1.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index 0a8194130..e50e8e51a 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -282,7 +282,7 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation case '''NumgradEpsilon''' epsilon = str2double(options_list{2*(o-1)+2}); otherwise - error(['csminwel: Unknown option (' options_list{2*(o-1)+1} ')!']) + warning(['csminwel: Unknown option (' options_list{2*(o-1)+1} ')!']) end o = o + 1; end From 51be957fb60eda3db4dcc9a1de8fbdb909c12ec4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Fri, 4 Oct 2013 16:12:14 +0200 Subject: [PATCH 57/80] Changed the organization of the options for gmhmaxlik (mode_compute=6) so that options can be set using the optim option of the estimation command. Added an option (targeted acceptance rate). --- matlab/dynare_estimation_1.m | 89 +++++++++++++++++++++++++--------- matlab/global_initialization.m | 1 + matlab/gmhmaxlik.m | 2 +- 3 files changed, 69 insertions(+), 23 deletions(-) diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index e50e8e51a..49dfecfc8 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -331,41 +331,85 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation parameter_names = bayestopt_.name; save([M_.fname '_mode.mat'],'xparam1','hh','gg','fval','invhess','parameter_names'); case 6 + % Set default options + gmhmaxlikOptions = options_.gmhmaxlik; + if ~isempty(hh); + gmhmaxlikOptions.varinit = 'previous'; + else + gmhmaxlikOptions.varinit = 'prior'; + end + if isfield(options_,'optim_opt') + options_list = strsplit(options_.optim_opt,','); + number_of_options = length(options_list)/2; + o = 1; + while o<=number_of_options + switch strtrim(options_list{2*(o-1)+1}) + case '''NumberOfMh''' + gmhmaxlikOptions.iterations = str2num(options_list{2*(o-1)+2}); + case '''ncov-mh''' + gmhmaxlikOptions.number = str2num(options_list{2*(o-1)+2}); + case '''nscale''' + gmhmaxlikOptions.nscale = str2double(options_list{2*(o-1)+2}); + case '''nclimb''' + gmhmaxlikOptions.nclimb = str2num(options_list{2*(o-1)+2}); + case '''InitialCovarianceMatrix''' + switch eval(options_list{2*(o-1)+2}) + case 'previous' + if isempty(hh) + error('gmhmaxlik: No previous estimate of the Hessian matrix available!') + else + gmhmaxlikOptions.varinit = 'previous' + end + case {'prior', 'identity'} + gmhmaxlikOptions.varinit = eval(options_list{2*(o-1)+2}); + otherwise + error('gmhmaxlik: Unknown value for option ''InitialCovarianceMatrix''!') + end + case '''AcceptanceRateTarget''' + gmhmaxlikOptions.target = str2num(options_list{2*(o-1)+2}); + if gmhmaxlikOptions.target>1 || gmhmaxlikOptions.target1 + if gmhmaxlikOptions.iterations>1 flag = ''; else flag = 'LastCall'; end [xparam1,PostVar,Scale,PostMean] = ... - gmhmaxlik(objective_function,xparam1,[lb ub],options_.gmhmaxlik,Scale,flag,MeanPar,CovJump,dataset_,options_,M_,estim_params_,bayestopt_,oo_); + gmhmaxlik(objective_function,xparam1,[lb ub],gmhmaxlikOptions,Scale,flag,MeanPar,CovJump,dataset_,options_,M_,estim_params_,bayestopt_,oo_); fval = feval(objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_); options_.mh_jscale = Scale; mouvement = max(max(abs(PostVar-OldPostVar))); @@ -378,14 +422,14 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation OldMode = fval; else OldPostVar = PostVar; - if i Date: Fri, 4 Oct 2013 16:12:43 +0200 Subject: [PATCH 58/80] Document the options of gmhmaxlik (mode_compute=6). --- doc/dynare.texi | 27 ++++++++++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/doc/dynare.texi b/doc/dynare.texi index bea289596..2b4fbf05e 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -4535,7 +4535,32 @@ Initial approximation for the inverse of the Hessian matrix of the posterior ker @end table -@item +@item 6 +Available options are: + +@table @code + +@item 'NumberOfMh' +Number of MCMC run sequentially. Default: @code{3} + +@item 'ncov-mh' +Number of iterations used for updating the covariance matrix of the jumping distribution. Default: @code{20000} + +@item 'nscale-mh' +Maximum number of iterations used for adjusting the scale parameter of the jumping distribution. @code{200000} + +@item 'nclimb' +Number of iterations in the last MCMC (climbing mode). + +@item 'InitialCovarianceMatrix' +Initial covariance matrix of the jumping distribution. Default is @code{'previous'} if option @code{mode_file} is used, @code{'prior'} otherwise. + +@item 'AcceptanceRateTarget' +A real number between zero and one. The scale parameter of the jumping distribution is adjusted so that the effective acceptance rate matches the value of option @code{'AcceptanceRateTarget'}. Default: @code{1.0/3.0} + +@end table + + @end table From c737f35ca7d74fc60cb18ee75a7211fcaff46677 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Fri, 4 Oct 2013 16:17:29 +0200 Subject: [PATCH 59/80] remove error message because odd numbers of apostrophes are valid for transposing matrices in Matlab syntax --- preprocessor/Statement.cc | 6 ------ 1 file changed, 6 deletions(-) diff --git a/preprocessor/Statement.cc b/preprocessor/Statement.cc index 4017d9959..38da6e366 100644 --- a/preprocessor/Statement.cc +++ b/preprocessor/Statement.cc @@ -93,12 +93,6 @@ NativeStatement::computingPass() else apostrophes.push_back(idx); - if (apostrophes.size() % 2) - { - cerr << "ERROR: A Matlab Statement has an odd number of apostrophes: " << native_statement << endl; - exit(EXIT_FAILURE); - } - bool skip = false; string newstr = ""; sregex date_expr = sregex::compile( "-?[0-9]+[Mm](1[0-2]|[1-9])|-?[0-9]+[Qq][1-4]|-?[0-9]+[Ww]([1-4][0-9]|5[0-2]|[1-9])" ); From 878ce60f5157bd0c0dcd112a77029db36fdc0a75 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Sun, 6 Oct 2013 20:22:20 +0200 Subject: [PATCH 60/80] Take care of NaN or Inf in residuals or endogenous values Also adds more debugging information. Closes #491 --- matlab/sim1.m | 53 +++++++++++++++++++++++++++++++++------------------ 1 file changed, 34 insertions(+), 19 deletions(-) diff --git a/matlab/sim1.m b/matlab/sim1.m index ffe622eb6..f3ca28205 100644 --- a/matlab/sim1.m +++ b/matlab/sim1.m @@ -70,9 +70,7 @@ i_upd = ny+(1:periods*ny); Y = endo_simul(:); disp (['-----------------------------------------------------']) ; -disp (['MODEL SIMULATION :']) ; -fprintf('\n') ; - +fprintf('MODEL SIMULATION:\n'); model_dynamic = str2func([M_.fname,'_dynamic']); z = Y(find(lead_lag_incidence')); @@ -115,18 +113,15 @@ for iter = 1:options_.maxit_ end err = max(abs(res)); - + if options_.debug + fprintf('\nLargest absolute residual at iteration %d: %10.3f\n',iter,err); + if any(isnan(res)) || any(isinf(res)) || any(isnan(Y)) || any(isinf(Y)) + fprintf('\nWARNING: NaN or Inf detected in the residuals or endogenous variables.\n'); + end + skipline() + end if err < options_.dynatol.f stop = 1 ; - fprintf('\n') ; - disp([' Total time of simulation :' num2str(etime(clock,h1))]) ; - fprintf('\n') ; - disp([' Convergency obtained.']) ; - fprintf('\n') ; - oo_.deterministic_simulation.status = 1;% Convergency obtained. - oo_.deterministic_simulation.error = err; - oo_.deterministic_simulation.iterations = iter; - oo_.endo_simul = reshape(Y,ny,periods+2); break end @@ -137,16 +132,36 @@ for iter = 1:options_.maxit_ end -if ~stop - fprintf('\n') ; - disp([' Total time of simulation :' num2str(etime(clock,h1))]) ; - fprintf('\n') ; - disp(['WARNING : maximum number of iterations is reached (modify options_.maxit_).']) ; - fprintf('\n') ; +if stop + if any(isnan(res)) || any(isinf(res)) || any(isnan(Y)) || any(isinf(Y)) + oo_.deterministic_simulation.status = 0;% NaN or Inf occurred + oo_.deterministic_simulation.error = err; + oo_.deterministic_simulation.iterations = iter; + oo_.endo_simul = reshape(Y,ny,periods+2); + skipline(); + fprintf('\nSimulation terminated after %d iterations.\n',iter); + fprintf('Total time of simulation : %10.3f\n',etime(clock,h1)); + error('Simulation terminated with NaN or Inf in the residuals or endogenous variables. There is most likely something wrong with your model.'); + else + skipline(); + fprintf('\nSimulation concluded successfully after %d iterations.\n',iter); + fprintf('Total time of simulation : %10.3f\n',etime(clock,h1)); + fprintf('Convergency obtained.\n'); + oo_.deterministic_simulation.status = 1;% Convergency obtained. + oo_.deterministic_simulation.error = err; + oo_.deterministic_simulation.iterations = iter; + oo_.endo_simul = reshape(Y,ny,periods+2); + end +elseif ~stop + skipline(); + fprintf('\nSimulation terminated after %d iterations.\n',iter); + fprintf('Total time of simulation : %10.3f\n',etime(clock,h1)); + fprintf('WARNING : maximum number of iterations is reached (modify options_.maxit_).\n') ; oo_.deterministic_simulation.status = 0;% more iterations are needed. oo_.deterministic_simulation.error = err; %oo_.deterministic_simulation.errors = c/abs(err) oo_.deterministic_simulation.iterations = options_.maxit_; end disp (['-----------------------------------------------------']) ; +skipline(); From 5f483c7d2c8a2070fa7bcd28c6ab156bcb814594 Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Mon, 7 Oct 2013 17:16:47 +0200 Subject: [PATCH 61/80] the field 'time' was not updated when adding two time series --- matlab/@dynSeries/plus.m | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/matlab/@dynSeries/plus.m b/matlab/@dynSeries/plus.m index 969374b5c..8df7178c6 100644 --- a/matlab/@dynSeries/plus.m +++ b/matlab/@dynSeries/plus.m @@ -75,7 +75,7 @@ if ~isequal(B.vobs,C.vobs) && ~(isequal(B.vobs,1) || isequal(C.vobs,1)) else if B.vobs>C.vobs idB = 1:B.vobs; - idC = ones(1:B.vobs); + idC = ones(1,B.vobs); elseif B.vobs Date: Tue, 8 Oct 2013 11:02:46 +0200 Subject: [PATCH 62/80] Do not print H0. --- matlab/dynare_estimation_1.m | 1 - 1 file changed, 1 deletion(-) diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index 49dfecfc8..dac8a0f4d 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -294,7 +294,6 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation analytic_grad=[]; end % Call csminwell. - H0 [fval,xparam1,grad,hessian_csminwel,itct,fcount,retcodehat] = ... csminwel1(objective_function, xparam1, H0, analytic_grad, crit, nit, numgrad, epsilon, dataset_, options_, M_, estim_params_, bayestopt_, oo_); % Disp value at the mode. From 842277447cd5b0da67fc788be182bd3bb8bc6ded Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Tue, 8 Oct 2013 11:06:23 +0200 Subject: [PATCH 63/80] Added comment. --- matlab/simplex_optimization_routine.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/simplex_optimization_routine.m b/matlab/simplex_optimization_routine.m index b06135c62..98aa783eb 100644 --- a/matlab/simplex_optimization_routine.m +++ b/matlab/simplex_optimization_routine.m @@ -145,7 +145,7 @@ else end % Set delta parameter. -if isfield(options,'delta_parameter') +if isfield(options,'delta_parameter')% Size of the simplex delta = options.delta_parameter; else delta = 0.05; From f5c9621ca95fbe2e18d0a3c959254471385d0c97 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Tue, 8 Oct 2013 11:11:35 +0200 Subject: [PATCH 64/80] Changed the definition of the maximum number of function evaluations in dynare'es iplementation of simplex algorithm (mode_compute=8). --- matlab/global_initialization.m | 4 +++- matlab/simplex_optimization_routine.m | 11 +++++------ 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index 6eaece758..b9a9bd51b 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -443,7 +443,9 @@ options_.homotopy_steps = 1; options_.homotopy_force_continue = 0; % Simplex optimization routine (variation on Nelder Mead algorithm). -options_.simplex = []; +simplex.maxfcallfactor = 500; +simplex.maxfcall = []; +options_.simplex = simplex; % CMAES optimization routine. cmaes.SaveVariables='off'; diff --git a/matlab/simplex_optimization_routine.m b/matlab/simplex_optimization_routine.m index 98aa783eb..78b3fcfe6 100644 --- a/matlab/simplex_optimization_routine.m +++ b/matlab/simplex_optimization_routine.m @@ -41,6 +41,11 @@ verbose = 2; % Set number of control variables. number_of_variables = length(x); +% get options. +if isempty(simplex.maxfcall) + max_func_calls = simplex.maxfcallfactor*number_of_variables +end + % Set tolerance parameter. if isfield(options,'tolerance') && isfield(options.tolerance,'x') x_tolerance = options.tolerance.x; @@ -61,12 +66,6 @@ if isfield(options,'maxiter') else max_iterations = 1000; end -% Set maximum number of iterations. -if isfield(options,'maxfcall') - max_func_calls = options.maxfcall; -else - max_func_calls = 500*number_of_variables; -end % Set reflection parameter. if isfield(options,'reflection_parameter') From b01aee4daf1c4c34d49e3b5d115403475522c4aa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Tue, 8 Oct 2013 12:53:19 +0200 Subject: [PATCH 65/80] Changed the default value for the maximum number of iterations (defined in the dynare'es implementation of the simplex optimization algorithm). --- matlab/simplex_optimization_routine.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/matlab/simplex_optimization_routine.m b/matlab/simplex_optimization_routine.m index 78b3fcfe6..5016ccb7b 100644 --- a/matlab/simplex_optimization_routine.m +++ b/matlab/simplex_optimization_routine.m @@ -64,7 +64,7 @@ end if isfield(options,'maxiter') max_iterations = options.maxiter; else - max_iterations = 1000; + max_iterations = 5000; end % Set reflection parameter. From 60e1d1b75e51722b7757088ba177739e7bf1c40a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Tue, 8 Oct 2013 12:55:11 +0200 Subject: [PATCH 66/80] Changed the organization of the options for the dynare's implementation of the simplex optimization algorithm. --- matlab/dynare_estimation_1.m | 28 ++++++++++++++++++++++++++-- matlab/global_initialization.m | 3 +++ 2 files changed, 29 insertions(+), 2 deletions(-) diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index dac8a0f4d..5fcab25c4 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -461,7 +461,6 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation elseif ~exist('OCTAVE_VERSION') && ~user_has_matlab_license('optimization_toolbox') error('Option mode_compute=7 requires the Optimization Toolbox') end - optim_options = optimset('display','iter','MaxFunEvals',1000000,'MaxIter',6000,'TolFun',1e-8,'TolX',1e-6); if isfield(options_,'optim_opt') eval(['optim_options = optimset(optim_options,' options_.optim_opt ');']); @@ -469,7 +468,32 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation [xparam1,fval,exitflag] = fminsearch(objective_function,xparam1,optim_options,dataset_,options_,M_,estim_params_,bayestopt_,oo_); case 8 % Dynare implementation of the simplex algorithm. - [xparam1,fval,exitflag] = simplex_optimization_routine(objective_function,xparam1,options_.simplex,dataset_,options_,M_,estim_params_,bayestopt_,oo_); + simplexOptions = options_.simplex; + if isfield(options_,'optim_opt') + options_list = strsplit(options_.optim_opt,','); + number_of_options = length(options_list)/2; + o = 1; + while o<=number_of_options + switch strtrim(options_list{2*(o-1)+1}) + case '''MaxIter''' + simplexOptions.maxiter = str2num(options_list{2*(o-1)+2}); + case '''TolFun''' + simplexOptions.tolerance.f = str2double(options_list{2*(o-1)+2}); + case '''TolX''' + simplexOptions.tolerance.x = str2double(options_list{2*(o-1)+2}); + case '''MaxFunEvals''' + simplexOptions.maxfcall = str2num(options_list{2*(o-1)+2}); + case '''MaxFunEvalFactor''' + simplexOptions.maxfcallfactor = str2num(options_list{2*(o-1)+2}); + case '''InitialSimplexSize''' + simplexOptions.delta_factor = str2double(options_list{2*(o-1)+2}); + otherwise + warning(['simplex: Unknown option (' options_list{2*(o-1)+1} ')!']) + end + o = o + 1; + end + end + [xparam1,fval,exitflag] = simplex_optimization_routine(objective_function,xparam1,simplexOptions,dataset_,options_,M_,estim_params_,bayestopt_,oo_); case 9 H0 = 1e-4*ones(nx,1); warning('off','CMAES:NonfinitenessRange'); diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index b9a9bd51b..a1046205b 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -443,6 +443,9 @@ options_.homotopy_steps = 1; options_.homotopy_force_continue = 0; % Simplex optimization routine (variation on Nelder Mead algorithm). +simplex.tolerance.x = 1e-4; +simplex.tolerance.f = 1e-4; +simplex.maxiter = 5000; simplex.maxfcallfactor = 500; simplex.maxfcall = []; options_.simplex = simplex; From 75fb6afa61ef2e94a917a27ed256b9b36ee418d7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Tue, 8 Oct 2013 15:15:58 +0200 Subject: [PATCH 67/80] Added documentation about the options of the simplex optimization routine (mode_compute=8). --- doc/dynare.texi | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/doc/dynare.texi b/doc/dynare.texi index 2b4fbf05e..775f266bf 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -4560,6 +4560,30 @@ A real number between zero and one. The scale parameter of the jumping distribut @end table +@item 8 +Available options are: + +@table @code + +@item 'MaxIter' +Maximum number of iterations. Default: @code{5000} + +@item 'MaxFunvEvals' +Maximum number of objective function evaluations. No default. + +@item 'MaxFunvEvalFactor' +Set @code{MaxFunvEvals} equal to @code{MaxFunvEvalFactor} times the number of estimated parameters. Default: @code{500}. + +@item 'TolFun' +Tolerance parameter (w.r.t the objective function). Default: @code{1e-4} + +@item 'TolX' +Tolerance parameter (w.r.t the instruments). Default: @code{1e-4} + +@item 'InitialSimplexSize' +Initial size of the simplex, expressed as percentage deviation from the provided initial guess in each direction. Default: @code{.05} + +@end table @end table From 36b69355f25b53f9c7553742493bfaaec2ecd8db Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Tue, 8 Oct 2013 15:17:18 +0200 Subject: [PATCH 68/80] Added a routine to copy options organized in struct into a cell. --- matlab/options2cell.m | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) create mode 100644 matlab/options2cell.m diff --git a/matlab/options2cell.m b/matlab/options2cell.m new file mode 100644 index 000000000..3ff8ef227 --- /dev/null +++ b/matlab/options2cell.m @@ -0,0 +1,36 @@ +function c = options2cell(o) + +% Converts an option structure as a cell of NAME and VALUE pairs. +% +% INPUTS +% o o matlab's structure holding a set of options (each field name is the name of an option and the associated content is the value of the option). +% +% OUTPUTS +% o c matlab's cell row array of the form {NAME1, VALUE1, NAME2, VALUE2, NAME3, VALUE3, ...}. + +% Copyright (C) 2013 Dynare Team. +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +s = fieldnames(o); +c = {}; +j = 1; + +for i=1:length(s) + c(j) = {s{i}}; + c(j+1) = {o.(s{i})}; + j = j+2; +end \ No newline at end of file From 334d9976d62965b736908cf9489c59fcb80594d9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Tue, 8 Oct 2013 15:18:14 +0200 Subject: [PATCH 69/80] Add interface to the main options of the simpsa optimization algorithm. --- matlab/dynare_estimation_1.m | 28 ++++++++++++++++++++++++++-- matlab/global_initialization.m | 15 +++++++++++++++ 2 files changed, 41 insertions(+), 2 deletions(-) diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index 5fcab25c4..544dc0654 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -502,8 +502,32 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation xparam1=BESTEVER.x; disp(sprintf('\n Objective function at mode: %f',fval)) case 10 - options = simpsaset('TOLX', options_.dynatol.x,'TOLFUN', options_.dynatol.f); - [xparam1, fval, exitflag] = simpsa(func2str(objective_function),xparam1,lb,ub,options,dataset_,options_,M_,estim_params_,bayestopt_,oo_); + simpsaOptions = options_.simpsa; + if isfield(options_,'optim_opt') + options_list = strsplit(options_.optim_opt,','); + number_of_options = length(options_list)/2; + o = 1; + while o<=number_of_options + switch strtrim(options_list{2*(o-1)+1}) + case '''MaxIter''' + simpsaOptions.MAX_ITER_TOTAL = str2num(options_list{2*(o-1)+2}); + case '''TolFun''' + simpsaOptions.TOLFUN = str2double(options_list{2*(o-1)+2}); + case '''TolX''' + simpsaOptions.TOLX = str2double(options_list{2*(o-1)+2}); + case '''EndTemparature''' + simpsaOptions.TEMP_END = str2double(options_list{2*(o-1)+2}); + case '''MaxFunEvals''' + simpsaOptions.MAX_FUN_EVALS = str2num(options_list{2*(o-1)+2}); + otherwise + warning(['simpsa: Unknown option (' options_list{2*(o-1)+1} ')!']) + end + o = o + 1; + end + end + simpsaOptionsList = options2cell(simpsaOptions); + simpsaOptions = simpsaset(simpsaOptionsList{:}); + [xparam1, fval, exitflag] = simpsa(func2str(objective_function),xparam1,lb,ub,simpsaOptions,dataset_,options_,M_,estim_params_,bayestopt_,oo_); case 11 options_.cova_compute = 0 ; [xparam1,stdh,lb_95,ub_95,med_param] = online_auxiliary_filter(xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_) ; diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index a1046205b..c456b52e0 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -459,6 +459,21 @@ cmaes.LogModulo='0'; cmaes.LogTime='0'; options_.cmaes = cmaes; +% simpsa optimization routine. +simpsa.TOLFUN = 1e-4; +simpsa.TOLX = 1e-4; +simpsa.TEMP_END = .1; +simpsa.COOL_RATE = 10; +simpsa.INITIAL_ACCEPTANCE_RATIO = .95; +simpsa.MIN_COOLING_FACTOR = .9; +simpsa.MAX_ITER_TEMP_FIRST = 50; +simpsa.MAX_ITER_TEMP_LAST = 2000; +simpsa.MAX_ITER_TEMP = 10; +simpsa.MAX_ITER_TOTAL = 5000; +simpsa.MAX_TIME = 2500; +simpsa.MAX_FUN_EVALS = 20000; +simpsa.DISPLAY = 'iter'; +options_.simpsa = simpsa; % prior analysis options_.prior_mc = 20000; From 5e63aa10da614176a5318add204455ea802c3169 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Tue, 8 Oct 2013 15:28:26 +0200 Subject: [PATCH 70/80] Added documentation about the options of the simpsa algorithm. --- doc/dynare.texi | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/doc/dynare.texi b/doc/dynare.texi index 775f266bf..761118d16 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -4585,6 +4585,27 @@ Initial size of the simplex, expressed as percentage deviation from the provided @end table +@item 10 +Available options are: + +@table @code + +@item 'MaxIter' +Maximum number of iterations. Default: @code{5000} + +@item 'MaxFunvEvals' +Maximum number of objective function evaluations. No default. + +@item 'TolFun' +Tolerance parameter (w.r.t the objective function). Default: @code{1e-4} + +@item 'TolX' +Tolerance parameter (w.r.t the instruments). Default: @code{1e-4} + +@item 'EndTemperature' +Terminal condition w.r.t the temperature. When the temperature reaches @code{EndTemperature}, the temperature is set to zero and the algorithm falls back into a standard simplex algorithm. Default: @code{.1} + +@end table @end table From 0f124042461e8beeff392ba1582bb442226ecd06 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Tue, 8 Oct 2013 15:53:55 +0200 Subject: [PATCH 71/80] Added interface for some options of cmaes. --- matlab/dynare_estimation_1.m | 27 ++++++++++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index 544dc0654..1ff2c6ec4 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -495,10 +495,35 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation end [xparam1,fval,exitflag] = simplex_optimization_routine(objective_function,xparam1,simplexOptions,dataset_,options_,M_,estim_params_,bayestopt_,oo_); case 9 + % Set defaults H0 = 1e-4*ones(nx,1); + cmaesOptions = options_.cmaes; + % Modify defaults + if isfield(options_,'optim_opt') + options_list = strsplit(options_.optim_opt,','); + number_of_options = length(options_list)/2; + o = 1; + while o<=number_of_options + switch strtrim(options_list{2*(o-1)+1}) + case '''MaxIter''' + cmaesOptions.MaxIter = str2num(options_list{2*(o-1)+2}); + case '''TolFun''' + cmaesOptions.TolFun = str2double(options_list{2*(o-1)+2}); + case '''TolX''' + cmaesOptions.TolX = str2double(options_list{2*(o-1)+2}); + case '''MaxFunEvals''' + cmaesOptions.MaxFunEvals = str2num(options_list{2*(o-1)+2}); + case '''H0''' + H0 = eval(eval(options_list{2*(o-1)+2})); + otherwise + warning(['cmaes: Unknown option (' options_list{2*(o-1)+1} ')!']) + end + o = o + 1; + end + end warning('off','CMAES:NonfinitenessRange'); warning('off','CMAES:InitialSigma'); - [x, fval, COUNTEVAL, STOPFLAG, OUT, BESTEVER] = cmaes(func2str(objective_function),xparam1,H0,options_.cmaes,dataset_,options_,M_,estim_params_,bayestopt_,oo_); + [x, fval, COUNTEVAL, STOPFLAG, OUT, BESTEVER] = cmaes(func2str(objective_function),xparam1,H0,cmaesOptions,dataset_,options_,M_,estim_params_,bayestopt_,oo_); xparam1=BESTEVER.x; disp(sprintf('\n Objective function at mode: %f',fval)) case 10 From 8d8407f906bc4520453084c01065aabb5cf21b5a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Tue, 8 Oct 2013 15:56:14 +0200 Subject: [PATCH 72/80] Changed default options for cmaes (reduced tolfun and tolx). --- matlab/global_initialization.m | 2 ++ 1 file changed, 2 insertions(+) diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index c456b52e0..154dee373 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -457,6 +457,8 @@ cmaes.WarnOnEqualFunctionValues='no'; cmaes.DispModulo='10'; cmaes.LogModulo='0'; cmaes.LogTime='0'; +cmaes.TolFun = 1e-7; +cmaes.TolX = 1e-7; options_.cmaes = cmaes; % simpsa optimization routine. From 6a250f894d5ce1d29606a53c3e64d7a51c114207 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Tue, 8 Oct 2013 16:18:00 +0200 Subject: [PATCH 73/80] Removed interface for H0 (cmaes). --- matlab/dynare_estimation_1.m | 2 -- 1 file changed, 2 deletions(-) diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index 1ff2c6ec4..e98cf6f9f 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -513,8 +513,6 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation cmaesOptions.TolX = str2double(options_list{2*(o-1)+2}); case '''MaxFunEvals''' cmaesOptions.MaxFunEvals = str2num(options_list{2*(o-1)+2}); - case '''H0''' - H0 = eval(eval(options_list{2*(o-1)+2})); otherwise warning(['cmaes: Unknown option (' options_list{2*(o-1)+1} ')!']) end From 1f884db8ea0b81cd481de3af18a51035ace753db Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Tue, 8 Oct 2013 16:18:54 +0200 Subject: [PATCH 74/80] If TolX is set to any negative number, let cmaes choose the value of TolX. --- matlab/dynare_estimation_1.m | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index e98cf6f9f..2c8783bfd 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -537,7 +537,12 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation case '''TolFun''' simpsaOptions.TOLFUN = str2double(options_list{2*(o-1)+2}); case '''TolX''' - simpsaOptions.TOLX = str2double(options_list{2*(o-1)+2}); + tolx = str2double(options_list{2*(o-1)+2}); + if tolx<0 + simpsaOptions = rmfield(simpsaOptions,'TOLX'); % Let cmaes choose the default. + else + simpsaOptions.TOLX = tolx; + end case '''EndTemparature''' simpsaOptions.TEMP_END = str2double(options_list{2*(o-1)+2}); case '''MaxFunEvals''' From 0157127d0508d0dd894066042392b325d208719a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Tue, 8 Oct 2013 16:26:49 +0200 Subject: [PATCH 75/80] Document some of the options for cmaes. --- doc/dynare.texi | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/doc/dynare.texi b/doc/dynare.texi index 761118d16..d3275cab3 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -4581,10 +4581,30 @@ Tolerance parameter (w.r.t the objective function). Default: @code{1e-4} Tolerance parameter (w.r.t the instruments). Default: @code{1e-4} @item 'InitialSimplexSize' + Initial size of the simplex, expressed as percentage deviation from the provided initial guess in each direction. Default: @code{.05} @end table +@item 9 +Available options are: + +@table @code + +@item 'MaxIter' +Maximum number of iterations. + +@item 'MaxFunEvals' +Maximum number of objective function evaluations. Default: @code{Inf}. + +@item 'TolFun' +Tolerance parameter (w.r.t the objective function). Default: @code{1e-7} + +@item 'TolX' +Tolerance parameter (w.r.t the instruments). Default: @code{1e-7} + +@end table + @item 10 Available options are: From 137662d5b22c20842f8d616103c82853197536cf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Tue, 8 Oct 2013 16:26:56 +0200 Subject: [PATCH 76/80] Fixed typo. --- doc/dynare.texi | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/dynare.texi b/doc/dynare.texi index d3275cab3..f6ad50fc7 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -4568,7 +4568,7 @@ Available options are: @item 'MaxIter' Maximum number of iterations. Default: @code{5000} -@item 'MaxFunvEvals' +@item 'MaxFunEvals' Maximum number of objective function evaluations. No default. @item 'MaxFunvEvalFactor' From 2967328bb0b4e238fa143b9d7418feb9075e72ef Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Tue, 8 Oct 2013 17:21:07 +0200 Subject: [PATCH 77/80] Changed behaviour of dynSeries objects. Let ts be a dynSeries object. The following syntaxes are equivalent: a = ts.lead a = ts(1) b = ts.lag(2) b = ts(-2) Advantage: If (some of) the variables used in the model block are known as dynSeries in matlab's workspace, then we can create new dynSeries objects with simple copy/pastes of the model's equations (because dynSeries objects understands leads and lags as Dynare's preprocessor). --- matlab/@dynSeries/subsref.m | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/matlab/@dynSeries/subsref.m b/matlab/@dynSeries/subsref.m index e55471b68..0edfca800 100644 --- a/matlab/@dynSeries/subsref.m +++ b/matlab/@dynSeries/subsref.m @@ -174,6 +174,16 @@ switch S(1).type % Do nothing. B = A; end + elseif isscalar(S(1).subs{1}) && isnumeric(S(1).subs{1}) && isint(S(1).subs{1}) + % Input is also interpreted as a backward/forward operator + if S(1).subs{1}>0 + B = feval('lead', A, S(1).subs{1}); + elseif S(1).subs{1}<0 + B = feval('lag', A, -S(1).subs{1}); + else + % Do nothing. + B = A; + end elseif isa(S(1).subs{1},'dynDates') % Extract a subsample using a dynDates object [junk,tdx] = intersect(A.time.time,S(1).subs{1}.time,'rows'); @@ -188,6 +198,8 @@ switch S(1).type B.time = A.time(tdx,:); elseif isvector(S(1).subs{1}) && all(isint(S(1).subs{1})) % Extract a subsample using a vector of integers (observation index). + % Note that this does not work if S(1).subs is an integer scalar... In which case S(1).subs is interpreted as a lead/lag operator (as in the Dynare syntax). + % To extract one observation, a dynDates with one element or a dynDate input must be used. if all(S(1).subs{1}>0) && all(S(1).subs{1}<=A.nobs) if size(A.data,2)>1 S(1).subs = [S(1).subs, ':']; @@ -203,6 +215,18 @@ switch S(1).type else error('dynSeries::subsref: Indices are out of bounds!') end + elseif isa(S(1).subs{1},'dynDate') + % Extract a subsample using a dynDates object + [junk,tdx] = intersect(A.time.time,S(1).subs{1}.time,'rows'); + B = dynSeries(); + B.data = A.data(tdx,:); + B.name = A.name; + B.tex = A.tex; + B.nobs = 1; + B.vobs = A.vobs; + B.freq = A.freq; + B.init = A.time(tdx,:); + B.time = A.time(tdx,:); else error('dynSeries::subsref: I have no idea of what you are trying to do!') end From 90d8efd147c38e09415d09df310c92827005d115 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Tue, 8 Oct 2013 17:21:17 +0200 Subject: [PATCH 78/80] Added unitary test. --- matlab/@dynSeries/subsref.m | 23 ++++++++++++++++++++++- 1 file changed, 22 insertions(+), 1 deletion(-) diff --git a/matlab/@dynSeries/subsref.m b/matlab/@dynSeries/subsref.m index 0edfca800..72ffd4adf 100644 --- a/matlab/@dynSeries/subsref.m +++ b/matlab/@dynSeries/subsref.m @@ -600,4 +600,25 @@ end %$ end %$ %$ T = all(t); -%@eof:13 \ No newline at end of file +%@eof:13 + +%@test:14 +%$ try +%$ data = transpose(0:1:50); +%$ ts = dynSeries(data,'1950Q1'); +%$ a = ts.lag; +%$ b = ts.lead; +%$ c = ts(-1); +%$ d = ts(1); +%$ t(1) = 1; +%$ catch +%$ t(1) = 0; +%$ end +%$ +%$ if t(1)>1 +%$ t(2) = (a==c); +%$ t(3) = (b==d); +%$ end +%$ +%$ T = all(t); +%@eof:14 \ No newline at end of file From c6602191887f8db65588bd3a1dd210468b3087c3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Tue, 8 Oct 2013 17:37:36 +0200 Subject: [PATCH 79/80] Fixes #493. --- matlab/@dynSeries/subsasgn.m | 1 + 1 file changed, 1 insertion(+) diff --git a/matlab/@dynSeries/subsasgn.m b/matlab/@dynSeries/subsasgn.m index a31f75391..31e849c3c 100644 --- a/matlab/@dynSeries/subsasgn.m +++ b/matlab/@dynSeries/subsasgn.m @@ -142,6 +142,7 @@ switch length(S) error('dynSeries::subsasgn: Dimension error! The number of variables on the left and right hand side must match.') end A.data(tdx,:) = B.data(tdy,:); + merge_dynSeries_objects = 0; elseif isnumeric(B) merge_dynSeries_objects = 0; if isequal(length(tdx),rows(B)) From 3dd62b37c6f4329e45835272a3b9cad776ccdbb5 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Fri, 4 Oct 2013 19:39:54 +0200 Subject: [PATCH 80/80] Add more explicit error message if mod-file cannot be located of file in different folder is called --- matlab/dynare.m | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/matlab/dynare.m b/matlab/dynare.m index baa87a788..793a33886 100644 --- a/matlab/dynare.m +++ b/matlab/dynare.m @@ -93,7 +93,13 @@ else end; d = dir(fname); if length(d) == 0 + fprintf('\nThe file %s could not be located in the "Current Folder". Check whether you typed in the correct filename\n',fname) + fprintf('and whether the file is really located in the "Current Folder".\n') error(['DYNARE: can''t open ' fname]) +elseif ~isempty(strfind(fname,'\')) || ~isempty(strfind(fname,'/')) + fprintf('\nIt seems you are trying to call a mod-file not located in the "Current Folder". This is not possible.\n') + fprintf('Please set your "Current Folder" to the folder where the mod-file is located.\n') + error(['DYNARE: can''t open ' fname, '. It seems to be located in a different folder (or has an invalid filename).']) end % pre-dynare-preprocessor-hook