Обобщенное распределение экстремальных значений

Распределение Generalized Extreme Value (GEV) объединяет распределения экстремальных значений типа I, типа II и типа III в одно семейство, чтобы обеспечить непрерывный диапазон возможных форм. Он параметризуется с параметрами местоположения и масштаба, mu и sigma и параметром формы k. Когда k < 0, GEV эквивалентен экстремальному значению типа III. Когда k > 0, GEV эквивалентен типу II. В пределе, когда k приближается к 0, GEV становится типом I.

x = linspace(-3,6,1000);
plot(x,gevpdf(x,-.5,1,0),'-', x,gevpdf(x,0,1,0),'-', x,gevpdf(x,.5,1,0),'-');
xlabel('(x-mu) / sigma');
ylabel('Probability Density');
legend({'k < 0, Type III' 'k = 0, Type I' 'k > 0, Type II'});

Обратите внимание, что для k < 0 или k > 0 плотность имеет нулевую вероятность выше или ниже, соответственно, верхней или нижней границы - (1/k). В пределе, когда k приближается к 0, GEV неограничен. Это можно суммировать как ограничение того, что 1 + k * (y-mu )/sigma должны быть положительными.

Моделирование максимального объема данных блока

GEV может быть конструктивно определен как ограничивающее распределение блоковых максимумов (или минимумов). То есть, если сгенерировать большое количество независимых случайных значений из одного распределения вероятностей и взять их максимальное значение, распределение этого максимума приблизительно равно GEV.

Исходное распределение определяет параметр k формы результирующего распределения GEV. Распределения, хвосты которых выпадают как полином, например t Стьюдента, приводят к положительному параметру формы. Распределения, хвосты которых уменьшаются экспоненциально, например нормаль, соответствуют нулевому параметру формы. Распределения с конечными хвостами, такие как бета, соответствуют отрицательному параметру формы.

Реальные заявки на GEV могут включать моделирование наибольшего дохода для акций в течение каждого месяца. Здесь мы будем моделировать данные, взяв максимум 25 значений из распределения Стьюдента с двумя степенями свободы. Смоделированные данные будут включать 75 максимальных значений случайных блоков.

rng(0,'twister');
y = max(trnd(2,25,75),[],1);

Подбор распределения по максимальной вероятности

Функция gevfit возвращает как оценки параметров максимального правдоподобия, так и (по умолчанию) 95% доверительные интервалы.

[paramEsts,paramCIs] = gevfit(y);

kMLE = paramEsts(1)        % Shape parameter
sigmaMLE = paramEsts(2)    % Scale parameter
muMLE = paramEsts(3)       % Location parameter

kMLE =

    0.4901


sigmaMLE =

    1.4856


muMLE =

    2.9710

kCI = paramCIs(:,1)
sigmaCI = paramCIs(:,2)
muCI = paramCIs(:,3)

kCI =

    0.2020
    0.7782


sigmaCI =

    1.1431
    1.9307


muCI =

    2.5599
    3.3821

Обратите внимание, что 95% доверительный интервал для k не включает нулевое значение. Распределение экстремальных значений типа I, по-видимому, не является хорошей моделью для этих данных. Это имеет смысл, потому что основное распределение для моделирования имело гораздо более тяжелые хвосты, чем нормальное, и распределение экстремальных значений типа II теоретически является правильным, поскольку размер блока становится большим.

В качестве альтернативы доверительным интервалам мы также можем вычислить аппроксимацию к асимптотической ковариационной матрице оценок параметров и извлечь из нее стандартные ошибки параметров.

[nll,acov] = gevlike(paramEsts,y);
paramSEs = sqrt(diag(acov))

paramSEs =

    0.1470
    0.1986
    0.2097

Визуальная проверка соответствия

Чтобы визуально оценить, насколько хорошо подходит, мы рассмотрим графики подогнанной функции плотности вероятности (PDF) и кумулятивной функции распределения (CDF).

Поддержка GEV зависит от значений параметров. В этом случае оценка для k является положительной, поэтому аппроксимированное распределение имеет нулевую вероятность ниже нижней границы.

lowerBnd = muMLE-sigmaMLE./kMLE;

Во-первых, мы построим масштабированную гистограмму данных, наложенную на PDF для установленной модели GEV. Эта гистограмма масштабируется таким образом, что высота полосы умножается на сумму их ширины до 1, чтобы сделать ее сопоставимой с PDF.

ymax = 1.1*max(y);
bins = floor(lowerBnd):ceil(ymax);
h = bar(bins,histc(y,bins)/length(y),'histc');
h.FaceColor = [.9 .9 .9];
ygrid = linspace(lowerBnd,ymax,100);
line(ygrid,gevpdf(ygrid,kMLE,sigmaMLE,muMLE));
xlabel('Block Maximum');
ylabel('Probability Density');
xlim([lowerBnd ymax]);

Мы также можем сравнить соответствие с данными с точки зрения кумулятивной вероятности, наложив эмпирический CDF и встроенный CDF.

