В этом примере показано, как смоделировать рыночный риск гипотетического портфеля глобальных индексов акций с помощью метода моделирования Монте-Карло с использованием студенческой копулы и экстремальной теории стоимости (EVT). Процесс сначала извлекает отфильтрованные остатки из каждой серии возврата с помощью асимметричной модели GARCH, затем конструирует выборочную краевую кумулятивную функцию распределения (CDF) каждого актива с использованием оценки ядра Гаусса для внутреннего пространства и обобщенной оценки распределения Парето (GPD) для верхнего и нижнего хвостов. Затем t-копула Стьюдента подгоняется к данным и используется для индукции корреляции между смоделированными остатками каждого актива. Наконец, моделирование оценивает стоимость рискованного (VaR) гипотетического глобального портфеля акций на протяжении одного месяца.
Обратите внимание, что это относительно продвинутый, комплексный пример, который предполагает некоторое знакомство с EVT и копулами. Для получения подробной информации об оценке обобщенных распределений Парето и моделирования копул см. Данные об хвосте моделирования с обобщенным распределением Парето и моделированием зависимых случайных величин с использованием копул в Toolbox™ статистики и машинного обучения. Подробнее о подходе, на котором основана большая часть этого примера, см. [5] и [6]).
Необработанные данные состоят из 2665 наблюдений значений дневного закрытия следующих репрезентативных индексов акций, охватывающих даты торгов 27-April-1993 14-July-2003:
Canada: TSX Composite (Ticker ^GSPTSE)
France: CAC 40 (Ticker ^FCHI)
Germany: DAX (Ticker ^GDAXI)
Japan: Nikkei 225 (Ticker ^N225)
UK: FTSE 100 (Ticker ^FTSE)
US: S&P 500 (Ticker ^GSPC)load Data_GlobalIdx1 % Import daily index closings
На следующем графике показаны относительные изменения цен по каждому индексу. Первоначальный уровень каждого индекса был нормализован до единицы для облегчения сопоставления относительных показателей, и никакие корректировки дивидендов явно не учитываются.
figure plot(dates, ret2price(price2ret(Data))) datetick('x') xlabel('Date') ylabel('Index Value') title ('Relative Daily Index Closings') legend(series, 'Location', 'NorthWest')
При подготовке к последующему моделированию преобразуйте уровень закрытия каждого индекса в суточные логарифмические возвращаемые значения (иногда называемые геометрическими или непрерывно скомпонованными возвращаемыми значениями).
returns = price2ret(Data); % Logarithmic returns T = size(returns,1); % # of returns (i.e., historical sample size)
Поскольку первый шаг в общем подходе к моделированию включает в себя повторное применение фильтрации GARCH и теории экстремальных значений для характеристики распределения каждой отдельной серии возврата индекса капитала, полезно изучить детали для конкретной страны. Можно изменить следующую строку кода на любое целое число в наборе {1,2,3,4,5,6}, чтобы проверить сведения для любого индекса.
index = 1; % 1 = Canada, 2 = France, 3 = Germany, 4 = Japan, 5 = UK, 6 = US figure plot(dates(2:end), returns(:,index)) datetick('x') xlabel('Date') ylabel('Return') title('Daily Logarithmic Returns')
Моделирование хвостов распределения с помощью GPD требует, чтобы наблюдения были приблизительно независимыми и одинаково распределенными (i.i.d.). Однако большинство серий финансовой отдачи демонстрируют некоторую степень автокорреляции и, что более важно, гетероскедастичности.
Например, выборочная функция автокорреляции (ACF) возвращений, связанных с выбранным индексом, обнаруживает некоторую мягкую последовательную корреляцию.
figure
autocorr(returns(:,index))
title('Sample ACF of Returns')
Однако выборка ACF возведенных в квадрат возвратов иллюстрирует степень стойкости в дисперсии и подразумевает, что моделирование GARCH может значительно обусловить данные, используемые в последующем процессе оценки хвостовой части.
figure
autocorr(returns(:,index).^2)
title('Sample ACF of Squared Returns')
Для создания серии i.i.d. наблюдения, подгоняют авторегрессивную модель первого порядка к условному среднему доходам каждого индекса капитала
и асимметричная модель GARCH для условной дисперсии
![$$\sigma^2_t = \kappa + \alpha\sigma^2_{t-1} + \phi\epsilon^2_{t-1} + \psi[\epsilon_{t-1}<0]\epsilon^2_{t-1}$$](../examples/econ/win64/Demo_RiskEVT_eq16675095131510669661.png)
Авторегрессивная модель первого порядка компенсирует автокорреляцию, в то время как модель GARCH компенсирует гетероскедастичность. В частности, последний термин включает асимметрию (леверидж) в дисперсию по логическому показателю, который принимает значение 1, если остаток предыдущей модели отрицательный, и 0 в противном случае (см. [3]).
Кроме того, стандартизированные остатки каждого индекса моделируются как стандартизированное распределение Стьюдента для компенсации жировых хвостов, часто связанных с доходностью справедливости. То есть
Следующий сегмент кода извлекает отфильтрованные остатки и условные отклонения из возвращений каждого индекса собственного капитала.
model = arima('AR', NaN, 'Distribution', 't', 'Variance', gjr(1,1)); nIndices = size(Data,2); % # of indices residuals = NaN(T, nIndices); % preallocate storage variances = NaN(T, nIndices); fit = cell(nIndices,1); options = optimoptions(@fmincon, 'Display', 'off', ... 'Diagnostics', 'off', 'Algorithm', 'sqp', 'TolCon', 1e-7); for i = 1:nIndices fit{i} = estimate(model, returns(:,i), 'Display', 'off', 'Options', options); [residuals(:,i), variances(:,i)] = infer(fit{i}, returns(:,i)); end
Для выбранного индекса сравните остатки модели и соответствующие условные стандартные отклонения, отфильтрованные из необработанных результатов. Нижний график ясно иллюстрирует изменение летучести (гетероскедастичности), присутствующей в отфильтрованных остатках.
figure subplot(2,1,1) plot(dates(2:end), residuals(:,index)) datetick('x') xlabel('Date') ylabel('Residual') title ('Filtered Residuals') subplot(2,1,2) plot(dates(2:end), sqrt(variances(:,index))) datetick('x') xlabel('Date') ylabel('Volatility') title ('Filtered Conditional Standard Deviations')
Фильтруя остатки модели из каждого возвращаемого ряда, стандартизируйте остатки по соответствующему условному стандартному отклонению. Эти стандартизированные остатки представляют нижележащее нулевое среднее, единичная дисперсия, i.i.d. серии, на которых основана оценка EVT образцов хвостов CDF.
residuals = residuals ./ sqrt(variances);
Чтобы завершить этот раздел, изучите ACF стандартизированных остатков и квадратичных стандартизированных остатков.
Сравнение ACF стандартизированных остатков с соответствующими ACF необработанных результатов показывает, что стандартизированные остатки теперь приблизительно i.i.d., тем самым гораздо больше поддаются последующей оценке хвоста.
figure autocorr(residuals(:,index)) title('Sample ACF of Standardized Residuals') figure autocorr(residuals(:,index).^2) title('Sample ACF of Squared Standardized Residuals')
Учитывая стандартизированный i.i.d. остатки от предыдущего шага, оценить эмпирический CDF каждого индекса с гауссовым ядром. Это сглаживает оценки CDF, устраняя лестничный рисунок несамотированных образцов CDF. Хотя непараметрические оценки CDF ядра хорошо подходят для внутренней части распределения, где обнаруживается большая часть данных, они имеют тенденцию плохо работать при применении к верхнему и нижнему хвостам. Чтобы лучше оценить хвосты распределения, примените EVT к тем остаткам, которые падают в каждом хвосте.
В частности, найти верхний и нижний пороги так, что 10% остатков зарезервировано для каждого хвоста. Затем подгоните величину, на которую эти экстремальные остатки в каждом хвосте попадают за соответствующий порог, к параметрическому GPD по максимальной вероятности. Этот подход часто называют распределением превышений или пиков по пороговому методу.
Учитывая превышения в каждом хвосте, оптимизируйте функцию отрицательного логарифмического правдоподобия для оценки параметров индекса хвоста (zeta) и масштаба (beta) GPD.
Следующий сегмент кода создает объекты типа paretotails, один такой объект для каждого возвращаемого индексом ряда. Эти объекты паретотейлов инкапсулируют оценки параметрического нижнего хвоста GP, непараметрического внутреннего пространства, сглаженного ядром, и параметрического верхнего хвоста GP для построения составного полупараметрического CDF для каждого индекса.
Полученный кусочно-распределительный объект позволяет осуществлять интерполяцию внутри CDF и экстраполяцию (оценку функции) в каждом хвосте. Экстраполяция очень желательна, позволяя оценивать квантили вне исторических данных, и является бесценной для приложений управления рисками.
Кроме того, хвостовые объекты Парето также обеспечивают методы для оценки CDF и обратной CDF (квантовой функции) и для запроса кумулятивных вероятностей и квантилей границ между каждым сегментом кусочного распределения.
nPoints = 200; % # of sampled points in each region of the CDF tailFraction = 0.1; % Decimal fraction of residuals allocated to each tail tails = cell(nIndices,1); % Cell array of Pareto tail objects for i = 1:nIndices tails{i} = paretotails(residuals(:,i), tailFraction, 1 - tailFraction, 'kernel'); end
Оценив три различных области составного полупараметрического эмпирического CDF, графически объедините и отобразите результат. Опять же, обратите внимание, что нижняя и верхняя хвостовые области, отображаемые соответственно красным и синим цветом, пригодны для экстраполяции, в то время как сглаженная ядром внутренняя часть, черная, пригодна для интерполяции.
Код ниже вызывает CDF и методы обратного CDF интересующего объекта Парето хвостов с данными, отличными от тех, на которых основана подгонка. В частности, указанные методы имеют доступ к подобранному состоянию и теперь вызываются для выбора и анализа конкретных областей кривой вероятности, действуя как мощный механизм фильтрации данных.
figure hold on grid on minProbability = cdf(tails{index}, (min(residuals(:,index)))); maxProbability = cdf(tails{index}, (max(residuals(:,index)))); pLowerTail = linspace(minProbability , tailFraction , nPoints); % sample lower tail pUpperTail = linspace(1 - tailFraction, maxProbability , nPoints); % sample upper tail pInterior = linspace(tailFraction , 1 - tailFraction, nPoints); % sample interior plot(icdf(tails{index}, pLowerTail), pLowerTail, 'red' , 'LineWidth', 2) plot(icdf(tails{index}, pInterior) , pInterior , 'black', 'LineWidth', 2) plot(icdf(tails{index}, pUpperTail), pUpperTail, 'blue' , 'LineWidth', 2) xlabel('Standardized Residual') ylabel('Probability') title('Empirical CDF') legend({'Pareto Lower Tail' 'Kernel Smoothed Interior' ... 'Pareto Upper Tail'}, 'Location', 'NorthWest')
Хотя предыдущий график иллюстрирует композитный CDF, поучительно исследовать вписывание GPD более подробно.
CDF распределения GP параметризуется как
для превышений (y), параметра индекса хвостовой части (zeta) и параметра шкалы (beta).
Для визуальной оценки соответствия GPD постройте график эмпирического CDF превышения остатков верхнего хвоста вместе с CDF, оборудованным GPD. Хотя используется только 10% стандартизированных остатков, соответствующее распределение точно соответствует данным о превышении, поэтому модель GPD, как представляется, является хорошим выбором.
figure
[P,Q] = boundary(tails{index}); % cumulative probabilities & quantiles at boundaries
y = sort(residuals(residuals(:,index) > Q(2), index) - Q(2)); % sort exceedances
plot(y, (cdf(tails{index}, y + Q(2)) - P(2))/P(1))
[F,x] = ecdf(y); % empirical CDF
hold on
stairs(x, F, 'r')
grid on
legend('Fitted Generalized Pareto CDF','Empirical CDF','Location','SouthEast');
xlabel('Exceedance')
ylabel('Probability')
title('Upper Tail of Standardized Residuals')
Учитывая стандартизированные остатки, теперь оцените параметр скалярных степеней свободы (DoF) и матрицу линейной корреляции (R) t-копулы, используя copulafit находится в панели инструментов статистики и машинного обучения. copulafit функция позволяет пользователям оценить параметры t-копулы двумя различными методами.
Метод по умолчанию - это формальный подход максимального правдоподобия, выполняемый в двухшаговом процессе.
В частности, хотя полное логарифмическое правдоподобие может быть максимизировано непосредственно, максимум целевой функции часто возникает в длинной, плоской, узкой долине в многомерном пространстве, и трудности сходимости могут возникать в больших измерениях. Чтобы обойти эти трудности, copulafit выполняет оценку максимального правдоподобия (MLE) в два этапа. Внутренний шаг максимизирует логарифмическое правдоподобие относительно матрицы линейной корреляции, заданной фиксированным значением для степеней свободы. Эта условная максимизация помещается в 1-D максимизацию относительно степеней свободы, таким образом максимизируя логарифмическую вероятность по всем параметрам. Функция, максимизированная на этом внешнем шаге, известна как логарифмическая вероятность профиля для степеней свободы.
Напротив, следующий сегмент кода использует альтернативу, которая аппроксимирует логарифмическое правдоподобие профиля для параметра степеней свободы для больших размеров выборки. Хотя этот метод часто значительно быстрее, чем MLE, его следует использовать с осторожностью, потому что оценки и доверительные пределы могут быть неточными для небольших или умеренных размеров выборки.
В частности, аппроксимация получается путем дифференцирования логарифмической функции правдоподобия относительно матрицы линейной корреляции, предполагая, что степени свободы являются фиксированной константой. Результирующее нелинейное уравнение затем решается итеративно для корреляционной матрицы. Эта итерация, в свою очередь, вложена в другую оптимизацию для оценки степеней свободы. Этот метод подобен по духу описанному выше методу логарифмического правдоподобия профиля, но не является истинным MLE, поскольку корреляционная матрица не сходится с условным MLE для заданных степеней свободы. Однако для больших размеров выборки оценки часто довольно близки к MLE (см. [1] и [7]).
Еще один способ оценивает матрицу линейной корреляции, сначала вычисляя матрицу ранговой корреляции выборки (тау Кендалла или ро Спирмена), затем преобразует ранговую корреляцию в линейную корреляцию с устойчивым синусоидальным преобразованием. Учитывая эту оценку матрицы линейной корреляции, логарифмическая функция правдоподобия затем максимизируется по отношению только к параметру степеней свободы. Этот метод также кажется мотивированным подходом логарифмического правдоподобия профиля, но не является методом MLE вообще. Однако этот способ не требует матричной инверсии и поэтому имеет преимущество, заключающееся в том, что он является численно стабильным при наличии матриц близкой к сингулярной корреляции (см. [8]).
Наконец, Nystrom и Skoglund [6] предлагают зарезервировать параметр степеней свободы в качестве пользовательского ввода моделирования, что позволяет пользователю субъективно индуцировать степень зависимости хвостовой части между активами. В частности, они рекомендуют относительно низкие степени свободы, где-то между 1 и 2, что позволяет более тщательно изучить вероятность экстремальных состояний суставов. Этот метод полезен для стресс-тестирования, в котором степень крайней взаимозависимости имеет решающее значение.
Сегмент кода, расположенный ниже, сначала преобразует стандартизированные остатки в однородные вариации с помощью полупараметрического эмпирического CDF, полученного выше, а затем подгоняет t-копулу к преобразованным данным. Когда однородные вариации преобразуются эмпирическим CDF каждого поля, метод калибровки часто известен как каноническое максимальное правдоподобие (CML).
U = zeros(size(residuals)); for i = 1:nIndices U(:,i) = cdf(tails{i}, residuals(:,i)); % transform margin to uniform end [R, DoF] = copulafit('t', U, 'Method', 'approximateml'); % fit the copula
Учитывая параметры t-копулы, теперь моделируют совместно зависимые доходности индекса собственного капитала, сначала моделируя соответствующие зависимые стандартизированные остатки.
Для этого сначала смоделируйте зависимые однородные вариации с помощью copularnd находится в панели инструментов статистики и машинного обучения.
Затем путем экстраполяции в хвосты ГП и интерполяции в сглаженный интерьер преобразуют однородные вариации в стандартизированные остатки посредством инверсии полупараметрического краевого CDF каждого индекса. Это дает смоделированные стандартизированные остатки, согласующиеся с остатками, полученными в вышеописанном процессе фильтрации AR (1) + GJR (1,1). Эти остатки независимы по времени, но зависимы в любой момент времени. Каждый столбец моделируемого массива стандартизированных остатков представляет i.i.d. одномерный стохастический процесс при рассмотрении изолированно, тогда как каждая строка имеет ранговую корреляцию, индуцированную копулой.
Следующий сегмент кода моделирует 2000 независимых случайных испытаний зависимых стандартизированных остатков индекса в течение месячного горизонта в 22 торговых дня.
s = RandStream.getGlobalStream(); reset(s) nTrials = 2000; % # of independent random trials horizon = 22; % VaR forecast horizon Z = zeros(horizon, nTrials, nIndices); % standardized residuals array U = copularnd('t', R, DoF, horizon * nTrials); % t copula simulation for j = 1:nIndices Z(:,:,j) = reshape(icdf(tails{j}, U(:,j)), horizon, nTrials); end
Использование моделируемых стандартизированных остатков в качестве i.i.d. процесс входного шума, повторное введение автокорреляции и гетероскедастичности, наблюдаемых в исходном индексе, возвращается через Econometrics Toolbox™ filter функция.
Обратите внимание, что filter принимает заданные пользователем стандартизированные возмущения, полученные из копулы, и моделирует несколько путей для одной индексной модели одновременно, что означает, что все пути выборки моделируются и сохраняются для каждого индекса последовательно.
Чтобы максимально использовать текущую информацию, укажите необходимые остатки, отклонения и возвращаемые значения модели предварительной выборки таким образом, чтобы каждый моделируемый путь развивался из общего начального состояния.
Y0 = returns(end,:); % presample returns Z0 = residuals(end,:); % presample standardized residuals V0 = variances(end,:); % presample variances simulatedReturns = zeros(horizon, nTrials, nIndices); for i = 1:nIndices simulatedReturns(:,:,i) = filter(fit{i}, Z(:,:,i), ... 'Y0', Y0(i), 'Z0', Z0(i), 'V0', V0(i)); end
Теперь измените форму моделируемого массива возвращений таким образом, чтобы каждая страница представляла одну пробную версию многомерной возвращаемой серии, а не несколько пробных версий одномерной возвращаемой серии.
simulatedReturns = permute(simulatedReturns, [1 3 2]);
Наконец, с учетом смоделированной доходности каждого индекса сформировать одинаково взвешенный портфель глобальных индексов, состоящий из отдельных индексов (глобальный индекс индексов стран). Поскольку мы работаем с ежедневной логарифмической доходностью, кумулятивная доходность за горизонт риска - это просто суммы доходности за каждый промежуточный период. Также следует отметить, что веса портфеля сохраняются фиксированными на протяжении всего горизонта риска и что моделирование игнорирует любые операционные затраты, необходимые для ребалансировки портфеля (предполагается, что ежедневный процесс ребалансировки является самофинансированием).
Следует отметить, что, несмотря на то, что моделируемые доходности являются логарифмическими (непрерывно скомпонованными), ряд возвратов портфеля строится путем преобразования сначала отдельных логарифмических доходностей в арифметические доходности (изменение цены, деленное на начальную цену), затем взвешивания отдельных арифметических возвратов для получения арифметических возвратов портфеля и, наконец, преобразования обратно в логарифмические доходности портфеля. При ежедневных данных и коротком горизонте VaR повторные преобразования мало чем отличаются, но в течение более длительных периодов времени различия могут быть значительными.
cumulativeReturns = zeros(nTrials, 1); weights = repmat(1/nIndices, nIndices, 1); % equally weighted portfolio for i = 1:nTrials cumulativeReturns(i) = sum(log(1 + (exp(simulatedReturns(:,:,i)) - 1) * weights)); end
Смоделировав доходность каждого индекса и сформировав глобальный портфель, сообщите о максимальных прибылях и убытках, а также о VaR на различных уровнях доверия за месячный горизонт риска. Также постройте график эмпирического CDF глобального портфеля.
VaR = 100 * quantile(cumulativeReturns, [0.10 0.05 0.01]'); disp(' ') fprintf('Maximum Simulated Loss: %8.4f%s\n' , -100*min(cumulativeReturns), '%') fprintf('Maximum Simulated Gain: %8.4f%s\n\n' , 100*max(cumulativeReturns), '%') fprintf(' Simulated 90%% VaR: %8.4f%s\n' , VaR(1), '%') fprintf(' Simulated 95%% VaR: %8.4f%s\n' , VaR(2), '%') fprintf(' Simulated 99%% VaR: %8.4f%s\n\n', VaR(3), '%') figure h = cdfplot(cumulativeReturns); h.Color = 'Red'; xlabel('Logarithmic Return') ylabel('Probability') title ('Simulated One-Month Global Portfolio Returns CDF')
Maximum Simulated Loss: 14.4161%
Maximum Simulated Gain: 12.4511%
Simulated 90% VaR: -4.7456%
Simulated 95% VaR: -6.3244%
Simulated 99% VaR: -10.6542%
[1] Буйе, Э., В. Дуррлеман, А. Никегбали, Г. Рибуле и Ронкалли, Т. «Копулас для финансов: руководство по чтению и некоторые приложения». Груп де Реч. Опер., Кредит Лионне, Париж, 2000.
[2] Эмбрехтс, П., А. Макнейл и Д. Штрауманн. «Корреляция и зависимость в управлении рисками: свойства и подводные камни». Управление рисками: ценность под угрозой и за ее пределами. Кембридж: Cambridge University Press, 1999, стр. 176-223.
[3] Glosten, L. R., Р. Джейгэннэзэн и Д. Э. Ранкл. «О связи между ожидаемой стоимостью и волатильностью номинальной избыточной доходности акций». Финансовый журнал. т. 48, № 5, 1993, с. 1779-1801.
[4] Макнейл, А. и Р. Фрэй. «Оценка показателя риска, связанного с хвостом, для гетероскедастического финансового временного ряда: подход с предельной стоимостью». Журнал эмпирических финансов. Том 7, 2000, стр. 271-300.
[5] Нистром, К. и Й. Скоглунд. «Одномерная теория экстремальных ценностей, GARCH и меры риска». Препринт, представлен в 2002 году.
[6] Нистром, К. и Й. Скоглунд. «Основа для управления рисками на основе сценариев». Препринт, представлен в 2002 году.
[7] Ронкалли, Т., А. Дуррлеман и А. Никегбали. «Какая Копула правильная?» Груп де Реч. Опер., Кредит Лионне, Париж, 2000.
[8] Машаль, Р. и А. Зееви. «Вне корреляции: экстремальные совместные движения между финансовыми активами». Колумбийский университет, Нью-Йорк, 2002 год.