Ниже приведен последовательный алгоритм построения портфеля с минимальным значением CVaR на один торговый день вперед при помощи копула-GARCH модели. Начнем с выгрузки данных. Для построения портфеля использутся торгуемые на бирже фонды SPDR S&P 500 ETF Trust (тикер – SPY) и iShares MSCI Hong Kong ETF (тикер – EWH). Котировки выгружаются из Yahoo Finance за период с 01.01.2015 по 01.11.2019.

Загрузка модулей

library(quantmod)
library(PMwR) 
library(rugarch)
library(VineCopula) 
library(copula)
library(FRAPO)
library(fPortfolio)
library(QRM)
library(plyr)
library(FinTS)
library(pastecs) 
library(tseries) 
library(LambertW)

Выгрузка данных и первичный анализ

Выгружаем котировки из Yahoo Finance

spy_data <- as.timeSeries(getSymbols("SPY", from = "2015-01-01", to = "2019-11-01"))
spy_close <- spy_data[,4]
ewh_data <- as.timeSeries(getSymbols("EWH", from = "2015-01-01", to = "2019-11-01"))
ewh_close <- ewh_data[,4]

Расчет логарифмических доходностей

spy_ret <- na.omit(diff(log(spy_close)))
ewh_ret <- na.omit(diff(log(ewh_close)))

Графики дневных логарифмических доходностей

plot(spy_ret, xlab = "Дни", ylab = "Логарифмическая доходность", main = "Дневные логарифмические доходности SPY")

plot(ewh_ret, xlab = "Дни", ylab = "Логарифмическая доходность", main = "Дневные логарифмические доходности EWH")

Находим числовые характеристики временных рядов

nrow(spy_ret) #длина временного ряда
## [1] 1216
mean(spy_ret) #средняя лог доходность
## [1] 0.0003204901
sd(spy_ret) #среднеквадратическое отклонение лог доходности
## [1] 0.008633183
max(spy_ret) #максимальное значение лог доходности
## [1] 0.04928991
min(spy_ret) #минимальное значение лог доходности
## [1] -0.04301906
skewness(spy_ret) #коэффициент асимметрии лог доходности
## [1] -0.5309839
kurtosis(spy_ret) #коэффициент эксцесса лог доходности
## [1] 6.779313
#заменить spy_ret на ewh_ret для аналогичных показателей EWH

Проводим статистические тесты

#тест Дики-Фуллера на стационарность
adf.test(spy_ret)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  spy_ret
## Dickey-Fuller = -11.285, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary
adf.test(ewh_ret)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  ewh_ret
## Dickey-Fuller = -10.623, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary
#LM-тест на наличие ARCH-эффекта
ArchTest(spy_ret, lags = 1)
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  spy_ret
## Chi-squared = 90.107, df = 1, p-value < 2.2e-16
ArchTest(ewh_ret, lags = 1)
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  ewh_ret
## Chi-squared = 50.615, df = 1, p-value = 1.124e-12
#Тест Колмогорова-Смирнова на нормальность распределения
ks.test(spy_ret, "pnorm", c(mean(spy_ret), sd(spy_ret)))
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  spy_ret
## D = 0.49025, p-value < 2.2e-16
## alternative hypothesis: two-sided
ks.test(ewh_ret, "pnorm", c(mean(ewh_ret), sd(ewh_ret)))
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  ewh_ret
## D = 0.48668, p-value < 2.2e-16
## alternative hypothesis: two-sided

Временные ряды являются стационарными, в них присутствует ARCH-эффект, а их распределение отлично от нормального. Все это говорит в пользу применения моделей условной гетероскедастичности. Перейдем к оценке модели GARCH(1,1) для каждого временного ряда.

Оценка одномерных моделей временных рядов

Оценим параметры модели GARCH(1,1) для логарифмической доходности SPY. В качестве распределения стандартизированных инноваций используем распределение Стьюдента.

garch_spy_spec <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1,1), external.regressors = NULL), 
                           mean.model = list(armaOrder = c(0,0), external.regressors = NULL),distribution.model = "std", start.pars = list(), fixed.pars = list())
