Изучение ежедневных закрытий данных Глобального индекса капитала

Необработанные данные состоят из 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 с условным отклонением

Авторегрессионная модель первого порядка компенсирует автокорреляцию, в то время как модель GARCH компенсирует гетероскедастичность. В частности, последний член включает асимметрию (рычаг) в отклонение булевым индикатором, который принимает значение 1, если невязка предыдущей модели отрицательна, и 0 в противном случае (см. [3]).

Кроме того, стандартизированные невязки каждого индекса моделируются как стандартизированное распределение Student's 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

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

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

Оценка полупараметрических CDF

Учитывая стандартизированный, i.i.d. невязки от предыдущего шага оценивают эмпирический CDF каждого индекса с Гауссовым ядром. Это сглаживает оценки CDF, устраняя диаграмму шаблона лестничной клетки не смоченных выборок CDF. Хотя непараметрические оценки CDF ядра хорошо подходят для внутреннего распределения, где большая часть данных найдена, они, как правило, плохо выполняются при применении к верхнему и нижнему хвостам. Чтобы лучше оценить хвосты распределения, примените EVT к тем невязкам, которые падают в каждом хвосте.

В частности, найдите верхний и нижний пороги, так что 10% невязок зарезервированы для каждого хвоста. Затем подгоняйте величину, на которую эти крайние невязки в каждом хвосте падают сверх связанного порога, к параметрическому GPD по максимальной вероятности. Этот подход часто упоминается как распределение превышений или peaks по пороговому методу.

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

Оценка подгонки GPD

Несмотря на то, что предыдущий график иллюстрирует составной 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')

Калибровка t-составляющей Копулы

Учитывая стандартизированные невязки, теперь оцените скалярный параметр степеней свободы (DoF) и линейную корреляционную матрицу (R) t-копулы, используя copulafit функция, найденная в наборе инструментов Statistics and Machine Learning Toolbox. The 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-тью Копулы

Учитывая параметры t-копулы, теперь моделируйте совместно зависимые возвраты индекса собственного капитала, сначала симулируя соответствующие зависимые стандартизированные невязки.

Для этого сначала моделируйте зависимые равномерные изменения с помощью copularnd функция, найденная в наборе инструментов Statistics and Machine Learning Toolbox.

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