From e71e89bb38270599a3e5749f3cc7efd336eeb125 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Thu, 26 Jan 2017 14:39:15 +0100 Subject: [PATCH 01/17] Added check on Octave version. --- matlab/dynare_config.m | 9 +++++++++ matlab/supported_octave_version.m | 20 ++++++++++++++++++++ 2 files changed, 29 insertions(+) create mode 100644 matlab/supported_octave_version.m diff --git a/matlab/dynare_config.m b/matlab/dynare_config.m index fc89f1539..c70e0a89f 100644 --- a/matlab/dynare_config.m +++ b/matlab/dynare_config.m @@ -91,6 +91,15 @@ if isoctave p{end+1} = '/missing/ordeig'; end +if isoctave && ~isequal(supported_octave_version(), version()) + skipline() + warning(['This version of Octave is not supported. Consider installing ' ... + 'version %s of Octave,\notherwise m files will be used instead ' ... + 'of precompiled mex files and some features, like solution\n' ... + 'of models approximated at third order, will not be available.'], supported_octave_version()) + skipline() +end + % ilu is missing in Octave < 4.0 if isoctave && octave_ver_less_than('4.0') p{end+1} = '/missing/ilu'; diff --git a/matlab/supported_octave_version.m b/matlab/supported_octave_version.m new file mode 100644 index 000000000..9f6f2d633 --- /dev/null +++ b/matlab/supported_octave_version.m @@ -0,0 +1,20 @@ +function v = supported_octave_version() + +% Copyright (C) 2017 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% Dynare is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see . + +v = '4.2.0'; \ No newline at end of file From 8d4abe4a9f45d50c8ee2bf0666a0944423b9fef9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Wed, 8 Feb 2017 12:24:13 +0100 Subject: [PATCH 02/17] Cosmetic changes. --- matlab/hessian.m | 77 ++++++++++++++++++++++++++---------------------- 1 file changed, 42 insertions(+), 35 deletions(-) diff --git a/matlab/hessian.m b/matlab/hessian.m index 13ddc2f44..4faa4c500 100644 --- a/matlab/hessian.m +++ b/matlab/hessian.m @@ -1,4 +1,4 @@ -function hessian_mat = hessian(func,x,gstep,varargin) % --*-- Unitary tests --*-- +function hessian_mat = hessian(func,x, gstep, varargin) % --*-- Unitary tests --*-- % Computes second order partial derivatives % @@ -26,7 +26,7 @@ function hessian_mat = hessian(func,x,gstep,varargin) % --*-- Unitary tests --*- % none % -% Copyright (C) 2001-2014 Dynare Team +% Copyright (C) 2001-2017 Dynare Team % % This file is part of Dynare. % @@ -46,48 +46,55 @@ function hessian_mat = hessian(func,x,gstep,varargin) % --*-- Unitary tests --*- if ~isa(func, 'function_handle') func = str2func(func); end -n=size(x,1); -h1=max(abs(x),sqrt(gstep(1))*ones(n,1))*eps^(1/6)*gstep(2); -h_1=h1; -xh1=x+h1; -h1=xh1-x; -xh1=x-h_1; -h_1=x-xh1; -xh1=x; -f0=feval(func,x,varargin{:}); -f1=zeros(size(f0,1),n); -f_1=f1; + +n = size(x,1); +h1 = max(abs(x), sqrt(gstep(1))*ones(n, 1))*eps^(1/6)*gstep(2); +h_1 = h1; +xh1 = x+h1; +h1 = xh1-x; +xh1 = x-h_1; +h_1 = x-xh1; +xh1 = x; +f0 = feval(func, x, varargin{:}); +f1 = zeros(size(f0, 1), n); +f_1 = f1; + for i=1:n %do step up - xh1(i)=x(i)+h1(i); - f1(:,i)=feval(func,xh1,varargin{:}); - %do step up - xh1(i)=x(i)-h_1(i); - f_1(:,i)=feval(func,xh1,varargin{:}); - xh1(i)=x(i);%reset parameter + xh1(i) = x(i)+h1(i); + f1(:,i) = feval(func, xh1, varargin{:}); + %do step down + xh1(i) = x(i)-h_1(i); + f_1(:,i) = feval(func, xh1, varargin{:}); + %reset parameter + xh1(i) = x(i); end -xh_1=xh1; -hessian_mat = zeros(size(f0,1),n*n); -temp=f1+f_1-f0*ones(1,n); %term f_(1,0)+f_(-1,0)-f_(0,0) used later + +xh_1 = xh1; +temp = f1+f_1-f0*ones(1, n); %term f_(1,0)+f_(-1,0)-f_(0,0) used later + +hessian_mat = zeros(size(f0,1), n*n); + for i=1:n - if i > 1 %fill symmetric part of Hessian based on previously computed results - k=[i:n:n*(i-1)]; - hessian_mat(:,(i-1)*n+1:(i-1)*n+i-1)=hessian_mat(:,k); + if i > 1 + %fill symmetric part of Hessian based on previously computed results + k = [i:n:n*(i-1)]; + hessian_mat(:,(i-1)*n+1:(i-1)*n+i-1) = hessian_mat(:,k); end - hessian_mat(:,(i-1)*n+i)=(f1(:,i)+f_1(:,i)-2*f0)./(h1(i)*h_1(i)); %formula 25.3.23 + hessian_mat(:,(i-1)*n+i) = (f1(:,i)+f_1(:,i)-2*f0)./(h1(i)*h_1(i)); %formula 25.3.23 for j=i+1:n %step in up direction - xh1(i)=x(i)+h1(i); - xh1(j)=x(j)+h_1(j); + xh1(i) = x(i)+h1(i); + xh1(j) = x(j)+h_1(j); %step in down direction - xh_1(i)=x(i)-h1(i); - xh_1(j)=x(j)-h_1(j); - hessian_mat(:,(i-1)*n+j)=-(-feval(func,xh1,varargin{:})-feval(func,xh_1,varargin{:})+temp(:,i)+temp(:,j))./(2*h1(i)*h_1(j)); %formula 25.3.27 + xh_1(i) = x(i)-h1(i); + xh_1(j) = x(j)-h_1(j); + hessian_mat(:,(i-1)*n+j) =-(-feval(func, xh1, varargin{:})-feval(func, xh_1, varargin{:})+temp(:,i)+temp(:,j))./(2*h1(i)*h_1(j)); %formula 25.3.27 %reset grid points - xh1(i)=x(i); - xh1(j)=x(j); - xh_1(i)=x(i); - xh_1(j)=x(j); + xh1(i) = x(i); + xh1(j) = x(j); + xh_1(i) = x(i); + xh_1(j) = x(j); end end From 35078e13ee885f3efbf88ff2f1c422f061b140d7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Wed, 8 Feb 2017 12:25:06 +0100 Subject: [PATCH 03/17] Completed header + cosmetic changes. --- .../optimization/penalty_objective_function.m | 25 +++++++++++++------ 1 file changed, 18 insertions(+), 7 deletions(-) diff --git a/matlab/optimization/penalty_objective_function.m b/matlab/optimization/penalty_objective_function.m index 639a248aa..c245800c2 100644 --- a/matlab/optimization/penalty_objective_function.m +++ b/matlab/optimization/penalty_objective_function.m @@ -1,7 +1,19 @@ -function [fval,exit_flag,arg1,arg2] = penalty_objective_function(x0,fcn,penalty,varargin) -%function [fval,exit_flag,arg1,arg2] = penalty_objective_function(x0,fcn,penalty,varargin) +function [fval, exit_flag, arg1, arg2] = penalty_objective_function(x, fcn, base_penalty, varargin) -% Copyright (C) 2016 Dynare Team +% Encapsulates an objective function to be minimized, adding a penalty if necessary. +% +% INPUTS +% - x [double] n*1 vector of instrument values. +% - fcn [fhandle] objective function. +% - base_penalty [double] scalar, base of the penality (typically the value of the objective at the previous iteration). +% - varagin [cell] additional parameters for fcn. +% +% OUTPUTS +% - fval [double] scalar, value of the objective function at x. +% - exit_flag [integer] scalar, flag returned by fcn (third output). +% - arg1, arg2 fourth and fifth output arguments of the objective function. + +% Copyright (C) 2016-2017 Dynare Team % % This file is part of Dynare. % @@ -18,9 +30,8 @@ function [fval,exit_flag,arg1,arg2] = penalty_objective_function(x0,fcn,penalty, % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . - [fval,info,exit_flag,arg1,arg2] = fcn(x0,varargin{:}); +[fval, info, exit_flag, arg1, arg2] = fcn(x, varargin{:}); - if info(1) ~= 0 - fval = penalty + info(4); - end +if info(1) + fval = base_penalty + info(4); end \ No newline at end of file From 5b73d7d59c6417dab93e8dd34b89d9aa55cf5471 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Wed, 8 Feb 2017 12:35:08 +0100 Subject: [PATCH 04/17] Added new optimization routine (particleswarm). Only available under Matlab if the Global Optimization Toolbox is installed. --- doc/dynare.texi | 6 ++- matlab/dsge_likelihood.m | 5 +- matlab/dsge_var_likelihood.m | 5 +- matlab/global_initialization.m | 17 ++++++ matlab/non_linear_dsge_likelihood.m | 8 +-- .../optimization/dynare_minimize_objective.m | 54 ++++++++++++++++++- tests/Makefile.am | 1 + tests/optimizers/fs2000_12.mod | 3 ++ 8 files changed, 90 insertions(+), 9 deletions(-) create mode 100644 tests/optimizers/fs2000_12.mod diff --git a/doc/dynare.texi b/doc/dynare.texi index 1ca3bc42e..bf63caf70 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -5410,6 +5410,10 @@ Uses the simpsa algorithm, based on the combination of the non-linear simplex an @item 11 This is not strictly speaking an optimization algorithm. The (estimated) parameters are treated as state variables and estimated jointly with the original state variables of the model using a nonlinear filter. The algorithm implemented in Dynare is described in @cite{Liu and West (2001)}. +@item 12 +Uses @code{particleswarm} optimization routine (available under MATLAB if +the Global Optimization Toolbox is installed; not available under Octave). + @item 101 Uses the SolveOpt algorithm for local nonlinear optimization problems proposed by @cite{Kuntsevich and Kappel (1997)}. @@ -5513,7 +5517,7 @@ A list of @var{NAME} and @var{VALUE} pairs. Can be used to set options for the o @table @code -@item 1, 3, 7 +@item 1, 3, 7, 12 Available options are given in the documentation of the MATLAB Optimization Toolbox or in Octave's documentation. @item 2 diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m index 8a4a7d6b6..d911adc6e 100644 --- a/matlab/dsge_likelihood.m +++ b/matlab/dsge_likelihood.m @@ -141,7 +141,10 @@ trend_coeff = []; exit_flag = 1; info = zeros(4,1); DLIK = []; -Hess = []; +Hess = []; + +% Ensure that xparam1 is a column vector. +xparam1 = xparam1(:); if DynareOptions.estimation_dll [fval,exit_flag,SteadyState,trend_coeff,info,params,H,Q] ... diff --git a/matlab/dsge_var_likelihood.m b/matlab/dsge_var_likelihood.m index 646120c5f..b97beef26 100644 --- a/matlab/dsge_var_likelihood.m +++ b/matlab/dsge_var_likelihood.m @@ -37,7 +37,7 @@ function [fval,info,exit_flag,grad,hess,SteadyState,trend_coeff,PHI_tilde,SIGMA_ % SPECIAL REQUIREMENTS % None. -% Copyright (C) 2006-2016 Dynare Team +% Copyright (C) 2006-2017 Dynare Team % % This file is part of Dynare. % @@ -68,6 +68,9 @@ iXX = []; prior = []; trend_coeff=[]; +% Ensure that xparam1 is a column vector. +xparam1 = xparam1(:); + % Initialization of of the index for parameter dsge_prior_weight in Model.params. if isempty(dsge_prior_weight_idx) dsge_prior_weight_idx = strmatch('dsge_prior_weight',Model.param_names); diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index e65ca4a62..aa1531bb4 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -645,6 +645,23 @@ options_.saopt.nt=10; options_.saopt.step_length_c=0.1; options_.saopt.initial_step_length=1; +% particleswarm (global optimization toolbox needed) +particleswarm.Display = 'iter'; +particleswarm.DisplayInterval = 1; +particleswarm.FunctionTolerance = 1e-6; +particleswarm.FunValCheck = 'on'; +particleswarm.HybridFcn = []; +particleswarm.InertiaRange = [0.1, 1.1]; +particleswarm.MaxIterations = 100000; +particleswarm.MaxStallIterations = 20; +particleswarm.MaxStallTime = Inf; +particleswarm.MaxTime = Inf; +particleswarm.MinNeighborsFraction = .25; +particleswarm.ObjectiveLimit = -Inf; +particleswarm.UseParallel = false; +particleswarm.UseVectorized = false; +options_.particleswarm = particleswarm; + % prior analysis options_.prior_mc = 20000; options_.prior_analysis_endo_var_list = []; diff --git a/matlab/non_linear_dsge_likelihood.m b/matlab/non_linear_dsge_likelihood.m index f8ac92e23..17d448927 100644 --- a/matlab/non_linear_dsge_likelihood.m +++ b/matlab/non_linear_dsge_likelihood.m @@ -107,7 +107,7 @@ function [fval,info,exit_flag,DLIK,Hess,ys,trend_coeff,Model,DynareOptions,Bayes %! @end deftypefn %@eod: -% Copyright (C) 2010-2016 Dynare Team +% Copyright (C) 2010-2017 Dynare Team % % This file is part of Dynare. % @@ -124,9 +124,6 @@ function [fval,info,exit_flag,DLIK,Hess,ys,trend_coeff,Model,DynareOptions,Bayes % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -% AUTHOR(S) stephane DOT adjemian AT univ DASH lemans DOT fr -% frederic DOT karame AT univ DASH lemans DOT fr - % Declaration of the penalty as a persistent variable. persistent init_flag persistent restrict_variables_idx observed_variables_idx state_variables_idx mf0 mf1 @@ -140,6 +137,9 @@ exit_flag = 1; DLIK = []; Hess = []; +% Ensure that xparam1 is a column vector. +xparam1 = xparam1(:); + % Issue an error if loglinear option is used. if DynareOptions.loglinear error('non_linear_dsge_likelihood: It is not possible to use a non linear filter with the option loglinear!') diff --git a/matlab/optimization/dynare_minimize_objective.m b/matlab/optimization/dynare_minimize_objective.m index 4724fc91a..0aebc6b50 100644 --- a/matlab/optimization/dynare_minimize_objective.m +++ b/matlab/optimization/dynare_minimize_objective.m @@ -401,8 +401,58 @@ switch minimizer_algorithm [LB, UB]=set_bounds_to_finite_values(bounds, options_.huge_number); [opt_par_values, fval, exitflag] = simpsa(func2str(objective_function),start_par_value,LB,UB,simpsaOptions,varargin{:}); case 11 - options_.cova_compute = 0 ; - [opt_par_values,stdh,lb_95,ub_95,med_param] = online_auxiliary_filter(start_par_value,varargin{:}) ; + options_.cova_compute = 0; + [opt_par_values, stdh, lb_95, ub_95, med_param] = online_auxiliary_filter(start_par_value, varargin{:}); + case 12 + [LB, UB] = set_bounds_to_finite_values(bounds, options_.huge_number); + tmp = transpose([fieldnames(options_.particleswarm), struct2cell(options_.particleswarm)]); + particleswarmOptions = optimoptions(@particleswarm); + particleswarmOptions = optimoptions(particleswarmOptions, tmp{:}); + if ~isempty(options_.optim_opt) + options_list = read_key_value_string(options_.optim_opt); + SupportedListOfOptions = {'CreationFcn', 'Display', 'DisplayInterval', 'FunctionTolerance', ... + 'FunValCheck', 'HybridFcn', 'InertiaRange', 'InitialSwarmMatrix', 'InitialSwarmSpan', ... + 'MaxIterations', 'MaxStallIterations', 'MaxStallTime', 'MaxTime', ... + 'MinNeighborsFraction', 'ObjectiveLimit', 'OutputFcn', 'PlotFcn', 'SelfAdjustmentWeight', ... + 'SocialAdjustmentWeight', 'SwarmSize', 'UseParallel', 'UseVectorized'}; + for i=1:rows(options_list) + if ismember(options_list{i,1}, SupportedListOfOptions) + particleswarmOptions = optimoptions(particleswarmOptions, options_list{i,1}, options_list{i,2}); + else + warning(['particleswarm: Unknown option (' options_list{i,1} ')!']) + end + end + end + % Get number of instruments. + numberofvariables = length(start_par_value); + % Set objective function. + objfun = @(x) objective_function(x, varargin{:}); + if ischar(particleswarmOptions.SwarmSize) + eval(['particleswarmOptions.SwarmSize = ' particleswarmOptions.SwarmSize ';']) + end + if isempty(particleswarmOptions.InitialSwarmMatrix) + particleswarmOptions.InitialSwarmMatrix = zeros(particleswarmOptions.SwarmSize, numberofvariables); + p = 1; + FVALS = zeros(particleswarmOptions.SwarmSize, 1); + while p<=particleswarmOptions.SwarmSize + candidate = rand(numberofvariables, 1).*(UB-LB)+LB; + [fval, info, exit_flag] = objfun(candidate); + if exit_flag + particleswarmOptions.InitialSwarmMatrix(p,:) = transpose(candidate); + FVALS(p) = fval; + p = p + 1; + end + end + end + % Set penalty to the worst value of the objective function. + TMP = [particleswarmOptions.InitialSwarmMatrix, FVALS]; + TMP = sortrows(TMP, length(start_par_value)+1); + penalty = TMP(end,end); + % Define penalized objective. + objfun = @(x) penalty_objective_function(x, objective_function, penalty, varargin{:}); + % Minimize the penalized objective (note that the penalty is not updated). + [opt_par_values, fval, exitflag, output] = particleswarm(objfun, length(start_par_value), LB, UB, particleswarmOptions); + opt_par_values = opt_par_values(:); case 101 solveoptoptions = options_.solveopt; if ~isempty(options_.optim_opt) diff --git a/tests/Makefile.am b/tests/Makefile.am index c2b935217..a85b41085 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -309,6 +309,7 @@ MODFILES = \ optimizers/fs2000_4_with_optim.mod \ optimizers/fs2000_5.mod \ optimizers/fs2000_7.mod \ + optimizers/fs2000_12.mod \ optimizers/fs2000_101.mod \ optimizers/fs2000_102.mod \ optimizers/fs2000_w.mod \ diff --git a/tests/optimizers/fs2000_12.mod b/tests/optimizers/fs2000_12.mod new file mode 100644 index 000000000..306d2aa80 --- /dev/null +++ b/tests/optimizers/fs2000_12.mod @@ -0,0 +1,3 @@ +@#include "fs2000.common.inc" + +estimation(mode_compute=12,order=1, datafile='../fs2000/fsdat_simul', nobs=192, mh_replic=0); From b14cf33c8ffe936c3c3ad79846335e918be7fb77 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Wed, 8 Feb 2017 12:36:55 +0100 Subject: [PATCH 05/17] Changed priors in tests/optimizers/fs2000_* integration tests. To avoid asymptote at zero on the autoregressive parameter rho. --- tests/optimizers/fs2000.common.inc | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/tests/optimizers/fs2000.common.inc b/tests/optimizers/fs2000.common.inc index 51c871bca..ee8510a3b 100644 --- a/tests/optimizers/fs2000.common.inc +++ b/tests/optimizers/fs2000.common.inc @@ -1,7 +1,7 @@ /* - * This file replicates the estimation of the cash in advance model described - * Frank Schorfheide (2000): "Loss function-based evaluation of DSGE models", - * Journal of Applied Econometrics, 15(6), 645-670. + * This file replicates (except for a small difference in the priors)the estimation of + * the cash in advance model described Frank Schorfheide (2000): "Loss function-based + * evaluation of DSGE models", Journal of Applied Econometrics, 15(6), 645-670. * * The data are in file "fs2000/fsdat_simul.m", and have been artificially generated. * They are therefore different from the original dataset used by Schorfheide. @@ -17,7 +17,7 @@ */ /* - * Copyright (C) 2004-2010 Dynare Team + * Copyright (C) 2004-2017 Dynare Team * * This file is part of Dynare. * @@ -106,7 +106,16 @@ alp, beta_pdf, 0.356, 0.02; bet, beta_pdf, 0.993, 0.002; gam, normal_pdf, 0.0085, 0.003; mst, normal_pdf, 1.0002, 0.007; -rho, beta_pdf, 0.129, 0.223; +/* + * The prior on the autoregressive parameter rho is different from the one considered in the original paper. + * We lowered the prior variance to ensure that the prior does not have an asymptote at zero, which causes + * troubles for some optimizers. For instance, mode_compute=12 (which relies on Mathworks' particleswarm + * routine) finds a posterior at zero for rho (with a much higher posterior density than the one reported + * by csminwell or other optimization routines). Because zero is on the boundary of the prior, this is an + * an issue when we have to compute the hessian matrix (which is not full rank) at the mode or if we try to + * use this solution as an initial guess for csminwell. + */ +rho, beta_pdf, 0.129, 0.100; psi, beta_pdf, 0.65, 0.05; del, beta_pdf, 0.01, 0.005; stderr e_a, inv_gamma_pdf, 0.035449, inf; From 9fbef0c1074280cff0c44f27d87d5fb58bbff29d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Wed, 8 Feb 2017 12:41:56 +0100 Subject: [PATCH 06/17] Removed penalty_hessian routine. + Code factorization. + Added an option for using the penalized objective when computing numerically the hessian at the mode. Previous behaviour (introduced with penalty_hessian routine) was to compute the hessian matrix at the mode with the penalized objective function (instead of the original objective function). This behaviour hides problematic situations, where the computed hessian (using the original objective) would not be full rank. For instance, if the estimation ends up with a parameter on (or very close to) the bounds of its possible values (which is often not a desirable outcome), the estimated posterior variance would be zero for this parameter (with the original objective) because the hessian is not finite in this direction, while the posterior variance would be positive if the penalized objective is used instead. But this estimate would not be reliable by construction of the penalty which is quite ad-hoc (more fundamentally I do not think that there exists any rational for approximating the covariance matrix with the inverse of the hessian matrix if the mode is on the border of the set of possible values). This commit restore the behaviour previous to 2446ab02ba4b3ed88c9c5021aced076078d96007. An option is available for computing the hessian with the penalized objective function. --- doc/dynare.texi | 10 +++ matlab/dynare_estimation_1.m | 12 ++-- matlab/global_initialization.m | 5 ++ matlab/optimization/penalty_hessian.m | 94 --------------------------- preprocessor/DynareBison.yy | 3 + preprocessor/DynareFlex.ll | 1 + 6 files changed, 27 insertions(+), 98 deletions(-) delete mode 100644 matlab/optimization/penalty_hessian.m diff --git a/doc/dynare.texi b/doc/dynare.texi index bf63caf70..0dc42dacf 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -6169,6 +6169,16 @@ This is the convergence criterion used in the fixed point Lyapunov solver. Its d @anchor{lyapunov_doubling_tol} This is the convergence criterion used in the doubling algorithm to solve the Lyapunov equation. Its default value is 1e-16. +@item use_penalized_objective_for_hessian +Use the penalized objective instead of the objective function to compute +numerically the hessian matrix at the mode. The penalties decrease the value of +the posterior density (or likelihood) when, for some perturbations, Dynare is +not able to solve the model (issues with steady state existence, Blanchard and +Kahn conditions, ...). In pratice, the penalized and original +objectives will only differ if the posterior mode is found to be near a region +where the model is ill-behaved. By default the original objective function is +used. + @item analytic_derivation Triggers estimation with analytic gradient. The final hessian is also computed analytically. Only works for stationary models without diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m index 6679f0cee..6077501ca 100644 --- a/matlab/dynare_estimation_1.m +++ b/matlab/dynare_estimation_1.m @@ -220,10 +220,14 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_); options_.analytic_derivation = ana_deriv_old; elseif ~isnumeric(options_.mode_compute) || ~(isequal(options_.mode_compute,5) && newratflag~=1), - % with flag==0, we force to use the hessian from outer - % product gradient of optimizer 5 - hh = reshape(penalty_hessian(objective_function,xparam1,fval, ... - options_.gstep,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_),nx,nx); + % with flag==0, we force to use the hessian from outer product gradient of optimizer 5 + if options_.hessian.use_penalized_objective + penalized_objective_function = str2func('penalty_objective_function'); + hh = hessian(penalized_objective_function, xparam1, options_.gstep, objective_function, fval, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_); + else + hh = hessian(objective_function, xparam1, options_.gstep, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_); + end + hh = reshape(hh, nx, nx); elseif isnumeric(options_.mode_compute) && isequal(options_.mode_compute,5) % other numerical hessian options available with optimizer 5 % diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index aa1531bb4..e8d1627e2 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -565,6 +565,11 @@ options_.homotopy_mode = 0; options_.homotopy_steps = 1; options_.homotopy_force_continue = 0; +% numerical hessian +hessian.use_penalized_objective = false; + +options_.hessian = hessian; + %csminwel optimization routine csminwel.tolerance.f=1e-7; csminwel.maxiter=1000; diff --git a/matlab/optimization/penalty_hessian.m b/matlab/optimization/penalty_hessian.m deleted file mode 100644 index 8da7bf6a9..000000000 --- a/matlab/optimization/penalty_hessian.m +++ /dev/null @@ -1,94 +0,0 @@ -function hessian_mat = penalty_hessian(func,x,penalty,gstep,varargin) % --*-- Unitary tests --*-- - -% Computes second order partial derivatives with penalty_objective_function -% -% INPUTS -% func [string] name of the function -% x [double] vector, the Hessian of "func" is evaluated at x. -% penalty [double] penalty base used if function fails -% gstep [double] scalar, size of epsilon. -% varargin [void] list of additional arguments for "func". -% -% OUTPUTS -% hessian_mat [double] Hessian matrix -% -% ALGORITHM -% Uses Abramowitz and Stegun (1965) formulas 25.3.23 -% \[ -% \frac{\partial^2 f_{0,0}}{\partial {x^2}} = \frac{1}{h^2}\left( f_{1,0} - 2f_{0,0} + f_{ - 1,0} \right) -% \] -% and 25.3.27 p. 884 -% -% \[ -% \frac{\partial ^2f_{0,0}}{\partial x\partial y} = \frac{-1}{2h^2}\left(f_{1,0} + f_{-1,0} + f_{0,1} + f_{0,-1} - 2f_{0,0} - f_{1,1} - f_{-1,-1} \right) -% \] -% -% SPECIAL REQUIREMENTS -% none -% - -% Copyright (C) 2001-2014 Dynare Team -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see . - -if ~isa(func, 'function_handle') - func = str2func(func); -end -n=size(x,1); -h1=max(abs(x),sqrt(gstep(1))*ones(n,1))*eps^(1/6)*gstep(2); -h_1=h1; -xh1=x+h1; -h1=xh1-x; -xh1=x-h_1; -h_1=x-xh1; -xh1=x; -f0=penalty_objective_function(x,func,penalty,varargin{:}); -f1=zeros(size(f0,1),n); -f_1=f1; -for i=1:n - %do step up - xh1(i)=x(i)+h1(i); - f1(:,i)=penalty_objective_function(xh1,func,penalty,varargin{:}); - %do step up - xh1(i)=x(i)-h_1(i); - f_1(:,i)=penalty_objective_function(xh1,func,penalty,varargin{:}); - xh1(i)=x(i);%reset parameter -end -xh_1=xh1; -hessian_mat = zeros(size(f0,1),n*n); -temp=f1+f_1-f0*ones(1,n); %term f_(1,0)+f_(-1,0)-f_(0,0) used later -for i=1:n - if i > 1 %fill symmetric part of Hessian based on previously computed results - k=[i:n:n*(i-1)]; - hessian_mat(:,(i-1)*n+1:(i-1)*n+i-1)=hessian_mat(:,k); - end - hessian_mat(:,(i-1)*n+i)=(f1(:,i)+f_1(:,i)-2*f0)./(h1(i)*h_1(i)); %formula 25.3.23 - for j=i+1:n - %step in up direction - xh1(i)=x(i)+h1(i); - xh1(j)=x(j)+h_1(j); - %step in down direction - xh_1(i)=x(i)-h1(i); - xh_1(j)=x(j)-h_1(j); - hessian_mat(:,(i-1)*n+j)=-(-penalty_objective_function(xh1,func,penalty,varargin{:})-penalty_objective_function(xh_1,func,penalty,varargin{:})+temp(:,i)+temp(:,j))./(2*h1(i)*h_1(j)); %formula 25.3.27 - %reset grid points - xh1(i)=x(i); - xh1(j)=x(j); - xh_1(i)=x(i); - xh_1(j)=x(j); - end -end - diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy index bbae98758..7f6d28834 100644 --- a/preprocessor/DynareBison.yy +++ b/preprocessor/DynareBison.yy @@ -113,6 +113,7 @@ class ParsingDriver; %token CPF_WEIGHTS AMISANOTRISTANI MURRAYJONESPARSLOW %token FILTER_ALGORITHM PROPOSAL_APPROXIMATION CUBATURE UNSCENTED MONTECARLO DISTRIBUTION_APPROXIMATION %token NAME +%token USE_PENALIZED_OBJECTIVE_FOR_HESSIAN %token NAN_CONSTANT NO_STATIC NOBS NOCONSTANT NODISPLAY NOCORR NODIAGNOSTIC NOFUNCTIONS NO_HOMOTOPY %token NOGRAPH POSTERIOR_NOGRAPH POSTERIOR_GRAPH NOMOMENTS NOPRINT NORMAL_PDF SAVE_DRAWS %token OBSERVATION_TRENDS OPTIM OPTIM_WEIGHTS ORDER OSR OSR_PARAMS MAX_DIM_COVA_GROUP ADVANCED OUTFILE OUTVARS OVERWRITE @@ -1825,6 +1826,7 @@ estimation_options : o_datafile | o_posterior_sampling_method | o_posterior_sampler_options | o_keep_kalman_algo_if_singularity_is_detected + | o_use_penalized_objective_for_hessian ; list_optim_option : QUOTED_STRING COMMA QUOTED_STRING @@ -3232,6 +3234,7 @@ o_mcmc_jumping_covariance : MCMC_JUMPING_COVARIANCE EQUAL HESSIAN | MCMC_JUMPING_COVARIANCE EQUAL filename { driver.option_str("MCMC_jumping_covariance", $3); } ; +o_use_penalized_objective_for_hessian : USE_PENALIZED_OBJECTIVE_FOR_HESSIAN { driver.option_num("hessian.use_penalized_objective","true"); }; o_irf_plot_threshold : IRF_PLOT_THRESHOLD EQUAL non_negative_number { driver.option_num("impulse_responses.plot_threshold", $3); }; o_dr_display_tol : DR_DISPLAY_TOL EQUAL non_negative_number { driver.option_num("dr_display_tol", $3); }; o_consider_all_endogenous : CONSIDER_ALL_ENDOGENOUS { driver.option_str("endo_vars_for_moment_computations_in_estimation", "all_endogenous_variables"); }; diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll index 8d00b2874..e1e0aa4c5 100644 --- a/preprocessor/DynareFlex.ll +++ b/preprocessor/DynareFlex.ll @@ -400,6 +400,7 @@ DATE -?[0-9]+([YyAa]|[Mm]([1-9]|1[0-2])|[Qq][1-4]|[Ww]([1-9]{1}|[1-4][0-9]|5[0-2 distribution_approximation {return token::DISTRIBUTION_APPROXIMATION;} proposal_distribution {return token::PROPOSAL_DISTRIBUTION;} no_posterior_kernel_density {return token::NO_POSTERIOR_KERNEL_DENSITY;} +use_penalized_objective_for_hessian {return token::USE_PENALIZED_OBJECTIVE_FOR_HESSIAN;} alpha { yylval->string_val = new string(yytext); From 593d62a7c14d7a2632e0d6295f5205de05a9e20e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Wed, 8 Feb 2017 14:29:34 +0100 Subject: [PATCH 07/17] Removed optimizers/fs2000_12.mod from testsuite... Because the Global Optimization Toolbox is not installed on dynbot@sedna. --- tests/Makefile.am | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/Makefile.am b/tests/Makefile.am index a85b41085..c2b935217 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -309,7 +309,6 @@ MODFILES = \ optimizers/fs2000_4_with_optim.mod \ optimizers/fs2000_5.mod \ optimizers/fs2000_7.mod \ - optimizers/fs2000_12.mod \ optimizers/fs2000_101.mod \ optimizers/fs2000_102.mod \ optimizers/fs2000_w.mod \ From edd2e980869856656e44ccc8848444a309c06bcb Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Thu, 9 Feb 2017 12:46:37 +0100 Subject: [PATCH 08/17] preprocessor: replace exit(1) with exit(EXIT_FAILURE) --- preprocessor/ComputingTasks.cc | 4 ++-- preprocessor/DynamicModel.cc | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc index f12c40867..ca83e94e8 100644 --- a/preprocessor/ComputingTasks.cc +++ b/preprocessor/ComputingTasks.cc @@ -1,5 +1,5 @@ /* - * Copyright (C) 2003-2016 Dynare Team + * Copyright (C) 2003-2017 Dynare Team * * This file is part of Dynare. * @@ -336,7 +336,7 @@ RamseyConstraintsStatement::writeOutput(ostream &output, const string &basename, break; default: cerr << "Ramsey constraints: this shouldn't happen." << endl; - exit(1); + exit(EXIT_FAILURE); } output << "', '"; it->expression->writeOutput(output); diff --git a/preprocessor/DynamicModel.cc b/preprocessor/DynamicModel.cc index 08dd2bc45..21b18d7a6 100644 --- a/preprocessor/DynamicModel.cc +++ b/preprocessor/DynamicModel.cc @@ -5012,7 +5012,7 @@ DynamicModel::writeFirstDerivativesC_csr(const string &basename, bool cuda) cons break; default: std::cerr << "This case shouldn't happen" << std::endl; - exit(1); + exit(EXIT_FAILURE); } derivative deriv(col_id + eq*cols_nbr,col_id,eq,it->second); D.push_back(deriv); From 7ced4b5580c0eba3846132bb2e44f3d1a322b9dc Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Thu, 9 Feb 2017 16:21:18 +0100 Subject: [PATCH 09/17] =?UTF-8?q?preprocessor:=20change=20error=20to=20war?= =?UTF-8?q?ning=20because=20we=20don=E2=80=99t=20end=20processing=20in=20t?= =?UTF-8?q?his=20situation;=20fix=20error=20message?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- preprocessor/ComputingTasks.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc index ca83e94e8..9843ba7d5 100644 --- a/preprocessor/ComputingTasks.cc +++ b/preprocessor/ComputingTasks.cc @@ -1058,7 +1058,7 @@ ObservationTrendsStatement::writeOutput(ostream &output, const string &basename, output << "';" << endl; } else - cout << "Error : Non-variable symbol used in TREND_COEFF: " << it->first << endl; + cerr << "Warning : Non-variable symbol used in observation_trends: " << it->first << endl; } } From fabca90345709e0e92deb96efb96a88891fc4398 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Thu, 9 Feb 2017 16:22:00 +0100 Subject: [PATCH 10/17] preprocessor: aesthetic fix --- preprocessor/ComputingTasks.cc | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc index 9843ba7d5..ae580372d 100644 --- a/preprocessor/ComputingTasks.cc +++ b/preprocessor/ComputingTasks.cc @@ -1044,10 +1044,7 @@ void ObservationTrendsStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const { output << "options_.trend_coeff = {};" << endl; - - trend_elements_t::const_iterator it; - - for (it = trend_elements.begin(); it != trend_elements.end(); it++) + for (trend_elements_t::const_iterator it = trend_elements.begin(); it != trend_elements.end(); it++) { SymbolType type = symbol_table.getType(it->first); if (type == eEndogenous) From 39855c95d2724f60af749a127afd28ff88a621a8 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Fri, 10 Feb 2017 11:07:45 +0100 Subject: [PATCH 11/17] preprocessor: aesthetic fix --- preprocessor/ComputingTasks.cc | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc index ae580372d..a05727d32 100644 --- a/preprocessor/ComputingTasks.cc +++ b/preprocessor/ComputingTasks.cc @@ -821,9 +821,7 @@ EstimatedParamsStatement::writeOutput(ostream &output, const string &basename, b << "estim_params_.corrn = [];" << endl << "estim_params_.param_vals = [];" << endl; - vector::const_iterator it; - - for (it = estim_params_list.begin(); it != estim_params_list.end(); it++) + for (vector::const_iterator it = estim_params_list.begin(); it != estim_params_list.end(); it++) { int symb_id = symbol_table.getTypeSpecificID(it->name) + 1; SymbolType symb_type = symbol_table.getType(it->name); @@ -892,9 +890,7 @@ EstimatedParamsInitStatement::writeOutput(ostream &output, const string &basenam if (use_calibration) output << "options_.use_calibration_initialization = 1;" << endl; - vector::const_iterator it; - - for (it = estim_params_list.begin(); it != estim_params_list.end(); it++) + for (vector::const_iterator it = estim_params_list.begin(); it != estim_params_list.end(); it++) { int symb_id = symbol_table.getTypeSpecificID(it->name) + 1; SymbolType symb_type = symbol_table.getType(it->name); @@ -955,9 +951,7 @@ EstimatedParamsBoundsStatement::EstimatedParamsBoundsStatement(const vector::const_iterator it; - - for (it = estim_params_list.begin(); it != estim_params_list.end(); it++) + for (vector::const_iterator it = estim_params_list.begin(); it != estim_params_list.end(); it++) { int symb_id = symbol_table.getTypeSpecificID(it->name) + 1; SymbolType symb_type = symbol_table.getType(it->name); From 8925e7ca1fe6b1f8db461b8a7e09650ed2b877b6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Mon, 13 Feb 2017 13:57:21 +0100 Subject: [PATCH 12/17] Cosmetic patch (remove some warnings about potentially uninitialized variables in bytecode). --- mex/build/bytecode.am | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mex/build/bytecode.am b/mex/build/bytecode.am index 11936ea5b..3fff91434 100644 --- a/mex/build/bytecode.am +++ b/mex/build/bytecode.am @@ -1,6 +1,6 @@ mex_PROGRAMS = bytecode -bytecode_CPPFLAGS = $(AM_CPPFLAGS) -I$(top_srcdir)/../../sources -I$(top_srcdir)/../../sources/bytecode -I$(top_srcdir)/../../../preprocessor +bytecode_CPPFLAGS = -Wno-maybe-uninitialized $(AM_CPPFLAGS) -I$(top_srcdir)/../../sources -I$(top_srcdir)/../../sources/bytecode -I$(top_srcdir)/../../../preprocessor TOPDIR = $(top_srcdir)/../../sources/bytecode From ffebf6b262039a2682ad2cd2f1bbbaaab192dcc1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?= Date: Wed, 15 Feb 2017 12:41:56 +0100 Subject: [PATCH 13/17] Check positive steady state for loglinear option... 1e-9 is too far from zero. Set 1e-16 as a threshold value for the positivity test. --- matlab/resol.m | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/matlab/resol.m b/matlab/resol.m index 2ad99cf2c..ec390ba0f 100644 --- a/matlab/resol.m +++ b/matlab/resol.m @@ -109,9 +109,11 @@ if info(1) end if options.loglinear + threshold = 1e-16; % Find variables with non positive steady state. Skip auxiliary % variables for lagges/leaded exogenous variables - idx = find(dr.ys(get_all_variables_but_lagged_leaded_exogenous(M)) < 1e-9); + idx = find(dr.ys(get_all_variables_but_lagged_leaded_exogenous(M)) Date: Wed, 15 Feb 2017 13:35:33 +0100 Subject: [PATCH 14/17] Removed printing of debugging info. Introduced in ffebf6b262039a2682ad2cd2f1bbbaaab192dcc1. --- matlab/resol.m | 1 - 1 file changed, 1 deletion(-) diff --git a/matlab/resol.m b/matlab/resol.m index ec390ba0f..31ce81ede 100644 --- a/matlab/resol.m +++ b/matlab/resol.m @@ -113,7 +113,6 @@ if options.loglinear % Find variables with non positive steady state. Skip auxiliary % variables for lagges/leaded exogenous variables idx = find(dr.ys(get_all_variables_but_lagged_leaded_exogenous(M)) Date: Fri, 24 Feb 2017 11:20:54 +0100 Subject: [PATCH 15/17] preprocessor: for consistency, use equations.size() instead of equation_number() in classes that have equations as a field --- preprocessor/DynamicModel.cc | 12 ++++++------ preprocessor/ModelTree.cc | 20 ++++++++++---------- preprocessor/StaticModel.cc | 8 ++++---- 3 files changed, 20 insertions(+), 20 deletions(-) diff --git a/preprocessor/DynamicModel.cc b/preprocessor/DynamicModel.cc index 21b18d7a6..ebf6b9b0e 100644 --- a/preprocessor/DynamicModel.cc +++ b/preprocessor/DynamicModel.cc @@ -3252,8 +3252,8 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative if (!no_tmp_terms) computeTemporaryTermsOrdered(); int k = 0; - equation_block = vector(equation_number()); - variable_block_lead_lag = vector< pair< int, pair< int, int> > >(equation_number()); + equation_block = vector(equations.size()); + variable_block_lead_lag = vector< pair< int, pair< int, int> > >(equations.size()); for (unsigned int i = 0; i < getNbBlocks(); i++) { for (unsigned int j = 0; j < getBlockSize(i); j++) @@ -4210,10 +4210,10 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con << "% from model file (.mod)" << endl << endl << model_local_vars_output.str() << model_output.str() - << "rp = zeros(" << equation_number() << ", " + << "rp = zeros(" << equations.size() << ", " << symbol_table.param_nbr() << ");" << endl << jacobian_output.str() - << "gp = zeros(" << equation_number() << ", " << dynJacobianColsNbr << ", " << symbol_table.param_nbr() << ");" << endl + << "gp = zeros(" << equations.size() << ", " << dynJacobianColsNbr << ", " << symbol_table.param_nbr() << ");" << endl << hessian_output.str() << "if nargout >= 3" << endl << "rpp = zeros(" << residuals_params_second_derivatives.size() << ",4);" << endl @@ -4238,10 +4238,10 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con << "ss_param_deriv, ss_param_2nd_deriv)" << endl << model_local_vars_output.str() << model_output.str() - << "rp = zeros(" << equation_number() << ", " + << "rp = zeros(" << equations.size() << ", " << symbol_table.param_nbr() << ");" << endl << jacobian_output.str() - << "gp = zeros(" << equation_number() << ", " << dynJacobianColsNbr << ", " << symbol_table.param_nbr() << ");" << endl + << "gp = zeros(" << equations.size() << ", " << dynJacobianColsNbr << ", " << symbol_table.param_nbr() << ");" << endl << hessian_output.str() << "rpp = zeros(" << residuals_params_second_derivatives.size() << ",4);" << endl << hessian1_output.str() diff --git a/preprocessor/ModelTree.cc b/preprocessor/ModelTree.cc index 70db28e01..1faac8ff6 100644 --- a/preprocessor/ModelTree.cc +++ b/preprocessor/ModelTree.cc @@ -36,7 +36,7 @@ using namespace MFS; bool ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, bool verbose) { - const int n = equation_number(); + const int n = equations.size(); assert(n == symbol_table.endo_nbr()); @@ -96,7 +96,7 @@ ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, boo #endif // Create the resulting map, by copying the n first elements of mate_map, and substracting n to them - endo2eq.resize(equation_number()); + endo2eq.resize(equations.size()); transform(mate_map.begin(), mate_map.begin() + n, endo2eq.begin(), bind2nd(minus(), n)); #ifdef DEBUG @@ -143,7 +143,7 @@ ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian cout << "Normalizing the model..." << endl; - int n = equation_number(); + int n = equations.size(); // compute the maximum value of each row of the contemporaneous Jacobian matrix //jacob_map normalized_contemporaneous_jacobian; @@ -329,7 +329,7 @@ ModelTree::writeRevXrefs(ostream &output, const map > &xrefmap, co void ModelTree::computeNormalizedEquations(multimap &endo2eqs) const { - for (int i = 0; i < equation_number(); i++) + for (int i = 0; i < equations.size(); i++) { VariableNode *lhs = dynamic_cast(equations[i]->get_arg1()); if (lhs == NULL) @@ -417,11 +417,11 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, jacob_m void ModelTree::computePrologueAndEpilogue(const jacob_map_t &static_jacobian_arg, vector &equation_reordered, vector &variable_reordered) { - vector eq2endo(equation_number(), 0); - equation_reordered.resize(equation_number()); - variable_reordered.resize(equation_number()); + vector eq2endo(equations.size(), 0); + equation_reordered.resize(equations.size()); + variable_reordered.resize(equations.size()); bool *IM; - int n = equation_number(); + int n = equations.size(); IM = (bool *) calloc(n*n, sizeof(bool)); int i = 0; for (vector::const_iterator it = endo2eq.begin(); it != endo2eq.end(); it++, i++) @@ -1697,7 +1697,7 @@ ModelTree::addEquation(expr_t eq, int lineno) void ModelTree::addEquation(expr_t eq, int lineno, vector > &eq_tags) { - int n = equation_number(); + int n = equations.size(); for (size_t i = 0; i < eq_tags.size(); i++) equation_tags.push_back(make_pair(n, eq_tags[i])); addEquation(eq, lineno); @@ -1741,7 +1741,7 @@ ModelTree::addNonstationaryVariables(vector nonstationary_vars, bool log_de void ModelTree::initializeVariablesAndEquations() { - for (int j = 0; j < equation_number(); j++) + for (int j = 0; j < equations.size(); j++) { equation_reordered.push_back(j); variable_reordered.push_back(j); diff --git a/preprocessor/StaticModel.cc b/preprocessor/StaticModel.cc index a7652f9b8..2cdb2f0a7 100644 --- a/preprocessor/StaticModel.cc +++ b/preprocessor/StaticModel.cc @@ -2368,10 +2368,10 @@ StaticModel::writeParamsDerivativesFile(const string &basename, bool julia) cons << "% from model file (.mod)" << endl << endl << model_local_vars_output.str() << model_output.str() - << "rp = zeros(" << equation_number() << ", " + << "rp = zeros(" << equations.size() << ", " << symbol_table.param_nbr() << ");" << endl << jacobian_output.str() - << "gp = zeros(" << equation_number() << ", " << symbol_table.endo_nbr() << ", " + << "gp = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ", " << symbol_table.param_nbr() << ");" << endl << hessian_output.str() << "if nargout >= 3" << endl @@ -2396,10 +2396,10 @@ StaticModel::writeParamsDerivativesFile(const string &basename, bool julia) cons << "function params_derivs(y, x, params)" << endl << model_local_vars_output.str() << model_output.str() - << "rp = zeros(" << equation_number() << ", " + << "rp = zeros(" << equations.size() << ", " << symbol_table.param_nbr() << ");" << endl << jacobian_output.str() - << "gp = zeros(" << equation_number() << ", " << symbol_table.endo_nbr() << ", " + << "gp = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ", " << symbol_table.param_nbr() << ");" << endl << hessian_output.str() << "rpp = zeros(" << residuals_params_second_derivatives.size() << ",4);" << endl From 2d25858886caf58ff79c2d0977584d12d4093442 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Wed, 1 Mar 2017 12:20:55 +0100 Subject: [PATCH 16/17] preprocessor: fix bug that caused temporary terms not to be printed for the parameter derivatives of the static model. closes #1397 --- preprocessor/StaticModel.hh | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/preprocessor/StaticModel.hh b/preprocessor/StaticModel.hh index bf6cad6f3..8535a7140 100644 --- a/preprocessor/StaticModel.hh +++ b/preprocessor/StaticModel.hh @@ -1,5 +1,5 @@ /* - * Copyright (C) 2003-2016 Dynare Team + * Copyright (C) 2003-2017 Dynare Team * * This file is part of Dynare. * @@ -30,9 +30,6 @@ using namespace std; class StaticModel : public ModelTree { private: - //! Temporary terms for the file containing parameters dervicatives - temporary_terms_t params_derivs_temporary_terms; - //! global temporary terms for block decomposed models vector > v_temporary_terms; From 33b38ab643c5dc33bf17af464f0aa8ca2c62d949 Mon Sep 17 00:00:00 2001 From: Houtan Bastani Date: Wed, 1 Mar 2017 12:44:31 +0100 Subject: [PATCH 17/17] preprocessor: write error messages to cerr instead of cout, replace \n with endl --- preprocessor/DynamicModel.cc | 6 +++--- preprocessor/ModelTree.cc | 2 +- preprocessor/StaticModel.cc | 6 +++--- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/preprocessor/DynamicModel.cc b/preprocessor/DynamicModel.cc index ebf6b9b0e..f3c90c257 100644 --- a/preprocessor/DynamicModel.cc +++ b/preprocessor/DynamicModel.cc @@ -800,7 +800,7 @@ DynamicModel::writeModelEquationsCode(string &file_name, const string &bin_basen code_file.open(main_name.c_str(), ios::out | ios::binary | ios::ate); if (!code_file.is_open()) { - cout << "Error : Can't open file \"" << main_name << "\" for writing\n"; + cerr << "Error : Can't open file \"" << main_name << "\" for writing" << endl; exit(EXIT_FAILURE); } @@ -1076,7 +1076,7 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin code_file.open(main_name.c_str(), ios::out | ios::binary | ios::ate); if (!code_file.is_open()) { - cout << "Error : Can't open file \"" << main_name << "\" for writing\n"; + cerr << "Error : Can't open file \"" << main_name << "\" for writing" << endl; exit(EXIT_FAILURE); } //Temporary variables declaration @@ -1768,7 +1768,7 @@ DynamicModel::Write_Inf_To_Bin_File_Block(const string &dynamic_basename, const SaveCode.open((bin_basename + "_dynamic.bin").c_str(), ios::out | ios::binary); if (!SaveCode.is_open()) { - cout << "Error : Can't open file \"" << bin_basename << "_dynamic.bin\" for writing\n"; + cerr << "Error : Can't open file \"" << bin_basename << "_dynamic.bin\" for writing" << endl; exit(EXIT_FAILURE); } u_count_int = 0; diff --git a/preprocessor/ModelTree.cc b/preprocessor/ModelTree.cc index 1faac8ff6..f295423e3 100644 --- a/preprocessor/ModelTree.cc +++ b/preprocessor/ModelTree.cc @@ -1594,7 +1594,7 @@ ModelTree::Write_Inf_To_Bin_File(const string &basename, SaveCode.open(bin_basename.c_str(), ios::out | ios::binary); if (!SaveCode.is_open()) { - cout << "Error : Can't open file \"" << bin_basename << "\" for writing\n"; + cerr << "Error : Can't open file \"" << bin_basename << "\" for writing" << endl; exit(EXIT_FAILURE); } u_count_int = 0; diff --git a/preprocessor/StaticModel.cc b/preprocessor/StaticModel.cc index 2cdb2f0a7..1d1760274 100644 --- a/preprocessor/StaticModel.cc +++ b/preprocessor/StaticModel.cc @@ -412,7 +412,7 @@ StaticModel::writeModelEquationsCode(const string file_name, const string bin_ba code_file.open(main_name.c_str(), ios::out | ios::binary | ios::ate); if (!code_file.is_open()) { - cout << "Error : Can't open file \"" << main_name << "\" for writing\n"; + cerr << "Error : Can't open file \"" << main_name << "\" for writing" << endl; exit(EXIT_FAILURE); } int count_u; @@ -596,7 +596,7 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string code_file.open(main_name.c_str(), ios::out | ios::binary | ios::ate); if (!code_file.is_open()) { - cout << "Error : Can't open file \"" << main_name << "\" for writing\n"; + cerr << "Error : Can't open file \"" << main_name << "\" for writing" << endl; exit(EXIT_FAILURE); } //Temporary variables declaration @@ -990,7 +990,7 @@ StaticModel::Write_Inf_To_Bin_File_Block(const string &static_basename, const st SaveCode.open((bin_basename + "_static.bin").c_str(), ios::out | ios::binary); if (!SaveCode.is_open()) { - cout << "Error : Can't open file \"" << bin_basename << "_static.bin\" for writing\n"; + cerr << "Error : Can't open file \"" << bin_basename << "_static.bin\" for writing" << endl; exit(EXIT_FAILURE); } u_count_int = 0;