Исследуйте ежедневные закрытия глобальных данных о фондовом индексе

Необработанные данные состоят из 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 в противном случае (см. [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')

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

Учитывая стандартизированный, 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')

Оцените подгонку GPD

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

Калибруйте t Связку

Учитывая стандартизированные остаточные значения, теперь оцените скалярный параметр степеней свободы (степень свободы) и матрица (R) линейной корреляции t связки с помощью copulafit функция найдена в Statistics and Machine Learning Toolbox. copulafit функция позволяет пользователям оценивать параметры t связки двумя отличными методами.

Метод по умолчанию является формальным подходом наибольшего правдоподобия, выполняемым в двухступенчатом процессе.

А именно, несмотря на то, что полная логарифмическая правдоподобность могла быть максимизирована непосредственно, максимум целевой функции часто происходит в длинном, плоском, узком овраге на многомерном пробеле, и трудности со сходимостью могут возникнуть в высоких размерностях. Обойти эти трудности, copulafit выполняет оценку наибольшего правдоподобия (MLE) на двух шагах. Внутренний шаг максимизирует логарифмическую правдоподобность относительно матрицы линейной корреляции, учитывая фиксированное значение для степеней свободы. Та условная максимизация помещается в рамках 1D максимизации относительно степеней свободы, таким образом максимизируя логарифмическую правдоподобность по всем параметрам. Функция, максимизируемая на этом внешнем шаге, известна как логарифмическую правдоподобность профиля для степеней свободы.

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

А именно, приближение выведено путем дифференциации функции логарифмической правдоподобности относительно матрицы линейной корреляции, предположения, что степени свободы являются фиксированной постоянной. Получившееся нелинейное уравнение затем решено итеративно для корреляционной матрицы. Эта итерация, в свою очередь, вкладывается в рамках другой оптимизации, чтобы оценить степени свободы. Этот метод похож в духе на метод логарифмической правдоподобности профиля выше, но не является истинным MLE, в котором корреляционная матрица не сходится к условному MLE для данного степени свободы. Однако для размеров большой выборки оценки часто вполне близко к MLE (см. [1] и [7]).

Еще один метод оценивает матрицу линейной корреляции первым вычислением демонстрационной матрицы порядковой корреляции (tau Кендалла или ро Копьеносца), затем преобразует порядковую корреляцию в линейную корреляцию с устойчивым преобразованием синуса. Учитывая эту оценку матрицы линейной корреляции, функция логарифмической правдоподобности затем максимизируется относительно одного только параметра степеней свободы. Этот метод также кажется мотивированным подходом логарифмической правдоподобности профиля, но не является методом MLE вообще. Однако этот метод не требует матричной инверсии и поэтому имеет преимущество того, чтобы быть численно устойчивым в присутствии близко-к-сингулярному корреляционных матриц (см. [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 Связкой

Учитывая параметры 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:  14.4161%
Maximum Simulated Gain:  12.4511%

     Simulated 90% VaR:  -4.7456%
     Simulated 95% VaR:  -6.3244%
     Simulated 99% VaR: -10.6542%