diff --git a/doc/dynare.texi b/doc/dynare.texi index 2462b1d78..86c058224 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: @@ -4212,12 +4218,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 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}. +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 @@ -4358,8 +4371,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 @@ -4388,13 +4401,14 @@ 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 @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) @@ -4444,6 +4458,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 @@ -4487,11 +4505,135 @@ 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 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 + +@item 8 +Available options are: + +@table @code + +@item 'MaxIter' +Maximum number of iterations. Default: @code{5000} + +@item 'MaxFunEvals' +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 + +@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: + +@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 + +@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 @@ -4758,6 +4900,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} @@ -5048,6 +5202,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{}; @@ -8352,11 +8556,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, 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}. @@ -8381,13 +8590,17 @@ 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 @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, tableSubSectionHeader, zerotol Adds a @code{Series} to a @code{Graph} or a @code{Table}. @optionshead @table @code @@ -8416,6 +8629,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 @@ -8441,6 +8658,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 @@ -8622,6 +8844,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 @@ -8644,6 +8869,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'' @@ -8688,6 +8917,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 diff --git a/license.txt b/license.txt index 2dcfdcf07..c9be4ebec 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 @@ -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/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/@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.vobs0 + 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 @@ -576,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 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.m b/matlab/dynare.m index 06d81f699..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 @@ -108,6 +114,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') 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/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index 7e5071ded..2c8783bfd 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 @@ -243,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 ');']); @@ -258,21 +256,48 @@ 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 + % 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 '''InitialInverseHessian''' + H0 = eval(eval(options_list{2*(o-1)+2})); + case '''TolFun''' + 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 + warning(['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. + [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; @@ -305,41 +330,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))); @@ -352,14 +421,14 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation OldMode = fval; else OldPostVar = PostVar; - if i0 ) && 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/dynare_estimation_init.m b/matlab/dynare_estimation_init.m index 793b96c8e..1eced4f1e 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) @@ -253,6 +253,18 @@ if ~isempty(estim_params_) end error(['Initial value(s) of ', disp_string ,' are outside parameter bounds. Potentially, you should set prior_trunc=0. If you used the mode_file-option, check whether your mode-file is consistent with the priors.']) end + inadmissible_inverse_gamma_values=find(bayestopt_.pshape==4 & xparam1 == 0); + if ~isempty(inadmissible_inverse_gamma_values) + for ii=1:length(inadmissible_inverse_gamma_values) + inadmissible_inverse_gamma_par_names{ii,1}=get_the_name(inadmissible_inverse_gamma_values(ii),0,M_,estim_params_,options_); + end + disp_string=[inadmissible_inverse_gamma_par_names{1,:}]; + for ii=2:size(inadmissible_inverse_gamma_par_names,1) + disp_string=[disp_string,', ',inadmissible_inverse_gamma_par_names{ii,:}]; + end + error(['Initial value(s) of ', disp_string ,' is zero. This is not allowed when using an inverse gamma prior.\n']) + end + lb = bounds(:,1); ub = bounds(:,2); bayestopt_.lb = lb; 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 8154c5bac..31d925dcd 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -108,6 +108,7 @@ gmhmaxlik.iterations = 3; gmhmaxlik.number = 20000; gmhmaxlik.nclimb = 200000; gmhmaxlik.nscale = 200000; +gmhmaxlik.target = 1/3; % Target for the acceptance rate. options_.gmhmaxlik = gmhmaxlik; % Request user input. @@ -445,7 +446,12 @@ options_.homotopy_steps = 1; options_.homotopy_force_continue = 0; % Simplex optimization routine (variation on Nelder Mead algorithm). -options_.simplex = []; +simplex.tolerance.x = 1e-4; +simplex.tolerance.f = 1e-4; +simplex.maxiter = 5000; +simplex.maxfcallfactor = 500; +simplex.maxfcall = []; +options_.simplex = simplex; % CMAES optimization routine. cmaes.SaveVariables='off'; @@ -454,8 +460,25 @@ 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. +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; @@ -580,6 +603,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() diff --git a/matlab/gmhmaxlik.m b/matlab/gmhmaxlik.m index b57f0593b..79c76e087 100644 --- a/matlab/gmhmaxlik.m +++ b/matlab/gmhmaxlik.m @@ -84,7 +84,7 @@ npar = length(xparam1); NumberOfIterations = options.number; MaxNumberOfTuningSimulations = options.nscale; MaxNumberOfClimbingSimulations = options.nclimb; -AcceptanceTarget = 1/3; +AcceptanceTarget = options.target; CovJump = VarCov; ModePar = xparam1; 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 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 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..ce30bf417 100644 --- a/matlab/reports/@series/series.m +++ b/matlab/reports/@series/series.m @@ -49,8 +49,11 @@ o.tableNegColor = 'red'; o.tablePosColor = 'blue'; o.tableMarkerLimit = 1e-4; +o.tableSubSectionHeader = ''; o.tableAlignRight = false; +o.tableRowColor = 'white'; + o.zerotol = 1e-6; if nargin == 1 diff --git a/matlab/reports/@series/write.m b/matlab/reports/@series/write.m index e0f0ffcae..eb2b6f277 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 @@ -37,11 +39,15 @@ 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'); +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'); @@ -51,6 +57,17 @@ 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 ~isempty(o.tableSubSectionHeader) + fprintf(fid, '%s', o.tableSubSectionHeader); + for i=1:size(dates)+length(yrsForAvgs) + fprintf(fid, '&'); + end + fprintf(fid, '\\\\%%\n'); + return; +end if o.tableAlignRight fprintf(fid, '\\multicolumn{1}{r}{'); end @@ -61,7 +78,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); @@ -77,5 +94,26 @@ for i=1:size(data,1) fprintf(fid, ']'); end end -fprintf(fid, ' \\\\\n\n'); + +% 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 55cdde174..eb742409c 100644 --- a/matlab/reports/@table/table.m +++ b/matlab/reports/@table/table.m @@ -39,6 +39,9 @@ o.titleSize = 'large'; o.showHlines = false; o.showVlines = false; o.vlineAfter = ''; +o.vlineAfterEndOfPeriod = false; + +o.annualAverages = false; o.data = ''; o.seriesToUse = ''; @@ -72,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'); @@ -87,6 +91,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..9e881778d 100644 --- a/matlab/reports/@table/write.m +++ b/matlab/reports/@table/write.m @@ -55,13 +55,34 @@ 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 +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'); @@ -72,10 +93,7 @@ end fprintf(fid, '\\toprule%%\n'); % Column Headers -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) @@ -94,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 @@ -114,17 +131,17 @@ switch dates.freq otherwise error('@table.write: invalid dynSeries frequency'); end -fprintf(fid, '\\\\[-10pt]%%\n'); -for i=1:ndates - fprintf(fid, ' & \\hrulefill'); +for i=1:length(yrsForAvgs) + fprintf(fid, ' & %d', years(i)); end -fprintf(fid, '\\\\%%\n'); +fprintf(fid, '\\\\[-2pt]%%\n'); +fprintf(fid, '\\hline%%\n'); 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 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(); diff --git a/matlab/simplex_optimization_routine.m b/matlab/simplex_optimization_routine.m index b06135c62..5016ccb7b 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; @@ -59,13 +64,7 @@ end if isfield(options,'maxiter') max_iterations = 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; + max_iterations = 5000; end % Set reflection parameter. @@ -145,7 +144,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; diff --git a/matlab/simpsa.m b/matlab/simpsa.m new file mode 100644 index 000000000..260607a24 --- /dev/null +++ b/matlab/simpsa.m @@ -0,0 +1,530 @@ +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); + +% 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',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',20000,... % maximum number of function evaluations + 'TOLX',1e-6,... % maximum difference between best and worst function evaluation in simplex + '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 + +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,transpose(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,transpose(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 = transpose(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..c69338918 --- /dev/null +++ b/matlab/simpsaget.m @@ -0,0 +1,124 @@ +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 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 . + + +% 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 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 diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc index 129bfff9e..458a1db26 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; @@ -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); diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy index 621a092b5..856632730 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 %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 @@ -2389,6 +2391,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 b045719c5..19829085d 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;} @@ -688,7 +690,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]))) { +-?[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; } 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..38da6e366 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,76 @@ 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( "\\s*\%.*" ); + 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); + + 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])" ); + string format( "dynDate('$&')" ); + size_t lastidx = 0; + for (size_t i = 0; i < apostrophes.size(); i++) + if (apostrophes[i] == 0) + skip = true; + else + 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; + + 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; +} + +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 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