本文是时间序列分析课程的作业，基于R、Rnw和Latex进行编写。
   GARCH代码实现主要参考自《经济与金融计量方法：原理、应用案例及R语言实现》和对应包的官方文档，代码进一步整合，但每次执行时可能需要较长的时间，建议执行完后将结果导出成excel。如果本文存在问题，随时欢迎交流~
 
 
一、数据来源
 
  沪深300指数，是由沪深证券交易所于 2005 年 4 月 8 日联合发布的反映沪深 300 指数编制目标和运行状况的金融指标，并能够作为投资业绩的评价标准，为指数化投资和指数衍生产品创新提供基础条件。因此，本次数据来源于网易财经，研究的数据集对象是沪深 300 指数（股票代码为000300），此次分析选取了沪深 300 指数2000 年1月到2019年12月的工作日收盘价格数据。
 
二、数据分析
 
（一）时序图
 
  为了分析数据的波动情况，对其进行对数化和差分得到对数收益率，下图为沪深 300 指数的收盘价时序图和对数收益率时序图。
 
 
 
（二）平稳性检验
 
  Augmented Dickey-Fuller Test (ADF) 是 DF 检验的拓展形式，可以对存在高阶滞后的序列进行单位根检验，原假设存在单位根，即序列不平稳。本文使用adf.test()进行单位根检验，检验结果如下所示，p 值远小于 0.01 说明拒绝原假设，即序列是平稳的。
 
 
三、模型建立
 
（一）均值模型
 
1. 均值模型的识别
 
  序列平稳后，使用auto.arima()对序列自动识别均值模型。
 
识别出来的模型为ARMA(4, 4)。经过模型识别后，对模型$ARMA(4,4)$进行参数显著性检验。检验结果发现部分参数不显著，采用建立疏系数的均值模型，将不显著的参数强制为0。
 
采用疏系数的$ARMA$模型后的参数显著性检验结果来看，非零参数结果都显著。（注：这一步存有疑问，后续依旧是采用$ARMA(4,4)$的非疏系数模型进行拟合）
 
 
（二）方差模型建立
 
1. ARCH效应检验
 
  上述建立模型后，对残差进行ARCH效应检验。Ljung-Box统计量$Q(m)$对残差序列进行自相关检验。原假设是序列不存在自相关，在残差的平方序列中可以检验条件异方差。
 
使用MTS包中的archTest()进行检验
 检验结果显示，滞后10阶和滞后20阶的残差序列存在自相关，因此拒绝原假设，残差序列存在ARCH效应。
## [1] "m = 10"
## Q(m) of squared series(LM test):
## Test statistic: 1043.197 p-value: 0
## Rank-based Test:
## Test statistic: 849.5804 p-value: 0
## [2] "m = 20"
## Q(m) of squared series(LM test):
## Test statistic: 1705.592 p-value: 0
## Rank-based Test:
## Test statistic: 1565.588 p-value: 0
 
 
2. 标准GARCH模型建立
 
  上述ARCH效应表明，条件方差是依赖于过去值。因此可以考虑GARCH模型对方差方程进行参数估计。
 
使用tseries包中的garch()函数进行拟合标准GARCH模型。
 从结果上看，拟合出来的参数都显著，Box-Ljung test结果中的P值大于显著性，因此可以认为模型的残差无序列相关，说明该模型拟合效果较好。但实际上，其中Jarque Bera Test用于对回归残差的正态性进行检验，Shapiro - Wilk Normality Test也可以用于正态性检验，原假设都是是残差序列服从正态分布，检验结果表明，残差序列是不服从正态分布，因此可以对模型进行优化，考虑其他GARCH模型。
## Call:
## garch(x = r.data, order = c(1, 1))
##
## Model:
## GARCH(1,1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.91948 -0.53035 0.04107 0.57992 5.60582
##
## Coefficient(s):
## Estimate Std. Error t value Pr(>|t|)
## a0 0.011014 0.002491 4.422 9.78e-06 ***
## a1 0.063407 0.004076 15.556 < 2e-16 ***
## b1 0.934953 0.003839 243.530 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Diagnostic Tests:
## Jarque Bera Test
##
## data: Residuals
## X-squared = 659.37, df = 2, p-value < 2.2e-16
##
##
## Box-Ljung test
##
## data: Squared.Residuals
## X-squared = 1.5224, df = 1, p-value = 0.2173
 
 
四、模型优化
 
（一）模型拟合
 
  上述的模型考虑的是$ARIMA(0,0,0)$的标准$GARCH(1,1)$模型，即均值模型的参数均设置为0。而从均值模型的分析来看，可以拟合$ARIMA(4,0,4)$的均值模型与$GARCH(1,1)$的方差模型。而上述的正态性检验结果表明，残差的分布不适用标准正态分布，应该考虑其他类型的分布。
 
使用$rugarch$包对多个模型进行拟合。
对于均值模型，考虑不带截距项的$ARIMA(4,0,4)$。
对于方差模型，阶数设定1阶ARCH和1阶GARCH，考虑标准GARCH（$sGARCH$）、指数GARCH（$eGARCH$）、$GJR-GARCH$、渐近幂ARCH（$APARCH$）、门限GARCH（$TGARCH$）、非线性非对称GARCH（$NAGARCH$）六类模型。
对于残差分布类型，考虑标准正态分布（$norm$）、标准t分布（$std$）、偏t分布（$sstd$）、广义误差分布（$ged$）和Johnson’SU分布（$jsu$）五类分布。
 
 
（二）模型比较
 
  从拟合的模型结果来看，非对称类型的指数GARCH模型在最大似然估计值达到最大，同时AIC和BIC都能达到最小，说明指数GARCH模型比较好。从拟合的模型残差分布来看，非正态分布的AIC和BIC都明显低于正态分布，说明残差是服从重尾类型的分布。
 
选取指数GARCH模型，对比不同分布的参数显著性、检验结果以及模型效果，这里分布考虑正态分布、t分布、广义误差分布与三种分布对应的偏态分布以及Johnson’SU分布，结果如下所示。
 
 
 
 
 
 
 
  从模型拟合的系数显著性表可以看出，拟合的系数基本都显著，仅有少数参数不显著。从系数稳定性的个别检验和联合检验可以看出，在5%显著性水平下正态分布、偏正态分布、广义误差分布和偏广义误差分布都接收参数是稳定的原假设。从符号偏误检验的结果来看，指数GARCH模型正负残差的受到冲击的差异不明显，说明该非对称模型有效消除了杠杆效应。从调整皮尔逊拟合优度检验的结果来看，原假设是残差分布与理论分布没有差异，结果表明广义误差分布和偏广义误差分布是没有拒绝原假设，说明这两种分布与模型适配较好。
 
（三）模型选择
 
  从模型拟合的信息准则表来看，偏广义误差分布的LLH达到最大，而BIC和HQ都达到最小，说明$ARMA(4,4)-eGARCH(1,1)-SGED$模型最优。拟合$ARMA(4,4)-EGARCH(1,1)-SGED$模型，理论模型和模型参数估计及显著性如下所示。
 $$\begin{array}{rl}{X}_t& =\sum _{i=1}^p{\phi }_i{X}_{t-i}-\sum _{j=1}^q{\vartheta }_j{\epsilon }_{t-j}+{\epsilon }_t\\ z_t& =\frac{{\epsilon }_t}{{\sigma }_t}=\frac{X_t}{{\sigma }_t}\\ ln{\sigma }_t^2& =\omega +\alpha z_{t-1}+\beta ln{\sigma }_{t-1}^2+\gamma (\frac{|{\epsilon }_{t-1}|}{{\sigma }_{t-1}}-E|z_{t-1}|)\\ 其中，& {\epsilon }_t\sim SGED(0,1,shape,skew)\end{array}$$
 
 
五、结论
 
  本文通过对沪深300指数的波动性分析发现，我国股票市场有两段时间出现较大的波动。第一次波动出现在2008年前后，这段期间为全球金融危机，明显可以看出收盘价时序图出现明显的陡峭的波峰，持续时间较长；第二次波动出现在2015年前后，该阶段是由于杠杆资金的加入和政策收紧，形成短暂的牛市和熊市，波动程度不亚于2008年金融危机，但持续时间比较短。
   本文使用多个GARCH模型进行比对，发现非对称模型$ARIMA-EGARCH$模型与沪深300指数的对数收益率有较高的匹配度，同时偏广义误差分布与理论分布较为接近，说明沪深300指数的波动性呈现的是尖峰厚尾、非对称性特征，说明我国股票市场存在杠杆效应。
 
六、实现代码
 
# Created by: Enguanei
# Created on: 2021/5/7

# 导入包
library(knitr)
library(pedquant)
library(zoo)
library(imputeTS)
library(tseries)
library(forecast)
library(MTS)
library(rugarch)
library(dplyr)

# 爬取网易财经的数据
datt <- md_stock(symbol = '000300.ss',
                 from = "2002-01-01",
                 to = "2019-12-31", source = "163")

