stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'master' into maxit

time-shift
Stéphane Adjemian (Charybdis) 2013-10-09 12:19:17 +02:00
parent d369b8ad0e 40a4678afd
commit 8997ae8a70
35 changed files with 2014 additions and 147 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

277
doc/dynare.texi
View File

@ -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

10
license.txt
View File

@ -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

3
m4/ax_matlab_version.m4
View File

@ -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"
;;

3
matlab/@dynSeries/plus.m
View File

@ -75,7 +75,7 @@ if ~isequal(B.vobs,C.vobs) && ~(isequal(B.vobs,1) || isequal(C.vobs,1))
else
if B.vobs>C.vobs
idB = 1:B.vobs;
idC = ones(1:B.vobs);
idC = ones(1,B.vobs);
elseif B.vobs<C.vobs
idB = ones(1,C.vobs);
idC = 1:C.vobs;
@ -116,6 +116,7 @@ for i=1:A.vobs
A.tex(i) = {['(' B.tex{idB(i)} '+' C.tex{idC(i)} ')']};
end
A.data = bsxfun(@plus,B.data,C.data);
A.time = A.init:A.init+A.nobs;
%@test:1
%$ % Define a datasets.

88
matlab/@dynSeries/subsasgn.m
View File

@ -97,8 +97,25 @@ switch length(S)
end
end
end
case '.' % Single variable selection.
if ~isequal(S(1).subs,B.name)
case '.'
if isequal(S(1).subs,'init') && isa(B,'dynDate')
% Overwrite the init member...
A.init = B;
% ... and update freq and time members.
A.freq = A.init.freq;
A.time = A.init:(A.init+(A.nobs-1));
return
elseif isequal(S(1).subs,'time') && isa(B,'dynDates')
% Overwrite the time member...
A.time = B;
% ... and update the freq and init members.
A.init = B(1);
A.freq = A.init.freq;
return
elseif ismember(S(1).subs,{'freq','nobs','vobs','data','name','tex'})
error(['dynSeries::subsasgn: You cannot overwrite ' S(1).subs ' member!'])
elseif ~isequal(S(1).subs,B.name)
% Single variable selection.
if ~isequal(S(1).subs,B.name{1})
% Rename a variable.
id = strmatch(S(1).subs,A.name);
@ -125,6 +142,7 @@ switch length(S)
error('dynSeries::subsasgn: Dimension error! The number of variables on the left and right hand side must match.')
end
A.data(tdx,:) = B.data(tdy,:);
merge_dynSeries_objects = 0;
elseif isnumeric(B)
merge_dynSeries_objects = 0;
if isequal(length(tdx),rows(B))
@ -691,4 +709,68 @@ end
%$ t(7) = dyn_assert(ts1.data,[[A(1:2,1); ones(5,1); A(8:end,1)], [A(1:2,2); ones(5,1); A(8:end,2)], A(:,3)],1e-15);
%$ end
%$ T = all(t);
%@eof:18
%@eof:18
%@test:19
%$ % Define a datasets.
%$ A = rand(40,3);
%$
%$ % Instantiate two dynSeries object.
%$ ts1 = dynSeries(A,'1950Q1',{'A1';'A2';'A3'},[]);
%$
%$ % Instantiate a dynDate object.
%$ dd = dynDate('1952Q1');
%$
%$ % modify first object.
%$ try
%$ ts1.init = dd;
%$ t(1) = 1;
%$ catch
%$ t(1) = 0;
%$ end
%$
%$ % Instantiate a time series object.
%$ if t(1)
%$ t(2) = dyn_assert(ts1.vobs,3);
%$ t(3) = dyn_assert(ts1.nobs,40);
%$ t(4) = dyn_assert(ts1.name{2},'A2');
%$ t(5) = dyn_assert(ts1.name{1},'A1');
%$ t(6) = dyn_assert(ts1.name{3},'A3');
%$ t(7) = dyn_assert(ts1.data,A,1e-15);
%$ t(8) = dyn_assert(isequal(ts1.init,dd),1);
%$ t(9) = dyn_assert(isequal(ts1.time(1),dd),1);
%$ end
%$ T = all(t);
%@eof:19
%@test:20
%$ % Define a datasets.
%$ A = rand(40,3);
%$
%$ % Instantiate two dynSeries object.
%$ ts1 = dynSeries(A,'1950Q1',{'A1';'A2';'A3'},[]);
%$
%$ % Instantiate a dynDate object.
%$ dd = dynDate('1952Q1');
%$
%$ % modify first object.
%$ try
%$ ts1.time = dd:(dd+(ts1.nobs-1));
%$ t(1) = 1;
%$ catch
%$ t(1) = 0;
%$ end
%$
%$ % Instantiate a time series object.
%$ if t(1)
%$ t(2) = dyn_assert(ts1.vobs,3);
%$ t(3) = dyn_assert(ts1.nobs,40);
%$ t(4) = dyn_assert(ts1.name{2},'A2');
%$ t(5) = dyn_assert(ts1.name{1},'A1');
%$ t(6) = dyn_assert(ts1.name{3},'A3');
%$ t(7) = dyn_assert(ts1.data,A,1e-15);
%$ t(8) = dyn_assert(isequal(ts1.init,dd),1);
%$ t(9) = dyn_assert(isequal(ts1.time(1),dd),1);
%$ end
%$ T = all(t);
%@eof:20

