From 27886047cf1c889dab59ba6c56316aa56c56dbbf Mon Sep 17 00:00:00 2001 From: adjemian Date: Mon, 7 Jul 2008 12:58:13 +0000 Subject: [PATCH] * Changed numgrad.m. If fcn(x+h) is not well defined (for instance cost_flag=0 because of BK conditions) we compute the gradient using (fcn(x)-fcn(x-h))/h instead of the original formula (fcn(x+h)-fcn(x))/h. * Added function numgrad3.m. This new function uses a three point formula to compute the gradient. By default Dynare uses the two point formula. By setting options_.gradient_method=3 Dynare switches to the three point formula. Note that an input argument has been added to csminwel. * Added automatic detection of needed constant in dsge-var model (may be problematic is some cases) git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1935 ac1d8469-bf42-47a9-8791-bf33cf982152 --- matlab/csminwel.m | 79 ++++++++++++++++------------------ matlab/dynare_estimation.m | 19 +++++--- matlab/global_initialization.m | 5 ++- matlab/numgrad.m | 26 +++++++---- matlab/numgrad3.m | 45 +++++++++++++++++++ matlab/osr1.m | 2 +- 6 files changed, 116 insertions(+), 60 deletions(-) create mode 100644 matlab/numgrad3.m diff --git a/matlab/csminwel.m b/matlab/csminwel.m index 44bb756ef..6d0edd845 100644 --- a/matlab/csminwel.m +++ b/matlab/csminwel.m @@ -1,5 +1,5 @@ -function [fh,xh,gh,H,itct,fcount,retcodeh] = csminwel(fcn,x0,H0,grad,crit,nit,varargin) -%[fhat,xhat,ghat,Hhat,itct,fcount,retcodehat] = csminwel(fcn,x0,H0,grad,crit,nit,varargin) +function [fh,xh,gh,H,itct,fcount,retcodeh] = csminwel(fcn,x0,H0,grad,crit,nit,method,varargin) +%[fhat,xhat,ghat,Hhat,itct,fcount,retcodehat] = csminwel(fcn,x0,H0,grad,crit,nit,method,varargin) % fcn: string naming the objective function to be minimized % x0: initial value of the parameter vector % H0: initial value for the inverse Hessian. Must be positive definite. @@ -9,6 +9,7 @@ function [fh,xh,gh,H,itct,fcount,retcodeh] = csminwel(fcn,x0,H0,grad,crit,nit,va % crit: Convergence criterion. Iteration will cease when it proves impossible to improve the % function value by more than crit. % nit: Maximum number of iterations. +% method: integer scalar, 2 or 3 points formula. % varargin: A list of optional length of additional parameters that get handed off to fcn each % time it is called. % Note that if the program ends abnormally, it is possible to retrieve the current x, @@ -37,25 +38,21 @@ snit=100; % tailstr=[ ',P' num2str(i) tailstr]; % stailstr=[' P' num2str(i) stailstr]; %end -f0 = feval(fcn,x0,varargin{:}); -%ARGLIST -%f0 = feval(fcn,x0,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13); -% disp('first fcn in csminwel.m ----------------') % Jinill on 9/5/95 -if f0 > 1e50, disp('Bad initial parameter.'), return, end + +[f0,cost_flag] = feval(fcn,x0,varargin{:}); +if ~cost_flag + disp('Bad initial parameter.') + return +end + if NumGrad - if length(grad)==0 - [g badg] = numgrad(fcn,x0, varargin{:}); - %ARGLIST - %[g badg] = numgrad(fcn,x0,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13); - else - badg=any(find(grad==0)); - g=grad; - end - %numgrad(fcn,x0,P1,P2,P3,P4); + if method==2 + [g,badg] = numgrad(fcn,x0, varargin{:}); + elseif method==3 + [g,badg] = numgrad3(fcn,x0, varargin{:}); + end else - [g badg] = feval(grad,x0,varargin{:}); - %ARGLIST - %[g badg] = feval(grad,x0,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13); + [g,badg] = feval(grad,x0,varargin{:}); end retcode3=101; x=x0; @@ -93,15 +90,13 @@ while ~done wall1=1; badg1=1; else if NumGrad - [g1 badg1] = numgrad(fcn, x1,varargin{:}); - %ARGLIST - %[g1 badg1] = numgrad(fcn, x1,P1,P2,P3,P4,P5,P6,P7,P8,P9,... - % P10,P11,P12,P13); + if method == 2 + [g1 badg1] = numgrad(fcn, x1,varargin{:}); + elseif method == 3 + [g1 badg1] = numgrad3(fcn, x1,varargin{:}); + end else [g1 badg1] = feval(grad,x1,varargin{:}); - %ARGLIST - %[g1 badg1] = feval(grad, x1,P1,P2,P3,P4,P5,P6,P7,P8,P9,... - % P10,P11,P12,P13); end wall1=badg1; % g1 @@ -126,15 +121,13 @@ while ~done wall2=1; badg2=1; else if NumGrad - [g2 badg2] = numgrad(fcn, x2,varargin{:}); - %ARGLIST - %[g2 badg2] = numgrad(fcn, x2,P1,P2,P3,P4,P5,P6,P7,P8,... - % P9,P10,P11,P12,P13); + if method==2 + [g2 badg2] = numgrad(fcn, x2,varargin{:}); + elseif method==3 + [g2 badg2] = numgrad3(fcn, x2,varargin{:}); + end else [g2 badg2] = feval(grad,x2,varargin{:}); - %ARGLIST - %[g2 badg2] = feval(grad,x2,P1,P2,P3,P4,P5,P6,P7,P8,... - % P9,P10,P11,P12,P13); end wall2=badg2; % g2 @@ -160,15 +153,13 @@ while ~done wall3=1; badg3=1; else if NumGrad - [g3 badg3] = numgrad(fcn, x3,varargin{:}); - %ARGLIST - %[g3 badg3] = numgrad(fcn, x3,P1,P2,P3,P4,P5,P6,P7,P8,... - % P9,P10,P11,P12,P13); + if method==2 + [g3 badg3] = numgrad(fcn, x3,varargin{:}); + elseif method==3 + [g3 badg3] = numgrad3(fcn, x3,varargin{:}); + end else [g3 badg3] = feval(grad,x3,varargin{:}); - %ARGLIST - %[g3 badg3] = feval(grad,x3,P1,P2,P3,P4,P5,P6,P7,P8,... - % P9,P10,P11,P12,P13); end wall3=badg3; % g3 @@ -224,7 +215,11 @@ while ~done end if nogh if NumGrad - [gh badgh] = numgrad(fcn,xh,varargin{:}); + if method==2 + [gh,badgh] = numgrad(fcn,xh,varargin{:}); + elseif method==3 + [gh,badgh] = numgrad3(fcn,xh,varargin{:}); + end else [gh badgh] = feval(grad, xh,varargin{:}); end @@ -274,4 +269,4 @@ while ~done badg=badgh; end % what about making an m-file of 10 lines including numgrad.m -% since it appears three times in csminwel.m +% since it appears three times in csminwel.m \ No newline at end of file diff --git a/matlab/dynare_estimation.m b/matlab/dynare_estimation.m index ed50be843..7c8745af2 100644 --- a/matlab/dynare_estimation.m +++ b/matlab/dynare_estimation.m @@ -1,5 +1,4 @@ function dynare_estimation(var_list_) - % function dynare_estimation(var_list_) % runs the estimation of the model % @@ -17,8 +16,7 @@ function dynare_estimation(var_list_) global M_ options_ oo_ estim_params_ -global bayestopt_ dsge_prior_weight - +global bayestopt_ options_.varlist = var_list_; options_.lgyidx2varobs = zeros(size(M_.endo_names,1),1); @@ -259,9 +257,16 @@ else% if the steady state file is not provided. [dd,info] = resol(oo_.steady_state,0); oo_.steady_state = dd.ys; clear('dd'); end +if all(abs(oo_.steady_state(bayestopt_.mfys))<1e-9) + disp('no constant') + options_.noconstant = 1; +else + options_.noconstant = 0; +end + %% compute sample moments if needed (bvar-dsge) -if options_.bvar_dsge~isempty(strmatch('dsge_prior_weight',M_.param_names)) +if options_.bvar_dsge if options_.noconstant evalin('base',... ['[mYY,mXY,mYX,mXX,Ydata,Xdata] = ' ... @@ -337,7 +342,7 @@ if options_.mode_compute > 0 & options_.posterior_mode_estimation if isfield(options_,'optim_opt') eval(['optim_options = optimset(optim_options,' options_.optim_opt ');']); end - if isempty(strmatch('dsge_prior_weight',M_.param_names)) + if ~options_.bvar_dsge [xparam1,fval,exitflag] = fminunc(fh,xparam1,optim_options,gend,data); else [xparam1,fval,exitflag] = fminunc(fh,xparam1,optim_options,gend); @@ -349,12 +354,12 @@ if options_.mode_compute > 0 & options_.posterior_mode_estimation verbose = 2; if ~options_.bvar_dsge [fval,xparam1,grad,hessian_csminwel,itct,fcount,retcodehat] = ... - csminwel('DsgeLikelihood',xparam1,H0,[],crit,nit,gend,data); + csminwel('DsgeLikelihood',xparam1,H0,[],crit,nit,options_.gradient_method,gend,data); disp(sprintf('Objective function at mode: %f',fval)) disp(sprintf('Objective function at mode: %f',DsgeLikelihood(xparam1,gend,data))) else [fval,xparam1,grad,hessian_csminwel,itct,fcount,retcodehat] = ... - csminwel('DsgeVarLikelihood',xparam1,H0,[],crit,nit,gend); + csminwel('DsgeVarLikelihood',xparam1,H0,[],crit,nit,options_.gradient_method,gend); disp(sprintf('Objective function at mode: %f',fval)) disp(sprintf('Objective function at mode: %f',DsgeVarLikelihood(xparam1,gend))) end diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index a89860f89..3e27eb0e3 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -50,8 +50,8 @@ function global_initialization() options_.varlag = 4; % Optimization algorithm [6] gmhmaxlik - options_.Opt6Iter = 3; - options_.Opt6Numb = 100000; + options_.Opt6Iter = 2; + options_.Opt6Numb = 5000; % Graphics options_.graphics.nrows = 3; @@ -148,6 +148,7 @@ function global_initialization() options_.subdraws = []; options_.unit_root_vars = []; options_.use_mh_covariance_matrix = 0; + options_.gradient_method = 2; % Misc options_.conf_sig = 0.6; diff --git a/matlab/numgrad.m b/matlab/numgrad.m index 03457df2f..e6076017b 100644 --- a/matlab/numgrad.m +++ b/matlab/numgrad.m @@ -1,7 +1,7 @@ function [g, badg] = numgrad(fcn,x,varargin) % function [g badg] = numgrad(fcn,xvarargin) % -delta = 1e-6; +delta = 1e-6; %delta=1e-2; n=length(x); tvec=delta*eye(n); @@ -15,29 +15,40 @@ g=zeros(n,1); %end %f0 = eval([fcn '(x' tailstr]); % Is there a way not to do this? %---------------------------------------------------------------^yes -f0 = eval([fcn '(x,varargin{:})']); +[f0,cost_flag] = feval(fcn, x, varargin{:}); +%f0 = eval([fcn '(x,varargin{:})']); % disp(' first fcn in numgrad.m ------------------') %home % disp('numgrad.m is working. ----') % Jiinil on 9/5/95 % sizex=size(x),sizetvec=size(tvec),x, % Jinill on 9/6/95 badg=0; +goog=1;% stepan 07/07/2008 +scale=1; % stepan 07/07/2008 for i=1:n - scale=1; % originally 1 % i,tveci=tvec(:,i)% ,plus=x+scale*tvec(:,i) % Jinill Kim on 9/6/95 if size(x,1)>size(x,2) tvecv=tvec(i,:); else tvecv=tvec(:,i); end - g0 = (eval([fcn '(x+scale*tvecv'', varargin{:})']) - f0) ... - /(scale*delta); + [fh,cost_flag] = feval(fcn, x+scale*transpose(tvecv), varargin{:});% stepan 07/07/2008 + if cost_flag% stepan 07/07/2008 + g0 = (fh - f0) / (scale*delta); + else + [fh,cost_flag] = feval(fcn, x-scale*transpose(tvecv), varargin{:}); + if cost_flag + g0 = (f0-fh) / (scale*delta); + else + goog=0; + end + end % disp(' fcn in the i=1:n loop of numgrad.m ------------------')% Jinill 9/6/95 % disp(' and i is') % Jinill % i % Jinill % fprintf('Gradient w.r.t. %3d: %10g\n',i,g0) %see below Jinill 9/6/95 % -------------------------- special code to essentially quit here % absg0=abs(g0) % Jinill on 9/6/95 - if abs(g0)< 1e15 + if goog && abs(g0)< 1e15 % stepan 07/07/2008 g(i)=g0; % disp('good gradient') % Jinill Kim else @@ -94,5 +105,4 @@ end % end %end %save g.dat g x f0 -%eval(['save g g x f0 ' stailstr]); - +%eval(['save g g x f0 ' stailstr]); \ No newline at end of file diff --git a/matlab/numgrad3.m b/matlab/numgrad3.m new file mode 100644 index 000000000..e74c7eb5e --- /dev/null +++ b/matlab/numgrad3.m @@ -0,0 +1,45 @@ +function [g, badg] = numgrad3(fcn,x,varargin) +% Computes the gradient of the objective function fcn using a three point +% formula if possible. +% +% Adapted from Sims' numgrad routine. +% +% part of DYNARE, copyright Dynare Team (2008) +% Gnu Public License. + +delta = 1e-6; +n=length(x); +tvec=delta*eye(n); +g=zeros(n,1); + +[f0,cost_flag] = feval(fcn, x, varargin{:}); + +badg=0; +goog=1; +scale=1; +for i=1:n + if size(x,1)>size(x,2) + tvecv=tvec(i,:); + else + tvecv=tvec(:,i); + end + [f1,cost_flag1] = feval(fcn, x+scale*transpose(tvecv), varargin{:}); + [f2,cost_flag2] = feval(fcn, x-scale*transpose(tvecv), varargin{:}); + if cost_flag1 && cost_flag2 + g0 = (f1 - f2) / (2*scale*delta); + elseif cost_flag1==1 && cost_flag2==0 + g0 = (f1-f0) / (scale*delta); + elseif cost_flag1==0 && cost_flag2==1 + g0 = (f0-f2) / (scale*delta); + else + goog=0; + end + + if goog && abs(g0)< 1e15 + g(i)=g0; + else + disp('bad gradient ------------------------') + g(i)=0; + badg=1; + end +end \ No newline at end of file diff --git a/matlab/osr1.m b/matlab/osr1.m index ae4ebe79c..d1da93095 100644 --- a/matlab/osr1.m +++ b/matlab/osr1.m @@ -45,7 +45,7 @@ function osr1(i_params,i_var,weights) nit = 1000; verbose = 2; - [f,p]=csminwel('osr_obj',t0,H0,[],crit,nit,i_params,... + [f,p]=csminwel('osr_obj',t0,H0,[],crit,nit,options_.gradient_method,i_params,... inv_order_var(i_var),weights(i_var,i_var)); % options = optimset('fminunc');