# 查看数据
head(datt$`000300.ss`$close)
head(datt$`000300.ss`$date)

# 保存数据
write.csv(datt, "D:\\Rproject\\hs300.csv")

# 导入数据
raw_data <- read.csv("D:\\Rproject\\hs300.csv", header = TRUE)

# 转成时间序列
my_date <- as.Date(raw_data$X000300.ss.date)
my_data <- raw_data$X000300.ss.close
data.ts <- zoo(my_data, my_date)

# 查看是否存在缺失值
sum(is.na(my_data))

# 转成时间序列（不包括日期）和对数差分化处理
r.data <- diff(log(ts(my_data))) * 100

# 画出时序图
par(mfrow = c(2, 1))
plot(data.ts, xlab = 'Year', ylab = 'Close.Price',
     main = 'Close of CSI300')
plot(r.data, ylab = 'Log.Return',
     main = 'Log.Return of CSI300')

# 平稳性检验
show(adf.test(r.data))

# 拟合arma均值模型
md1 <- auto.arima(r.data)
md1

# 系数显著性检验
t <- md1$coef / sqrt(diag(md1$var.coef))
p_1 <- 2 * (1 - pnorm(abs(t)))
kable(rbind(p = p_1), caption = "Result of Coef. Test")

# 拟合疏系数arma
mean_md_1 <- Arima(r.data,
                   order = c(4, 0, 4),
                   fixed = c(NA, 0, 0, NA, NA, 0, 0, NA, 0))

# 疏系数模型的系数显著性检验
t <- mean_md_1$coef / sqrt(diag(mean_md_1$var.coef))
p_2 <- 2 * (1 - pnorm(abs(t)))
kable(rbind(p = p_2[p_2 != 1]), caption = "Result of Coef. Test")

# ARCH效应检验（原假设是没有ARCH效应）
print(paste0('m = 10'))
archTest(mean_md_1$res, lag = 10)
print(paste0('m = 20'))
archTest(mean_md_1$res, lag = 20)

# 使用tseries拟合sGARCH（只能拟合方差模型，没有考虑均值模型）
out <- garch(r.data, order = c(1, 1))
summary(out)

# 模型优化
# 模型选择
md_name <- c("sGARCH", "eGARCH", "gjrGARCH",
             "apARCH", "TGARCH", "NAGARCH")
md_dist <- c("norm", "std", "sstd", "ged", "jsu")

# 设定函数
model_train <- function(data_, m_name, m_dist) {
  res <- NULL
  for (i in m_name) {
    sub_model <- NULL
    main_model <- i
    if (i == "TGARCH" | i == "NAGARCH") {
      sub_model <- i
      main_model <- "fGARCH"
    }

    for (j in m_dist) {
      model_name <- paste0('ARMA(4,4)', '-',
                           i, '-', j,
                           sep = "")
      if (main_model == "fGARCH") {
        model_name <- paste0('ARMA(4,4)', '-',
                             sub_model, '-', j,
                             sep = "")
      }
      print(model_name)
      # 模型设定
      mean.spec <- list(armaOrder = c(4, 4), include.mean = F,
                        archm = F, archpow = 1, arfima = F,
                        external.regressors = NULL)
      var.spec <- list(model = main_model, garchOrder = c(1, 1),
                       submodel = sub_model,
                       external.regressors = NULL,
                       variance.targeting = F)
      dist.spec <- j
      my_spec <- ugarchspec(mean.model = mean.spec,
                            variance.model = var.spec,
                            distribution.model = dist.spec)
      # 模型拟合
      my_fit <- ugarchfit(data = data_, spec = my_spec)
      res <- rbind(res, c(list(ModelName = model_name),
                          list(LogL = likelihood(my_fit)),
                          list(AIC = infocriteria(my_fit)[1]),
                          list(BIC = infocriteria(my_fit)[2]),
                          list(RMSE = sqrt(mean(residuals(my_fit)^2))),
                          list(md = my_fit)))
    }
  }
  return(res)
}

# 模型拟合
my_model <- model_train(r.data, md_name, md_dist)

# 保存模型结果
write.table(my_model[, 1:5],
            "D:\\Rproject\\result.txt",
            sep = " ")

# 读取模型结果
res_table <- read.table("D:\\Rproject\\result.txt",
                        header = TRUE,
                        sep = " ")
kable(res_table, caption = "Result of GARCH Model", align = 'l')

# 选择模型再拟合
md_name <- c("eGARCH")
md_dist <- c("norm", "snorm", "std", "sstd", "ged", "sged", "jsu")
new_model <- model_train(r.data, md_name, md_dist)

# 信息准则表
my_info_cri <- NULL
for (i in 1:dim(new_model)[1]) {
  new_aic <- round(infocriteria(new_model[i, 6]$md)[1], digits = 4)
  new_bic <- round(infocriteria(new_model[i, 6]$md)[2], digits = 4)
  new_sib <- round(infocriteria(new_model[i, 6]$md)[3], digits = 4)
  new_hq <- round(infocriteria(new_model[i, 6]$md)[4], digits = 4)
  new_llk <- round(likelihood(new_model[i, 6]$md), digits = 4)
  my_info_cri <- rbind(my_info_cri, c(new_model[i, 1],
                                      Akaike = new_aic,
                                      Bayes = new_bic,
                                      Shibata = new_sib,
                                      HQ = new_hq,
                                      LLH = new_llk))
}


# 模型系数显著性表
my_cof_table <- NULL
my_cof_names <- NULL
for (i in 1:dim(new_model)[1]) {
  my_cof_names <- rbind(my_cof_names, new_model[i, 1]$ModelName)
  my_df <- data.frame(t(round(new_model[i, 6]$md@fit$robust.matcoef[, 4],
                              digits = 4)))
  if (i == 1) {
    my_cof_table <- data.frame(my_df)
  }
  else {
    my_cof_table <- full_join(my_cof_table, my_df)
  }
}
row.names(my_cof_table) <- my_cof_names[, 1]
my_cof_table <- t(my_cof_table)


# 参数稳定性个别检验表
my_nyblom_table <- NULL
my_nyblom_name <- NULL
for (i in 1:dim(new_model)[1]) {
  my_nyblom_name <- rbind(my_nyblom_name,
                          new_model[i, 1]$ModelName)
  my_df <- data.frame(t(round(nyblom(new_model[i, 6]$md)$IndividualStat,
                              digits = 4)))
  if (i == 1) {
    my_nyblom_table <- data.frame(my_df)
  }
  else {
    my_nyblom_table <- full_join(my_nyblom_table, my_df)
  }
}
row.names(my_nyblom_table) <- my_nyblom_name[, 1]
my_IC <- NULL
for (i in 1:dim(new_model)[1]) {
  my_IC <- cbind(my_IC, nyblom(new_model[1, 6]$md)$IndividualCritical)
}
my_nyblom_table <- rbind(t(my_nyblom_table), my_IC)

# 参数稳定性联合检验表
my_nyblom_table2 <- NULL
my_nyblom_name2 <- NULL
for (i in 1:dim(new_model)[1]) {
  my_nyblom_name2 <- rbind(my_nyblom_name2,
                           new_model[i, 1]$ModelName)
  df1 <- data.frame(JoinStat = round(nyblom(new_model[i,6]$md)$JointStat,
                                     digits = 4))
  df2 <- data.frame(mm = round(nyblom(new_model[i, 6]$md)$JointCritical,
                               digits = 4))
  temp_table <- cbind(df1, t(df2))
  if (i == 1) {
    my_nyblom_table2 <- temp_table
  }
  else {
    my_nyblom_table2 <- full_join(my_nyblom_table2, temp_table)
  }
}
row.names(my_nyblom_table2) <- my_nyblom_name2

# 模型分布拟合度p值检验表
my_gof_table <- NULL
my_gof_name <- NULL
group_name <- c(20, 30, 40, 50)
for (i in 1:dim(new_model)[1]) {
  my_gof_name <- rbind(my_gof_name,
                       new_model[i, 1]$ModelName)
  df1 <- round(gof(new_model[i, 6]$md,
                   groups = group_name)[1:4, 3],
               digits = 4)
  if (i == 1) {
    my_gof_table <- df1
  }
  else {
    my_gof_table <- rbind(my_gof_table, df1)
  }
}
row.names(my_gof_table) <- my_gof_name
colnames(my_gof_table) <- c(20, 30, 40, 50)

# 符号偏误显著性检验表
my_sign_table <- NULL
my_sign_name <- NULL
row_sign_name <- row.names(signbias(new_model[1, 6]$md))
for (i in 1:dim(new_model)[1]) {
  my_sign_name <- rbind(my_sign_name,
                        new_model[i, 1]$ModelName)
  df1 <- round(signbias(new_model[i, 6]$md)$prob, digits = 4)
  if (i == 1) {
    my_sign_table <- df1
  }
  else {
    my_sign_table <- rbind(my_sign_table, df1)
  }
}
row.names(my_sign_table) <- my_sign_name
colnames(my_sign_table) <- row_sign_name

# 展示表格
kable(my_cof_table,
      caption = "P-value Table of Coefficients")
