gpd.m

极值理论中各种函数及图像的程序。matlab实现。

function res=gpd(data,threshold,nextremes,information),
%Fits a generalized Pareto model to excesses over a high threshold
%
%      USAGE: res=gpd(data,threshold,nextremes,information)
%
%!Either threshold or nextremes should be defined. The undefined one should be entered as '[]'
%
%       data: Data vector
%  threshold: Excesses over this value will be fitted a model
%  nextremes: Implies a threshold value that number of observations remaining above is nextremes
%information: Default is 'observed'. Can be entered as 'expected' also. Determines whether
%             standard errors will be calculated with observed or expected information
%
%        res: Fitted distribution
%
%             res.par_ests: Estimated parameters. 1X2 vector: 
%                                                            1st element: xi
%                                                            2nd element: beta
%               res.funval: Value of the negative log likelihood
%           res.terminated: Termination condition. 1 if successfully terminated
%              res.details: Details of the nonlinear minimization process of the negative
%                           likelihood
%               res.varcov: Variance-covariance matrix of the parameters
%              res.par_ses: Standard deviations of the parameters of the distribution
%                 res.data: Elements that are exceeding the threshold.


warning off
n=length(data);
if (isempty(threshold)&isempty(nextremes))
    disp('Enter either a threshold or the number of upper extremes')
    return
end
if (~isempty(threshold)&~isempty(nextremes))
    disp('Enter Either a threshold or the number of upper extremes')
    return
end

if (~isempty(nextremes))
    threshold=findthresh(data,nextremes);
end
if nargin<4,
    information='observed';
end
    

exceedances=data(data>threshold);
excess=exceedances-threshold;
xbar=mean(excess);
s2=var(excess);
% xi0=-0.5*xbar*(((xbar^2)/s2)-1);
xi0=-0.5*(((xbar^2)/s2)-1); % Correction by Andrea Colombo May 9, 2005
beta0=0.5*xbar*(((xbar^2)/s2)+1);

theta=[xi0,beta0];

opts=optimset('MaxFunEvals',5000,'MaxIter',1000,'TolX',1e-6,'TolFun',1e-6,'Display','off');
xi=theta(1);
beta=theta(2);
cond1 = beta <= 0;
cond2 = ((xi <= 0) & (max(excess) > ( - beta/xi)));
if (cond1 | cond2),
    theta(1)=1;
    theta(2)=1;
end
[res.par_ests,res.funval,res.terminated,res.details] = fminsearch('negloglikgpd',theta,opts,excess);
[res.par_ests,res.funval,res.terminated,res.details] = fminunc('negloglikgpd',res.par_ests,opts,excess);
if strcmp(information,'observed'),
    res.varcov=hessigpd('negloglikgpd',res.par_ests,excess);
    res.par_ses=sqrt(diag(res.varcov))';
elseif strcmp(information,'expected'),
    one = (1 + res.par_ests(1))^2/length(excess);
    two = (2 * (1 + res.par_ests(1)) * res.par_ests(2)^2)/length(excess);
	 cov =  - ((1 + res.par_ests(1)) * res.par_ests(2))/length(excess);
	 res.varcov = [one,cov;cov,two];
    res.par_ses=sqrt(diag(res.varcov))';
else
    
    
    disp('WARNING 4th input should be either observed or expected');
    return
end

res.threshold=threshold;
res.data=exceedances;
res.p_less_thresh=1-length(excess)/n;


warning on