From 0769fd1cce9f8564cc746aaf7cbc812c50e18a73 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Wed, 13 May 2015 20:58:47 +0200 Subject: [PATCH 1/8] Add verbosity and SaveFile options to some optimizers Also moves additional optimizer-related files to corresponding subfolder Closes #894 --- matlab/global_initialization.m | 6 + matlab/{ => optimization}/bfgsi1.m | 20 +-- matlab/optimization/csminit1.m | 19 ++- matlab/optimization/csminwel1.m | 69 +++++----- .../optimization/dynare_minimize_objective.m | 51 +++++-- matlab/optimization/mr_gstep.m | 15 ++- matlab/{ => optimization}/newrat.m | 125 ++++++++++-------- .../simplex_optimization_routine.m | 4 +- matlab/utilities/general/disp_verbose.m | 25 ++++ tests/optimizers/optimizer_function_wrapper.m | 5 +- 10 files changed, 217 insertions(+), 122 deletions(-) rename matlab/{ => optimization}/bfgsi1.m (69%) rename matlab/{ => optimization}/newrat.m (65%) create mode 100644 matlab/utilities/general/disp_verbose.m diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index b98417452..bb7b1f5a3 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -487,6 +487,9 @@ options_.homotopy_force_continue = 0; %csminwel optimization routine csminwel.tolerance.f=1e-7; csminwel.maxiter=1000; +csminwel.verbosity=1; +csminwel.Save_files=1; + options_.csminwel=csminwel; %newrat optimization routine @@ -494,6 +497,9 @@ newrat.hess=1; % dynare numerical hessian newrat.tolerance.f=1e-5; newrat.tolerance.f_analytic=1e-7; newrat.maxiter=1000; +newrat.verbosity=1; +newrat.Save_files=1; + options_.newrat=newrat; % Simplex optimization routine (variation on Nelder Mead algorithm). diff --git a/matlab/bfgsi1.m b/matlab/optimization/bfgsi1.m similarity index 69% rename from matlab/bfgsi1.m rename to matlab/optimization/bfgsi1.m index c3b87ca49..696ee410b 100644 --- a/matlab/bfgsi1.m +++ b/matlab/optimization/bfgsi1.m @@ -1,13 +1,15 @@ -function H = bfgsi1(H0,dg,dx) -% H = bfgsi1(H0,dg,dx) +function H = bfgsi1(H0,dg,dx,Verbose,Save_files) +% H = bfgsi1(H0,dg,dx,Verbose,Save_files) % Update Inverse Hessian matrix % % Inputs: % H0 [npar by npar] initial inverse Hessian matrix % dg [npar by 1] previous change in gradient % dx [npar by 1] previous change in x; +% Verbose [scalar] Indicator for silent mode +% Save_files [scalar] Indicator whether to save files % -% 6/8/93 version that updates inverse hessian instead of hessian +% 6/8/93 version that updates inverse Hessian instead of Hessian % itself. % % Original file downloaded from: @@ -42,10 +44,12 @@ dgdx = dg'*dx; if (abs(dgdx) >1e-12) H = H0 + (1+(dg'*Hdg)/dgdx)*(dx*dx')/dgdx - (dx*Hdg'+Hdg*dx')/dgdx; else - disp('bfgs update failed.') - disp(['|dg| = ' num2str(sqrt(dg'*dg)) '|dx| = ' num2str(sqrt(dx'*dx))]); - disp(['dg''*dx = ' num2str(dgdx)]) - disp(['|H*dg| = ' num2str(Hdg'*Hdg)]) + disp_verbose('bfgs update failed.',Verbose) + disp_verbose(['|dg| = ' num2str(sqrt(dg'*dg)) '|dx| = ' num2str(sqrt(dx'*dx))],Verbose); + disp_verbose(['dg''*dx = ' num2str(dgdx)],Verbose) + disp_verbose(['|H*dg| = ' num2str(Hdg'*Hdg)],Verbose) H=H0; end -save('H.dat','H') +if Save_files + save('H.dat','H') +end diff --git a/matlab/optimization/csminit1.m b/matlab/optimization/csminit1.m index b6574b005..428c704c6 100644 --- a/matlab/optimization/csminit1.m +++ b/matlab/optimization/csminit1.m @@ -1,4 +1,4 @@ -function [fhat,xhat,fcount,retcode] = csminit1(fcn,x0,f0,g0,badg,H0,varargin) +function [fhat,xhat,fcount,retcode] = csminit1(fcn,x0,f0,g0,badg,H0,Verbose,varargin) % [fhat,xhat,fcount,retcode] = csminit1(fcn,x0,f0,g0,badg,H0,varargin) % % Inputs: @@ -101,7 +101,7 @@ else % toc dxnorm = norm(dx); if dxnorm > 1e12 - disp('Near-singular H problem.') + disp_verbose('Near-singular H problem.',Verbose) dx = dx*FCHANGE/dxnorm; end dfhat = dx'*g0; @@ -115,10 +115,10 @@ else dx = dx - (ANGLE*dxnorm/gnorm+dfhat/(gnorm*gnorm))*g; dfhat = dx'*g; dxnorm = norm(dx); - disp(sprintf('Correct for low angle: %g',a)) + disp_verbose(sprintf('Correct for low angle: %g',a),Verbose) end end - disp(sprintf('Predicted improvement: %18.9f',-dfhat/2)) + disp_verbose(sprintf('Predicted improvement: %18.9f',-dfhat/2),Verbose) % % Have OK dx, now adjust length of step (lambda) until min and % max improvement rate criteria are met. @@ -141,7 +141,7 @@ else %ARGLIST %f = feval(fcn,dxtest,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13); % f = feval(fcn,x0+dx*lambda,P1,P2,P3,P4,P5,P6,P7,P8); - disp(sprintf('lambda = %10.5g; f = %20.7f',lambda,f )) + disp_verbose(sprintf('lambda = %10.5g; f = %20.7f',lambda,f ),Verbose) %debug %disp(sprintf('Improvement too great? f0-f: %g, criterion: %g',f0-f,-(1-THETA)*dfhat*lambda)) if f 0) && (f0 <= fhat) % try going against gradient, which may be inaccurate - lambda = -lambda*factor^6 + if Verbose + lambda = -lambda*factor^6 + else + lambda = -lambda*factor^6; + end else if lambda < 0 retcode = 6; @@ -222,4 +226,5 @@ else end end end -disp(sprintf('Norm of dx %10.5g', dxnorm)) + +disp_verbose(sprintf('Norm of dx %10.5g', dxnorm),Verbose) diff --git a/matlab/optimization/csminwel1.m b/matlab/optimization/csminwel1.m index dcc579f50..5c641ba9e 100644 --- a/matlab/optimization/csminwel1.m +++ b/matlab/optimization/csminwel1.m @@ -1,4 +1,4 @@ -function [fh,xh,gh,H,itct,fcount,retcodeh] = csminwel1(fcn,x0,H0,grad,crit,nit,method,epsilon,varargin) +function [fh,xh,gh,H,itct,fcount,retcodeh] = csminwel1(fcn,x0,H0,grad,crit,nit,method,epsilon,Verbose,Save_files,varargin) %[fhat,xhat,ghat,Hhat,itct,fcount,retcodeh] = csminwel1(fcn,x0,H0,grad,crit,nit,method,epsilon,varargin) % Inputs: % fcn: [string] string naming the objective function to be minimized @@ -65,7 +65,6 @@ fh = []; xh = []; [nx,no]=size(x0); nx=max(nx,no); -Verbose=1; NumGrad= isempty(grad); done=0; itct=0; @@ -87,7 +86,7 @@ retcodeh = []; [f0,junk1,junk2,cost_flag] = feval(fcn,x0,varargin{:}); if ~cost_flag - disp('Bad initial parameter.') + disp_verbose('Bad initial parameter.',Verbose) return end @@ -110,15 +109,13 @@ while ~done g1=[]; g2=[]; g3=[]; %addition fj. 7/6/94 for control - if Verbose - disp('-----------------') - disp(sprintf('f at the beginning of new iteration, %20.10f',f)) - end + disp_verbose('-----------------',Verbose) + disp_verbose(sprintf('f at the beginning of new iteration, %20.10f',f),Verbose) %-----------Comment out this line if the x vector is long---------------- - % disp([sprintf('x = ') sprintf('%15.8g %15.8g %15.8g %15.8g\n',x)]); + % disp_verbose([sprintf('x = ') sprintf('%15.8g %15.8g %15.8g %15.8g\n',x)]); %------------------------- itct=itct+1; - [f1, x1, fc, retcode1] = csminit1(fcn,x,f,g,badg,H,varargin{:}); + [f1, x1, fc, retcode1] = csminit1(fcn,x,f,g,badg,H,Verbose,varargin{:}); fcount = fcount+fc; % erased on 8/4/94 % if (retcode == 1) || (abs(f1-f) < crit) @@ -142,7 +139,9 @@ while ~done end wall1=badg1; % g1 - save g1.mat g1 x1 f1 varargin; + if Save_files + save g1.mat g1 x1 f1 varargin; + end end if wall1 % && (~done) by Jinill % Bad gradient or back and forth on step length. Possibly at @@ -150,10 +149,8 @@ while ~done % %fcliff=fh;xcliff=xh; Hcliff=H+diag(diag(H).*rand(nx,1)); - if Verbose - disp('Cliff. Perturbing search direction.') - end - [f2, x2, fc, retcode2] = csminit1(fcn,x,f,g,badg,Hcliff,varargin{:}); + disp_verbose('Cliff. Perturbing search direction.',Verbose) + [f2, x2, fc, retcode2] = csminit1(fcn,x,f,g,badg,Hcliff,Verbose,varargin{:}); fcount = fcount+fc; % put by Jinill if f2 < f if retcode2==2 || retcode2==4 @@ -169,11 +166,15 @@ while ~done end wall2=badg2; % g2 - badg2 - save g2.mat g2 x2 f2 varargin + if Verbose + badg2 + end + if Save_files + save g2.mat g2 x2 f2 varargin + end end if wall2 - disp('Cliff again. Try traversing') + disp_verbose('Cliff again. Try traversing',Verbose) if norm(x2-x1) < 1e-13 f3=f; x3=x; badg3=1;retcode3=101; else @@ -181,7 +182,7 @@ while ~done if(size(x0,2)>1) gcliff=gcliff'; end - [f3, x3, fc, retcode3] = csminit1(fcn,x,f,gcliff,0,eye(nx),varargin{:}); + [f3, x3, fc, retcode3] = csminit1(fcn,x,f,gcliff,0,eye(nx),Verbose,varargin{:}); fcount = fcount+fc; % put by Jinill if retcode3==2 || retcode3==4 wall3=1; @@ -197,7 +198,9 @@ while ~done end wall3=badg3; % g3 - save g3.mat g3 x3 f3 varargin; + if Save_files + save g3.mat g3 x3 f3 varargin; + end end end else @@ -225,7 +228,7 @@ while ~done fh=f1;xh=x1;gh=g1;badgh=badg1;retcodeh=retcode1; else [fh,ih] = min([f1,f2,f3]); - %disp(sprintf('ih = %d',ih)) + %disp_verbose(sprintf('ih = %d',ih)) %eval(['xh=x' num2str(ih) ';']) switch ih case 1 @@ -259,18 +262,16 @@ while ~done %end of picking stuck = (abs(fh-f) < crit); if (~badg) && (~badgh) && (~stuck) - H = bfgsi1(H,gh-g,xh-x); - end - if Verbose - disp('----') - disp(sprintf('Improvement on iteration %d = %18.9f',itct,f-fh)) + H = bfgsi1(H,gh-g,xh-x,Verbose,Save_files); end + disp_verbose('----',Verbose) + disp_verbose(sprintf('Improvement on iteration %d = %18.9f',itct,f-fh),Verbose) % if Verbose if itct > nit - disp('iteration count termination') + disp_verbose('iteration count termination',Verbose) done = 1; elseif stuck - disp('improvement < crit termination') + disp_verbose('improvement < crit termination',Verbose) done = 1; end rc=retcodeh; @@ -278,19 +279,19 @@ while ~done if rc ==0 %do nothing, just a normal step elseif rc == 1 - disp('zero gradient') + disp_verbose('zero gradient',Verbose) elseif rc == 6 - disp('smallest step still improving too slow, reversed gradient') + disp_verbose('smallest step still improving too slow, reversed gradient',Verbose) elseif rc == 5 - disp('largest step still improving too fast') + disp_verbose('largest step still improving too fast',Verbose) elseif (rc == 4) || (rc==2) - disp('back and forth on step length never finished') + disp_verbose('back and forth on step length never finished',Verbose) elseif rc == 3 - disp('smallest step still improving too slow') + disp_verbose('smallest step still improving too slow',Verbose) elseif rc == 7 - disp('warning: possible inaccuracy in H matrix') + disp_verbose('warning: possible inaccuracy in H matrix',Verbose) else - error('Unaccounted Case, please contact the developers') + error('Unaccounted Case, please contact the developers',Verbose) end end diff --git a/matlab/optimization/dynare_minimize_objective.m b/matlab/optimization/dynare_minimize_objective.m index bf2ed53ad..368003421 100644 --- a/matlab/optimization/dynare_minimize_objective.m +++ b/matlab/optimization/dynare_minimize_objective.m @@ -112,13 +112,15 @@ switch minimizer_algorithm end npar=length(start_par_value); [LB, UB]=set_bounds_to_finite_values(bounds, options_.huge_number); - fprintf('\nNumber of parameters= %d, initial temperatur= %4.3f \n', npar,sa_options.initial_temperature); - fprintf('rt= %4.3f; TolFun= %4.3f; ns= %4.3f;\n',sa_options.rt,sa_options.TolFun,sa_options.ns); - fprintf('nt= %4.3f; neps= %4.3f; MaxIter= %d\n',sa_options.nt,sa_options.neps,sa_options.MaxIter); - fprintf('Initial step length(vm): %4.3f; step_length_c: %4.3f\n', sa_options.initial_step_length,sa_options.step_length_c); - fprintf('%-20s %-6s %-6s %-6s\n','Name:', 'LB;','Start;','UB;'); - for pariter=1:npar - fprintf('%-20s %6.4f; %6.4f; %6.4f;\n',parameter_names{pariter}, LB(pariter),start_par_value(pariter),UB(pariter)); + if sa_options.verbosity + fprintf('\nNumber of parameters= %d, initial temperatur= %4.3f \n', npar,sa_options.initial_temperature); + fprintf('rt= %4.3f; TolFun= %4.3f; ns= %4.3f;\n',sa_options.rt,sa_options.TolFun,sa_options.ns); + fprintf('nt= %4.3f; neps= %4.3f; MaxIter= %d\n',sa_options.nt,sa_options.neps,sa_options.MaxIter); + fprintf('Initial step length(vm): %4.3f; step_length_c: %4.3f\n', sa_options.initial_step_length,sa_options.step_length_c); + fprintf('%-20s %-6s %-6s %-6s\n','Name:', 'LB;','Start;','UB;'); + for pariter=1:npar + fprintf('%-20s %6.4f; %6.4f; %6.4f;\n',parameter_names{pariter}, LB(pariter),start_par_value(pariter),UB(pariter)); + end end sa_options.initial_step_length= sa_options.initial_step_length*ones(npar,1); %bring step length to correct vector size sa_options.step_length_c= sa_options.step_length_c*ones(npar,1); %bring step_length_c to correct vector size @@ -152,6 +154,8 @@ switch minimizer_algorithm nit = options_.csminwel.maxiter; numgrad = options_.gradient_method; epsilon = options_.gradient_epsilon; + Verbose = options_.csminwel.verbosity; + Save_files = options_.csminwel.Save_files; % Change some options. if ~isempty(options_.optim_opt) options_list = read_key_value_string(options_.optim_opt); @@ -167,6 +171,10 @@ switch minimizer_algorithm numgrad = options_list{i,2}; case 'NumgradEpsilon' epsilon = options_list{i,2}; + case 'verbosity' + Verbose = options_list{i,2}; + case 'SaveFiles' + Save_files = options_list{i,2}; otherwise warning(['csminwel: Unknown option (' options_list{i,1} ')!']) end @@ -180,7 +188,7 @@ switch minimizer_algorithm end % Call csminwell. [fval,opt_par_values,grad,inverse_hessian_mat,itct,fcount,exitflag] = ... - csminwel1(objective_function, start_par_value, H0, analytic_grad, crit, nit, numgrad, epsilon, varargin{:}); + csminwel1(objective_function, start_par_value, H0, analytic_grad, crit, nit, numgrad, epsilon, Verbose, Save_files, varargin{:}); hessian_mat=inv(inverse_hessian_mat); case 5 if options_.analytic_derivation==-1 %set outside as code for use of analytic derivation @@ -193,6 +201,8 @@ switch minimizer_algorithm newratflag = options_.newrat.hess; %default end nit=options_.newrat.maxiter; + Verbose = options_.newrat.verbosity; + Save_files = options_.newrat.Save_files; if ~isempty(options_.optim_opt) options_list = read_key_value_string(options_.optim_opt); for i=1:rows(options_list) @@ -208,12 +218,16 @@ switch minimizer_algorithm end case 'TolFun' crit = options_list{i,2}; + case 'verbosity' + Verbose = options_list{i,2}; + case 'SaveFiles' + Save_files = options_list{i,2}; otherwise warning(['newrat: Unknown option (' options_list{i,1} ')!']) end end end - [opt_par_values,hessian_mat,gg,fval,invhess] = newrat(objective_function,start_par_value,analytic_grad,crit,nit,0,varargin{:}); + [opt_par_values,hessian_mat,gg,fval,invhess] = newrat(objective_function,start_par_value,analytic_grad,crit,nit,0,Verbose, Save_files,varargin{:}); %hessian_mat is the plain outer product gradient Hessian case 6 [opt_par_values, hessian_mat, Scale, fval] = gmhmaxlik(objective_function, start_par_value, ... @@ -255,6 +269,8 @@ switch minimizer_algorithm simplexOptions.maxfcallfactor = options_list{i,2}; case 'InitialSimplexSize' simplexOptions.delta_factor = options_list{i,2}; + case 'verbosity' + simplexOptions.verbose = options_list{i,2}; otherwise warning(['simplex: Unknown option (' options_list{i,1} ')!']) end @@ -278,6 +294,17 @@ switch minimizer_algorithm cmaesOptions.TolX = options_list{i,2}; case 'MaxFunEvals' cmaesOptions.MaxFunEvals = options_list{i,2}; + case 'verbosity' + if options_list{i,2}==0 + cmaesOptions.DispFinal = 'off'; % display messages like initial and final message'; + cmaesOptions.DispModulo = '0'; % [0:Inf], disp messages after every i-th iteration'; + end + case 'SaveFiles' + if options_list{i,2}==0 + cmaesOptions.SaveVariables='off'; + cmaesOptions.LogModulo = '0'; % [0:Inf] if >1 record data less frequently after gen=100'; + cmaesOptions.LogTime = '0'; % [0:100] max. percentage of time for recording data'; + end otherwise warning(['cmaes: Unknown option (' options_list{i,1} ')!']) end @@ -309,6 +336,12 @@ switch minimizer_algorithm simpsaOptions.TEMP_END = options_list{i,2}; case 'MaxFunEvals' simpsaOptions.MAX_FUN_EVALS = options_list{i,2}; + case 'verbosity' + if options_list{i,2} == 0 + simpsaOptions.DISPLAY = 'none'; + else + simpsaOptions.DISPLAY = 'iter'; + end otherwise warning(['simpsa: Unknown option (' options_list{i,1} ')!']) end diff --git a/matlab/optimization/mr_gstep.m b/matlab/optimization/mr_gstep.m index 82c8d70f4..51bfb09c4 100644 --- a/matlab/optimization/mr_gstep.m +++ b/matlab/optimization/mr_gstep.m @@ -1,4 +1,4 @@ -function [f0, x, ig] = mr_gstep(h1,x,func0,htol0,varargin) +function [f0, x, ig] = mr_gstep(h1,x,func0,htol0,Verbose,Save_files,varargin) % function [f0, x, ig] = mr_gstep(h1,x,func0,htol0,varargin) % % Gibbs type step in optimisation @@ -69,14 +69,19 @@ while i htol - [f0 x fc retcode] = csminit1(func0,x,f0,gg,0,diag(hh),varargin{:}); + [f0 x fc retcode] = csminit1(func0,x,f0,gg,0,diag(hh),Verbose,varargin{:}); ig(i)=1; - fprintf(['Done for param %s = %8.4f\n'],varargin{6}.name{i},x(i)) + if Verbose + fprintf(['Done for param %s = %8.4f\n'],varargin{6}.name{i},x(i)) + end end xh1=x; end + if Save_files + save gstep.mat x h1 f0 + end +end +if Save_files save gstep.mat x h1 f0 end -save gstep.mat x h1 f0 - diff --git a/matlab/newrat.m b/matlab/optimization/newrat.m similarity index 65% rename from matlab/newrat.m rename to matlab/optimization/newrat.m index 424db0c1c..52676db9b 100644 --- a/matlab/newrat.m +++ b/matlab/optimization/newrat.m @@ -1,4 +1,4 @@ -function [xparam1, hh, gg, fval, igg] = newrat(func0, x, analytic_derivation, ftol0, nit, flagg, varargin) +function [xparam1, hh, gg, fval, igg] = newrat(func0, x, analytic_derivation, ftol0, nit, flagg, Verbose, Save_files, varargin) % [xparam1, hh, gg, fval, igg] = newrat(func0, x, hh, gg, igg, ftol0, nit, flagg, varargin) % % Optimiser with outer product gradient and with sequences of univariate steps @@ -49,6 +49,7 @@ function [xparam1, hh, gg, fval, igg] = newrat(func0, x, analytic_derivation, ft global objective_function_penalty_base + icount=0; nx=length(x); xparam1=x; @@ -95,14 +96,19 @@ else h1=[]; end H = igg; -disp(['Gradient norm ',num2str(norm(gg))]) +disp_verbose(['Gradient norm ',num2str(norm(gg))],Verbose) ee=eig(hh); -disp(['Minimum Hessian eigenvalue ',num2str(min(ee))]) -disp(['Maximum Hessian eigenvalue ',num2str(max(ee))]) +disp_verbose(['Minimum Hessian eigenvalue ',num2str(min(ee))],Verbose) +disp_verbose(['Maximum Hessian eigenvalue ',num2str(max(ee))],Verbose) g=gg; check=0; -if max(eig(hh))<0, disp('Negative definite Hessian! Local maximum!'), pause, end, -save m1.mat x hh g hhg igg fval0 +if max(eig(hh))<0 + disp_verbose('Negative definite Hessian! Local maximum!',Verbose) + pause +end +if Save_files + save m1.mat x hh g hhg igg fval0 +end igrad=1; igibbs=1; @@ -116,13 +122,13 @@ while norm(gg)>gtol && check==0 && jit1 - disp('Gradient step!!') + disp_verbose('Gradient step!!',Verbose) else igrad=0; end @@ -139,27 +145,27 @@ while norm(gg)>gtol && check==0 && jit=ftol - disp('Diagonal Hessian successful') + disp_verbose('Diagonal Hessian successful',Verbose) end fval=fval2; end if (fval0(icount)-fval)=ftol - disp('Gradient direction successful') + disp_verbose('Gradient direction successful',Verbose) end fval=fval3; end @@ -167,7 +173,7 @@ while norm(gg)>gtol && check==0 && jitgtol && check==0 && jit1.e-12 && analytic_derivation==0, try - save m1.mat x fval0 nig -append + if Save_files + save m1.mat x fval0 nig -append + end catch - save m1.mat x fval0 nig + if Save_files + save m1.mat x fval0 nig + end end [dum, gg, htol0, igg, hhg, h1]=mr_hessian(0,xparam1,func0,flagit,htol,varargin{:}); if isempty(dum), @@ -220,12 +230,12 @@ while norm(gg)>gtol && check==0 && jithtol htol=htol0; skipline() - disp('Numerical noise in the likelihood') - disp('Tolerance has to be relaxed') + disp_verbose('Numerical noise in the likelihood',Verbose) + disp_verbose('Tolerance has to be relaxed',Verbose) skipline() end if ~outer_product_gradient, - H = bfgsi1(H,gg-g(:,icount),xparam1-x(:,icount)); + H = bfgsi1(H,gg-g(:,icount),xparam1-x(:,icount),Verbose,Save_files); hh=inv(H); hhg=hh; else @@ -246,39 +256,44 @@ while norm(gg)>gtol && check==0 && jitftol0 skipline() - disp('Numerical noise in the likelihood') - disp('Tolerance had to be relaxed') + disp_verbose('Numerical noise in the likelihood',Verbose) + disp_verbose('Tolerance had to be relaxed',Verbose) skipline() end if jit==nit skipline() - disp('Maximum number of iterations reached') + disp_verbose('Maximum number of iterations reached',Verbose) skipline() end if norm(gg)<=gtol - disp(['Estimation ended:']) - disp(['Gradient norm < ', num2str(gtol)]) + disp_verbose(['Estimation ended:'],Verbose) + disp_verbose(['Gradient norm < ', num2str(gtol)],Verbose) end if check==1, - disp(['Estimation successful.']) + disp_verbose(['Estimation successful.'],Verbose) end -return \ No newline at end of file +return diff --git a/matlab/optimization/simplex_optimization_routine.m b/matlab/optimization/simplex_optimization_routine.m index 33843cf0d..58b3ff1b0 100644 --- a/matlab/optimization/simplex_optimization_routine.m +++ b/matlab/optimization/simplex_optimization_routine.m @@ -507,12 +507,12 @@ fval = fv(1); exitflag = 1; if func_count>= max_func_calls - disp('Maximum number of objective function calls has been exceeded!') + disp_verbose('Maximum number of objective function calls has been exceeded!',verbose) exitflag = 0; end if iter_count>= max_iterations - disp('Maximum number of iterations has been exceeded!') + disp_verbose('Maximum number of iterations has been exceeded!',verbose) exitflag = 0; end diff --git a/matlab/utilities/general/disp_verbose.m b/matlab/utilities/general/disp_verbose.m new file mode 100644 index 000000000..09c09e4d7 --- /dev/null +++ b/matlab/utilities/general/disp_verbose.m @@ -0,0 +1,25 @@ +function disp_verbose(input_string,Verbose) +% function disp_verbose(input_string,Verbose) +% Prints input_string unless Verbose=0 is requested +% +% Copyright (C) 2015 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 Verbose + disp(input_string) +end \ No newline at end of file diff --git a/tests/optimizers/optimizer_function_wrapper.m b/tests/optimizers/optimizer_function_wrapper.m index a7992bab3..262f3d86d 100644 --- a/tests/optimizers/optimizer_function_wrapper.m +++ b/tests/optimizers/optimizer_function_wrapper.m @@ -9,8 +9,9 @@ crit = 1e-7; numgrad = 2; epsilon = 1e-6; analytic_grad=[]; - +Verbose=1; +Save_files=1; %call optimizer [fval,opt_par_values,grad,hessian_mat,itct,fcount,exitflag] = ... - csminwel1(objective_function_handle, start_par_value, H0, analytic_grad, crit, nit, numgrad, epsilon, varargin{:}); + csminwel1(objective_function_handle, start_par_value, H0, analytic_grad, crit, nit, numgrad, epsilon, Verbose,Save_files, varargin{:}); end \ No newline at end of file From 3bf13dd53b51ad6fb7b47f0cb83def5fed02ead7 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Sat, 2 May 2015 14:21:13 +0200 Subject: [PATCH 2/8] Add functionality for TaRB --- matlab/TaRB_metropolis_hastings_core.m | 293 +++++++++++++++++ matlab/TaRB_optimizer_wrapper.m | 40 +++ matlab/chol_SE.m | 309 ++++++++++++++++++ matlab/global_initialization.m | 7 + .../optimization/dynare_minimize_objective.m | 39 +++ matlab/random_walk_metropolis_hastings.m | 17 +- 6 files changed, 702 insertions(+), 3 deletions(-) create mode 100644 matlab/TaRB_metropolis_hastings_core.m create mode 100644 matlab/TaRB_optimizer_wrapper.m create mode 100644 matlab/chol_SE.m diff --git a/matlab/TaRB_metropolis_hastings_core.m b/matlab/TaRB_metropolis_hastings_core.m new file mode 100644 index 000000000..28716c427 --- /dev/null +++ b/matlab/TaRB_metropolis_hastings_core.m @@ -0,0 +1,293 @@ +function myoutput = TaRB_metropolis_hastings_core(myinputs,fblck,nblck,whoiam, ThisMatlab) +% function myoutput = TaRB_metropolis_hastings_core(myinputs,fblck,nblck,whoiam, ThisMatlab) +% Contains the most computationally intensive portion of code in +% random_walk_metropolis_hastings (the 'for xxx = fblck:nblck' loop) using the TaRB algorithm. +% The branches in that 'for' +% cycle are completely independent to be suitable for parallel execution. +% +% INPUTS +% o myimput [struc] The mandatory variables for local/remote +% parallel computing obtained from random_walk_metropolis_hastings.m +% function. +% o fblck and nblck [integer] The Metropolis-Hastings chains. +% o whoiam [integer] In concurrent programming a modality to refer to the different threads running in parallel is needed. +% The integer whoaim is the integer that +% allows us to distinguish between them. Then it is the index number of this CPU among all CPUs in the +% cluster. +% o ThisMatlab [integer] Allows us to distinguish between the +% 'main' Matlab, the slave Matlab worker, local Matlab, remote Matlab, +% ... Then it is the index number of this slave machine in the cluster. +% OUTPUTS +% o myoutput [struc] +% If executed without parallel, this is the original output of 'for b = +% fblck:nblck'. Otherwise, it's a portion of it computed on a specific core or +% remote machine. In this case: +% record; +% irun; +% NewFile; +% OutputFileName +% +% ALGORITHM +% Portion of Tailored Randomized Block Metropolis-Hastings proposed in +% Chib/Ramamurthy (2010): Tailored randomized block MCMC methods with +% application to DSGE models, Journal of Econometrics 155, pp. 19-38 +% +% This implementation differs from the originally proposed one in the +% treatment of non-positive definite Hessians. Here we +% - use the Jordan decomposition +% +% SPECIAL REQUIREMENTS. +% None. +% +% PARALLEL CONTEXT +% See the comments in the random_walk_metropolis_hastings.m funtion. + + +% Copyright (C) 2006-2015 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 . +global objective_function_penalty_base; + +if nargin<4, + whoiam=0; +end + +% reshape 'myinputs' for local computation. +% In order to avoid confusion in the name space, the instruction struct2local(myinputs) is replaced by: + +TargetFun=myinputs.TargetFun; +ProposalFun=myinputs.ProposalFun; +xparam1=myinputs.xparam1; +mh_bounds=myinputs.mh_bounds; +last_draw=myinputs.ix2; +last_posterior=myinputs.ilogpo2; +fline=myinputs.fline; +npar=myinputs.npar; +nruns=myinputs.nruns; +NewFile=myinputs.NewFile; +MAX_nruns=myinputs.MAX_nruns; +d=myinputs.d; +InitSizeArray=myinputs.InitSizeArray; +record=myinputs.record; +dataset_ = myinputs.dataset_; +dataset_info = myinputs.dataset_info; +bayestopt_ = myinputs.bayestopt_; +estim_params_ = myinputs.estim_params_; +options_ = myinputs.options_; +M_ = myinputs.M_; +oo_ = myinputs.oo_; + +% Necessary only for remote computing! +if whoiam + % initialize persistent variables in priordens() + priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7, bayestopt_.p3,bayestopt_.p4,1); +end + +MetropolisFolder = CheckPath('metropolis',M_.dname); +ModelName = M_.fname; +BaseName = [MetropolisFolder filesep ModelName]; + +options_.lik_algo = 1; +OpenOldFile = ones(nblck,1); + +% +% Now I run the (nblck-fblck+1) Metropolis-Hastings chains +% + +block_iter=0; + +for curr_chain = fblck:nblck, + block_iter=block_iter+1; + try + % This will not work if the master uses a random number generator not + % available in the slave (different Matlab version or + % Matlab/Octave cluster). Therefore the trap. + % + % Set the random number generator type (the seed is useless but needed by the function) + set_dynare_seed(options_.DynareRandomStreams.algo, options_.DynareRandomStreams.seed); + % Set the state of the RNG + set_dynare_random_generator_state(record.InitialSeeds(curr_chain).Unifor, record.InitialSeeds(curr_chain).Normal); + catch + % If the state set by master is incompatible with the slave, we only reseed + set_dynare_seed(options_.DynareRandomStreams.seed+curr_chain); + end + if (options_.load_mh_file~=0) && (fline(curr_chain)>1) && OpenOldFile(curr_chain) %load previous draws and likelihood + load([BaseName '_mh' int2str(NewFile(curr_chain)) '_blck' int2str(curr_chain) '.mat']) + x2 = [x2;zeros(InitSizeArray(curr_chain)-fline(curr_chain)+1,npar)]; + logpo2 = [logpo2;zeros(InitSizeArray(curr_chain)-fline(curr_chain)+1,1)]; + OpenOldFile(curr_chain) = 0; + else + x2 = zeros(InitSizeArray(curr_chain),npar); + logpo2 = zeros(InitSizeArray(curr_chain),1); + end + %Prepare waiting bars + if whoiam + prc0=(curr_chain-fblck)/(nblck-fblck+1)*(isoctave || options_.console_mode); + hh = dyn_waitbar({prc0,whoiam,options_.parallel(ThisMatlab)},['MH (' int2str(curr_chain) '/' int2str(options_.mh_nblck) ')...']); + else + hh = dyn_waitbar(0,['Metropolis-Hastings (' int2str(curr_chain) '/' int2str(options_.mh_nblck) ')...']); + set(hh,'Name','Metropolis-Hastings'); + end + accepted_draws_this_chain = 0; + accepted_draws_this_file = 0; + blocked_draws_counter_this_chain=0; + blocked_draws_counter_this_chain_this_file=0; + draw_index_current_file = fline(curr_chain); %get location of first draw in current block + draw_iter = 1; + while draw_iter <= nruns(curr_chain) + + %% randomize indices for blocking in this iteration + indices=randperm(npar)'; + blocks=[1; (1+cumsum((rand(length(indices)-1,1)>(1-options_.TaRB.new_block_probability))))]; + nblocks=blocks(end,1); %get number of blocks this iteration + current_draw=last_draw(curr_chain,:)'; %get starting point for current draw for updating + + for block_iter=1:nblocks + blocked_draws_counter_this_chain=blocked_draws_counter_this_chain+1; + blocked_draws_counter_this_chain_this_file=blocked_draws_counter_this_chain_this_file+1; + nxopt=length(indices(blocks==block_iter,1)); %get size of current block + par_start_current_block=current_draw(indices(blocks==block_iter,1)); + [xopt_current_block, fval, exitflag, hess_mat_optimizer, options_, Scale] = dynare_minimize_objective(@TaRB_optimizer_wrapper,par_start_current_block,options_.TaRB.mode_compute,options_,[mh_bounds.lb(indices(blocks==block_iter,1),1) mh_bounds.ub(indices(blocks==block_iter,1),1)],bayestopt_.name,bayestopt_,[],... + current_draw,indices(blocks==block_iter,1),TargetFun,...% inputs for wrapper + dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_); %inputs for objective + objective_function_penalty_base=Inf; %reset penalty that may have been changed by optimizer + %% covariance for proposal density + hessian_mat = reshape(hessian('TaRB_optimizer_wrapper',xopt_current_block, ... + options_.gstep,... + current_draw,indices(blocks==block_iter,1),TargetFun,...% inputs for wrapper + dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_),nxopt,nxopt); + + if any(any(isnan(hessian_mat))) || any(any(isinf(hessian_mat))) + inverse_hessian_mat=eye(nxopt)*1e-4; %use diagonal + else + inverse_hessian_mat=inv(hessian_mat); %get inverse Hessian + if any(any((isnan(inverse_hessian_mat)))) || any(any((isinf(inverse_hessian_mat)))) + inverse_hessian_mat=eye(nxopt)*1e-4; %use diagonal + end + end + [proposal_covariance_Cholesky_decomposition_upper,negeigenvalues]=cholcov(inverse_hessian_mat,0); + %if not positive definite, use generalized Cholesky if + %Eskow/Schnabel + if negeigenvalues~=0 + proposal_covariance_Cholesky_decomposition_upper=chol_SE(inverse_hessian_mat,0); + end + proposal_covariance_Cholesky_decomposition_upper=proposal_covariance_Cholesky_decomposition_upper*diag(bayestopt_.jscale(indices(blocks==block_iter,1),:)); + %get proposal draw + if strcmpi(ProposalFun,'rand_multivariate_normal') + n = nxopt; + elseif strcmpi(ProposalFun,'rand_multivariate_student') + n = options_.student_degrees_of_freedom; + end + + proposed_par = feval(ProposalFun, xopt_current_block', proposal_covariance_Cholesky_decomposition_upper, n); + % chech whether draw is valid and compute posterior + if all( proposed_par(:) > mh_bounds.lb(indices(blocks==block_iter,1),:) ) && all( proposed_par(:) < mh_bounds.ub(indices(blocks==block_iter,1),:) ) + try + logpost = - feval('TaRB_optimizer_wrapper', proposed_par(:),... + current_draw,indices(blocks==block_iter,1),TargetFun,...% inputs for wrapper + dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_); + catch + logpost = -inf; + end + else + logpost = -inf; + end + %get ratio of proposal densities, required because proposal depends + %on current mode via Hessian and is thus not symmetric anymore + if strcmpi(ProposalFun,'rand_multivariate_normal') + proposal_density_proposed_move_forward=multivariate_normal_pdf(proposed_par,xopt_current_block',proposal_covariance_Cholesky_decomposition_upper,n); + proposal_density_proposed_move_backward=multivariate_normal_pdf(par_start_current_block',xopt_current_block',proposal_covariance_Cholesky_decomposition_upper,n); + elseif strcmpi(ProposalFun,'rand_multivariate_student') + proposal_density_proposed_move_forward=multivariate_student_pdf(proposed_par,xopt_current_block',proposal_covariance_Cholesky_decomposition_upper,n); + proposal_density_proposed_move_backward=multivariate_student_pdf(par_start_current_block',xopt_current_block',proposal_covariance_Cholesky_decomposition_upper,n); + end + accprob=logpost-last_posterior(curr_chain)+ log(proposal_density_proposed_move_backward)-log(proposal_density_proposed_move_forward); %Formula (6), Chib/Ramamurthy + + if (logpost > -inf) && (log(rand) < accprob) + current_draw(indices(blocks==block_iter,1))=proposed_par; + last_posterior(curr_chain)=logpost; + accepted_draws_this_chain =accepted_draws_this_chain +1; + accepted_draws_this_file = accepted_draws_this_file + 1; + else %no updating + %do nothing, keep old value + end + end + %save draws and update stored last values + x2(draw_index_current_file,:) = current_draw; + last_draw(curr_chain,:) = current_draw; + %save posterior after full run through all blocks + logpo2(draw_index_current_file) = last_posterior(curr_chain); + + prtfrc = draw_iter/nruns(curr_chain); + if (mod(draw_iter, 3)==0 && ~whoiam) || (mod(draw_iter,50)==0 && whoiam) + dyn_waitbar(prtfrc,hh,[ 'MH (' int2str(curr_chain) '/' int2str(options_.mh_nblck) ') ' sprintf('Current acceptance ratio %4.3f', accepted_draws_this_chain/blocked_draws_counter_this_chain)]); + end + if (draw_index_current_file == InitSizeArray(curr_chain)) || (draw_iter == nruns(curr_chain)) % Now I save the simulations, either because the current file is full or the chain is done + save([BaseName '_mh' int2str(NewFile(curr_chain)) '_blck' int2str(curr_chain) '.mat'],'x2','logpo2'); + fidlog = fopen([MetropolisFolder '/metropolis.log'],'a'); + fprintf(fidlog,['\n']); + fprintf(fidlog,['%% Mh' int2str(NewFile(curr_chain)) 'Blck' int2str(curr_chain) ' (' datestr(now,0) ')\n']); + fprintf(fidlog,' \n'); + fprintf(fidlog,[' Number of simulations.: ' int2str(length(logpo2)) '\n']); + fprintf(fidlog,[' Acceptance ratio......: ' num2str(accepted_draws_this_file /blocked_draws_counter_this_chain_this_file) '\n']); + fprintf(fidlog,[' Posterior mean........:\n']); + for i=1:length(x2(1,:)) + fprintf(fidlog,[' params:' int2str(i) ': ' num2str(mean(x2(:,i))) '\n']); + end + fprintf(fidlog,[' log2po:' num2str(mean(logpo2)) '\n']); + fprintf(fidlog,[' Minimum value.........:\n']); + for i=1:length(x2(1,:)) + fprintf(fidlog,[' params:' int2str(i) ': ' num2str(min(x2(:,i))) '\n']); + end + fprintf(fidlog,[' log2po:' num2str(min(logpo2)) '\n']); + fprintf(fidlog,[' Maximum value.........:\n']); + for i=1:length(x2(1,:)) + fprintf(fidlog,[' params:' int2str(i) ': ' num2str(max(x2(:,i))) '\n']); + end + fprintf(fidlog,[' log2po:' num2str(max(logpo2)) '\n']); + fprintf(fidlog,' \n'); + fclose(fidlog); + %reset counters; + accepted_draws_this_file = 0; + blocked_draws_counter_this_chain_this_file=0; + if draw_iter == nruns(curr_chain) % I record the last draw... + record.LastParameters(curr_chain,:) = x2(end,:); + record.LastLogPost(curr_chain) = logpo2(end); + end + % size of next file in chain curr_chain + InitSizeArray(curr_chain) = min(nruns(curr_chain)-draw_iter,MAX_nruns); + % initialization of next file if necessary + if InitSizeArray(curr_chain) + x2 = zeros(InitSizeArray(curr_chain),npar); + logpo2 = zeros(InitSizeArray(curr_chain),1); + NewFile(curr_chain) = NewFile(curr_chain) + 1; + draw_index_current_file = 0; + end + end + draw_iter=draw_iter+1; + draw_index_current_file = draw_index_current_file + 1; + end% End of the simulations for one mh-block. + record.AcceptanceRatio(curr_chain) = accepted_draws_this_chain/blocked_draws_counter_this_chain; + dyn_waitbar_close(hh); + [record.LastSeeds(curr_chain).Unifor, record.LastSeeds(curr_chain).Normal] = get_dynare_random_generator_state(); + OutputFileName(block_iter,:) = {[MetropolisFolder,filesep], [ModelName '_mh*_blck' int2str(curr_chain) '.mat']}; +end% End of the loop over the mh-blocks. + + +myoutput.record = record; +myoutput.irun = draw_index_current_file; +myoutput.NewFile = NewFile; +myoutput.OutputFileName = OutputFileName; \ No newline at end of file diff --git a/matlab/TaRB_optimizer_wrapper.m b/matlab/TaRB_optimizer_wrapper.m new file mode 100644 index 000000000..9c7ed4ccc --- /dev/null +++ b/matlab/TaRB_optimizer_wrapper.m @@ -0,0 +1,40 @@ +function [fval,DLIK,Hess,exit_flag] = TaRB_optimizer_wrapper(optpar,par_vector,parameterindices,TargetFun,varargin) +% function [fval,DLIK,Hess,exit_flag] = TaRB_optimizer_wrapper(optpar,par_vector,parameterindices,TargetFun,varargin) +% Wrapper function for target function used in TaRB algorithm; reassembles +% full parameter vector before calling target function +% +% INPUTS +% o optpar [double] (p_opt*1) vector of subset of parameters to be considered +% o par_vector [double] (p*1) full vector of parameters +% o parameterindices [double] (p_opt*1) index of optpar entries in +% par_vector +% o TargetFun [char] string specifying the name of the objective +% function (posterior kernel). +% o varargin [structure] other inputs of target function +% +% OUTPUTS +% o fval [scalar] value of (minus) the likelihood. +% o DLIK [double] (p*1) score vector of the likelihood. +% o Hess [double] (p*p) asymptotic Hessian matrix. +% o exit_flag [scalar] equal to zero if the routine return with a penalty (one otherwise). +% +% Copyright (C) 2015 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 . + +par_vector(parameterindices,:)=optpar; %reassemble parameter +[fval,DLIK,Hess,exit_flag] = feval(TargetFun,par_vector,varargin{:}); %call target function + diff --git a/matlab/chol_SE.m b/matlab/chol_SE.m new file mode 100644 index 000000000..0ec5a900b --- /dev/null +++ b/matlab/chol_SE.m @@ -0,0 +1,309 @@ +function [R,indef, E, P]=chol_SE(A,pivoting) +% [R,indef, E, P]=chol_SE(A,pivoting) +% Performs a (modified) Cholesky factorization of the form +% +% P'*A*P + E = R'*R +% +% As detailed in Schnabel/Eskow (1990), the factorization has 2 phases: +% Phase 1: While A can still be positive definite, pivot on the maximum diagonal element. +% Check that the standard Cholesky update would result in a positive diagonal +% at the current iteration. If so, continue with the normal Cholesky update. +% Otherwise switch to phase 2. +% If A is safely positive definite, stage 1 is never left, resulting in +% the standard Cholesky decomposition. +% +% Phase 2: Pivot on the minimum of the negatives of the lower Gershgorin bound +% estimates. To prevent negative diagonals, compute the amount to add to the +% pivot element and add this. Then, do the Cholesky update and update the estimates of the +% Gershgorin bounds. +% +% Notes: +% - During factorization, L=R' is stored in the lower triangle of the original matrix A, +% miminizing the memory requirements +% - Conforming with the original Schnabel/Eskow (1990) algorithm: +% - at each iteration the updated Gershgorin bounds are estimated instead of recomputed, +% reducing the computational requirements from o(n^3) to o (n^2) +% - For the last 2 by 2 submatrix, an eigenvalue-based decomposition is used +% - While pivoting is not necessary, it improves the size of E, the add-on on to the diagonal. But this comes at +% the cost of introduding a permutation. +% +% +% Inputs +% A [n*n] Matrix to be factorized +% pivoting [scalar] dummy whether pivoting is used +% +% Outputs +% R [n*n] originally stored in lower triangular portion of A, including the main diagonal +% +% E [n*1] Elements added to the diagonal of A +% P [n*1] record of how the rows and columns of the matrix were permuted while +% performing the decomposition +% +% REFERENCES: +% This implementation is based on +% +% Robert B. Schnabel and Elizabeth Eskow. 1990. "A New Modified Cholesky +% Factorization," SIAM Journal of Scientific Statistical Computating, +% 11, 6: 1136-58. +% +% Elizabeth Eskow and Robert B. Schnabel 1991. "Algorithm 695 - Software for a New Modified Cholesky +% Factorization," ACM Transactions on Mathematical Software, Vol 17, No 3: 306-312 +% +% +% Author: Johannes Pfeifer based on Eskow/Schnabel (1991) +% +% Copyright (C) 2015 Johannes Pfeifer +% Copyright (C) 2015 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 sum(sum(abs(A-A'))) > 0 + error('A is not symmetric') +end + +if nargin==1 + pivoting=0; +end + + +n=size(A,1); + +tau1=eps^(1/3); %tolerance parameter for determining when to switch phase 2 +tau2=eps^(1/3); %tolerance used for determining the maximum condition number of the final 2X2 submatrix. + +phase1 = 1; +delta = 0; + +P=1:n; +g=zeros(n,1); +E=zeros(n,1); + +% Find the maximum magnitude of the diagonal elements. If any diagonal element is negative, then phase1 is false. +gammma=max(diag(A)); +if any(diag(A)) < 0 + phase1 = 0; +end + +taugam = tau1*gammma; +% If not in phase1, then calculate the initial Gershgorin bounds needed for the start of phase2. +if ~phase1 + g=gersh_nested(A,1,n); +end + +% check for n=1 +if n==1 + delta = tau2*abs(A(1,1)) - A(1,1); + if delta > 0 + E(1) = delta; + end + if A(1,1) == 0 + E(1) = tau2; + end + A(1,1)=sqrt(A(1,1)+E(1)); +end + +for j = 1:n-1 + % PHASE 1 + if phase1 + if pivoting==1 + % Find index of maximum diagonal element A(i,i) where i>=j + [tmp,imaxd] = max(diag(A(j:n,j:n))); + imaxd=imaxd+j-1; + % Pivot to the top the row and column with the max diag + if (imaxd ~= j) + % Swap row j with row of max diag + for ii = 1:j-1 + temp = A(j,ii); + A(j,ii) = A(imaxd,ii); + A(imaxd,ii) = temp; + end + % Swap colj and row maxdiag between j and maxdiag + for ii = j+1:imaxd-1 + temp = A(ii,j); + A(ii,j) = A(imaxd,ii); + A(imaxd,ii) = temp; + end + % Swap column j with column of max diag + for ii = imaxd+1:n + temp = A(ii,j); + A(ii,j) = A(ii,imaxd); + A(ii,imaxd) = temp; + end + % Swap diag elements + temp = A(j,j); + A(j,j) = A(imaxd,imaxd); + A(imaxd,imaxd) = temp; + % Swap elements of the permutation vector + itemp = P(j); + P(j) = P(imaxd); + P(imaxd) = itemp; + end + end + % check to see whether the normal Cholesky update for this + % iteration would result in a positive diagonal, + % and if not then switch to phase 2. + jp1 = j+1; + if A(j,j)>0 + if min((diag(A(j+1:n,j+1:n)) - A(j+1:n,j).^2/A(j,j))) < taugam %test whether stage 2 is required + phase1=0; + end + else + phase1 = 0; + end + if phase1 + % Do the normal cholesky update if still in phase 1 + A(j,j) = sqrt(A(j,j)); + tempjj = A(j,j); + for ii = jp1: n + A(ii,j) = A(ii,j)/tempjj; + end + for ii=jp1:n + temp=A(ii,j); + for k = jp1:ii + A(ii,k) = A(ii,k)-(temp * A(k,j)); + end + end + if j == n-1 + A(n,n)=sqrt(A(n,n)); + end + else + % Calculate the negatives of the lower Gershgorin bounds + g=gersh_nested(A,j,n); + end + end + + % PHASE 2 + if ~phase1 + if j ~= n-1 + if pivoting + % Find the minimum negative Gershgorin bound + [tmp,iming] = min(g(j:n)); + iming=iming+j-1; + % Pivot to the top the row and column with the + % minimum negative Gershgorin bound + if iming ~= j + % Swap row j with row of min Gershgorin bound + for ii = 1:j-1 + temp = A(j,ii); + A(j,ii) = A(iming,ii); + A(iming,ii) = temp; + end + + % Swap col j with row iming from j to iming + for ii = j+1:iming-1 + temp = A(ii,j); + A(ii,j) = A(iming,ii); + A(iming,ii) = temp; + end + % Swap column j with column of min Gershgorin bound + for ii = iming+1:n + temp = A(ii,j); + A(ii,j) = A(ii,iming); + A(ii,iming) = temp; + end + % Swap diagonal elements + temp = A(j,j); + A(j,j) = A(iming,iming); + A(iming,iming) = temp; + % Swap elements of the permutation vector + itemp = P(j); + P(j) = P(iming); + P(iming) = itemp; + % Swap elements of the negative Gershgorin bounds vector + temp = g(j); + g(j) = g(iming); + g(iming) = temp; + end + end + % Calculate delta and add to the diagonal. delta=max{0,-A(j,j) + max{normj,taugam},delta_previous} + % where normj=sum of |A(i,j)|,for i=1,n, delta_previous is the delta computed at the previous iter and taugam is tau1*gamma. + + normj=sum(abs(A(j+1:n,j))); + + delta = max([0;delta;-A(j,j)+normj;-A(j,j)+taugam]); % get adjustment based on formula on bottom of p. 309 of Eskow/Schnabel (1991) + + E(j) = delta; + A(j,j) = A(j,j) + E(j); + % Update the Gershgorin bound estimates (note: g(i) is the negative of the Gershgorin lower bound.) + if A(j,j) ~= normj + temp = (normj/A(j,j)) - 1; + for ii = j+1:n + g(ii) = g(ii) + abs(A(ii,j)) * temp; + end + end + % Do the cholesky update + A(j,j) = sqrt(A(j,j)); + tempjj = A(j,j); + for ii = j+1:n + A(ii,j) = A(ii,j) / tempjj; + end + for ii = j+1:n + temp = A(ii,j); + for k = j+1: ii + A(ii,k) = A(ii,k) - (temp * A(k,j)); + end + end + else + % Find eigenvalues of final 2 by 2 submatrix + % Find delta such that: + % 1. the l2 condition number of the final 2X2 submatrix + delta*I <= tau2 + % 2. delta >= previous delta, + % 3. min(eigvals) + delta >= tau2 * gamma, where min(eigvals) is the smallest eigenvalue of the final 2X2 submatrix + + A(n-1,n)=A(n,n-1); %set value above diagonal for computation of eigenvalues + eigvals = eig(A(n-1:n,n-1:n)); + delta = max([ 0 ; delta ; -min(eigvals)+tau2*max((max(eigvals)-min(eigvals))/(1-tau1),gammma) ]); %Formula 5.3.2 of Schnabel/Eskow (1990) + + if delta > 0 + A(n-1,n-1) = A(n-1,n-1) + delta; + A(n,n) = A(n,n) + delta; + E(n-1) = delta; + E(n) = delta; + end + + % Final update + A(n-1,n-1) = sqrt(A(n-1,n-1)); + A(n,n-1) = A(n,n-1)/A(n-1,n-1); + A(n,n) = A(n,n) - (A(n,n-1)^2); + A(n,n) = sqrt(A(n,n)); + end + end +end + +R=(tril(A))'; +indef=~phase1; +Pprod=zeros(n,n); +Pprod(sub2ind([n,n],P,1:n))=1; +P=Pprod; +end + +function g=gersh_nested(A,j,n) + +g=zeros(n,1); +for ii = j:n + if ii == 1; + sum_up_to_i = 0; + else + sum_up_to_i = sum(abs(A(ii,j:(ii-1)))); + end; + if ii == n; + sum_after_i = 0; + else + sum_after_i = sum(abs(A((ii+1):n,ii))); + end; + g(ii) = sum_up_to_i + sum_after_i- A(ii,ii); +end +end + diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index bb7b1f5a3..70303d3aa 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -416,7 +416,14 @@ options_.recursive_estimation_restart = 0; options_.MCMC_jumping_covariance='hessian'; options_.use_calibration_initialization = 0; options_.endo_vars_for_moment_computations_in_estimation=[]; + +% Tailored Random Block Metropolis-Hastings options_.TaRB.use_TaRB = 0; +options_.TaRB.mode_compute=4; +options_.TaRB.new_block_probability=0.25; %probability that next parameter belongs to new block + +% Run optimizer silently +options_.silent_optimizer=0; % Prior restrictions options_.prior_restrictions.status = 0; diff --git a/matlab/optimization/dynare_minimize_objective.m b/matlab/optimization/dynare_minimize_objective.m index 368003421..f570a1a3e 100644 --- a/matlab/optimization/dynare_minimize_objective.m +++ b/matlab/optimization/dynare_minimize_objective.m @@ -73,6 +73,9 @@ switch minimizer_algorithm if ~isempty(options_.optim_opt) eval(['optim_options = optimset(optim_options,' options_.optim_opt ');']); end + if options_.silent_optimizer + optim_options = optimset(optim_options,'display','off'); + end if options_.analytic_derivation, optim_options = optimset(optim_options,'GradObj','on','TolX',1e-7); end @@ -110,6 +113,9 @@ switch minimizer_algorithm end end end + if options_.silent_optimizer + sa_options.verbosity = 0; + end npar=length(start_par_value); [LB, UB]=set_bounds_to_finite_values(bounds, options_.huge_number); if sa_options.verbosity @@ -140,6 +146,9 @@ switch minimizer_algorithm if options_.analytic_derivation, optim_options = optimset(optim_options,'GradObj','on'); end + if options_.silent_optimizer + optim_options = optimset(optim_options,'display','off'); + end if ~isoctave [opt_par_values,fval,exitflag] = fminunc(objective_function,start_par_value,optim_options,varargin{:}); else @@ -180,6 +189,10 @@ switch minimizer_algorithm end end end + if options_.silent_optimizer + Save_files = 0; + Verbose = 0; + end % Set flag for analytical gradient. if options_.analytic_derivation analytic_grad=1; @@ -227,6 +240,10 @@ switch minimizer_algorithm end end end + if options_.silent_optimizer + Save_files = 0; + Verbose = 0; + end [opt_par_values,hessian_mat,gg,fval,invhess] = newrat(objective_function,start_par_value,analytic_grad,crit,nit,0,Verbose, Save_files,varargin{:}); %hessian_mat is the plain outer product gradient Hessian case 6 @@ -243,6 +260,9 @@ switch minimizer_algorithm if ~isempty(options_.optim_opt) eval(['optim_options = optimset(optim_options,' options_.optim_opt ');']); end + if options_.silent_optimizer + optim_options = optimset(optim_options,'display','off'); + end if ~isoctave [opt_par_values,fval,exitflag] = fminsearch(objective_function,start_par_value,optim_options,varargin{:}); else @@ -276,6 +296,9 @@ switch minimizer_algorithm end end end + if options_.silent_optimizer + simplexOptions.verbose = options_list{i,2}; + end [opt_par_values,fval,exitflag] = simplex_optimization_routine(objective_function,start_par_value,simplexOptions,parameter_names,varargin{:}); case 9 % Set defaults @@ -310,6 +333,13 @@ switch minimizer_algorithm end end end + if options_.silent_optimizer + cmaesOptions.DispFinal = 'off'; % display messages like initial and final message'; + cmaesOptions.DispModulo = '0'; % [0:Inf], disp messages after every i-th iteration'; + cmaesOptions.SaveVariables='off'; + cmaesOptions.LogModulo = '0'; % [0:Inf] if >1 record data less frequently after gen=100'; + cmaesOptions.LogTime = '0'; % [0:100] max. percentage of time for recording data'; + end warning('off','CMAES:NonfinitenessRange'); warning('off','CMAES:InitialSigma'); [x, fval, COUNTEVAL, STOPFLAG, OUT, BESTEVER] = cmaes(func2str(objective_function),start_par_value,H0,cmaesOptions,varargin{:}); @@ -347,6 +377,9 @@ switch minimizer_algorithm end end end + if options_.silent_optimizer + simpsaOptions.DISPLAY = 'none'; + end simpsaOptionsList = options2cell(simpsaOptions); simpsaOptions = simpsaset(simpsaOptionsList{:}); [LB, UB]=set_bounds_to_finite_values(bounds, options_.huge_number); @@ -377,6 +410,9 @@ switch minimizer_algorithm end end end + if options_.silent_optimizer + solveoptoptions.verbosity = 0; + end [opt_par_values,fval]=solvopt(start_par_value,objective_function,[],[],[],solveoptoptions,varargin{:}); case 102 if isoctave @@ -389,6 +425,9 @@ switch minimizer_algorithm if ~isempty(options_.optim_opt) eval(['optim_options = saoptimset(optim_options,' options_.optim_opt ');']); end + if options_.silent_optimizer + optim_options = optimset(optim_options,'display','off'); + end func = @(x)objective_function(x,varargin{:}); [opt_par_values,fval,exitflag,output] = simulannealbnd(func,start_par_value,bounds(:,1),bounds(:,2),optim_options); otherwise diff --git a/matlab/random_walk_metropolis_hastings.m b/matlab/random_walk_metropolis_hastings.m index 2092922ac..c5ff16c03 100644 --- a/matlab/random_walk_metropolis_hastings.m +++ b/matlab/random_walk_metropolis_hastings.m @@ -87,6 +87,10 @@ load_last_mh_history_file(MetropolisFolder, ModelName); % on many cores). The mandatory variables for local/remote parallel % computing are stored in the localVars struct. +if options_.TaRB.use_TaRB + options_.silent_optimizer=1; %locally set optimizer to silent mode +end + localVars = struct('TargetFun', TargetFun, ... 'ProposalFun', ProposalFun, ... 'xparam1', xparam1, ... @@ -117,9 +121,12 @@ localVars = struct('TargetFun', TargetFun, ... % single chain compute Random walk Metropolis-Hastings algorithm sequentially. if isnumeric(options_.parallel) || (nblck-fblck)==0, - fout = random_walk_metropolis_hastings_core(localVars, fblck, nblck, 0); + if options_.TaRB.use_TaRB + fout = TaRB_metropolis_hastings_core(localVars, fblck, nblck, 0); + else + fout = random_walk_metropolis_hastings_core(localVars, fblck, nblck, 0); + end record = fout.record; - % Parallel in Local or remote machine. else % Global variables for parallel routines. @@ -138,7 +145,11 @@ else end % from where to get back results % NamFileOutput(1,:) = {[M_.dname,'/metropolis/'],'*.*'}; - [fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'random_walk_metropolis_hastings_core', localVars, globalVars, options_.parallel_info); + if options_.TaRB.use_TaRB + [fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'TaRB_metropolis_hastings_core', localVars, globalVars, options_.parallel_info); + else + [fout, nBlockPerCPU, totCPU] = masterParallel(options_.parallel, fblck, nblck,NamFileInput,'random_walk_metropolis_hastings_core', localVars, globalVars, options_.parallel_info); + end for j=1:totCPU, offset = sum(nBlockPerCPU(1:j-1))+fblck-1; record.LastLogPost(offset+1:sum(nBlockPerCPU(1:j)))=fout(j).record.LastLogPost(offset+1:sum(nBlockPerCPU(1:j))); From 55d44f098383e194cc205169f1c20fa820fca605 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Thu, 14 May 2015 08:27:28 +0200 Subject: [PATCH 3/8] Add preprocessor option silent_optimizer --- doc/dynare.texi | 8 ++++++++ preprocessor/DynareBison.yy | 5 ++++- preprocessor/DynareFlex.ll | 1 + 3 files changed, 13 insertions(+), 1 deletion(-) diff --git a/doc/dynare.texi b/doc/dynare.texi index 19d1c4e1d..f0617f7a9 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -4949,6 +4949,11 @@ value of that function as the posterior mode. @noindent Default value is @code{4}. +@item silent_optimizer +@anchor{silent_optimizer} +Instructs Dynare to run mode computing/optimization silently without displaying results or +saving files in between. Useful when running loops. + @item mcmc_jumping_covariance = hessian|prior_variance|identity_matrix|@var{FILENAME} Tells Dynare which covariance to use for the proposal density of the MCMC sampler. @code{mcmc_jumping_covariance} can be one of the following: @@ -6666,6 +6671,9 @@ Specifies the optimizer for minimizing the objective function. The same solvers @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 (i.e. on the value of option @ref{opt_algo}). @xref{optim}. +@item silent_optimizer +@pxref{silent_optimizer} + @item huge_number = @var{DOUBLE} Value for replacing the infinite bounds on parameters by finite numbers. Used by some optimizers for numerical reasons (@pxref{huge_number}). Users need to make sure that the optimal parameters are not larger than this value. Default: @code{1e7} diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy index 808fa6a76..58dfb5618 100644 --- a/preprocessor/DynareBison.yy +++ b/preprocessor/DynareBison.yy @@ -83,7 +83,7 @@ class ParsingDriver; } %token AIM_SOLVER ANALYTIC_DERIVATION AR AUTOCORR TARB_MODE_COMPUTE -%token BAYESIAN_IRF BETA_PDF BLOCK USE_CALIBRATION USE_TARB TARB_NEW_BLOCK_PROBABILITY +%token BAYESIAN_IRF BETA_PDF BLOCK USE_CALIBRATION USE_TARB TARB_NEW_BLOCK_PROBABILITY SILENT_OPTIMIZER %token BVAR_DENSITY BVAR_FORECAST NODECOMPOSITION DR_DISPLAY_TOL HUGE_NUMBER %token BVAR_PRIOR_DECAY BVAR_PRIOR_FLAT BVAR_PRIOR_LAMBDA TARB_OPTIM %token BVAR_PRIOR_MU BVAR_PRIOR_OMEGA BVAR_PRIOR_TAU BVAR_PRIOR_TRAIN @@ -1722,6 +1722,7 @@ estimation_options : o_datafile | o_tarb_mode_compute | o_tarb_new_block_probability | o_tarb_optim + | o_silent_optimizer | o_proposal_distribution | o_student_degrees_of_freedom ; @@ -1791,6 +1792,7 @@ osr_options : stoch_simul_primary_options | o_opt_algo | o_optim | o_huge_number + | o_silent_optimizer ; osr : OSR ';' @@ -2917,6 +2919,7 @@ o_equations : EQUATIONS EQUAL vec_int o_use_tarb : USE_TARB { driver.option_num("TaRB.use_TaRB", "1"); }; o_tarb_mode_compute : TARB_MODE_COMPUTE EQUAL INT_NUMBER { driver.option_num("TaRB.mode_compute", $3); }; o_tarb_new_block_probability : TARB_NEW_BLOCK_PROBABILITY EQUAL non_negative_number {driver.option_num("TaRB.new_block_probability",$3); }; +o_silent_optimizer : SILENT_OPTIMIZER { driver.option_num("silent_optimizer", "1"); }; o_instruments : INSTRUMENTS EQUAL '(' symbol_list ')' {driver.option_symbol_list("instruments"); }; o_ext_func_name : EXT_FUNC_NAME EQUAL filename { driver.external_function_option("name", $3); }; diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll index 43466da55..3687d4c25 100644 --- a/preprocessor/DynareFlex.ll +++ b/preprocessor/DynareFlex.ll @@ -576,6 +576,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 tarb_mode_compute {return token::TARB_MODE_COMPUTE;} tarb_new_block_probability {return token::TARB_NEW_BLOCK_PROBABILITY;} tarb_optim {return token::TARB_OPTIM;} +silent_optimizer {return token::SILENT_OPTIMIZER;} lmmcp {return token::LMMCP;} occbin {return token::OCCBIN;} From aa4b07610ec152ff47e2604b026f57bc2f0eb55d Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Thu, 14 May 2015 14:26:18 +0200 Subject: [PATCH 4/8] Add unit test for TaRB MH-algorithm --- tests/Makefile.am | 1 + tests/estimation/TaRB/fs2000_tarb.mod | 122 ++++++++++++++++++++++++++ 2 files changed, 123 insertions(+) create mode 100644 tests/estimation/TaRB/fs2000_tarb.mod diff --git a/tests/Makefile.am b/tests/Makefile.am index 1eefe3098..a73323ec0 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -4,6 +4,7 @@ MODFILES = \ estimation/fs2000_calibrated_covariance.mod \ estimation/MH_recover/fs2000_recover.mod \ estimation/t_proposal/fs2000_student.mod \ + estimation/TaRB/fs2000_tarb.mod \ gsa/ls2003.mod \ gsa/ls2003a.mod \ gsa/cod_ML_morris/cod_ML_morris.mod \ diff --git a/tests/estimation/TaRB/fs2000_tarb.mod b/tests/estimation/TaRB/fs2000_tarb.mod new file mode 100644 index 000000000..a3941d5a6 --- /dev/null +++ b/tests/estimation/TaRB/fs2000_tarb.mod @@ -0,0 +1,122 @@ +/* + * 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. + * + * The data are in file "fsdat_simul.m", and have been artificially generated. + * They are therefore different from the original dataset used by Schorfheide. + * + * The equations are taken from J. Nason and T. Cogley (1994): "Testing the + * implications of long-run neutrality for monetary business cycle models", + * Journal of Applied Econometrics, 9, S37-S70. + * Note that there is an initial minus sign missing in equation (A1), p. S63. + * + * This implementation was written by Michel Juillard. Please note that the + * following copyright notice only applies to this Dynare implementation of the + * model. + */ + +/* + * Copyright (C) 2004-2010 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 . + */ + +var m P c e W R k d n l gy_obs gp_obs y dA; +varexo e_a e_m; + +parameters alp bet gam mst rho psi del; + +alp = 0.33; +bet = 0.99; +gam = 0.003; +mst = 1.011; +rho = 0.7; +psi = 0.787; +del = 0.02; + +model; +dA = exp(gam+e_a); +log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m; +-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0; +W = l/n; +-(psi/(1-psi))*(c*P/(1-n))+l/n = 0; +R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W; +1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0; +c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1); +P*c = m; +m-1+d = l; +e = exp(e_a); +y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a)); +gy_obs = dA*y/y(-1); +gp_obs = (P/P(-1))*m(-1)/dA; +end; + +shocks; +var e_a; stderr 0.014; +var e_m; stderr 0.005; +end; + +steady_state_model; + dA = exp(gam); + gst = 1/dA; + m = mst; + khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1)); + xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1); + nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp ); + n = xist/(nust+xist); + P = xist + nust; + k = khst*n; + + l = psi*mst*n/( (1-psi)*(1-n) ); + c = mst/P; + d = l - mst + 1; + y = k^alp*n^(1-alp)*gst^alp; + R = mst/bet; + W = l/n; + ist = y-c; + q = 1 - d; + + e = 1; + + gp_obs = m/dA; + gy_obs = dA; +end; + +steady; + +check; + +estimated_params; +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; +psi, beta_pdf, 0.65, 0.05; +del, beta_pdf, 0.01, 0.005; +stderr e_a, inv_gamma_pdf, 0.035449, inf; +stderr e_m, inv_gamma_pdf, 0.008862, inf; +end; + +varobs gp_obs gy_obs; + +options_.silent_optimizer=1; +estimation(order=1, datafile='../fsdat_simul',nobs=192, loglinear, mh_replic=30, mh_nblocks=2, mh_jscale=0.5,mode_compute=4, +use_TaRB, +tarb_mode_compute=4, +tarb_new_block_probability=0.3 +); From c2f07c720f5b9cafd849dd35aacffaff5c9a132c Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Thu, 14 May 2015 14:26:54 +0200 Subject: [PATCH 5/8] Add check whether tarb_new_block_probability is between 0 and 1 --- matlab/initial_estimation_checks.m | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/matlab/initial_estimation_checks.m b/matlab/initial_estimation_checks.m index 48a0b4749..f3b9ce06c 100644 --- a/matlab/initial_estimation_checks.m +++ b/matlab/initial_estimation_checks.m @@ -54,6 +54,10 @@ if DynareOptions.student_degrees_of_freedom <= 0 error('initial_estimation_checks:: the student_degrees_of_freedom takes a positive integer argument'); end +if DynareOptions.TaRB.use_TaRB && (DynareOptions.TaRB.new_block_probability<0 || DynareOptions.TaRB.new_block_probability>1) + error(['initial_estimation_checks:: The tarb_new_block_probability must be between 0 and 1!']) +end + old_steady_params=Model.params; %save initial parameters for check if steady state changes param values % % check if steady state solves static model (except if diffuse_filter == 1) From 35040725bb39ded0bf9cb85e9cc4ed25f17951ee Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Thu, 14 May 2015 14:27:32 +0200 Subject: [PATCH 6/8] Document TaRB algorithm --- doc/dynare.texi | 47 +++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 41 insertions(+), 6 deletions(-) diff --git a/doc/dynare.texi b/doc/dynare.texi index f0617f7a9..1dae09c08 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -5255,6 +5255,32 @@ deprecated and will be removed in a future release of Dynare. @anchor{dsge_varlag} The number of lags used to estimate a DSGE-VAR model. Default: @code{4}. + +@item use_tarb +Instructs Dynare to use the Tailored randomized block (TaRB) Metropolis-Hastings algorithm +proposed by @cite{Chib and Ramamurthy 2010} instead of the standard Random-Walk Metropolis-Hastings. +In this algorithm, at each iteration the estimated parameters are randomly assigned to different +blocks. For each of these blocks a mode-finding step is conducted. The inverse Hessian at this mode +is then used as the covariance of the proposal density for a Random-Walk Metropolis-Hastings step. +If the numerical Hessian is not positive definite, the generalized Cholesky decomposition of +@cite{Schnabel and Eskow 1990} is used, but without pivoting. The TaRB-MH algorithm massively reduces +the autocorrelation in the MH draws and thus reduces the number of draws required to +representatively sample from the posterior. However, this comes at a computational costs as the +algorithm takes more time to run. + +@item tarb_new_block_probability = @var{DOUBLE} +Specifies the probability of the next parameter belonging to a new block when the random blocking in the TaRB +Metropolis-Hastings algorithm is conducted. The higher this number, the smaller is the average block size and the +more random blocks are formed during each parameter sweep. Default: @code{0.25}. + +@item tarb_mode_compute = @var{INTEGER} +Specifies the mode-finder run in every iteration for every block of the +TaRB Metropolis-Hastings algorithm. @pxref{mode_compute}. Default: @code{4}. + +@item tarb_optim = @var{INTEGER} +Specifies the options for the mode-finder used in the TaRB +Metropolis-Hastings algorithm. @pxref{optim}. + @item moments_varendo @anchor{moments_varendo} Triggers the computation of the posterior distribution of the theoretical moments of the endogenous @@ -12560,6 +12586,17 @@ 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 +Chib, Siddhartha and Srikanth Ramamurthy (2010): +``Tailored randomized block MCMC methods with application to DSGE models,'' +@i{Journal of Econometrics}, 155, 19--38 + +@item +Christiano, Lawrence J., Mathias Trabandt and Karl Walentin (2011): +``Introducing financial frictions and unemployment into a small open +economy model,'' @i{Journal of Economic Dynamics and Control}, 35(12), +1999--2041 + @item Collard, Fabrice (2001): ``Stochastic simulations with Dynare: A practical guide'' @@ -12579,12 +12616,6 @@ Corona, Angelo, M. Marchesi, Claudio Martini, and Sandro Ridella (1987): ``Minimizing multimodal functions of continuous variables with the ``simulated annealing'' algorithm'', @i{ACM Transactions on Mathematical Software}, 13(3), 262--280 -@item -Christiano, Lawrence J., Mathias Trabandt and Karl Walentin (2011): -``Introducing financial frictions and unemployment into a small open -economy model,'' @i{Journal of Economic Dynamics and Control}, 35(12), -1999--2041 - @item Del Negro, Marco and Franck Schorfheide (2004): ``Priors from General Equilibrium Models for VARS'', @i{International Economic Review}, 45(2), 643--673 @@ -12724,6 +12755,10 @@ General Equilibrium Models Using a Second-Order Approximation to the Policy Function,'' @i{Journal of Economic Dynamics and Control}, 28(4), 755--775 +@item +Schnabel, Robert B. and Elizabeth Eskow (1990): ``A new modified Cholesky algorithm,'' +@i{SIAM Journal of Scientific and Statistical Computing}, 11, 1136--1158 + @item Sims, Christopher A., Daniel F. Waggoner and Tao Zha (2008): ``Methods for inference in large multiple-equation Markov-switching models,'' From dec37f372a3a0d9eeac6fb8a3537abfad03e2b18 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Thu, 14 May 2015 14:37:41 +0200 Subject: [PATCH 7/8] Document options for silent mode-finding --- doc/dynare.texi | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/doc/dynare.texi b/doc/dynare.texi index 1dae09c08..3f25ae366 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -5079,6 +5079,12 @@ Size of the perturbation used to compute numerically the gradient of the objecti @item 'TolFun' Stopping criteria. Default: @code{1e-7} +@item 'verbosity' +Controls verbosity of display during optimization. Set to 0 to set to silent. Default: @code{1} + +@item 'SaveFiles' +Controls saving of intermediate results during optimization. Set to 0 to shut off saving. Default: @code{1} + @end table @item 5 @@ -5098,6 +5104,12 @@ Maximum number of iterations. Default: @code{1000} @item 'TolFun' Stopping criteria. Default: @code{1e-5} for numerical derivatives @code{1e-7} for analytic derivatives. +@item 'verbosity' +Controls verbosity of display during optimization. Set to 0 to set to silent. Default: @code{1} + +@item 'SaveFiles' +Controls saving of intermediate results during optimization. Set to 0 to shut off saving. Default: @code{1} + @end table @item 6 @@ -5148,6 +5160,9 @@ 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 'verbosity' +Controls verbosity of display during optimization. Set to 0 to set to silent. Default: @code{1} + @end table @item 9 @@ -5167,6 +5182,12 @@ 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} +@item 'verbosity' +Controls verbosity of display during optimization. Set to 0 to set to silent. Default: @code{1} + +@item 'SaveFiles' +Controls saving of intermediate results during optimization. Set to 0 to shut off saving. Default: @code{1} + @end table @item 10 @@ -5189,6 +5210,9 @@ 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 'verbosity' +Controls verbosity of display during optimization. Set to 0 to set to silent. Default: @code{1} + @end table @item 101 @@ -5211,6 +5235,9 @@ Tolerance parameter (w.r.t the objective function). Default: @code{1e-6} @item 'TolX' Tolerance parameter (w.r.t the instruments). Default: @code{1e-6} +@item 'verbosity' +Controls verbosity of display during optimization. Set to 0 to set to silent. Default: @code{1} + @end table @item 102 From f0167d8cd5e70a1a895de55d8e9feeebd9d28894 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Fri, 15 May 2015 15:22:34 +0200 Subject: [PATCH 8/8] Cosmetic fixes to manual --- doc/dynare.texi | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/doc/dynare.texi b/doc/dynare.texi index 3f25ae366..3a2bcaf5e 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -4797,7 +4797,7 @@ Default value is @code{1}. For advanced use only. For internal use and testing only. @item conf_sig = @var{DOUBLE} -Confidence interval used for classical forecasting after estimation. See @xref{conf_sig}. +Confidence interval used for classical forecasting after estimation. @xref{conf_sig}. @item mh_conf_sig = @var{DOUBLE} @anchor{mh_conf_sig} @@ -5285,12 +5285,12 @@ model. Default: @code{4}. @item use_tarb Instructs Dynare to use the Tailored randomized block (TaRB) Metropolis-Hastings algorithm -proposed by @cite{Chib and Ramamurthy 2010} instead of the standard Random-Walk Metropolis-Hastings. +proposed by @cite{Chib and Ramamurthy (2010)} instead of the standard Random-Walk Metropolis-Hastings. In this algorithm, at each iteration the estimated parameters are randomly assigned to different blocks. For each of these blocks a mode-finding step is conducted. The inverse Hessian at this mode is then used as the covariance of the proposal density for a Random-Walk Metropolis-Hastings step. If the numerical Hessian is not positive definite, the generalized Cholesky decomposition of -@cite{Schnabel and Eskow 1990} is used, but without pivoting. The TaRB-MH algorithm massively reduces +@cite{Schnabel and Eskow (1990)} is used, but without pivoting. The TaRB-MH algorithm massively reduces the autocorrelation in the MH draws and thus reduces the number of draws required to representatively sample from the posterior. However, this comes at a computational costs as the algorithm takes more time to run. @@ -5302,11 +5302,11 @@ more random blocks are formed during each parameter sweep. Default: @code{0.25}. @item tarb_mode_compute = @var{INTEGER} Specifies the mode-finder run in every iteration for every block of the -TaRB Metropolis-Hastings algorithm. @pxref{mode_compute}. Default: @code{4}. +TaRB Metropolis-Hastings algorithm. @xref{mode_compute}. Default: @code{4}. @item tarb_optim = @var{INTEGER} Specifies the options for the mode-finder used in the TaRB -Metropolis-Hastings algorithm. @pxref{optim}. +Metropolis-Hastings algorithm. @xref{optim}. @item moments_varendo @anchor{moments_varendo} Triggers the computation of the posterior @@ -9427,7 +9427,7 @@ Note that creating the configuration file is not enough in order to trigger parallelization of the computations: you also need to specify the @code{parallel} option to the @code{dynare} command. For more details, and for other options related to the parallelization engine, -see @pxref{Dynare invocation}. +@pxref{Dynare invocation}. You also need to verify that the following requirements are met by your cluster (which is composed of a master and of one or more