Обобщенное распределение Парето

Обобщенный парето (GP) - это скошенное вправо распределение, параметризованное параметром формы k и параметром масштаба sigma. k также известен как параметр «хвостовой индекс» и может быть положительным, нулевым или отрицательным.

x = linspace(0,10,1000);
plot(x,gppdf(x,-.4,1),'-', x,gppdf(x,0,1),'-', x,gppdf(x,2,1),'-');
xlabel('x / sigma');
ylabel('Probability density');
legend({'k < 0' 'k = 0' 'k > 0'});

Обратите внимание, что для k < 0 вероятность GP выше верхнего предела - (1/k). Для k > = 0 GP не имеет верхнего предела. Кроме того, GP часто используется в сочетании с третьим пороговым параметром, который смещает нижний предел от нуля. Нам не нужна здесь такая общность.

Распределение ГП является обобщением как экспоненциального распределения (k = 0), так и распределения Парето (k > 0). GP включает эти два распределения в большем семействе, так что возможен непрерывный диапазон форм.

Моделирование данных превышения

Распределение ГП может быть определено конструктивно с точки зрения превышений. Начиная с вероятностного распределения, правый хвост которого падает до нуля, такого как нормаль, мы можем отсчитывать случайные значения независимо от этого распределения. Если мы фиксируем пороговое значение, выбрасываем все значения, которые ниже порогового значения, и вычитаем пороговое значение из значений, которые не выбрасываются, результат называется превышениями. Распределение превышений составляет приблизительно ГП. Аналогично, мы можем установить порог в левом хвосте распределения и игнорировать все значения выше этого порога. Порог должен быть достаточно далеко в хвосте исходного распределения, чтобы аппроксимация была разумной.

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

Реальные приложения для распределения ГП включают моделирование экстремальных показателей доходности фондовых рынков и моделирование экстремальных наводнений. В этом примере мы будем использовать смоделированные данные, созданные из распределения Стьюдента с 5 степенями свободы. Мы возьмем самые большие 5% из 2000 наблюдений из распределения t, а затем вычтем 95% квантиль, чтобы получить превышения.

rng(3,'twister');
x = trnd(5,2000,1);
q = quantile(x,.95);
y = x(x>q) - q;
n = numel(y)

n =

   100

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

Распределение GP определяется для 0 < sigma и -Inf < k < Inf. Однако интерпретация результатов оценки максимального правдоподобия проблематична, когда k < -1/2. К счастью, эти случаи соответствуют подгонке хвостов из распределений, таких как бета или треугольные, и поэтому не будут представлять проблемы здесь.

paramEsts = gpfit(y);
kHat      = paramEsts(1)   % Tail index parameter
sigmaHat  = paramEsts(2)   % Scale parameter

kHat =

    0.0987


sigmaHat =

    0.7156

Как можно ожидать, поскольку смоделированные данные были сгенерированы с использованием распределения t, оценка k является положительной.

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

Чтобы визуально оценить, насколько хорошо подходит, мы построим масштабированную гистограмму данных хвоста, наложенную на функцию плотности ГП, которую мы оценили. Гистограмма масштабируется таким образом, что высота полосы умножает сумму ее ширины на 1.

bins = 0:.25:7;
h = bar(bins,histc(y,bins)/(length(y)*.25),'histc');
h.FaceColor = [.9 .9 .9];
ygrid = linspace(0,1.1*max(y),100);
line(ygrid,gppdf(ygrid,kHat,sigmaHat));
xlim([0,6]);
xlabel('Exceedance');
ylabel('Probability Density');

Мы использовали довольно небольшую ширину бункера, так что в гистограмме много шума. Тем не менее, подогнанная плотность соответствует форме данных, и поэтому модель GP кажется хорошим выбором.

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

[F,yi] = ecdf(y);
plot(yi,gpcdf(yi,kHat,sigmaHat),'-');
hold on;
stairs(yi,F,'r');
hold off;
legend('Fitted Generalized Pareto CDF','Empirical CDF','location','southeast');

Вычисление стандартных ошибок для оценок параметров

Чтобы количественно определить точность оценок, мы будем использовать стандартные ошибки, вычисленные из асимптотической ковариационной матрицы оценщиков максимального правдоподобия. Функция gplike вычисляет в качестве второго выхода численное приближение к этой ковариационной матрице. В качестве альтернативы мы могли бы позвонить gpfit с двумя выходными аргументами, и он бы вернул доверительные интервалы для параметров.

[nll,acov] = gplike(paramEsts, y);
stdErr = sqrt(diag(acov))

stdErr =

    0.1158
    0.1093

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

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

Интерпретация стандартных ошибок обычно предполагает, что, если одна и та же подгонка может быть повторена много раз для данных, полученных из одного и того же источника, оценки максимального правдоподобия параметров будут приблизительно следовать нормальному распределению. Например, доверительные интервалы часто основаны на этом предположении.

Однако это нормальное приближение может быть хорошим, а может и нет. Чтобы оценить, насколько это хорошо в этом примере, мы можем использовать моделирование начальной загрузки. Мы создадим 1000 наборов данных репликации путем повторной выборки из данных, подберем распределение GP для каждого и сохраним все оценки репликации.

replEsts = bootstrp(1000,@gpfit,y);

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

subplot(2,1,1);
hist(replEsts(:,1));
title('Bootstrap estimates of k');
subplot(2,1,2);
hist(replEsts(:,2));
title('Bootstrap estimates of sigma');

Использование преобразования параметров

Гистограмма оценок начальной загрузки для k кажется лишь немного асимметричной, в то время как для оценок sigma определенно выглядит скошенной вправо. Обычным средством устранения этого перекоса является оценка параметра и его стандартной ошибки на логарифмической шкале, где нормальная аппроксимация может быть более разумной. График Q-Q является лучшим способом оценки нормальности, чем гистограмма, потому что ненормальность проявляется как точки, которые приблизительно не следуют за прямой линией. Давайте проверим, подходит ли преобразование журнала для sigma.

subplot(1,2,1);
qqplot(replEsts(:,1));
title('Bootstrap estimates of k');
subplot(1,2,2);
qqplot(log(replEsts(:,2)));
title('Bootstrap estimates of log(sigma)');

Оценки начальной загрузки для k и log (sigma) оказываются приемлемо близкими к нормальности. График Q-Q для оценок сигмы, на незарегистрированной шкале, подтвердит перекос, который мы уже видели на гистограмме. Таким образом, было бы более разумно построить доверительный интервал для сигмы, сначала вычисляя один для логарифма (сигма) в предположении нормальности, а затем возведя в степень, чтобы преобразовать этот интервал обратно в исходный масштаб для сигмы.

На самом деле, это именно то, что функция gpfit делает за кадром.

[paramEsts,paramCI] = gpfit(y);

kHat
kCI  = paramCI(:,1)

kHat =

    0.0987


kCI =

   -0.1283
    0.3258

sigmaHat
sigmaCI  = paramCI(:,2)

sigmaHat =

    0.7156


sigmaCI =

    0.5305
    0.9654

Обратите внимание, что в то время как 95% доверительный интервал для k является симметричным относительно оценки максимального правдоподобия, доверительный интервал для sigma не является. Это потому, что он был создан путем преобразования симметричного CI для log (sigma).