标签：Pp 06 方差 模型 SV 笔记 课程 GARCH xp
本文是课程《数据科学与金融计算》第6章的学习笔记，主要介绍GARCH类、SV类模型和高频波动模型，用于知识点总结和代码练习，Q&A为问题及解决方案。
往期回顾：
| 博文 | 内容 |
|---|---|
| 【R】【课程笔记】01 R软件基础知识 | 数据类型、数据结构、运算、绘图等 |
| 【R】【课程笔记】02+03 基于R软件的计算 | 聚类分析、因子分析、神经网络、支持向量机等 |
| 【R】【课程笔记】04+05 数据预处理+收益率计算 | 金融数据处理、收益率、R与C++等 |
| 【R】【课程笔记】06 金融波动模型 | GARCH、SV、高频波动模型等 |
| 【R】【课程笔记】07 分位数回归与VaR（ES）计算 | VaR、ES、极值模型等 |
| 【R】【课程笔记】08 金融投资组合决策分析 | 均值-方差模型、均值-VaR模型、均值-CVaR模型等 |
目录
第六章 金融数据整理与预处理
- 课程链接：理论篇《金融风险管理》02
6.1 GARCH类模型
1、ARCH模型
（1）模型
（2）ARCH效应检验
LM检验（原假设：不存在ARCH效应αi=0）：
（3）估计方法：极大似然估计
2、GARCH模型
（1）模型
（2）模型定阶
q：残差序列（平方），PACF
p：AIC 、BIC准则
（3）模型估计：MLE
（4）模型检验
①标准化残差
②Ljung-Box统计量：均值方程与方差方程设定的正确性。
Q-Q图：分布假设的正确性。
（5）模型预测
两边取期望：
用递推方式对条件方差进行求解
案例：恒生指数 GARCH模型
函数：chartSeries()： 画出价格与交易的时序图
library(fGarch)
source("ArchTest.R")
# (1) 读数据
library(quantmod) # 加载包
getSymbols('^HSI', from='2000-01-01', to='2020-05-08') # 从Yahoo网站下载恒生指数日价格数据
dim(HSI) # 数据规模
names(HSI) # 数据变量名称
chartSeries(HSI, theme='white') # 画出价格与交易的时序图
# HSI <- read.table('HSI.txt') # 或者从硬盘中读取恒生指数日价格数据
# HSI <- as.xts(HSI) # 将数据格式转化为xts格式
# (2) 计算收益率数据
ptd.HSI <- HSI$HSI.Adjusted # 提取日收盘价信息
rtd.HSI <- diff(log(ptd.HSI))*100 # 计算日对数收益
rtd.HSI <- rtd.HSI[-1,] # 删除一期缺失值
plot(rtd.HSI) # 画出日收益序列的时序图
ptm.HSI <- to.monthly(HSI)$HSI.Adjusted # 提取月收盘价信息
rtm.HSI <- diff(log(ptm.HSI))*100 # 计算月对数收益
rtm.HSI <- rtm.HSI[-1,] # 删除一期缺失值
plot(rtm.HSI) # 画出月收益序列的时序图
detach(package:quantmod)
# (3) ARCH效应检验
# rtm.HSI <- as.numeric(rtm.HSI)
ind.outsample <- sub(' ','',substr(index(rtm.HSI), 4, 8)) %in% '2013' # 设置样本外下标：2013年为样本外
ind.insample <- !ind.outsample # 设置样本内下标：其余为样本内
rtm.insample <- rtm.HSI[ind.insample]
rtm.outsample <- rtm.HSI[ind.outsample]
Box.test(rtm.insample, lag=12, type='Ljung-Box')
# 月收益序列不存在自相关
Box.test(rtm.insample^2, lag=12, type='Ljung-Box') # 平方月收益序列存在自相关
ArchTest(x=rtm.insample, lags=12) # 存在显著的ARCH效应
# (4) 模型定阶
epst <- rtm.insample - mean(rtm.insample) # 均值调整对数收益
par(mfrow=c(1,2))
acf(as.numeric(epst)^2, lag.max=20, main='平方序列')
pacf(as.numeric(epst)^2, lag.max=20, main='平方序列')
# (5) 建立GARCH类模型
GARCH.model_1 <- garchFit(~garch(1,1), data=rtm.insample, trace=FALSE) # GARCH(1,1)-N模型
GARCH.model_2 <- garchFit(~garch(2,1), data=rtm.insample, trace=FALSE) # GARCH(1,2)-N模型
GARCH.model_3 <- garchFit(~garch(1,1), data=rtm.insample, cond.dist='std', trace=FALSE) # GARCH(1,1)-t模型
GARCH.model_4 <- garchFit(~garch(1,1), data=rtm.insample, cond.dist='sstd', trace=FALSE) # GARCH(1,1)-st模型
GARCH.model_5 <- garchFit(~garch(1,1), data=rtm.insample, cond.dist='ged', trace=FALSE) # GARCH(1,1)-GED模型
GARCH.model_6 <- garchFit(~garch(1,1), data=rtm.insample, cond.dist='sged', trace=FALSE) # GARCH(1,1)-SGED模型
summary(GARCH.model_1)
summary(GARCH.model_3)
plot(GARCH.model_1) # 请键入相应数字获取信息
# (6) 提取GARCH类模型信息
vol_1 <- fBasics::volatility(GARCH.model_1) # 提取GARCH(1,1)-N模型得到的波动率估计
sres_1 <- residuals(GARCH.model_1, standardize=TRUE) # 提取GARCH(1,1)-N模型得到的标准化残差
vol_1.ts <- ts(vol_1, frequency=12, start=c(1990, 1))
sres_1.ts <- ts(sres_1, frequency=12, start=c(1990, 1))
par(mfcol=c(2,1))
plot(vol_1.ts, xlab='年', ylab='波动率')
plot(sres_1.ts, xlab='年', ylab='标准化残差')
# (7) 模型检验
par(mfrow=c(2,2))
acf(sres_1, lag=24)
pacf(sres_1, lag=24)
acf(sres_1^2, lag=24)
pacf(sres_1^2, lag=24)
par(mfrow=c(1,1))
qqnorm(sres_1)
qqline(sres_1)
# (8) 模型预测
pred.model_1 <- predict(GARCH.model_1, n.ahead = 11, trace = FALSE, mse = 'cond', plot=FALSE)
pred.model_2 <- predict(GARCH.model_2, n.ahead = 11, trace = FALSE, mse = 'cond', plot=FALSE)
pred.model_3 <- predict(GARCH.model_3, n.ahead = 11, trace = FALSE, mse = 'cond', plot=FALSE)
pred.model_4 <- predict(GARCH.model_4, n.ahead = 11, trace = FALSE, mse = 'cond', plot=FALSE)
pred.model_5 <- predict(GARCH.model_5, n.ahead = 11, trace = FALSE, mse = 'cond', plot=FALSE)
pred.model_6 <- predict(GARCH.model_6, n.ahead = 11, trace = FALSE, mse = 'cond', plot=FALSE)
predVol_1 <- pred.model_1$standardDeviation
predVol_2 <- pred.model_2$standardDeviation
predVol_3 <- pred.model_3$standardDeviation
predVol_4 <- pred.model_4$standardDeviation
predVol_5 <- pred.model_5$standardDeviation
predVol_6 <- pred.model_6$standardDeviation
et <- abs(rtm.outsample - mean(rtm.outsample))
rtd.HSI.2013 <- rtd.HSI['2013']
rv <- sqrt(aggregate(rtd.HSI.2013^2, by=substr(index(rtd.HSI.2013), 1, 7), sum))
predVol <- round(rbind(predVol_1,predVol_2,predVol_3,predVol_4,predVol_5,predVol_6,
as.numeric(et), as.numeric(rv)), digits=3)
colnames(predVol) <- 1:11
rownames(predVol) <- c('GARCH(1,1)-N模型','GARCH(1,2)-N模型','GARCH(1,1)-t模型','GARCH(1,1)-st模型',
'GARCH(1,1)-GED模型','GARCH(1,1)-SGED模型','残差绝对值', '已实现波动')
print(predVol)
# (9) 模型选择
cor(t(predVol))
3、GARCH模型拓展
（1）GARCH-M(p,q)模型
（2）IGARCH模型（单整）
（3）EGARCH模型（指数）
（4）TGARCH模型
（5）APARCH模型
library(rugarch)
# (1) GARCH-M模型
GARCHM.spec <- ugarchspec(variance.model=list(model='fGARCH', garchOrder=c(1,1), submodel='GARCH'),
mean.model=list(armaOrder=c(0,0), include.mean=TRUE, archm=TRUE),
distribution.model='norm')
GARCHM.fit <- ugarchfit(GARCHM.spec, data=rtm.insample)
#（2）EGARCH模型
EGARCH.spec <- ugarchspec(variance.model=list(model='eGARCH', garchOrder=c(1,1)),
mean.model=list(armaOrder=c(0,0), include.mean=TRUE),
distribution.model='norm')
EGARCH.fit <- ugarchfit(EGARCH.spec, data=rtm.insample,solver="gosolnp") #solver="hybrid"
EGARCH.fit
#（3）IGARCH模型
IGARCH.spec <- ugarchspec(variance.model=list(model='iGARCH', garchOrder=c(1,1)),
mean.model=list(armaOrder=c(0,0), include.mean=TRUE),
distribution.model='norm')
IGARCH.fit <- ugarchfit(IGARCH.spec, data=rtm.insample)
IGARCH.fit
# (4) TGARCH模型
TGARCH.spec <- ugarchspec(variance.model=list(model='fGARCH', garchOrder=c(1,1), submodel='TGARCH'),
mean.model=list(armaOrder=c(0,0), include.mean=TRUE, archm=FALSE),
distribution.model='norm')
TGARCH.fit <- ugarchfit(TGARCH.spec, data=rtm.insample)
# (5) APARCH模型
APARCH.model_1 <- garchFit(~1+aparch(1,1), data=rtm.insample, trace=FALSE) # GARCH(1,1)-N模型
summary(APARCH.model_1)
APARCH.model_2 <- garchFit(~1+aparch(1,1), data=rtm.insample, delta=2, trace=FALSE) # GARCH(1,1)-N模型
summary(APARCH.model_2)
# 6. 例6-6：案例分析, MGARCH模型
library(rmgarch)
library(xlsx)
# (1) DCC-GARCH模型
StockExchange <- read.xlsx(file='Data.xlsx', sheetIndex='StockExchange')
rstock <- diff(log(StockExchange$stock)) # 计算对数收益率
rexchange <- diff(log(StockExchange$exchange)) # 计算对数收益率
r.StockExchange <- data.frame(rstock=rstock, rexchange=rexchange, row.names=StockExchange$date[-1])
r.StockExchange <- as.xts(r.StockExchange) # 转化为时间序列结构
par(mfrow=c(2,1)) # 时间序列图形
plot(r.StockExchange$rstock, ylab='收益率',main='沪深300',lty=1,cex.lab=0.8,cex.main=0.8)
plot(r.StockExchange$rexchange, ylab='收益率',main='汇率',lty=1,cex.lab=0.8,cex.main=0.8)
garch11.spec <- ugarchspec(mean.model = list(armaOrder = c(1,1),include.mean =TRUE),
variance.model = list(garchOrder = c(1,1),model = "eGARCH"),
distribution.model = "sstd") # 设定GARCH和DCC过程
dcc.garch11.spec <- dccspec(uspec = multispec(replicate(2,garch11.spec)),VAR=FALSE,dccOrder = c(1,1),distribution="mvt")
dcc.fit = dccfit(dcc.garch11.spec, data=r.StockExchange) # DCC估计
dcc.fit # 显示DCC估计结果
plot(dcc.fit) # DCC结果输出，按照要求键入一个数字（1-5）
dcc.fcst = dccforecast(dcc.fit, n.ahead=20) # 外推预测
plot(dcc.fcst) # 显示外推预测，按照要求键入一个数字（1-5）
多元GARCH模型
6.2 SV类模型
一、SV模型（随机效应）
随机波动模型（SV）认为波动率由潜在的不可观测的随机过程所决定，即在波动率方程中引入一个新的随机变量，该变量可能服从马尔科夫过程、随机游走或其他。
- “看不见的积分积掉，平均化”
- “似然函数是观测数据的概率分布”
1、基本模型
补充，给定随机变量的分布，求期望和方差。
n=10000
x=rnorm(n)
mean(log(x^2))
var(log(x^2))
2、参数估计（QML方法，转化为标准的线性状态空间形式）
3、参数估计（证明过程见笔记）
在这个线性状态空间模型中，可以使用Kalman滤波实现模型参数估计，包含三个步骤。
第一步：初始化
第二步：序列迭代
第三步：样本内拟合
拟极大似然（不断变化更新的）：“条件密度可以求和加起来”
4、波动预测（样本外预测与平滑预测）
（1）向前l步预测
（2）向后的平滑预测
- 注：卡尔曼滤波（Kalman filtering）[1]是一种利用线性系统状态方程，通过系统输入输出观测数据，对系统状态进行最优估计的算法。由于观测数据中包括系统中的噪声和干扰的影响，所以最优估计也可看作是滤波过程。
案例：SV模型
建立上证综指收益序列的SV模型，随机误差项服从混合正态分布
# 0. initializing
rm(list=ls())
# 1. read data
# (1) get price
library(quantmod)
getSymbols('^SSEC', from='2008-01-01', to='2013-12-31') # downloads daily price data of APPLE,please wait
dim(SSEC) # dimensions of the data
names(SSEC) # names of the data
detach(package:quantmod)
# (2) comput return of stocks
price.SSEC <- SSEC$SSEC.Adjusted
r.SSEC <- 100*diff(log(price.SSEC)) # 计算日收益率
r <- as.numeric(r.SSEC[-1]) # 移除第一个NA值
# (3) logrithm square return
y <- log(r^2 + 1e-4) # 加上一个极小的数使得log里面的数保持正值
sum(is.na(y))
y=na.omit(y)
导入函数文件：source('sub-06.R')
# 2. load filter function
source('sub-06.R')
#滤波方法，写出似然
# 1. SVfilter through Kalman filter
SVfilter <- function (y, phi0, phi1, sQ, alpha, mu1, sR1, mu2, sR2)
{
y = as.numeric(y)
num = length(y)
Q = sQ^2
R1 = sR1^2
R2 = sR2^2
xp = matrix(0, num, 1)
Pp = matrix(0, num, 1)
xp[1] = 0 #初始值
Pp[1] = phi1^2 + Q #初始方差
pi1 = 1/2
pi2 = 1/2
pit1 = 1/2
pit2 = 1/2
like = 0
# 返回似然值
for (i in 2:num) {
sig1 = Pp[i - 1] + R1 #两组数据，给定t-1，方差
sig2 = Pp[i - 1] + R2
k1 = phi1 * Pp[i - 1]/sig1
k2 = phi1 * Pp[i - 1]/sig2
e1 = y[i] - xp[i - 1] - mu1 - alpha #观测数据减均值
e2 = y[i] - xp[i - 1] - mu2 - alpha
den1 = (1/sqrt(sig1)) * exp(-0.5 * e1^2/sig1) #密度函数
den2 = (1/sqrt(sig2)) * exp(-0.5 * e2^2/sig2)
denom = pi1 * den1 + pi2 * den2 #联合密度
pit1[i] = pi1 * den1/denom
pit2[i] = pi2 * den2/denom
xp[i] = phi0 + phi1 * xp[i-1] + pit1[i]*k1*e1 + pit2[i]*k2*e2 # 更新
Pp[i] = (phi1^2)*Pp[i-1] + Q - pit1[i]*(k1^2)*sig1 - pit2[i]*(k2^2)*sig2
like = like - 0.5 * log(pit1[i]*den1 + pit2[i]*den2) # 写出密度函数
# like = like - 0.5 * log(pi1 * den1 + pi0 * den0)
}
list(xp = xp, Pp = Pp, like = like)
}
# 2. innovations Likelihood
Linn <- function(para, y){
phi0<-para[1]; phi1<-para[2]; sQ<-para[3]; alpha<-para[4];
mu1<-para[5]; sR1<-para[6]; mu2<-para[7]; sR2<-para[8]
sv <- SVfilter(y, phi0, phi1, sQ, alpha, mu1, sR1, mu2, sR2)
return(sv$like)
}
# 3. SV modeling
# (1) Initial Parameters
phi0<-0; phi1<-.95; sQ<-.2; alpha <- mean(y); mu1<- 0; sR1<-2; mu2<- 0; sR2<-2
init.par <- c(phi0, phi1, sQ, alpha, mu1, sR1, mu2, sR2)
# (2) Innovations Likelihood
# see sub-06.R
# (3) Estimation
est <- optim(init.par, Linn, y=y, method="BFGS", hessian=TRUE, control=list(trace=1,REPORT=1)) #参数极大化，极小化负的极大似然
SE <- sqrt(diag(solve(est$hessian))) #取方差对角元开根号
t_test <- est$par/SE
prob <- 2*pt(abs(t_test), df=length(y)-1, lower.tail=FALSE)
u <- round(cbind(estimates=est$par, SE, t_test, prob), digits=4)
rownames(u)<-c("phi0","phi1","sQ","alpha","mu1","sigv1","mu2","sigv2")
print(u)
# 4. compare distribution
# (1) filters at the estimated parameters
phi0<-est$par[1]; phi1<-est$par[2]; sQ<-est$par[3]; alpha<-est$par[4];
mu1<-est$par[5]; sR1<-est$par[6]; mu2<-est$par[7]; sR2<-est$par[8]
sv <- SVfilter(y, phi0, phi1, sQ, alpha, mu1, sR1, mu2, sR2) #随机波动模型
# (2) densities plot (f is chi-sq, fm is fitted mixture)
x <- seq(-15, 6, by=.01)
f <- exp(-.5*(exp(x)-x))/(sqrt(2*pi)) #双指数分布
f1 <- exp(-.5*(x-mu1)^2/sR1^2)/(sR1*sqrt(2*pi))
f2 <- exp(-.5*(x-mu2)^2/sR2^2)/(sR2*sqrt(2*pi))
fm <- (f1+f2)/2 #混合正态
plot(density(y), ylim=c(0, 0.25), main='', xlab='x', ylab='密度', lty=1, lwd=2, col='black')
lines(x, f, lty=2, lwd=2, col='blue')
lines(x, fm, lty=3,lwd=1, col='red')
legend('topleft', legend=c('y的核密度', '卡方自然对数', '拟合的混合正态'),
lty=c(1,2,3), lwd=c(2,1,1), col=c('black', 'blue','red'))
# 5. volatility plot
Time <- 1:100
par(mfrow=c(3,1))
plot(Time, y[Time], type='l', main='SSEC平方收益自然对数)', xlab='时间', ylab='y')
plot(Time, sv$xp[Time], type='l', main='预测波动的自然对数',
ylim=c(-3.0, 2.5), xlab='时间', ylab='对数波动')
lines(Time, sv$xp[Time]+2*sqrt(sv$Pp[Time]), lty='dashed')
lines(Time, sv$xp[Time]-2*sqrt(sv$Pp[Time]), lty='dashed')
plot(Time, exp(sv$xp[Time]), type='l', main='预测波动', xlab='时间', ylab='波动率')
par(mfrow=c(1,1))
案例：多元SV模型
# 0. initializing
rm(list=ls())
# 1. load package
library(stochvol)
# 2. numerical simualtion for univariate SV process
# (1) simulate a highly persistent SV process
set.seed(12345)
sim <- svsim(500, mu = -5, phi = 0.95, sigma = 0.25)
# (2) 有放回地抽取5000个，有先验信息
draws <- svsample(sim$y, draws=5000, burnin=500, priormu=c(-8, 1),priorphi=c(10, 1.2), priorsigma=0.2)
# (3) 估计量，抽完之后看后验
print(sim)
summary(draws)
# (4) predict 20 days ahead
fore <- predict(draws, 20)
plot(draws, forecast = fore)
# plot(draws, pages=1, all.terms=TRUE, forecast = fore)
volplot(draws, forecast = 10)
# (5) re-plot with different quantiles
newquants <- c(0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99)
draws <- updatesummary(draws, quantiles = newquants)
plot(draws, forecast = 20, showobs = FALSE, col = seq(along = newquants),
forecastlty = 3, showprior = FALSE)
volplot(draws, forecast = 10)
6.3 高频波动模型
一、概念
1、金融高频数据及其特征：非等间隔、离散价格、日历效应、多重交易
2、主要变量：
收盘价、高频收益率
3、描述统计：样本均值、标准差（具有相依性）
4、“已实现”方差：“对数收益是增量，相依性为0”
如果日内收益不相关，则：
积分方差
二、模型（详见笔记）
1、HAR-RV
2、HAR-RV-J（跳）
3、HEAVY（对方差、已实现测度建模）
4、ACD模型
用最大似然估计法
假定指数分布（具有无记忆性）
5、ACD模型的扩展
（1）分布角度的扩展
WACD模型、GGACD模型
（2）函数关系角度的扩展
LACD模型、TACD模型
案例：ACD模型
# 0. initializing
rm(list=ls())
# 1. read data
caterp <- read.table('taq-cat-t-jan04t082010.txt', header=T)
head(caterp)
(TT <- nrow(caterp)) # sample size
# 2. comput price change and durations
source('sub-06.R')
hf.series <- compRtnDura(da=caterp, plot.it=TRUE)
print(hf.series$len.total)
print(hf.series$len.norm)
print(hf.series$len.duration)
print(hf.series$ratio)
layout(matrix(c(1,2,3,3), 2, 2, byrow=FALSE), c(3,2), c(2,2), respect=TRUE)
layout.show(3)
par(mar=c(4, 4, 2, 2))
plot(hf.series$prs, xlab='时间', ylab='价格', type='l')
plot(hf.series$prc, xlab='时间', ylab='价格变化', type='l')
hist(hf.series$prc, nclass = 100,xlab='价格变化', ylab='频数', main='价格变化直方图')
# 3. extract trading information
extrTrad(da=caterp, date='2010-01-04', from='09:30', to='09:45', plot.ACF=TRUE)
# 4. comput returns
# (1) comput return series for different frequency
ints <- c(1:5, 10, 15, 20, 25, 30) # include 10 intervals
rtns <- list()
ntrads <- list()
for (i in seq_along(ints)){
cat('the current interval is', ints[i], '\n')
tmp <- hfanal(da=caterp, int=ints[i], basic=1)
rtns[[i]] <- tmp$returns
ntrads[[i]] <- tmp$ntrad
}
# (2) comput descriptive statistics for 1-min return series
par(mfrow=c(3,1))
par(mar=c(5.5, 5.5, 1.5, 2.5))
plot(rtns[[1]], xlab='时间', ylab='对数收益', type='l')
hist(rtns[[1]], xlab='对数收益', ylab='频数', main='')
acf(rtns[[1]], xlab='滞后期', lag.max=300, main='')
par(mfrow=c(1,1))
# # 5. comput realized volatility
rv <- NULL
for (i in seq_along(ints)){
cat('the current interval is', ints[i], '\n')
tmp <- hfanal(da=caterp, int=ints[i], basic=1)
rv <- cbind(rv, tmp$realized)
}
colnames(rv) <- paste('int=', ints, sep='')
matplot(1:nrow(rv), rv, xlab='time', ylab='realized volatility', type='l')
legend('topleft', legend=paste('int=', ints, sep=''), lty=1:10,col=1:10)
# 6. implement ACD modeling
# (1) prepare data
sec <- 3600*caterp$hour+60*caterp$minute+caterp$second # 转换为秒
ist <- 3600*9+30*60; end <- 3600*16
lunch <- 3600*12
idx <- c(1:length(sec))[sec < ist] # before market opens
jdx <- c(1:length(sec))[sec > end] # after market closes
sec.norm <- sec[-c(idx,jdx)] # normal trading hours only.
dt <- diff(sec.norm)
kdx <- c(1:length(dt))[dt > 0] # Positive durations only
ti <- sec.norm[2:length(sec.norm)]
dt <- dt[kdx]
ti <- ti[kdx]
st <- 3600*6.5
f1 <- (ti-lunch)/st # trading time
# (2) fit a simple time series model
m1 <- lm(log(dt)~f1+I(f1^2))
summary(m1)
names(m1)
fit <- m1$fitted.values
adjdt <- dt/exp(fit)
plot(adjdt[1:1200], type='l', xlab='时间', ylab='调整持续期')
# (3) ACD modeling
m2 <- acd(adjdt,order=c(1,1), cond.dist="exp")
m3 <- acd(adjdt,order=c(1,1), cond.dist="weibull")
m4 <- acd(adjdt,order=c(1,2), cond.dist="weibull")
m5 <- acd(adjdt,order=c(1,2), cond.dist="gamma")
names(m5)
ep2 <- m2$epsilon
ep3 <- m3$epsilon
ep4 <- m4$epsilon
ep5 <- m5$epsilon
par(mfrow=c(5,1))
par(mar=c(4,4,3,2))
acf(adjdt, xlab='滞后期', main='调整序列', ylim=c(-0.05,0.25))
acf(ep2, xlab='滞后期', main=expression(epsilon[2]), ylim=c(-0.05,0.25))
acf(ep3, xlab='滞后期', main=expression(epsilon[3]), ylim=c(-0.05,0.25))
acf(ep4, xlab='滞后期', main=expression(epsilon[4]), ylim=c(-0.05,0.25))
acf(ep5, xlab='滞后期', main=expression(epsilon[5]), ylim=c(-0.05,0.25))
par(mfrow=c(1,1))
Box.test(ep5, lag=5, type='Ljung')
Box.test(ep5, lag=10, type='Ljung')
Box.test(ep5, lag=15, type='Ljung')
Box.test(ep5, lag=20, type='Ljung')
Box.test(ep5^2, lag=5, type='Ljung')
Box.test(ep5^2, lag=10, type='Ljung')
Box.test(ep5^2, lag=15, type='Ljung')
Box.test(ep5^2, lag=20, type='Ljung')
案例：高频“已实现”方差
# 0. initializing
#setwd('F:/programe/book/R with application to financial quantitive analysis/CH-06')
#rm(list=ls())
# 1. clean HF data
library(highfrequency)
library(xts)
data(sample_tdataraw)
dim(sample_tdataraw)
head(sample_tdataraw)
args(tradesCleanup)
tdata_afterfirstcleaning <- tradesCleanup(tdataraw=sample_tdataraw, exchanges="N")
names(tdata_afterfirstcleaning)
print(tdata_afterfirstcleaning$report)
barplot(tdata_afterfirstcleaning$report)
dim(tdata_afterfirstcleaning$tdata)
# 2. aggregate HF data
# (1) load sample price data
data("sample_tdata")
ts <- sample_tdata$PRICE
head(ts)
# (2) aggregate previous tick to the 30-second sampling frequency
args(aggregatets)
tsagg30sec = aggregatets(ts, on="seconds", k=30) # 按秒进行整合，每30s做平均
head(tsagg30sec, 20)
# (3) aggregate previous tick to the 5-minute sampling frequency
tsagg5min <- aggregatets(ts, on="minutes", k=5)
head(tsagg5min, 20)
# 3. comput descriptive statistics for realized returns
source('sub-06.R')
data(sample_real5minprices)
# 12-19是时间，取出字符串，！=不等于的意思
pdata <- sample_real5minprices[substr(time(sample_real5minprices), 12, 19)!='00:00:00']
length(pdata)
mins <- c(5, 10, 15, 20, 25, 30)
rtn <- list()
stats <- matrix(NA, nrow=length(mins), ncol=7)
for (i in seq_along(mins)){
rtn[[i]] <- HFrtn(pdata=pdata, on="minutes", k=mins[i])
stats[i,] <- stat(as.vector(rtn[[i]]))
}
#对数差分，加起来，相当于第一个数减最后一个
head(rtn[[1]])
head(rtn[[2]])
rownames(stats) <- paste('min=', mins, sep='')
colnames(stats) <- c('mean', 'sd', 'skew', 'kurt', 'median', 'min', 'max')
print(stats)
par(mfrow=c(2,2))
par(mar=c(4,4,2,2))
plot(mins, stats[,'mean'], xlab='分钟', ylab='均值', type='o', pch='*')
plot(mins, stats[,'sd'], xlab='分钟', ylab='标准差', type='o', pch='*')
plot(mins, stats[,'skew'], xlab='分钟', ylab='偏度', type='o', pch='*')
plot(mins, stats[,'kurt'], xlab='分钟', ylab='峰度', type='o', pch='*')
par(mfrow=c(1,1))
# 4. comput (weighted) realized volatility
# 权重可以根据实际情况修改
days <- levels(factor(substr(index(pdata), 1, 10)))
rv <- matrix(NA, nrow=length(days), ncol=length(mins))
wrv <- matrix(NA, nrow=length(days), ncol=length(mins))
for (i in seq_along(mins)){
ind.day <- substr(index(rtn[[i]]), 1, 10)
rv[,i] <- aggregate(rtn[[i]], ind.day, sumsq)
ind.time <- substr(index(rtn[[i]]), 12, 19)
lambda <- aggregate(rtn[[i]], ind.time, sumsq)/sum(rtn[[i]]^2)
ind.time <- index(lambda)[lambda!=0]
ind <- substr(index(rtn[[i]]), 12, 19) %in% ind.time
rtn.f <- rtn[[i]][ind]
ind.day <- substr(index(rtn.f), 1, 10)
N <- length(ind.time)
lambda <- as.numeric(lambda)[lambda!=0]
weight <- 1/(N*lambda)
weight <- weight/sum(weight)*N
wrv[,i] <- aggregate(rtn.f, ind.day, wsumsq, w=weight)
}
rownames(rv) <- days
colnames(rv) <- paste('min=', mins, sep='')
matplot(1:nrow(rv), sqrt(rv), xlab='时间', ylab='已实现波动', type='l')
legend('topright', legend=paste('min=', mins, sep=''), lty=1:6,col=1:6)
rownames(wrv) <- days
colnames(wrv) <- paste('min=', mins, sep='')
matplot(1:nrow(wrv), sqrt(wrv), xlab='时间', ylab='加权已实现波动', type='l')
legend('topright', legend=paste('min=', mins, sep=''), lty=1:6,col=1:6)
# 5. modeling realized volatility
# (1) HARRV model
data(realized_library) # Get sample daily Realized Volatility data
DJI_RV <- read.csv("DJIRV.txt") #realized_library$Dow.Jones.Industrials.Realized.Variance; # Select DJI
DJI_RV <- xts(DJI_RV[,2], order.by = as.Date(DJI_RV[,1]))
DJI_RV <- DJI_RV[!is.na(DJI_RV)] # Remove NA's
DJI_RV <- DJI_RV['2008']
args(harModel)
HARRV <- harModel(data=DJI_RV, periods=c(1,5,22), RVest=c("rCov"), type="HARRV",h=1)
class(HARRV)
summary(HARRV)
plot(HARRV)
matplot(predict(HARRV, interval='confidence'), type='l', xlab='index', ylab='realized volatility')
# (2) HEAVY model
#returns <- realized_library$Dow.Jones.Industrials.Returns # realized return
#rk <- realized_library$Dow.Jones.Industrials.Realized.Kernel # realized measure
returns <- read.csv("DJIReturn.txt")
rk <- read.csv("DJIRK.txt")
returns <- xts(returns[,2], order.by = as.Date(returns[,1]))
rk <- xts(rk[,2], order.by = as.Date(rk[,1]))
returns <- returns[!is.na(rk)]
rk <- rk[!is.na(rk)] # remove NA's
data <- cbind(returns^2, rk ) # make data matrix with squared returns and realized measures
backcast <- matrix(c(var(returns),mean(rk)), ncol=1)
args(heavyModel)
startvalues <- c(0.004,0.02,0.44,0.41,0.74,0.56) # initial values
HEAVY <- heavyModel(data=as.matrix(data,ncol=2), compconst=FALSE, startingvalues=startvalues, backcast=backcast)
names(HEAVY)
print(HEAVY$estparams)
plot(sqrt(HEAVY$condvar)[,1], xlab='时间', ylab='波动率', ylim=c(-0.01,0.1),main='', lty=1)
lines(sqrt(HEAVY$condvar)[,2],col=2)
lines(sqrt(rk)[,1], lty=2,col=3)
legend('topleft', legend=c('可观测已实现核', 'HEAVY条件波动'), lty=c(1,2))
标签：Pp,06,方差,模型,SV,笔记,课程,GARCH,xp 来源： https://blog.csdn.net/weixin_46285914/article/details/106046732
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。