47
matlab/@dynSeries/subsref.m
View File

@ -174,6 +174,16 @@ switch S(1).type
% Do nothing.
B = A;
end
elseif isscalar(S(1).subs{1}) && isnumeric(S(1).subs{1}) && isint(S(1).subs{1})
% Input is also interpreted as a backward/forward operator
if S(1).subs{1}>0
B = feval('lead', A, S(1).subs{1});
elseif S(1).subs{1}<0
B = feval('lag', A, -S(1).subs{1});
else
% Do nothing.
B = A;
end
elseif isa(S(1).subs{1},'dynDates')
% Extract a subsample using a dynDates object
[junk,tdx] = intersect(A.time.time,S(1).subs{1}.time,'rows');
@ -188,6 +198,8 @@ switch S(1).type
B.time = A.time(tdx,:);
elseif isvector(S(1).subs{1}) && all(isint(S(1).subs{1}))
% Extract a subsample using a vector of integers (observation index).
% Note that this does not work if S(1).subs is an integer scalar... In which case S(1).subs is interpreted as a lead/lag operator (as in the Dynare syntax).
% To extract one observation, a dynDates with one element or a dynDate input must be used.
if all(S(1).subs{1}>0) && all(S(1).subs{1}<=A.nobs)
if size(A.data,2)>1
S(1).subs = [S(1).subs, ':'];
@ -203,6 +215,18 @@ switch S(1).type
else
error('dynSeries::subsref: Indices are out of bounds!')
end
elseif isa(S(1).subs{1},'dynDate')
% Extract a subsample using a dynDates object
[junk,tdx] = intersect(A.time.time,S(1).subs{1}.time,'rows');
B = dynSeries();
B.data = A.data(tdx,:);
B.name = A.name;
B.tex = A.tex;
B.nobs = 1;
B.vobs = A.vobs;
B.freq = A.freq;
B.init = A.time(tdx,:);
B.time = A.time(tdx,:);
else
error('dynSeries::subsref: I have no idea of what you are trying to do!')
end
@ -576,4 +600,25 @@ end
%$ end
%$
%$ T = all(t);
%@eof:13
%@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

59
matlab/McMCDiagnostics.m
View File

@ -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

10
matlab/dynare.m
View File

@ -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')

4
matlab/dynare_config.m
View File

@ -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

235
matlab/dynare_estimation_1.m
View File