garch_spy <- ugarchfit(spec = garch_spy_spec, data = spy_ret, solver.control = list(trace=0))
print(coef(garch_spy))
##           mu        omega       alpha1        beta1        shape 
## 7.232677e-04 2.381682e-06 2.004519e-01 7.876103e-01 4.971209e+00

Аналогичным образом оценим параметры модели GARCH(1,1) для лоагрифмической доходности EWH.

garch_ewh_spec <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1,1), external.regressors = NULL), 
                            mean.model = list(armaOrder = c(0,0), external.regressors = NULL), distribution.model = "std", start.pars = list(), fixed.pars = list())
garch_ewh <- ugarchfit(spec = garch_ewh_spec, data = ewh_ret, solver.control = list(trace=0))
print(coef(garch_ewh))
##           mu        omega       alpha1        beta1        shape 
## 6.388528e-04 4.399946e-06 1.147307e-01 8.521113e-01 6.718380e+00

Проверим качество модели GARCH(1,1) для каждого временного ряда. В стандартизированных остатках модели должен отсутствовать ARCH-эффект и автокорреляция, а теоретическое распределение стандартизированных инноваций должно совпадать с эмпирическим распределением остатков модели.

#вычисляем стандартизированные остатки и квадраты стандартизированных остатков моделей
spy_resid <- as.timeSeries(residuals(garch_spy)/sigma(garch_spy))
spy_resid_sq <- spy_resid**2
ewh_resid <- as.timeSeries(residuals(garch_ewh)/sigma(garch_ewh))
ewh_resid_sq <- ewh_resid**2
#Коррелограмма стандартизированных остатков SPY
acfPlot(spy_resid)

#Коррелограмма стандартизированных остатков EWH
acfPlot(ewh_resid)

Можно видеть, что автокорреляция в стандартизированных остатках отсутствует. Проведем LM-тест на наличие ARCH-эффекта в стандартизированных остатках модели.

ArchTest(spy_resid, lags = 10)
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  spy_resid
## Chi-squared = 7.9596, df = 10, p-value = 0.6328
ArchTest(ewh_resid, lags = 10)
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  ewh_resid
## Chi-squared = 6.9516, df = 10, p-value = 0.73

ARCH-эффект в стандартизированных остатках отсутствует. Проведем тест Колмогорова-Смирнова для определения того, есть ли значимые различия между предполагаемым теоретическим распределением инноваций и эмпирическим распределением ошибок.

ks.test.t(spy_resid)
## 
##  One-sample Kolmogorov-Smirnov test student-t with df=4.97,
##  location=-0.02, scale=0.77
## 
## data:  spy_resid
## D = 0.021262, p-value = 0.6416
## alternative hypothesis: two-sided
ks.test.t(ewh_resid)
## 
##  One-sample Kolmogorov-Smirnov test student-t with df=6.8,
##  location=-0.03, scale=0.84
## 
## data:  ewh_resid
## D = 0.016632, p-value = 0.8896
## alternative hypothesis: two-sided

Распределение Стьюдента является адекватным выбором распределения стандартизированных инноваций GARCH(1,1) для каждого временного ряда. По итогу проведенных тестов можно сделать вывод, что модель GARCH(1,1) является адекватной моделью для каждого временного ряда.

Моделирование совместного распределения временных рядов при помощи копулы

Стандартизированные остатки одномерных моделей при помощи метода обратной функции трансформируются в шкалу [0,1]

spy_resid_01 <- pobs(spy_resid)
ewh_resid_01 <- pobs(ewh_resid)
plot(spy_resid_01,ewh_resid_01, col = "blue", xlab = "Трансформированные стандартизированные остатки SPY",
     ylab = "Трансформированные стандартизированные остатки EWH",
     main = "Диаграмма рассеивания трансформированных стандартизированных остатков")

Копула, наиболее удачно описывающая зависимость временных рядов, выбирается на основе наименьшего значения критерия Акаике

selectedCopula <- BiCopSelect(spy_resid_01, ewh_resid_01, familyset = 1:5)
print(selectedCopula)
## Bivariate copula: t (par = 0.61, par2 = 10.4, tau = 0.42)