[F,yi] = ecdf(y);
plot(ygrid,gevcdf(ygrid,kMLE,sigmaMLE,muMLE),'-');
hold on;
stairs(yi,F,'r');
hold off;
xlabel('Block Maximum');
ylabel('Cumulative Probability');
legend('Fitted Generalized Extreme Value CDF','Empirical CDF','location','southeast');
xlim([lowerBnd ymax]);

Оценка квантилей модели

Хотя оценки параметров могут быть важны сами по себе, квантиль подогнанной модели GEV часто представляет собой количество, представляющее интерес для анализа данных максимумов блоков.

Например, уровень Rm возврата определяется как максимальное значение блока, которое должно быть превышено только один раз в m блоках. Это всего лишь квантиль (1-1/m). Мы можем подключить оценки параметров максимального правдоподобия к обратной CDF для оценки Rm для m = 10.

R10MLE = gevinv(1-1./10,kMLE,sigmaMLE,muMLE)

R10MLE =

    9.0724

Мы могли бы вычислить доверительные пределы для R10, используя асимптотические приближения, но они могут быть недействительными. Вместо этого мы будем использовать метод на основе правдоподобия для вычисления доверительных пределов. Этот способ часто дает более точные результаты, чем один, на основе матрицы оцененной ковариации оценок параметров.

Учитывая любой набор значений для параметров mu, sigma и k, мы можем вычислить логарифмическое правдоподобие - например, MLE являются значениями параметров, которые максимизируют логарифмическое правдоподобие GEV. По мере удаления значений параметров от MLE их логарифмическое правдоподобие обычно становится значительно меньше максимального. Если посмотреть на набор значений параметров, которые дают логарифмическое правдоподобие больше указанного критического значения, это сложная область в пространстве параметров. Однако для подходящего критического значения это доверительная область для параметров модели. Область содержит значения параметров, которые «совместимы с данными». Критическое значение, которое определяет регион, основано на приближении хи-квадрат, и мы будем использовать 95% в качестве нашего доверительного уровня. (Обратите внимание, что мы фактически будем работать с негативом логарифмического правдоподобия.)

nllCritVal = gevlike([kMLE,sigmaMLE,muMLE],y) + .5*chi2inv(.95,1)

nllCritVal =

  170.3044

Для любого набора значений параметров mu, sigma и k можно вычислить R10. Поэтому мы можем найти наименьшее значение R10, достигнутое в критической области пространства параметров, где отрицательное логарифмическое правдоподобие больше критического значения. Это наименьшее значение является более низким доверительным пределом на основе правдоподобия для R10.

Это трудно визуализировать во всех трех измерениях параметров, но как мысленный эксперимент, мы можем зафиксировать параметр формы, k, мы можем видеть, как процедура будет работать над двумя оставшимися параметрами, сигма и mu.

sigmaGrid = linspace(.8, 2.25, 110);
muGrid = linspace(2.4, 3.6);
nllGrid = zeros(length(sigmaGrid),length(muGrid));
R10Grid = zeros(length(sigmaGrid),length(muGrid));
for i = 1:size(nllGrid,1)
    for j = 1:size(nllGrid,2)
        nllGrid(i,j) = gevlike([kMLE,sigmaGrid(i),muGrid(j)],y);
        R10Grid(i,j) = gevinv(1-1./10,kMLE,sigmaGrid(i),muGrid(j));
    end
end
nllGrid(nllGrid>gevlike([kMLE,sigmaMLE,muMLE],y)+6) = NaN;
contour(muGrid,sigmaGrid,R10Grid,6.14:.64:12.14,'LineColor','r');
hold on
contour(muGrid,sigmaGrid,R10Grid,[7.42 11.26],'LineWidth',2,'LineColor','r');
contour(muGrid,sigmaGrid,nllGrid,[168.7 169.1 169.6 170.3:1:173.3],'LineColor','b');
contour(muGrid,sigmaGrid,nllGrid,[nllCritVal nllCritVal],'LineWidth',2,'LineColor','b');
hold off
axis([2.4 3.6 .8 2.25]);
xlabel('mu');
ylabel('sigma');

Синие контуры представляют логарифмическую поверхность правдоподобия, а жирный синий контур является границей критической области. Красные контуры представляют поверхность для R10 - большие значения находятся вверху справа, внизу слева внизу. Контуры являются прямыми, поскольку для фиксированных k Rm является линейной функцией сигма и mu. Полужирные красные контуры - это самые низкие и самые высокие значения R10, которые попадают в критическую область. В полном пространстве трехмерных параметров контуры логарифмического правдоподобия будут эллипсоидальными, а контуры R10 - поверхностями.

Поиск нижнего доверительного предела для R10 является задачей оптимизации с нелинейными ограничениями неравенства, и поэтому мы будем использовать функцию fmincon из Toolbox™ Оптимизация. Нам нужно найти наименьшее значение R10, и поэтому цель, которая должна быть минимизирована, сама R10, равна обратной CDF, оцененной для p = 1-1/m. Мы создадим оберточную функцию, которая вычисляет Rm специально для m = 10.

CIobjfun = @(params) gevinv(1-1./10,params(1),params(2),params(3));