kable(my_nyblom_table,
      caption = "P-value Table of Nyblom Individual Stability Test")
kable(my_nyblom_table2,
      caption = "P-value Table of Nyblom Joint Stability Test")
kable(my_sign_table,
      caption = "P-value Table of Sign Bias Test")
kable(my_gof_table,
      caption = "P-value Table of Goodness-of-fit")
kable(my_info_cri,
      caption = "Table of Information Criteria")

# 最终模型系数表
kable(new_model[6, 6]$md@fit$robust.matcoef,
      caption = "Table of Final Model Coef.")
 
七、参考资料
 
[1] Ruey S. Tsay, 李洪成, 尚秀芬,等. 金融数据分析导论[M]. 机械工业出版社, 2013.
 [2] 何宗武, 马卫锋. 经济与金融计量方法：原理、应用案例及R语言实现[M]. 机械工业出版社, 2019
 [3] 张东旭. 基于ARMA-GARCH模型族的上证指数收益率波动的实证分析[D]. 清华大学.

 
本文将说明金融数学中的R 语言优化投资组合，因子模型的实现和使用。
 
具有单一市场因素的宏观经济因素模型
 
我们将从一个包含单个已知因子(即市场指数)的简单示例开始。该模型为
 
其中显式因子ft为S＆P 500指数。我们将做一个简单的最小二乘(LS)回归来估计截距α和加载β：
 

 
大多数代码行用于准备数据，而不是执行因子建模。让我们开始准备数据：
 
# 设置开始结束日期和股票名称列表
 
begin_date 
 
end_date 
 
# 从YahooFinance下载数据
 
data_set 
 
for (stock_index in 1:length(stock_namelist))
 
data_set 
 
from = begin_date, to = end_date,
 
head(data_set)
 
#> AAPL AMD ADI ABBV AEZS A APD AA CF
 
#> 2016-01-04 98.74225 2.77 49.99239 49.46063 4.40 39.35598 107.89010 23.00764 35.13227
 
#> 2016-01-05 96.26781 2.75 49.62508 49.25457 4.21 39.22057 105.96097 21.96506 34.03059
 
#> 2016-01-06 94.38389 2.51 47.51298 49.26315 3.64 39.39467 103.38042 20.40121 31.08988
 
#> 2016-01-07 90.40047 2.28 46.30082 49.11721 3.29 37.72138 99.91463 19.59558 29.61520
 
#> 2016-01-08 90.87848 2.14 45.89677 47.77789 3.29 37.32482 99.39687 19.12169 29.33761
 
#> 2016-01-11 92.35001 2.34 46.98954 46.25827 3.13 36.69613 99.78938 18.95583 28.14919
 
head(SP500_index)
 
#> index
 
#> 2016-01-04 2012.66
 
#> 2016-01-05 2016.71
 
#> 2016-01-06 1990.26
 
#> 2016-01-07 1943.09
 
#> 2016-01-08 1922.03
 
#> 2016-01-11 1923.67
 
plot(SP500_index)
 

 
# 计算股票和SP500指数的对数收益率作为显式因子
 
X 
 
N 
 
T 
 
现在我们准备进行因子模型拟合。LS拟合很容易在R中实现，如下所示：
 

 
beta 
 
alpha 
 
sigma2 
 
print(alpha)
 
#> index
 
#> AAPL 0.0003999086
 
#> AMD 0.0013825599
 
#> ADI 0.0003609968
 
#> ABBV 0.0006684632
 
#> AEZS -0.0022091301
 
#> A 0.0002810616
 
#> APD 0.0001786375
 
#> AA 0.0006429140
 
#> CF -0.0006029705
 
print(beta)
 
#> index
 
#> AAPL 1.0957919
 
#> AMD 2.1738304
 
#> ADI 1.2683047
 
#> ABBV 0.9022748
 
#> AEZS 1.7115761
 
#> A 1.3277212
 
#> APD 1.0239453
 
#> AA 1.8593524
 
#> CF 1.5702493
 
或者，我们可以使用矩阵表示法进行拟合，我们定义和扩展因子。然后最小化
 
t(X) %*% F_ %*% solve(t(F_) %*% F_)
 
#> alpha beta
 
#> AAPL 0.0003999086 1.0957919
 
#> AMD 0.0013825599 2.1738304
 
#> ADI 0.0003609968 1.2683047
 
#> ABBV 0.0006684632 0.9022748
 
#> AEZS -0.0022091301 1.7115761
 
#> A 0.0002810616 1.3277212
 
#> APD 0.0001786375 1.0239453
 
#> AA 0.0006429140 1.8593524
 
#> CF -0.0006029705 1.5702493
 
E 
 
另外，我们可以简单地使用R为我们完成工作：
 
cbind(alpha = factor_model$alpha, beta = factor_model$beta)
 
#> alpha index
 
#> AAPL 0.0003999086 1.0957919
 
#> AMD 0.0013825599 2.1738304
 
#> ADI 0.0003609968 1.2683047
 
#> ABBV 0.0006684632 0.9022748
 
#> AEZS -0.0022091301 1.7115761
 
#> A 0.0002810616 1.3277212
 
#> APD 0.0001786375 1.0239453
 
#> AA 0.0006429140 1.8593524
 
#> CF -0.0006029705 1.5702493
 
可视化协方差矩阵
 
有趣的是，可视化对数收益率[算术处理误差]以及残差Ψ的估计协方差矩阵。让我们从对数收益率的协方差矩阵开始：
 
main = "单因子模型对数收益的协方差矩阵")
 

 
我们可以观察到所有股票都是高度相关的，这是市场因素的影响。为了检查股票相关关系，我们绘制相关图：
 
plot(cov2cor(Psi),
 
main = "残差协方差矩阵")
 

 
cbind(stock_namelist, sector_namelist) # 股票的行业
 
#> stock_namelist sector_namelist
 
#> [1,] "AAPL" "Information Technology"
 
#> [2,] "AMD" "Information Technology"
 
#> [3,] "ADI" "Information Technology"
 
#> [4,] "ABBV" "Health Care"
 
#> [5,] "AEZS" "Health Care"
 
#> [6,] "A" "Health Care"
 
#> [7,] "APD" "Materials"
 
#> [8,] "AA" "Materials"
 
#> [9,] "CF" "Materials"
 
有趣的是，我们可以观察到对Ψ执行的自动聚类可以正确识别股票的行业。
 
评估投资资金
 
在此示例中，我们将基于因子模型评估几种投资基金的绩效。我们将标准普尔500指数作为明确的市场因素，并假设无风险收益为零 rf = 0。特别是，我们考虑六种交易所买卖基金(ETF)：
 
我们首先加载数据：
 
# 设置开始结束日期和股票名称列表
 
begin_date 
 
end_date 
 
# 从YahooFinance下载数据
 
data_set 
 
for (stock_index in 1:length(stock_namelist))
 
data_set 
 
head(data_set)
 
#> SPY XIVH SPHB SPLV USMV JKD
 
#> 2016-10-03 203.6610 29.400 31.38322 38.55683 42.88382 119.8765
 
#> 2016-10-04 202.6228 30.160 31.29729 38.10687 42.46553 119.4081
 
#> 2016-10-05 203.5195 30.160 31.89880 38.02249 42.37048 119.9421
 
#> 2016-10-06 203.6610 30.160 31.83196 38.08813 42.39899 120.0826
 
#> 2016-10-07 202.9626 30.670 31.58372 37.98500 42.35146 119.8296
 
#> 2016-10-10 204.0197 31.394 31.87970 38.18187 42.56060 120.5978
 
head(SP500_index)
 
#> index
 
#> 2016-10-03 2161.20
 
#> 2016-10-04 2150.49
 
#> 2016-10-05 2159.73
 
#> 2016-10-06 2160.77
 
#> 2016-10-07 2153.74
 
#> 2016-10-10 2163.66
 
# 计算股票和SP500指数的对数收益率作为显式因子
 
X 
 
N 
 
T 
 
现在我们可以计算所有ETF的alpha和beta：
 
#> alpha beta
 
#> SPY 7.142225e-05 1.0071424
 
#> XIVH 1.810392e-03 2.4971086
 
#> SPHB -2.422107e-04 1.5613533
 
#> SPLV 1.070918e-04 0.6777149
 
#> USMV 1.166177e-04 0.6511667
 
#> JKD 2.569578e-04 0.8883843
 
现在可以进行一些观察：
 
SPY是S＆P 500的ETF，如预期的那样，其alpha值几乎为零，beta值几乎为1： α= 7.142211×10-5和 β= 1.0071423。
 
XIVH是具有高alpha值的ETF，计算出的alpha值是ETF中最高的(高1-2个数量级)： α= 1.810392×10-3。
 
SPHB是一种ETF，据推测具有很高的beta，而计算出的beta却是最高的，但不是最高的：β= 1.5613531。有趣的是，计算出的alpha为负，因此，该ETF应谨慎。
 
SPLV是降低波动性的ETF，实际上，计算得出的beta偏低：β= 0.6777072。
 
