quant.m

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

function out=quant(data,p,models,first,last,ci,opt),
%Creates a plot showing how the estimate of a high quantile in the tail of a dataset
%based on the GPD approximation varies with threshold or number of extremes
%
% USAGE: out=quant(data,p,models,first,last,ci,opt),
%
%  data: Data vector
%     p: Desired probability for quantile estimate (default=0.99)
%models: Number of consecutive gpd models to be fitted (default=30)
% first: Lowest number of exceedances to be considered (default=15)
%  last: Maximum number of exceedances to be considered (default=500). It can also
%        be entered as a percentage of the data such as '%5'.
%
%    ci: Probability for confidence band  (default=0.95)
%   opt: Option for plotting with respect to exceedances or thresholds. If 0 (default)
%        plots with respect to number of exceedances if 1 plots with respect to increasing
%        threshold  
%
%   out: Output matrix
%        1st column: Thresholds
%        2nd column: Estimates
%        3rd column: Number of exceedances
%        4th column: Lower Confidence Band
%        5th column: Upper Confidence Band

data=surecol(data);
if nargin<2,
    p=0.99;models=30;first=15;last=500;ci=0.95;opt=0;
elseif nargin<3,
   if isempty(p),
      p=0.99;
   end
   models=30;first=15;last=500;ci=0.95;opt=0;
elseif nargin<4,
   if isempty(p),
      p=0.99;
   end
   if isempty(models),
      models=30;
   end
   first=15;last=500;ci=0.95;opt=0;
elseif nargin<5,
   if isempty(p),
      p=0.99;
   end
   if isempty(models),
      models=30;
   end
   if isempty(first),
      first=15;
   end
   last=500;ci=0.95;opt=0;
elseif nargin<6
   if isempty(p),
      p=0.99;
   end
   if isempty(models),
      models=30;
   end
   if isempty(first),
      first=15;
   end
   if isempty(last),
      last=500;
   end
   
   ci=0.95;opt=0;
elseif nargin<7,
   if isempty(p),
      p=0.99;
   end
   if isempty(models),
      models=30;
   end
   if isempty(first),
      first=15;
   end
   if isempty(last),
      last=500;
   end
   if isempty(ci),
      ci=0.95;
   end
   
   opt=0;
elseif nargin==7,
    opt=opt;
else
    disp('Wrong number of inputs');
    return
end

n=length(data);

if nargin>4,
   if ~isempty(last),
      if isstr(last)==1,
			perc=str2num(last(2:end))/100;
			last=floor(perc*n);
   	else
      	last=last;
   	end
   end
end

qq=norminv(1-(1-ci)/2);
exceed=fix(linspace(min(last,n),first,models));
if p<1-min(exceed/n),
    disp(sprintf('\nGraph may look strange \n'));
    disp(['Suggestion 1: Increase ''p'' above ' , num2str(1-min(exceed)/n)]);
    disp(['Suggestion 2: Increase ''first'' above ' , num2str(ceil(length(data)*(1-p)))]);
end

mat=[];
for i=1:length(exceed),
     mat=[mat;gpd_dummy(exceed(i),data)];
end
mat=mat';
thresh=mat(1,:);
xihat=mat(2,:);
betahat=mat(3,:);
lambda=length(data)./exceed;
a=lambda*(1-p);
qest=thresh+betahat.*gfunc(a,xihat);
l=qest;u=qest;
yrange=[min(qest) max(qest)];
xivar=mat(4,:);betavar=mat(5,:);covar=mat(6,:);scaling=thresh;
betahat=betahat./scaling;
covar=covar./scaling;
betavar=betavar./scaling.^2;
term1=betavar.*(gfunc(a,xihat)).^2;
term2=xivar.*(betahat.^2).*(gfunc_deriv(a,xihat)).^2;
term3=2*covar.*betavar.*gfunc(a,xihat).*gfunc_deriv(a,xihat);
qvar=term1+term2+term3;
if(min(qvar) < 0)
    disp(sprintf('\nConditioning problems lead to estimated negative \nquantile variance'));
    return
end
qse=scaling.*sqrt(qvar);
u=qest+qse*qq;
l=qest-qse*qq;
out=[thresh;qest;exceed;l;u]';

if opt==1,
    index=thresh;
    
elseif opt==0,
    index=exceed;
else
    disp('opt should be 0 or 1');
end


plot(index,qest);


hold on
plot(index,u,'r:');
plot(index,l,'r:');
hold off

ylabel([num2str(p) ' Quantile (CI = ' num2str(ci) ')']);
if opt==1,
    xlabel('Threshold');
end
if opt==0;
    xlabel('Exceedances')
end

if p<1-min(exceed/n),
    disp(sprintf('\nGraph may look strange \n'));
    disp(['Suggestion 1: Increase ''p'' above ' , num2str(1-min(exceed)/n)]);
    disp(['Suggestion 2: Increase ''first'' above ' , num2str(ceil(length(data)*(1-p)))]);
end


function c=gfunc_deriv(a,xihat)
c=(-(a.^(-xihat)-1)./xihat-a.^(-xihat).*log(a))./xihat;

function c=gfunc(a, xihat),
c=(a.^( - xihat) - 1)./xihat;

function c=gpd_dummy(nex,data);
out=gpd(data,[],nex,'expected');
c=[out.threshold,out.par_ests(1),out.par_ests(2),out.varcov(1,1),out.varcov(2,2),out.varcov(1,2)];