@ -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.target<eps
error('gmhmaxlik: The value of option AcceptanceRateTarget should be a double between 0 and 1!')
end
otherwise
Warning(['gmhmaxlik: Unknown option (' options_list{2*(o-1)+1} ')!'])
end
o = o + 1;
end
end
% Evaluate the objective function.
fval = feval(objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
OldMode = fval;
if ~exist('MeanPar','var')
MeanPar = xparam1;
end
if ~isempty(hh)
switch gmhmaxlikOptions.varinit
case 'previous'
CovJump = inv(hh);
else% The covariance matrix is initialized with the prior
case 'prior'
% The covariance matrix is initialized with the prior
% covariance (a diagonal matrix) %%Except for infinite variances ;-)
varinit = 'prior';
if strcmpi(varinit,'prior')
stdev = bayestopt_.p2;
indx = find(isinf(stdev));
stdev(indx) = ones(length(indx),1)*sqrt(10);
vars = stdev.^2;
CovJump = diag(vars);
elseif strcmpi(varinit,'eye')
vars = ones(length(bayestopt_.p2),1)*0.1;
CovJump = diag(vars);
else
disp('gmhmaxlik :: Error!')
return
end
stdev = bayestopt_.p2;
indx = find(isinf(stdev));
stdev(indx) = ones(length(indx),1)*sqrt(10);
vars = stdev.^2;
CovJump = diag(vars);
case 'identity'
vars = ones(length(bayestopt_.p2),1)*0.1;
CovJump = diag(vars);
otherwise
error('gmhmaxlik: This is a bug! Please contact the developers.')
end
OldPostVar = CovJump;
Scale = options_.mh_jscale;
for i=1:options_.gmhmaxlik.iterations
for i=1:gmhmaxlikOptions.iterations
if i == 1
if options_.gmhmaxlik.iterations>1
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 i<options_.gmhmaxlik.iterations
if i<gmhmaxlikOptions.iterations
flag = '';
else
flag = 'LastCall';
end
[xparam1,PostVar,Scale,PostMean] = ...
gmhmaxlik(objective_function,xparam1,[lb ub],...
options_.gmhmaxlik,Scale,flag,PostMean,PostVar,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
gmhmaxlikOptions,Scale,flag,PostMean,PostVar,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)));
@ -373,7 +442,8 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
OldMode = fval;
end
hh = inv(PostVar);
save([M_.fname '_mode.mat'],'xparam1','hh');
parameter_names = bayestopt_.name;
save([M_.fname '_mode.mat'],'xparam1','hh','parameter_names');
save([M_.fname '_optimal_mh_scale_parameter.mat'],'Scale');
bayestopt_.jscale = ones(length(xparam1),1)*Scale;
end
@ -391,7 +461,6 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
elseif ~exist('OCTAVE_VERSION') && ~user_has_matlab_license('optimization_toolbox')
error('Option mode_compute=7 requires the Optimization Toolbox')
end
optim_options = optimset('display','iter','MaxFunEvals',1000000,'MaxIter',6000,'TolFun',1e-8,'TolX',1e-6);
if isfield(options_,'optim_opt')
eval(['optim_options = optimset(optim_options,' options_.optim_opt ');']);
@ -399,15 +468,95 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
[xparam1,fval,exitflag] = fminsearch(objective_function,xparam1,optim_options,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
case 8
% Dynare implementation of the simplex algorithm.
[xparam1,fval,exitflag] = simplex_optimization_routine(objective_function,xparam1,options_.simplex,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
simplexOptions = options_.simplex;
if isfield(options_,'optim_opt')
options_list = strsplit(options_.optim_opt,',');
number_of_options = length(options_list)/2;
o = 1;
while o<=number_of_options
switch strtrim(options_list{2*(o-1)+1})
case '''MaxIter'''
simplexOptions.maxiter = str2num(options_list{2*(o-1)+2});
case '''TolFun'''
simplexOptions.tolerance.f = str2double(options_list{2*(o-1)+2});
case '''TolX'''
simplexOptions.tolerance.x = str2double(options_list{2*(o-1)+2});
case '''MaxFunEvals'''
simplexOptions.maxfcall = str2num(options_list{2*(o-1)+2});
case '''MaxFunEvalFactor'''
simplexOptions.maxfcallfactor = str2num(options_list{2*(o-1)+2});
case '''InitialSimplexSize'''
simplexOptions.delta_factor = str2double(options_list{2*(o-1)+2});
otherwise
warning(['simplex: Unknown option (' options_list{2*(o-1)+1} ')!'])
end
o = o + 1;
end
end
[xparam1,fval,exitflag] = simplex_optimization_routine(objective_function,xparam1,simplexOptions,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
case 9
% Set defaults
H0 = 1e-4*ones(nx,1);
cmaesOptions = options_.cmaes;
% Modify defaults
if isfield(options_,'optim_opt')
options_list = strsplit(options_.optim_opt,',');
number_of_options = length(options_list)/2;
o = 1;
while o<=number_of_options
switch strtrim(options_list{2*(o-1)+1})
case '''MaxIter'''
cmaesOptions.MaxIter = str2num(options_list{2*(o-1)+2});
case '''TolFun'''
cmaesOptions.TolFun = str2double(options_list{2*(o-1)+2});
case '''TolX'''
cmaesOptions.TolX = str2double(options_list{2*(o-1)+2});
case '''MaxFunEvals'''
cmaesOptions.MaxFunEvals = str2num(options_list{2*(o-1)+2});
otherwise
warning(['cmaes: Unknown option (' options_list{2*(o-1)+1} ')!'])
end
o = o + 1;
end
end
warning('off','CMAES:NonfinitenessRange');
warning('off','CMAES:InitialSigma');
[x, fval, COUNTEVAL, STOPFLAG, OUT, BESTEVER] = cmaes(func2str(objective_function),xparam1,H0,options_.cmaes,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
[x, fval, COUNTEVAL, STOPFLAG, OUT, BESTEVER] = cmaes(func2str(objective_function),xparam1,H0,cmaesOptions,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
xparam1=BESTEVER.x;
disp(sprintf('\n Objective function at mode: %f',fval))
case 10
case 10
simpsaOptions = options_.simpsa;
if isfield(options_,'optim_opt')
options_list = strsplit(options_.optim_opt,',');
number_of_options = length(options_list)/2;
o = 1;
while o<=number_of_options
switch strtrim(options_list{2*(o-1)+1})
case '''MaxIter'''
simpsaOptions.MAX_ITER_TOTAL = str2num(options_list{2*(o-1)+2});
case '''TolFun'''
simpsaOptions.TOLFUN = str2double(options_list{2*(o-1)+2});
case '''TolX'''
tolx = str2double(options_list{2*(o-1)+2});
if tolx<0
simpsaOptions = rmfield(simpsaOptions,'TOLX'); % Let cmaes choose the default.
else
simpsaOptions.TOLX = tolx;
end
case '''EndTemparature'''
simpsaOptions.TEMP_END = str2double(options_list{2*(o-1)+2});
case '''MaxFunEvals'''
simpsaOptions.MAX_FUN_EVALS = str2num(options_list{2*(o-1)+2});
otherwise
warning(['simpsa: Unknown option (' options_list{2*(o-1)+1} ')!'])
end
o = o + 1;
end
end
simpsaOptionsList = options2cell(simpsaOptions);
simpsaOptions = simpsaset(simpsaOptionsList{:});
[xparam1, fval, exitflag] = simpsa(func2str(objective_function),xparam1,lb,ub,simpsaOptions,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
case 11
options_.cova_compute = 0 ;
[xparam1,stdh,lb_95,ub_95,med_param] = online_auxiliary_filter(xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_) ;
case 101
@ -502,7 +651,7 @@ if ~options_.mh_posterior_mode_estimation && options_.cova_compute
params_at_bound=find(xparam1==ub | xparam1==lb);
if ~isempty(params_at_bound)
for ii=1:length(params_at_bound)
params_at_bound_name{ii,1}=get_the_name(ii,0,M_,estim_params_,options_);
params_at_bound_name{ii,1}=get_the_name(params_at_bound(ii),0,M_,estim_params_,options_);
end
disp_string=[params_at_bound_name{1,:}];
for ii=2:size(params_at_bound_name,1)
@ -611,8 +760,8 @@ if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ...
CutSample(M_, options_, estim_params_);
return
else
if ~options_.nodiagnostic && options_.mh_replic > 2000 && options_.mh_nblck > 1
McMCDiagnostics(options_, estim_params_, M_);
if ~options_.nodiagnostic && options_.mh_replic > 2000
oo_= McMCDiagnostics(options_, estim_params_, M_,oo_);
end
%% Here i discard first half of the draws:
CutSample(M_, options_, estim_params_);

14
matlab/dynare_estimation_init.m
View File

@ -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;

74
matlab/geweke_chi2_test.m Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
%
% ------------------------------------------------
% 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;

109
matlab/geweke_moments.m Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
% 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

29
matlab/global_initialization.m
View File

@ -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()

2
matlab/gmhmaxlik.m
View File

@ -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;

37
matlab/missing/stats/chi2cdf.m Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
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

36
matlab/options2cell.m Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
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

3
matlab/reports/@report/write.m
View File

@ -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');

3
matlab/reports/@series/series.m
View File

@ -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

50
matlab/reports/@series/write.m
View File

@ -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

6
matlab/reports/@table/table.m
View File

@ -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)), ...

47
matlab/reports/@table/write.m
View File

@ -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

53
matlab/sim1.m
View File

@ -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();

15
matlab/simplex_optimization_routine.m
View File

@ -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;

530
matlab/simpsa.m Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
% 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

124
matlab/simpsaget.m Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
% 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

172
matlab/simpsaset.m Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
% 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

6
matlab/stoch_simul.m
View File

@ -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

14
preprocessor/ComputingTasks.cc
View File

@ -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);

8
preprocessor/DynareBison.yy
View File

@ -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 <string_val> 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"); };

4
preprocessor/DynareFlex.ll
View File

@ -222,6 +222,8 @@ string eofbuff;
<DYNARE_STATEMENT>presample {return token::PRESAMPLE;}
<DYNARE_STATEMENT>lik_algo {return token::LIK_ALGO;}
<DYNARE_STATEMENT>lik_init {return token::LIK_INIT;}
<DYNARE_STATEMENT>taper_steps {return token::TAPER_STEPS;}
<DYNARE_STATEMENT>geweke_interval {return token::GEWEKE_INTERVAL;}
<DYNARE_STATEMENT>graph {return token::GRAPH;}
<DYNARE_STATEMENT>nograph {return token::NOGRAPH;}
<DYNARE_STATEMENT>nodisplay {return token::NODISPLAY;}
@ -688,7 +690,7 @@ string eofbuff;
return token::INT_NUMBER;
}
<DYNARE_STATEMENT,DYNARE_BLOCK>([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]))) {
<DYNARE_STATEMENT,DYNARE_BLOCK>-?[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;
}

3
preprocessor/ParsingDriver.cc
View File

@ -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

68
preprocessor/Statement.cc
View File

@ -18,6 +18,7 @@
*/
#include "Statement.hh"
#include <boost/xpressive/xpressive.hpp>
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<string::const_iterator> results;
if (regex_match(native_statement, results, comment_expr))
return;
// Otherwise, look at the line and consider substituting date
size_t idx = -1;
vector<size_t> 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
{

6
preprocessor/Statement.hh
View File

@ -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

12
windows/dynare.nsi
View File

@ -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