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

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

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 часто используется в сочетании с третьим, пороговым параметром, который смещает нижний предел от нуля. Нам здесь такая общность не понадобится.

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

Симуляция данных о превышении

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

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

Реальные приложения для распределения GP включают моделирование экстремумов возвратов фондовых рынков и моделирование экстремальных наводнений. В этом примере мы будем использовать моделируемые данные, сгенерированные из распределения Student's 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 положительна.

Визуальная проверка подгонки

Чтобы визуально оценить, насколько хороша подгонка, мы построим масштабированную гистограмму хвостовых данных, покрытую функцией плотности GP, которую мы оценили. Гистограмма масштабируется так, чтобы высоты стержней умножали их ширину на 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 верна, и что у нас есть достаточно данных для асимптотического приближения к ковариационной матрице, чтобы удержать.

Проверка асимптотического допущения нормальности

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

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

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

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

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

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

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

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)');

Оценки bootstrap для 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).