Commit 71700dff authored by Stéphane Adjemian's avatar Stéphane Adjemian
Browse files

Merge pull request #942 from JohannesPfeifer/TaRB_integration

Integrate the TaRB-algorithm into Dynare
parents 46263910 f0167d8c
......@@ -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}
......@@ -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:
......@@ -5074,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
......@@ -5093,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
......@@ -5143,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
......@@ -5162,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
......@@ -5184,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
......@@ -5206,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
......@@ -5250,6 +5282,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. @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. @xref{optim}.
@item moments_varendo
@anchor{moments_varendo} Triggers the computation of the posterior
distribution of the theoretical moments of the endogenous
......@@ -6666,6 +6724,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}
......@@ -9366,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
......@@ -12552,6 +12613,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''
......@@ -12571,12 +12643,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
......@@ -12716,6 +12782,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,''
......
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 <http://www.gnu.org/licenses/>.
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
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 <http://www.gnu.org/licenses/>.
par_vector(parameterindices,:)=optpar; %reassemble parameter
[fval,DLIK,Hess,exit_flag] = feval(TargetFun,par_vector,varargin{:}); %call target function
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 <http://www.gnu.org/licenses/>.
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