Этот пример показывает, как смоделировать риск рынка гипотетического глобального портфеля фондового индекса с методом симуляции Монте-Карло с помощью t связки и Теории экстремума (EVT) Студента. Процесс сначала извлекает отфильтрованные невязки от каждого, возвращают ряд с асимметричной моделью GARCH, затем создает демонстрационную крайнюю кумулятивную функцию распределения (CDF) каждого актива с помощью Гауссовой оценки ядра для внутренней части и оценки обобщенного распределения Парето (GPD) для верхних и более низких хвостов. T связка Студента является затем подходящей к данным и используемая, чтобы вызвать корреляцию между моделируемыми невязками каждого актива. Наконец, симуляция оценивает Подверженный риску значения (VaR) гипотетического глобального портфеля акции по одному горизонту месяца.
Обратите внимание на то, что это - относительно усовершенствованный, всесторонний пример, который принимает некоторое знакомство с EVT и связками. Для получения дополнительной информации относительно оценки обобщенных дистрибутивов Парето и симуляции связки, смотрите Данные о Хвосте Моделирования с Обобщенным Распределением Парето (Statistics and Machine Learning Toolbox) и Симуляция Зависимых Случайных переменных Используя Связки (Statistics and Machine Learning Toolbox) в Statistics and Machine Learning Toolbox™. Для получения дополнительной информации относительно подхода, на котором базируется большая часть этого примера, смотрите ссылки [5] и [6] Nystrom и Skoglund в библиографии.
Необработанные данные состоят из 2 665 наблюдений за ежедневными заключительными значениями следующих представительных фондовых индексов, охватывающих торговые даты 27 апреля 1993 до 14 июля 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).. Однако большинство финансовых рядов возврата показывает определенную степень автокорреляции и, что еще более важно, heteroskedasticity.
Например, демонстрационная автокорреляционная функция (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 к условному отклонению
Авторегрессивная модель первого порядка компенсирует автокорреляцию, в то время как модель GARCH компенсирует heteroskedasticity. В частности, последний срок включает асимметрию (рычаги) в отклонение булевым индикатором, который принимает значение 1, если предшествующая образцовая невязка отрицательна и 0 в противном случае (см. Glosten, Jagannathan, & Runkle [3]).
Кроме того, стандартизированные невязки каждого индекса моделируются, когда t распределение стандартизированного Студента, чтобы компенсировать толстые хвосты, часто сопоставляемые с акцией, возвращается. Это
Следующий сегмент кода извлекает отфильтрованные невязки и условные отклонения от возвратов каждого фондового индекса.
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
Для выбранного индекса сравните образцовые невязки, и соответствующие условные стандартные отклонения, отфильтрованные от сырых данных, возвращается. Более низкий график ясно иллюстрирует изменение в энергозависимости (heteroskedasticity) существующий в отфильтрованных невязках.
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);
Чтобы завершить этот раздел, исследуйте ACFs стандартизированных невязок и стандартизированных невязок в квадрате.
Сравнение ACFs стандартизированных невязок к соответствующему ACFs сырых данных возвращается, показывает, что стандартизированные невязки теперь приблизительно 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, устраняя шаблон лестницы не сглаживавшего демонстрационного CDFS. Несмотря на то, что непараметрическое ядро, оценки CDF хорошо подходят для внутренней части распределения, где большинство данных найден, они имеют тенденцию выполнять плохо, когда применено верхние и более низкие хвосты. Чтобы лучше оценить хвосты распределения, примените EVT к тем невязкам, которые падают в каждом хвосте.
А именно, найдите верхние и более низкие пороги таким образом, что 10% невязок резервируются для каждого хвоста. Затем соответствуйте сумме, которой те экстремальные невязки в каждом хвосте падают вне связанного порога на параметрический GPD наибольшим правдоподобием. Этот подход часто упоминается как распределение exceedances или peaks по пороговому методу.
Учитывая exceedances в каждом хвосте, оптимизируйте отрицательную логарифмическую функцию правдоподобия, чтобы оценить индекс хвоста (дзэта) и шкала (бета) параметры GPD.
Следующий сегмент кода создает объекты типа paretotails
, один такой объект для каждого индекса возвращает ряд. Эти объекты 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 параметризован как
для exceedances (y), индексный параметр хвоста (дзэта) и масштабный коэффициент (бета).
Чтобы визуально оценить подгонку GPD, постройте эмпирический CDF верхнего хвоста exceedances невязок наряду с CDF, адаптированным GPD. Несмотря на то, что только 10% стандартизированных невязок используются, подходящее распределение тесно следует за exceedance данными, таким образом, модель 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')
Учитывая стандартизированные невязки, теперь оцените скалярный параметр степеней свободы (степень свободы) и матрица (R) линейной корреляции t связки с помощью функции copulafit
, найденной в Statistics and Machine Learning Toolbox. Функция copulafit
позволяет пользователям оценивать параметры t связки двумя отличными методами.
Метод по умолчанию является формальным подходом наибольшего правдоподобия, выполняемым в двухступенчатом процессе.
А именно, несмотря на то, что полная логарифмическая вероятность могла быть максимизирована непосредственно, максимум целевой функции часто происходит в длинной, плоской, узкой долине на многомерном пробеле, и трудности со сходимостью могут возникнуть в высоких размерностях. Чтобы обойти эти трудности, copulafit
выполняет оценку наибольшего правдоподобия (MLE) на двух шагах. Внутренний шаг максимизирует логарифмическую вероятность относительно матрицы линейной корреляции, учитывая фиксированное значение для степеней свободы. Та условная максимизация помещается в рамках 1D максимизации относительно степеней свободы, таким образом максимизируя логарифмическую вероятность по всем параметрам. Функция, максимизируемая на этом внешнем шаге, известна как логарифмическую вероятность профиля для степеней свободы.
Напротив, следующий сегмент кода использует альтернативу, которая аппроксимирует логарифмическую вероятность профиля для параметра степеней свободы для размеров большой выборки. Несмотря на то, что этот метод часто значительно быстрее, чем MLE, он должен использоваться с осторожностью, потому что оценки и пределы достоверности не могут быть точными для маленьких или умеренных объемов выборки.
А именно, приближение выведено путем дифференциации логарифмической функции правдоподобия относительно матрицы линейной корреляции, предположения, что степени свободы являются фиксированной постоянной. Получившееся нелинейное уравнение затем решено итеративно для корреляционной матрицы. Эта итерация, в свою очередь, вкладывается в рамках другой оптимизации, чтобы оценить степени свободы. Этот метод подобен в духе методу логарифмической вероятности профиля выше, но не является истинным MLE, в котором корреляционная матрица не сходится к условному MLE для данного степени свободы. Однако для размеров большой выборки оценки часто вполне близко к MLE (см. Bouye и др. [1], и Деррлмен, Nikeghbali, Ронкалли [7]).
Еще один метод оценивает матрицу линейной корреляции первым вычислением демонстрационной матрицы порядковой корреляции (tau Кендалла или ро Копьеносца), затем преобразовывает порядковую корреляцию в линейную корреляцию с устойчивым преобразованием синуса. Учитывая эту оценку матрицы линейной корреляции, логарифмическая функция правдоподобия затем максимизируется относительно одного только параметра степеней свободы. Этот метод также кажется мотивированным подходом логарифмической вероятности профиля, но не является методом MLE вообще. Однако этот метод не требует матричной инверсии и поэтому имеет преимущество того, чтобы быть численно стабильным в присутствии близко-к-сингулярному корреляционных матриц (см. Zeevi и Mashal [8]).
Наконец, Nystrom и Skoglund [6] предлагают резервировать параметр степеней свободы как заданный пользователями вход симуляции, таким образом позволяя пользователю субъективно вызвать степень зависимости хвоста между активами. В частности, они рекомендуют относительно низкие степени свободы, где-нибудь между 1 и 2, таким образом позволяя более близкое исследование вероятности объединенных экстремальных значений. Этот метод полезен для стресс-тестирования, в котором степень экстремального значения codependence имеет жизненное значение.
Сегмент кода ниже первого преобразовывает стандартизированные невязки к универсальным варьируемым величинам полупараметрическим эмпирическим 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
, найденной в Statistics and Machine Learning Toolbox.
Затем путем экстраполирования в хвосты GP и интерполяции в сглаживавшую внутреннюю часть, преобразуйте универсальные варьируемые величины к стандартизированным невязкам через инверсию полупараметрического крайнего CDF каждого индекса. Это производит моделируемые стандартизированные невязки, сопоставимые с полученными из AR (1) + GJR (1,1) процесс фильтрации выше. Эти невязки независимы вовремя, но зависимый в любом моменте времени. Каждый столбец моделируемого стандартизированного массива невязок представляет i.i.d. одномерный стохастический процесс, когда просматривается в изоляции, тогда как каждая строка совместно использует порядковую корреляцию, вызванную связкой.
Следующий сегмент кода моделирует 2 000 независимых случайных испытаний стандартизированных индексных невязок зависимого по одному горизонту месяца 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. вводят шумовой процесс, повторно вводят автокорреляцию, и heteroskedasticity, наблюдаемый в исходном индексе, возвращается через функцию filter
Econometrics Toolbox™.
Обратите внимание на то, что 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: 15.2681% Maximum Simulated Gain: 12.0603% Simulated 90% VaR: -4.8300% Simulated 95% VaR: -6.3838% Simulated 99% VaR: -10.7009%
Этот пример основан на следующих бумагах и статьях в журнале:
[1] Bouye, E., В. Деррлмен, А. Никегбали, Г. Рибулет и Ронкалли, T. "Связки для Финансов: Руководство Чтения и Некоторые Приложения". Groupe de Rech. Oper., Credit Lyonnais, Париж, 2000.
[2] Embrechts, P., А. Макнейл и Д. Штрауманн. "Корреляция и Зависимость в управлении рисками: Свойства и Ловушки". Управление рисками: Значение В опасности и Вне. Кембридж: Издательство Кембриджского университета, 1999, стр 176–223.
[3] Glosten, L. R. Р. Джейгэннэзэн и Д. Э. Ранкл. "На Отношении между Ожидаемым значением и Энергозависимостью Номинального Избыточного Возврата на Запасах". Журнал Финансов. Издание 48, 1993, стр 1779–1801.
[4] Макнейл, А. и Р. Фрэй. "Оценка Связанной с хвостом Меры по Риску для Финансовых Временных рядов Heteroscedastic: Подход Экстремума". Журнал Эмпирических Финансов. Издание 7, 2000, стр 271–300.
[5] Nystrom, K. и Дж. Скогланд. "Одномерная Теория Экстремума, GARCH и Меры Риска". Предварительно распечатайте, представленный 2002. http://gloria-mundi.com/UploadFile/2010-2/knjsug.pdf.
[6] Nystrom, K. и Дж. Скогланд. "Среда для Основанного на сценарии управления рисками". Предварительно распечатайте, представленный 2002. http://gloria-mundi.com/UploadFile/2010-2/knjs.pdf.
[7] Ронкалли, T., А. Деррлмен и А. Никегбали. "Какая Связка Является Правильной?" Groupe de Rech. Oper., Credit Lyonnais, Париж, 2000.
[8] Mashal, R. и А. Зиви. "Вне корреляции: экстремальные Co-перемещения между финансовыми активами". Колумбийский университет, Нью-Йорк, 2002.