Моделирование стандартизированных остатков из копулы Стьюдента

rho <- as.numeric(selectedCopula[2])
copula_dist <- mvdc(copula = tCopula(dim = 2, rho), margins = c("t", "t"), paramMargins = list(list(df = 4.97), list(df = 6.8)))
simulated_residuals <- rMvdc(10000, copula_dist)
plot(simulated_residuals, col = "blue", main = "Стандартизированные остатки, смоделированные из копулы Стьюдента", xlab = "Стандартизированные остатки SPY", ylab = "Стандартизированные остатки EWH")

Моделирование доходностей ETF при помощи копула-GARCH модели

При помощи копула-GARCH модели было смоделировано 10 000 доходностей каждого ETF на следующий торговый день.

sim_resid_spy <- matrix(simulated_residuals[,1], ncol = 10000)
sim_resid_ewh <- matrix(simulated_residuals[,2], ncol = 10000)
sim_result_spy <- ugarchsim(garch_spy, n.sim = 1, m.sim = 10000, 
                           custom.dist = list(name = "sample", distfit = sim_resid_spy))
sim_result_ewh <- ugarchsim(garch_ewh, n.sim = 1, m.sim = 10000, 
                           custom.dist = list(name = "sample", distfit = sim_resid_ewh))
sim_returns_spy <-as.timeSeries(matrix(fitted(sim_result_spy), ncol = 1))
sim_returns_ewh <-as.timeSeries(matrix(fitted(sim_result_ewh), ncol = 1))
portfolio_ret <- as.data.frame(matrix(c(sim_returns_spy, sim_returns_ewh), ncol = 2))
plot(portfolio_ret, col = "blue", xlab = "SPY", ylab = "EWH", main = "Смоделированные доходности ETF на 1 день вперед")

Построение портфеля с минимальным значением CVaR

На основе смоделированных из копула-GARCH модели доходностей был построен портфель с минимальным значением меры риска CVaR. Ниже представлена структура данного портфеля.

globminSpec <- portfolioSpec()
setType(globminSpec) <- "CVaR"
## Solver set to solveRquadprog
## setSolver: solveRglpk
setSolver(globminSpec) <- "solveRglpk.CVAR"
minriskPortfolio(data = as.timeSeries(portfolio_ret), spec = globminSpec, constraints = "LongOnly")
## 
## Title:
##  CVaR Minimum Risk Portfolio 
##  Estimator:         covEstimator 
##  Solver:            solveRglpk.CVAR 
##  Optimize:          minRisk 
##  Constraints:       LongOnly 
##  VaR Alpha:         0.05 
## 
## Portfolio Weights:
##     V1     V2 
## 0.2082 0.7918 
## 
## Covariance Risk Budgets:
##     V1     V2 
## 0.2073 0.7927 
## 
## Target Returns and Risks:
##   mean    Cov   CVaR    VaR 
## 0.0007 0.0130 0.0284 0.0199 
## 
## Description:
##  Wed Jun 02 12:54:47 2021 by user: ars99

Кроме того, портфель с минимальным CVaR был построен с использованием исторических логарифмических доходностей, который уже были получены ранее. Структура данного портфеля также представлена ниже.

portfolio_ret_hist <- as.timeSeries(matrix(c(spy_ret, ewh_ret), ncol = 2))
minriskPortfolio(data = portfolio_ret_hist, spec = globminSpec, constraints = "LongOnly")
## 
## Title:
##  CVaR Minimum Risk Portfolio 
##  Estimator:         covEstimator 
##  Solver:            solveRglpk.CVAR 
##  Optimize:          minRisk 
##  Constraints:       LongOnly 
##  VaR Alpha:         0.05 
## 
## Portfolio Weights:
##   SS.1   SS.2 
## 0.6327 0.3673 
## 
## Covariance Risk Budgets:
##   SS.1   SS.2 
## 0.5948 0.4052 
## 
## Target Returns and Risks:
##   mean    Cov   CVaR    VaR 
## 0.0002 0.0085 0.0211 0.0137 
## 
## Description:
##  Wed Jun 02 12:54:48 2021 by user: ars99