shape2.asv

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

function out=shape(data,models,first,last,ci,opt),
%Creates a plot showing how the estimate of shape varies with threshold or number of
%extremes
%
% USAGE: out=shape(data,models,first,last,ci,opt)
%
%  data: Data vector
%  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: Index (Number of Exceedances)
%        2nd column: Index (Threshold values)
%        3rd column: Corresponding Parmater Estimates
%        4th column: Lower Confidence Band
%        5th column: Upper Confidence Band
%
%
%   Modified: November 6, 2004
%   A bug (found by Andrea Colombo of bankaIntesa) in the output from the threshold option is fixed
%   
%  
%
data=surecol(data);
warning off
if nargin<2,
    models=30;first=15;last=500;ci=0.95;opt=0;
elseif nargin<3,
   if isempty(models),
      models=30;
   end
   first=15;last=500;ci=0.95;opt=0;
elseif nargin<4,
   if isempty(models),
      models=30;
   end
   if isempty(first),
      first=15;
   end
   last=500;ci=0.95;opt=0;
elseif nargin<5,
   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<6,
   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==6,
    opt=opt;
else
    disp('Wrong number of inputs');
    return
end

qq=0;
if ~isempty(ci),
    qq=norminv(1-(1-ci)/2);
end
n=length(data);
if nargin>3,
   if ~isempty(last),
      if isstr(last)==1,
			perc=str2num(last(2:end))/100;
			last=floor(perc*n);
   	else
      	last=last;
   	end
   end
end

x=fix(linspace(min(last,n),first,models));
mat=[];

for i=1:length(x),
    
    mat=[mat;gpd_dummy(x(i),data)];
    
end

mat=[mat,x']';
thresh=mat(1,:);
y=mat(2,:);
yrange=[min(y) max(y)];

if ~isempty(ci),
    u=y+mat(3,:)*qq;
    l=y-mat(3,:)*qq;
end

if opt==0,
index=x;
plot(index,y);
hold on
plot(index,u,'r:');
plot(index,l,'r:');
hold off
xlabel('Exceedances');
ylabel(['Shape(xi) (CI = ' num2str(ci) ')']);
elseif opt==1,
    index=x;
    index2=thresh;
    plot(index2,y);
    hold on
    plot(index2,u,'r:');
    plot(index2,l,'r:');
    hold off
    xlabel('Thresholds');
    ylabel(['Shape(xi) (CI = ' num2str(ci) ')']);
else
    disp('opt should be 0 or 1');
    return
end

warning on

out=[index',thresh',y',l',u'];

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