USMV还是降低波动性的ETF，实际上，计算出的beta是最低的：β= 0.6511671。
 
JKD显示出很好的折衷。
 
我们可以使用一些可视化：
 
barplot(rev(alpha), horiz = TRUE, main = "alph
 

 
我们还可以使用例如Sharpe比率，以更系统的比较不同的ETF。回顾一种资产和一个因素的因子模型
 
我们获得
 

 
夏普比率如下：
 

 
假设。因此，基于Sharpe比率对不同资产进行排名的一种方法是根据α/β比率对它们进行排名：
 
print(ranking)
 
#> alpha/beta SR alpha beta
 
#> XIVH 7.249952e-04 0.13919483 1.810392e-03 2.4971086
 
#> JKD 2.892417e-04 0.17682677 2.569578e-04 0.8883843
 
#> USMV 1.790904e-04 0.12280053 1.166177e-04 0.6511667
 
#> SPLV 1.580189e-04 0.10887903 1.070918e-04 0.6777149
 
#> SPY 7.091574e-05 0.14170591 7.142225e-05 1.0071424
 
#> SPHB -1.551287e-04 0.07401566 -2.422107e-04 1.5613533
 
可以看到：
 
就α/β而言，XIVH最佳(α最大)，而SPHB最差(α负)。
 
就夏普比率(更确切地说，是信息比率，因为我们忽略了无风险利率)而言，JDK是最好的，其次是SPY。这证实了大多数投资基金的表现不超过市场的观点。
 
显然，无论以哪种衡量标准，SPHB都是最差的：负α，负β比率和Sharpe比率。
 
JDK之所以能够取得最佳性能，是因为它的alpha值很好(尽管不是最好的)，而同时具有0.88的中等beta值。
 
XIVH和SPHB有大量不同的beta，因此在市场上具有极端敞口。
 
USMV在市场上的曝光率最小，有可接受的alpha值，并且其Sharpe比率接近第二和第三高的位置。Fama-French三因子模型
 
该示例将说明使用标准普尔500指数中的九种股票的Fama-French三因子模型。让我们从加载数据开始：
 
# 设置开始结束日期和股票名称列表
 
begin_date 
 
end_date 
 
# 从YahooFinance下载数据
 
data_set 
 
for (stock_index in 1:length(stock_namelist))
 
data_set 
 
# 下载Fama-French因子
 
head(fama_lib)
 
#> Mkt.RF SMB HML
 
#> 1926-07-01 0.10 -0.24 -0.28
 
#> 1926-07-02 0.45 -0.32 -0.08
 
#> 1926-07-06 0.17 0.27 -0.35
 
#> 1926-07-07 0.09 -0.59 0.03
 
#> 1926-07-08 0.21 -0.36 0.15
 
#> 1926-07-09 -0.71 0.44 0.56
 
tail(fama_lib)
 
#> Mkt.RF SMB HML
 
#> 2017-11-22 -0.05 0.10 -0.04
 
#> 2017-11-24 0.21 0.02 -0.44
 
#> 2017-11-27 -0.06 -0.36 0.03
 
#> 2017-11-28 1.06 0.38 0.84
 
#> 2017-11-29 0.02 0.04 1.45
 
#> 2017-11-30 0.82 -0.56 -0.50
 
# 计算股票的对数收益率和Fama-French因子
 
X 
 
N 
 
现在我们在矩阵F中具有三个因子，并希望拟合模型，其中现在的载荷是一个beta矩阵：。我们可以做最小二乘拟合，最小化。更方便地，我们定义和扩展因子 。然后可以将LS公式写为最小化
 
print(Gamma)
 
#> alpha b1 b2 b3
 
#> AAPL 1.437845e-04 0.9657612 -0.23339130 -0.49806858
 
#> AMD 6.181760e-04 1.4062105 0.80738336 -0.07240117
 
#> ADI -2.285017e-05 1.2124008 0.09025928 -0.20739271
 
#> ABBV 1.621380e-04 1.0582340 0.02833584 -0.72152627
 
#> AEZS -4.513235e-03 0.6989534 1.31318108 -0.25160182
 
#> A 1.146100e-05 1.2181429 0.10370898 -0.20487290
 
#> APD 6.281504e-05 1.0222936 -0.04394061 0.11060938
 
#> AA -4.587722e-05 1.3391852 0.62590136 0.99858692
 
#> CF -5.777426e-04 1.0387867 0.48430007 0.82014523
 
另外，我们可以使用R完成：
 
#> alpha Mkt.RF SMB HML
 
#> AAPL 1.437845e-04 0.9657612 -0.23339130 -0.49806858
 
#> AMD 6.181760e-04 1.4062105 0.80738336 -0.07240117
 
#> ADI -2.285017e-05 1.2124008 0.09025928 -0.20739271
 
#> ABBV 1.621380e-04 1.0582340 0.02833584 -0.72152627
 
#> AEZS -4.513235e-03 0.6989534 1.31318108 -0.25160182
 
#> A 1.146100e-05 1.2181429 0.10370898 -0.20487290
 
#> APD 6.281504e-05 1.0222936 -0.04394061 0.11060938
 
#> AA -4.587722e-05 1.3391852 0.62590136 0.99858692
 
#> CF -5.777426e-04 1.0387867 0.48430007 0.82014523统计因子模型
 
现在让我们考虑统计因子模型或隐式因子模型，其中因子和载荷均不可用。调用具有 K因子的模型 XT =α1T+ BFT + ET的主成分方法：
 
PCA：样本均值：矩阵：样本协方差矩阵：特征分解：估计：更新特征分解：重复步骤2-3，直到收敛为止。
 
#> alpha
 
#> AAPL 0.0007074564 0.0002732114 -0.004631647 -0.0044814226
 
#> AMD 0.0013722468 0.0045782146 -0.035202146 0.0114549515
 
#> ADI 0.0006533116 0.0004151904 -0.007379066 -0.0053058139
 
#> ABBV 0.0007787929 0.0017513359 -0.003967816 -0.0056000810
 
#> AEZS -0.0041576357 0.0769496344 0.002935950 0.0006249473
 
#> A 0.0006902482 0.0012690079 -0.005680162 -0.0061507654
 
#> APD 0.0006236565 0.0005442926 -0.004229364 -0.0057976394
 
#> AA 0.0006277163 0.0027405024 -0.009796620 -0.0149177957
 
#> CF -0.0000573028 0.0023108605 -0.007409061 -0.0153425661
 
同样，我们可以使用R完成工作：
 
#> alpha factor1 factor2 factor3
 
#> AAPL 0.0007074564 0.0002732114 -0.004631647 -0.0044814226
 
#> AMD 0.0013722468 0.0045782146 -0.035202146 0.0114549515
 
#> ADI 0.0006533116 0.0004151904 -0.007379066 -0.0053058139
 
#> ABBV 0.0007787929 0.0017513359 -0.003967816 -0.0056000810
 
#> AEZS -0.0041576357 0.0769496344 0.002935950 0.0006249473
 
#> A 0.0006902482 0.0012690079 -0.005680162 -0.0061507654
 
#> APD 0.0006236565 0.0005442926 -0.004229364 -0.0057976394
 
#> AA 0.0006277163 0.0027405024 -0.009796620 -0.0149177957
 
#> CF -0.0000573028 0.0023108605 -0.007409061 -0.0153425661通过不同因子模型进行协方差矩阵估计的最终比较
 
我们最终将比较以下不同的因子模型：
 
样本协方差矩阵
 
宏观经济一因素模型
 
基本的三因素Fama-French模型
 
统计因素模型
 
我们在训练阶段估计模型，然后将估计的协方差矩阵与测试阶段的样本协方差矩阵进行比较。估计误差将根据PRIAL(平均损失提高百分比)进行评估：
 

 
加载训练和测试集：
 
# 设置开始结束日期和股票名称列表
 
begin_date 
 
end_date 
 
# 准备股票数据
 
data_set 
 
for (stock_index in 1:length(stock_namelist))
 
data_set 
 
# Fama-French 因子
 
mydata 
 
# 准备指数
 
f_SP500 
 
# 将数据拆分为训练数据和测试数据
 
T_trn 
 
X_trn 
 
X_tst 
 
现在让我们用训练数据估算不同的因子模型：
 
# 样本协方差矩阵
 
Sigma_SCM 
 
# 单因素模型
 
Gamma 
 
E 
 
# Fama-French三因子模型
 
Sigma_FamaFrench 
 
# 统计单因子模型
 
while (norm(Sigma - Sigma_prev, "F")/norm(Sigma, "F") > 1e-3) {
 
B 
 
# 统计三因子模型
 
K 
 
while (norm(Sigma - Sigma_prev, "F")/norm(Sigma, "F") > 1e-3) {
 
B 
 
Psi 
 
Sigma_PCA3 
 
# 统计五因子模型
 
K 
 
eigSigma 
 
while (norm(Sigma - Sigma_prev, "F")/norm(Sigma, "F") > 1e-3) {
 
B 
 
Psi 
 
最后，让我们比较测试数据中的不同估计：
 
Sigma_true 
 
barplot(error, main = "协方差矩阵估计误差",
 

 
PRIAL 
 
barplot(PRIAL, main = "协方差矩阵估计的先验方法",
 

 
最终可以看到使用因子模型进行协方差矩阵估计会有所帮助。
 

 
最受欢迎的见解
 
1.用机器学习识别不断变化的股市状况—隐马尔科夫模型(HMM)的应用
 
2.R语言GARCH-DCC模型和DCC(MVT)建模估计
 
3.R语言实现 Copula 算法建模依赖性案例分析报告
 
4.R语言COPULAS和金融时间序列数据VaR分析
 
5.R语言多元COPULA GARCH 模型时间序列预测
 
6.用R语言实现神经网络预测股票实例
 
7.r语言预测波动率的实现：ARCH模型与HAR-RV模型
 
8.R语言如何做马尔科夫转换模型markov switching model
 
9.matlab使用Copula仿真优化市场风险

 
原文链接：http://tecdat.cn/?p=24182
 
 
原文出处：拓端数据部落公众号
 
 
概要
 
本文用 R 编程语言极值理论 (EVT) 以确定 10 只股票指数的风险价值（和条件 VaR）。使用 Anderson-Darling 检验对 10 只股票的组合数据进行正态性检验，并使用 Block Maxima 和 Peak-Over-Threshold 的 EVT 方法估计 VaR/CvaR。最后，使用条件异向性 (GARCH) 处理的广义自回归来预测未来 20 天后指数的未来值。本文将确定计算风险因素的不同方法对模型结果的影响。
 
极值理论（最初由Fisher、Tippett和Gnedenko提出）表明，独立同分布（iid）变量样本的分块最大值的分布会收敛到三个极值分布之一。
 
最近，统计学家对极端值建模的兴趣又有了新的变化。极限值分析已被证明在各种风险因素的案例中很有用。在1999年至2008年的金融市场动荡之后，极值分析获得了有效性，与之前的风险价值分析不同。极限值代表一个系统的极端波动。极限值分析提供了对极端事件的概率、规模和保护成本的关系进行建模的能力。
 

  
 
 
参考
 
https://arxiv.org/pdf/1310.3222.pdf
https://www.ma.utexas.edu/mp_arc/c/11/11-33.pdf
http://evt2013.weebly.com/uploads/1/2/6/9/12699923/penalva.pdf
Risk Measurement in Commodities Markets Using Conditional Extreme Value Theory
 

  
 
第 1a 部分 - 工作目录、所需的包和会话信息
 
为了开始分析，工作目录被设置为包含股票行情的文件夹。然后，安装所需的 R 编程语言包并包含在包库中。R 包包括极值理论函数、VaR 函数、时间序列分析、定量交易分析、回归分析、绘图和 html 格式的包。
 
 
library(ggplot2)
library(tseries)
library(vars)
library(evd)
library(POT)
library(rugarch)
 

  
 
第 1b 节 - 格式化专有数据
 
用于此分析的第一个文件是“Data_CSV.csv”。该文件包含在 DAX 证券交易所上市的 15 家公司的股票代码数据，以及 DAX 交易所的市场投资组合数据。从这个数据文件中选出了 10 家公司，这些公司最近十年的股价信息是从谷歌财经下载的。
 
 
 
 
第 1c 节 - 下载股票代码数据
 
股票价格数据下载并读入 R 编程环境。收益率是用“开盘价/收盘价 ”计算的，十家公司的数据合并在一个数据框中，（每家公司一列）。
 
结果数据帧的每一行代表记录股价的 10 年中的一个工作日。然后计算数据帧中每一行的均值。一列 10 年的日期被附加到数据框。还创建了仅包含行均值和日期信息的第二个数据框。
 
alDat <- cbind(retursDaa, returnDta_A,
                 retrnsata_Ss, reunsataDB,
                 retunsDta_H, reurnsDta_S, rtunsDaaA,
                 retrnsaa_senus,reursDtaAlnz,
                 reurnsData_ailer)
 

  
 
第 2a 节 - 探索性数据分析
 
创建一个数据框统计表，其中包含每列（或公司）的最小值、中值、平均值、最大值、标准偏差、1% 分位数、5% 分位数、95% 分位数、99% 分位数。分位数百分比适用于极值。还创建了所有收益率均值的时间序列图表。
 

  taeSs<- c(min(x), medan(x), man(x),
                    max(x), sd(x), quntile(x, .01),
                    quanile(x, .05), qunile(x, .95),
                    quatile(x, .99), lngth(x))


 
 
 
 

第 2b 节 - 10 只股票指数的 VaR 估计
 
all_va.2 <- VAR(lDvarts, p = 2, tpe= "cnst")

# 预测未来125天、250天和500天
aDFva100 <- pdc(alDva.c, n.aea = 100, ci = 0.9)
 
 
 
 
为了开始估算数据所隐含的未来事件，我们进行了初步的风险值估算。首先，所有行的平均值和日期信息的数据框架被转换为时间序列格式，然后从这个时间序列中计算出风险值。根据VaR计算对未来100天和500天的价值进行预测。在随后的预测图中，蓝色圆圈代表未来100天的数值，红色圆圈代表500天的回报值。
 
plot（ap0$t$Tme[1:1200],
     alF_ar.d.$fst[1:1200]）
 
 
 
第 2c 部分 - 估计期望shortfall(ES),条件VAR(CvaR) 10 股票指数
 
为便于比较，计算了10只股票指数数据的条件风险值（CvaR或估计亏损）。首先，利用数据的时间序列，找到最差的0.95%的跌幅的最大值。然后，通过 "高斯 "方法计算出估计亏损，这两种计算的结果都以表格形式呈现。
 
ES(s(lD1:2528, 2, rp=FAE]),p=0.95, mho="gausn")
 
 
 
第 2d 节 - 10 只股票指数的希尔Hill估计
 
 
由于假设10股指数数据为重尾分布，数据极少变化，所以采用Hill Estimation对尾指数进行参数估计。目的是验证 10 只股票数据是否为极值分布。Hill Estimation 生成的图证实了。
 
hil(orvtis, otio="x", trt=15, nd=45)
 
 
 
 
 
第 2e 节 - 正态分布的 Anderson-Darling 检验
 
Anderson-Darling 检验主要用于分布族，是分布非正态性的决定因素。在样本量较大的情况下（如在 10 股指数中），小于 0.05 的 P 值表明分布与正态性不同。这是极值分布的预期。使用 Anderson-Darling 检验发现的概率值为 3.7^-24，因此证实了非正态性。
 
 
 
第 2f 节 - 结果表
 
最后，给出了10个股票指数未来价值的估计结果表。3 个 VaR 估计值（和估计差额）的点估计值和范围被制成表格以比较。
 

VaRES[3,] <- c("ES", etFbl[1], 4)
                        eSFbe[2], estFtbl[3],
                        rond(eSab[4], 4))
 
 
 
 
 
第 3a 节 - 10 个股票指数的 EVT 分块最大值估计
 
极值理论中的 Block Maxima 方法是 EVT 分析的最基本方法。Block Maxima 包括将观察期划分为相同大小的不重叠的时期，并将注意力限制在每个时期的最大观察值上。创建的观察遵循吸引条件的域，近似于极值分布。然后将极值分布的参数统计方法应用于这些观察。
 
极值理论家开发了广义极值分布。GEV 包含一系列连续概率分布，即 Gumbel、Frechet 和 Weibull 分布（也称为 I、II 和 III 型极值分布）。
 
在以下 EVT Block Maxima 分析中，10 股指数数据拟合 GEV。绘制得到的分布。创建时间序列图以定位时间轴上的极端事件，从 2006 年到 2016 年。然后创建四个按 Block Maxima 数据顺序排列的图。最后，根据 gev() 函数创建 Block Maxima 分析参数表。
 
gev(ltMeans, x=0.8, m=0)


plt(alVF)
 

 
 
 
 

 
 
 
第 3b 节 - 分块最大值的 VaR 预测
 
为了从 Block Maxima 数据中创建风险价值 (VaR) 估计，将 10 股指数 GEV 数据转换为时间序列。VaR 估计是根据 GEV 时间序列数据进行的。未来值的预测（未来 100 天和 500 天）是从 VaR 数据推断出来的。在结果图中，蓝色圆圈表示未来 100 天的值，红色圆圈表示 500 天的收益率值。
 

# 预测未来500天
aGE500<- preit(aG_va.c, n.ad = 500, ci = 0.9)

plot(aGE500pd.500)
 

 
 
 
 
 
第 3c 部分 - 分块最大值的期望损失ES (CvaR)
 
10只股票指数GEV数据的条件风险值（"CvaR "或 "期望损失"）被计算。首先，利用数据的时间序列，找到最差的0.95%的缩水的最大值。然后，通过极端分布的 "修正 "方法来计算 "估计亏损"，这两种计算的结果都以表格形式呈现。
 
# 条件缩减是最差的0.95%缩减的平均值
ddGV <- xdrow(aEVts[,2])
# CvaR（预期亏损）估计值
CvaR(ts(alE), p=0.95, meho="miie")
 
 
 

  
 
第 3d 节 - 分块极大值的 Hill 估计
 
希尔估计（用于尾部指数的参数估计）验证 10 只股票的 GEV 数据是极值分布。
 
 
 
第 3e 节 - 正态分布的 Anderson-Darling 检验
 
Anderson-Darling 检验是确定大样本数量分布的非正态性的有力决定因素。如果 P 值小于 0.05，则分布与正态性不同。通过该测试发现了一个微小的概率值 3.7^-24。
 
 
 
 
第 3f 节 - 结果表
 
最后，给出了对 10 股指数 GEV 未来价值的估计结果表。3 个 GEV VaR 估计值（和 GEV 期望损失）的点估计值和范围制成表格比较。
 
G_t[3,] <- c("GEV ES",sFale[1],
                     sStble[2], SEble[3],
                     "NA")
GRst
 

 
 
第 3g 节 - 分块极大值的 100 天 GARCH 预测
 
通过将 Block Maxima GEV 分布（10 只股票的指数）拟合到 GARCH(1,1)（广义自回归条件异型）模型，对 Block Maxima EVT 数据进行预测。显示预测公式参数表。创建一个“自相关函数”(ACF) 图，显示随时间变化的重要事件。然后，显示拟合模型结果的一组图。创建对未来 20 天（股票指数表现）的预测。最后，20 天的预测显示在 2 个图中。
 
spec(aanc.ol = list(mel = 'eGARCH',
                                        garer= c(1, 1)),
                  dirion = 'sd')

# 用广义自回归条件异质性拟合模型
alimol = ugct(pec,allV, sovr = 'ybi')

cofale <- dtafe(cof(litol))
oeBal

plt(l.itodl)
 
 
 
 
 
 
 
 

  
 
第 4a 节 - 峰值超过阈值估计 - 10 个股票指数
 
在 EVT 中的峰值超过阈值方法中，选择超过某个高阈值的初始观测值。这些选定观测值的概率分布近似为广义帕累托分布。通过拟合广义帕累托分布来创建最大似然估计 (mle)。MLE 统计数据以表格形式呈现。然后通过 MLE 绘图以图形方式诊断所得估计值。
 
plot(Dseans, u.rg=c(0.3, 0.35))
 
 
 
 
 
 
 
 
第 4b 节 - POT 的 VaR 预测
 
POT 数据的风险价值 (VaR) 估计是通过将 10 个股票指数 MLE 数据转换为时间序列来创建的。VaR 估计是根据 MLE 时间序列数据进行的。未来值的预测（未来 100 天和 500 天）是从 MLE VaR 数据推断出来的。在结果图中，蓝色圆圈表示未来 100 天的值，红色圆圈表示 500 天的收益值。
 
VAR(merts, p = 2, tp = "cost")

# 预测未来125天、250天和500天
mle_r.pd <- prect(e.ar, n.ahad = 100, ci = 0.9)


plot(mea.prd)
 

 
 
 
 
 
第 4c 部分 - POT 的期望损失ES (CvaR) 预测
 
 
然后计算10只股票指数MLE数据的条件风险值（"CvaR "或 "期望损失ES"）。数据的时间序列被用来寻找最差的0.95%的跌幅的最大值。通过极端分布的 "修正 "方法，计算出 "期望损失ES"，两种计算的结果都以表格形式呈现。
 
# 最差的0.95%最大回撤的平均值
mdM <- maxdadw(mlvs[,2])

CvaR(ldaa), p=0.95, meto="mdii",
               pimeod = "comnen", weghts)
 

  
 
 
第 4d 节 - 峰值超过阈值的 Hill 估计
 
Hill 估计（用于尾部指数的参数估计）验证 10 只股票的 MLE 数据是一个极值分布。
 
 
 
第 4e 节 - 正态分布的 Anderson-Darling 检验
 
Anderson-Darling 检验是确定大样本数量分布的非正态性的有力决定因素。如果 P 值小于 0.05，则分布与正态性不同。此测试的结果 P 值为 3.7^-24。
 
 
 
 
第 4f 节 - 结果表
 
最后，给出了 10 个股票指数 MLE 未来价值的估计结果表。3 个 MLE VaR 估计值（和 MLE 期望损失ES）的点估计值和范围被制成表格来比较。
 
 
 
 
 
第 4g 节 - 峰值超过阈值的 100 天 GARCH 预测
 
通过将 MLE（10 只股票指数的最大似然估计）拟合到 GARCH(1,1)（广义自回归条件异型性）模型，对峰值超过阈值 EVT 数据进行预测。显示预测公式参数表。创建了一个“自相关函数”（ACF）图，显示了随时间变化的重要事件。然后，显示拟合模型结果的一组图。然后创建对接下来 20 天（股票指数表现）的预测。最后，20 天的预测（来自峰值超过阈值 EVT extimation）显示在 2 个图中。
 
fit(ec,ta, slvr = 'hybrid')



plot(pot.fite.ol)
 
 
 
 
 
 
 

  
 
第 5a 节 - 估计方法影响表
 
下表汇总了检验 极值分布的 10 个股票的四种方法的结果。第一列包含四种估计方法的名称。提供了 VaR、ES、mu统计量和 Anderson-Darling P 值的统计量。
 
c("VaR",
                     round(mean(cofets),4),
                     "NA", "NA", p.vau)
c("Block Maxm", round(mean(coffies),4),
                     MES, pr.ss[3],.vle)

c("POT", 
                     round(mean(cofies), 4),
                     MES, fitdaes, p.ale)

 
 

  
 
 
第 5b 节 - 结论
 
在对 10 家公司（在 DAX 证券交易所上市）的 10 年股票收益率进行检查后，确认将收益率百分比的变化表征为极值分布的有效性。四种分析方法的拟合值的所有 Anderson-Darling 检验都显示分布具有正态性或所有非极值值的概率不显着。这些方法在收益数据中的风险价值方面是一致的。Block Maxima 方法会产生 VaR 估计的轻微偏差。传统的 VaR 估计和 POT 估计产生相同的风险价值。与股票收益率数据的传统 CvaR 估计相比，这 2 种 EVT 方法预测的预期缺口较低。标准 QQ 图表明峰值超过阈值是最可靠的估计方法，
 
在对10家公司（在德国DAX证券交易所上市）10年的股票收益率进行检查后，证实了将收益率变化定性为极值分布的有效性。对四种分析方法的拟合值进行的所有安德森-达林测试显示，分布具有正态性或所有非极值的概率不大。这些方法在收益数据的风险值方面是一致的。分块最大值方法产生了一个风险值估计的偏差。传统的VaR估计和POT估计产生相同的风险值。相对于传统的股票收益率数据的CvaR估计，两种EVT方法预测的期望损失较低。标准Q-Q图表明，在10只股票的指数中，Peaks-Over-Threshold是最可靠的估计方法。
 

 
 
最受欢迎的见解
 
1.R语言基于ARMA-GARCH-VaR模型拟合和预测实证研究
 
2.R语言时变参数VAR随机模型
 
3.R语言估计时变VAR模型时间序列的实证研究
 
4.R语言基于ARMA-GARCH过程的VAR拟合和预测
 
5.GARCH（1,1），MA以及历史模拟法的VaR比较
 
6.R语言用向量自回归（VAR）进行经济数据脉冲响应
 
7.R语言实现向量自动回归VAR模型
 
8.R语言随机搜索变量选择SSVS估计贝叶斯向量自回归（BVAR）模型
 
9.R语言VAR模型的不同类型的脉冲响应分析

原文链接：http://tecdat.cn/?p=20360 
 
原文出处：拓端数据部落公众号
 
 
本文将说明金融数学中的R 语言优化投资组合，因子模型的实现和使用。
 
具有单一市场因素的宏观经济因素模型
 
我们将从一个包含单个已知因子（即市场指数）的简单示例开始。该模型为
 
 
其中显式因子ft为S＆P 500指数。我们将做一个简单的最小二乘（LS）回归来估计截距α和加载β：
 
 
大多数代码行用于准备数据，而不是执行因子建模。让我们开始准备数据：
 


# 设置开始结束日期和股票名称列表
begin_date <- "2016-01-01"
end_date <- "2017-12-31"


# 从YahooFinance下载数据
data_set <- xts()
for (stock_index in 1:length(stock_namelist))
  data_set <- cbind(data_set, Ad(getSymbols(stock_namelist[stock_index], 
                                            from = begin_date, to = end_date, 
head(data_set)
#>                AAPL  AMD      ADI     ABBV AEZS        A       APD       AA       CF
#> 2016-01-04 98.74225 2.77 49.99239 49.46063 4.40 39.35598 107.89010 23.00764 35.13227
#> 2016-01-05 96.26781 2.75 49.62508 49.25457 4.21 39.22057 105.96097 21.96506 34.03059
#> 2016-01-06 94.38389 2.51 47.51298 49.26315 3.64 39.39467 103.38042 20.40121 31.08988
#> 2016-01-07 90.40047 2.28 46.30082 49.11721 3.29 37.72138  99.91463 19.59558 29.61520
#> 2016-01-08 90.87848 2.14 45.89677 47.77789 3.29 37.32482  99.39687 19.12169 29.33761
#> 2016-01-11 92.35001 2.34 46.98954 46.25827 3.13 36.69613  99.78938 18.95583 28.14919

head(SP500_index)
#>              index
#> 2016-01-04 2012.66
#> 2016-01-05 2016.71
#> 2016-01-06 1990.26
#> 2016-01-07 1943.09
#> 2016-01-08 1922.03
#> 2016-01-11 1923.67
plot(SP500_index)
 
 

# 计算股票和SP500指数的对数收益率作为显式因子
X <- diff(log(data_set), na.pad = FALSE)
N <- ncol(X)  # 股票数量
T <- nrow(X)  # 天数
 
现在我们准备进行因子模型拟合。LS拟合很容易在R中实现，如下所示：
 
 

beta <- cov(X,f)/as.numeric(var(f))
alpha <- colMeans(X) - beta*colMeans(f)
sigma2 <- rep(NA, N)

print(alpha)
#>              index
#> AAPL  0.0003999086
#> AMD   0.0013825599
#> ADI   0.0003609968
#> ABBV  0.0006684632
#> AEZS -0.0022091301
#> A     0.0002810616
#> APD   0.0001786375
#> AA    0.0006429140
#> CF   -0.0006029705
print(beta)
#>          index
#> AAPL 1.0957919
#> AMD  2.1738304
#> ADI  1.2683047
#> ABBV 0.9022748
#> AEZS 1.7115761
#> A    1.3277212
#> APD  1.0239453
#> AA   1.8593524
#> CF   1.5702493
 
或者，我们可以使用矩阵表示法进行拟合，我们定义和扩展因子。然后最小化
 
 t(X) %*% F_ %*% solve(t(F_) %*% F_)  

#>              alpha      beta
#> AAPL  0.0003999086 1.0957919
#> AMD   0.0013825599 2.1738304
#> ADI   0.0003609968 1.2683047
#> ABBV  0.0006684632 0.9022748
#> AEZS -0.0022091301 1.7115761
#> A     0.0002810616 1.3277212
#> APD   0.0001786375 1.0239453
#> AA    0.0006429140 1.8593524
#> CF   -0.0006029705 1.5702493
E <- xts(t(t(X) - Gamma %*% t(F_)), index(X))  # 残差

 
另外，我们可以简单地使用R为我们完成工作：
 

cbind(alpha = factor_model$alpha, beta = factor_model$beta)
#>              alpha     index
#> AAPL  0.0003999086 1.0957919
#> AMD   0.0013825599 2.1738304
#> ADI   0.0003609968 1.2683047
#> ABBV  0.0006684632 0.9022748
#> AEZS -0.0022091301 1.7115761
#> A     0.0002810616 1.3277212
#> APD   0.0001786375 1.0239453
#> AA    0.0006429140 1.8593524
#> CF   -0.0006029705 1.5702493
 
可视化协方差矩阵
 
有趣的是，可视化对数收益率[算术处理误差]以及残差Ψ的估计协方差矩阵。让我们从对数收益率的协方差矩阵开始：
 

         main = "单因子模型对数收益的协方差矩阵")
 
 
我们可以观察到所有股票都是高度相关的，这是市场因素的影响。为了检查股票相关关系，我们绘制相关图：
 
plot(cov2cor(Psi),
         main = "残差协方差矩阵")
 
 

cbind(stock_namelist, sector_namelist)  # 股票的行业
#>       stock_namelist sector_namelist         
#>  [1,] "AAPL"         "Information Technology"
#>  [2,] "AMD"          "Information Technology"
#>  [3,] "ADI"          "Information Technology"
#>  [4,] "ABBV"         "Health Care"           
#>  [5,] "AEZS"         "Health Care"           
#>  [6,] "A"            "Health Care"           
#>  [7,] "APD"          "Materials"             
#>  [8,] "AA"           "Materials"             
#>  [9,] "CF"           "Materials"
 
有趣的是，我们可以观察到对Ψ执行的自动聚类可以正确识别股票的行业。
 
评估投资资金
 
在此示例中，我们将基于因子模型评估几种投资基金的绩效。我们将标准普尔500指数作为明确的市场因素，并假设无风险收益为零 rf = 0。特别是，我们考虑六种交易所买卖基金（ETF）：
 
我们首先加载数据：
 


# 设置开始结束日期和股票名称列表
begin_date <- "2016-10-01"
end_date <- "2017-06-30"

# 从YahooFinance下载数据
data_set <- xts()
for (stock_index in 1:length(stock_namelist))
  data_set <- cbind(data_set, Ad(getSymbols(stock_namelist[stock_index], 

head(data_set)
#>                 SPY   XIVH     SPHB     SPLV     USMV      JKD
#> 2016-10-03 203.6610 29.400 31.38322 38.55683 42.88382 119.8765
#> 2016-10-04 202.6228 30.160 31.29729 38.10687 42.46553 119.4081
#> 2016-10-05 203.5195 30.160 31.89880 38.02249 42.37048 119.9421
#> 2016-10-06 203.6610 30.160 31.83196 38.08813 42.39899 120.0826
#> 2016-10-07 202.9626 30.670 31.58372 37.98500 42.35146 119.8296
#> 2016-10-10 204.0197 31.394 31.87970 38.18187 42.56060 120.5978

head(SP500_index)
#>              index
#> 2016-10-03 2161.20
#> 2016-10-04 2150.49
#> 2016-10-05 2159.73
#> 2016-10-06 2160.77
#> 2016-10-07 2153.74
#> 2016-10-10 2163.66

# 计算股票和SP500指数的对数收益率作为显式因子
X <- diff(log(data_set), na.pad = FALSE)
N <- ncol(X)  # 股票数量
T <- nrow(X)  # 天数

 
现在我们可以计算所有ETF的alpha和beta：
 

#>              alpha      beta
#> SPY   7.142225e-05 1.0071424
#> XIVH  1.810392e-03 2.4971086
#> SPHB -2.422107e-04 1.5613533
#> SPLV  1.070918e-04 0.6777149
#> USMV  1.166177e-04 0.6511667
#> JKD   2.569578e-04 0.8883843
 
现在可以进行一些观察：
 
SPY是S＆P 500的ETF，如预期的那样，其alpha值几乎为零，beta值几乎为1： α= 7.142211×10-5和 β= 1.0071423。
XIVH是具有高alpha值的ETF，计算出的alpha值是ETF中最高的（高1-2个数量级）： α= 1.810392×10-3。
SPHB是一种ETF，据推测具有很高的beta，而计算出的beta却是最高的，但不是最高的：β= 1.5613531。有趣的是，计算出的alpha为负，因此，该ETF应谨慎。
SPLV是降低波动性的ETF，实际上，计算得出的beta偏低：β= 0.6777072。
USMV还是降低波动性的ETF，实际上，计算出的beta是最低的：β= 0.6511671。
JKD显示出很好的折衷。
我们可以使用一些可视化：
 

  barplot(rev(alpha), horiz = TRUE, main = "alph
 
 
我们还可以使用例如Sharpe比率，以更系统的比较不同的ETF。回顾一种资产和一个因素的因子模型
 
我们获得
 
 
夏普比率如下：
 
假设。因此，基于Sharpe比率对不同资产进行排名的一种方法是根据α/β比率对它们进行排名：
 

print(ranking)
#>         alpha/beta         SR         alpha      beta
#> XIVH  7.249952e-04 0.13919483  1.810392e-03 2.4971086
#> JKD   2.892417e-04 0.17682677  2.569578e-04 0.8883843
#> USMV  1.790904e-04 0.12280053  1.166177e-04 0.6511667
#> SPLV  1.580189e-04 0.10887903  1.070918e-04 0.6777149
#> SPY   7.091574e-05 0.14170591  7.142225e-05 1.0071424
#> SPHB -1.551287e-04 0.07401566 -2.422107e-04 1.5613533
 
可以看到：
 
就α/β而言，XIVH最佳（α最大），而SPHB最差（α负）。
就夏普比率（更确切地说，是信息比率，因为我们忽略了无风险利率）而言，JDK是最好的，其次是SPY。这证实了大多数投资基金的表现不超过市场的观点。
显然，无论以哪种衡量标准，SPHB都是最差的：负α，负β比率和Sharpe比率。
JDK之所以能够取得最佳性能，是因为它的alpha值很好（尽管不是最好的），而同时具有0.88的中等beta值。
XIVH和SPHB有大量不同的beta，因此在市场上具有极端敞口。
USMV在市场上的曝光率最小，有可接受的alpha值，并且其Sharpe比率接近第二和第三高的位置。
Fama-French三因子模型
 
该示例将说明使用标准普尔500指数中的九种股票的Fama-French三因子模型。让我们从加载数据开始：
 


# 设置开始结束日期和股票名称列表
begin_date <- "2013-01-01"
end_date <- "2017-08-31"

# 从YahooFinance下载数据
data_set <- xts()
for (stock_index in 1:length(stock_namelist))
  data_set <- cbind(data_set, Ad(getSymbols(stock_namelist[stock_index], 

# 下载Fama-French因子


head(fama_lib)
#>            Mkt.RF   SMB   HML
#> 1926-07-01   0.10 -0.24 -0.28
#> 1926-07-02   0.45 -0.32 -0.08
#> 1926-07-06   0.17  0.27 -0.35
#> 1926-07-07   0.09 -0.59  0.03
#> 1926-07-08   0.21 -0.36  0.15
#> 1926-07-09  -0.71  0.44  0.56
tail(fama_lib)
#>            Mkt.RF   SMB   HML
#> 2017-11-22  -0.05  0.10 -0.04
#> 2017-11-24   0.21  0.02 -0.44
#> 2017-11-27  -0.06 -0.36  0.03
#> 2017-11-28   1.06  0.38  0.84
#> 2017-11-29   0.02  0.04  1.45
#> 2017-11-30   0.82 -0.56 -0.50

# 计算股票的对数收益率和Fama-French因子
X <- diff(log(data_set), na.pad = FALSE)
N <- ncol(X)  #股票数量

 
现在我们在矩阵F中具有三个因子，并希望拟合模型，其中现在的载荷是一个beta矩阵：。我们可以做最小二乘拟合，最小化。更方便地，我们定义和扩展因子 。然后可以将LS公式写为最小化
 
print(Gamma)
#>              alpha        b1          b2          b3
#> AAPL  1.437845e-04 0.9657612 -0.23339130 -0.49806858
#> AMD   6.181760e-04 1.4062105  0.80738336 -0.07240117
#> ADI  -2.285017e-05 1.2124008  0.09025928 -0.20739271
#> ABBV  1.621380e-04 1.0582340  0.02833584 -0.72152627
#> AEZS -4.513235e-03 0.6989534  1.31318108 -0.25160182
#> A     1.146100e-05 1.2181429  0.10370898 -0.20487290
#> APD   6.281504e-05 1.0222936 -0.04394061  0.11060938
#> AA   -4.587722e-05 1.3391852  0.62590136  0.99858692
#> CF   -5.777426e-04 1.0387867  0.48430007  0.82014523
 
另外，我们可以使用R完成：
 
#>              alpha    Mkt.RF         SMB         HML
#> AAPL  1.437845e-04 0.9657612 -0.23339130 -0.49806858
#> AMD   6.181760e-04 1.4062105  0.80738336 -0.07240117
#> ADI  -2.285017e-05 1.2124008  0.09025928 -0.20739271
#> ABBV  1.621380e-04 1.0582340  0.02833584 -0.72152627
#> AEZS -4.513235e-03 0.6989534  1.31318108 -0.25160182
#> A     1.146100e-05 1.2181429  0.10370898 -0.20487290
#> APD   6.281504e-05 1.0222936 -0.04394061  0.11060938
#> AA   -4.587722e-05 1.3391852  0.62590136  0.99858692
#> CF   -5.777426e-04 1.0387867  0.48430007  0.82014523
 
统计因子模型
 
现在让我们考虑统计因子模型或隐式因子模型，其中因子和载荷均不可用。调用具有 K因子的模型 XT =α1T+ BFT + ET的主成分方法：
 
PCA： 
  
样本均值：
矩阵：
样本协方差矩阵：
特征分解：
估计： 
  
 
更新特征分解：
重复步骤2-3，直到收敛为止。
#>              alpha                                        
#> AAPL  0.0007074564 0.0002732114 -0.004631647 -0.0044814226
#> AMD   0.0013722468 0.0045782146 -0.035202146  0.0114549515
#> ADI   0.0006533116 0.0004151904 -0.007379066 -0.0053058139
#> ABBV  0.0007787929 0.0017513359 -0.003967816 -0.0056000810
#> AEZS -0.0041576357 0.0769496344  0.002935950  0.0006249473
#> A     0.0006902482 0.0012690079 -0.005680162 -0.0061507654
#> APD   0.0006236565 0.0005442926 -0.004229364 -0.0057976394
#> AA    0.0006277163 0.0027405024 -0.009796620 -0.0149177957
#> CF   -0.0000573028 0.0023108605 -0.007409061 -0.0153425661
 
同样，我们可以使用R完成工作：
 
#>              alpha      factor1      factor2       factor3
#> AAPL  0.0007074564 0.0002732114 -0.004631647 -0.0044814226
#> AMD   0.0013722468 0.0045782146 -0.035202146  0.0114549515
#> ADI   0.0006533116 0.0004151904 -0.007379066 -0.0053058139
#> ABBV  0.0007787929 0.0017513359 -0.003967816 -0.0056000810
#> AEZS -0.0041576357 0.0769496344  0.002935950  0.0006249473
#> A     0.0006902482 0.0012690079 -0.005680162 -0.0061507654
#> APD   0.0006236565 0.0005442926 -0.004229364 -0.0057976394
#> AA    0.0006277163 0.0027405024 -0.009796620 -0.0149177957
#> CF   -0.0000573028 0.0023108605 -0.007409061 -0.0153425661
 
通过不同因子模型进行协方差矩阵估计的最终比较
 
我们最终将比较以下不同的因子模型：
 
样本协方差矩阵
宏观经济一因素模型
基本的三因素Fama-French模型
统计因素模型
我们在训练阶段估计模型，然后将估计的协方差矩阵与测试阶段的样本协方差矩阵进行比较。估计误差将根据PRIAL（平均损失提高百分比）进行评估：
 
 
加载训练和测试集：
 
# 设置开始结束日期和股票名称列表
begin_date <- "2013-01-01"
end_date <- "2015-12-31"

# 准备股票数据
data_set <- xts()
for (stock_index in 1:length(stock_namelist))
  data_set <- cbind(data_set, Ad(getSymbols(stock_namelist[stock_index], 


#   Fama-French 因子
mydata <- mydata[-nrow(mydata), 


# 准备指数
f_SP500 <- diff(log(SP500_index), na.pad = FALSE)

# 将数据拆分为训练数据和测试数据
T_trn <- round(0.45*T)
X_trn <- X[1:T_trn, ]
X_tst <- X[(T_trn+1):T, ]
 
 现在让我们用训练数据估算不同的因子模型：
 

# 样本协方差矩阵
Sigma_SCM <- cov(X_trn)

# 单因素模型
Gamma <- t(solve(t(F_) %*% F_, t(F_) %*% X_trn))

E <- xts(t(t(X_trn) - Gamma %*% t(F_)), index(X_trn))

# Fama-French三因子模型

Sigma_FamaFrench <- B %*% cov(F_FamaFrench_trn) %*% t(B) + diag(diag(Psi))

# 统计单因子模型

while (norm(Sigma - Sigma_prev, "F")/norm(Sigma, "F") > 1e-3) {
  B <- eigSigma$vectors[, 1:K, drop = FALSE] %*% diag(sqrt(eigSigma$values[1:K]), K, K)



# 统计三因子模型
K <- 3

while (norm(Sigma - Sigma_prev, "F")/norm(Sigma, "F") > 1e-3) {
  B <- eigSigma$vectors[, 1:K] %*% diag(sqrt(eigSigma$values[1:K]), K, K)
  Psi <- diag(diag(Sigma - B %*% t(B)))

Sigma_PCA3 <- Sigma

# 统计五因子模型
K <- 5

eigSigma <- eigen(Sigma)
while (norm(Sigma - Sigma_prev, "F")/norm(Sigma, "F") > 1e-3) {
  B <- eigSigma$vectors[, 1:K] %*% diag(sqrt(eigSigma$values[1:K]), K, K)
  Psi <- diag(diag(Sigma - B %*% t(B)))
 
 
最后，让我们比较测试数据中的不同估计：
 

Sigma_true <- cov(X_tst)

barplot(error, main = "协方差矩阵估计误差", 
 
 
 

 

PRIAL <- 100*(ref - error^2)/ref

barplot(PRIAL, main = "协方差矩阵估计的先验方法", 
 
 
 
最终可以看到使用因子模型进行协方差矩阵估计会有所帮助。
 

 
最受欢迎的见解
 
1.用机器学习识别不断变化的股市状况—隐马尔科夫模型(HMM)的应用
 
2.R语言GARCH-DCC模型和DCC（MVT）建模估计
 
3.R语言实现 Copula 算法建模依赖性案例分析报告
 
4.R语言COPULAS和金融时间序列数据VaR分析
 
5.R语言多元COPULA GARCH 模型时间序列预测
 
6.用R语言实现神经网络预测股票实例
 
7.r语言预测波动率的实现：ARCH模型与HAR-RV模型
 
8.R语言如何做马尔科夫转换模型markov switching model
 
9.matlab使用Copula仿真优化市场风险