Для выполнения ограниченной оптимизации также потребуется функция, определяющая ограничение, то есть, чтобы отрицательное логарифмическое правдоподобие было меньше критического значения. Функция ограничения должна возвращать положительные значения при нарушении ограничения. Мы создадим анонимную функцию, используя смоделированные данные и критическое логарифмическое значение правдоподобия. Он также возвращает пустое значение, потому что мы не используем никаких ограничений равенства здесь.

CIconfun = @(params) deal(gevlike(params,y) - nllCritVal, []);

Наконец, мы называем fmincon, используя алгоритм активного набора для выполнения ограниченной оптимизации.

opts = optimset('Algorithm','active-set', 'Display','notify', 'MaxFunEvals',500, ...
                'RelLineSrchBnd',.1, 'RelLineSrchBndDuration',Inf);
[params,R10Lower,flag,output] = ...
    fmincon(CIobjfun,paramEsts,[],[],[],[],[],[],CIconfun,opts);

Feasible point with lower objective function value found.

Чтобы найти верхний предел достоверности вероятности для R10, мы просто обратим знак на целевой функции, чтобы найти наибольшее значение R10 в критической области, и вызовем fmincon второй раз.

CIobjfun = @(params) -gevinv(1-1./10,params(1),params(2),params(3));
[params,R10Upper,flag,output] = ...
    fmincon(CIobjfun,paramEsts,[],[],[],[],[],[],CIconfun,opts);
R10Upper = -R10Upper;

R10CI = [R10Lower, R10Upper]

R10CI =

    7.0841   13.4452

plot(ygrid,gevcdf(ygrid,kMLE,sigmaMLE,muMLE),'-');
hold on;
stairs(yi,F,'r');
plot(R10CI([1 1 1 1 2 2 2 2]), [.88 .92 NaN .9 .9 NaN .88 .92],'k-')
hold off;
xlabel('Block Maximum');
ylabel('Cumulative Probability');
legend('Fitted Generalized Extreme Value CDF','Empirical CDF', ...
       'R_{10} 95% CI','location','southeast');
xlim([lowerBnd ymax]);

Профиль правдоподобия для квантиля

Иногда только интервал не дает достаточной информации о оцениваемом количестве, и вместо этого требуется вероятность профиля. Чтобы найти логарифмический профиль правдоподобия для R10, мы зафиксируем возможное значение для R10, а затем максимизируем логарифмическое правдоподобие GEV с параметрами, ограниченными так, чтобы они соответствовали текущему значению R10. Это нелинейное ограничение равенства. Если мы сделаем это в диапазоне R10 значений, мы получим профиль правдоподобия.

Как и в случае доверительного интервала на основе правдоподобия, мы можем подумать о том, какой была бы эта процедура, если бы мы фиксировали k и работали над двумя оставшимися параметрами, сигмой и mu. Каждая красная линия контура на графике контура, показанном ранее, представляет фиксированное значение R10; оптимизация правдоподобия профиля состоит из шага вдоль одной линии контура R10, чтобы найти самый высокий контур логарифмического правдоподобия (синий).

В этом примере мы вычисляем вероятность профиля для R10 по значениям, включенным в доверительный интервал правдоподобия.

R10grid = linspace(R10CI(1)-.05*diff(R10CI), R10CI(2)+.05*diff(R10CI), 51);

Целевой функцией для оптимизации правдоподобия профиля является просто логарифмическое правдоподобие с использованием смоделированных данных.

PLobjfun = @(params) gevlike(params,y);

Использовать fmincon, нам понадобится функция, которая возвращает ненулевые значения при нарушении ограничения, то есть, когда параметры не согласуются с текущим значением R10. Для каждого значения R10 мы создадим анонимную функцию для конкретного значения рассматриваемого R10. Он также возвращает пустое значение, потому что мы не используем никаких ограничений неравенства здесь.

Наконец, мы позвоним fmincon при каждом значении R10 найти соответствующий ограниченный максимум логарифмического правдоподобия. Мы начнем с максимальной оценки вероятности R10 и будем работать в обоих направлениях.

Lprof = nan(size(R10grid));
params = paramEsts;
[dum,peak] = min(abs(R10grid-R10MLE));
for i = peak:1:length(R10grid)
    PLconfun = ...
        @(params) deal([], gevinv(1-1./10,params(1),params(2),params(3)) - R10grid(i));
    [params,Lprof(i),flag,output] = ...
        fmincon(PLobjfun,params,[],[],[],[],[],[],PLconfun,opts);
end
params = paramEsts;
for i = peak-1:-1:1
    PLconfun = ...
        @(params) deal([], gevinv(1-1./10,params(1),params(2),params(3)) - R10grid(i));
    [params,Lprof(i),flag,output] = ...
        fmincon(PLobjfun,params,[],[],[],[],[],[],PLconfun,opts);
end

plot(R10grid,-Lprof,'-', R10MLE,-gevlike(paramEsts,y),'ro', ...
     [R10grid(1), R10grid(end)],[-nllCritVal,-nllCritVal],'k--');
xlabel('R_{10}');
ylabel('Log-Likelihood');
legend('Profile likelihood','MLE','95% Conf. Limit');