[关闭]
@fanxy 2020-05-10T11:23:38.000000Z 字数 11589 阅读 4327

第十二讲 波动率模型及应用(III)

樊潇彦 复旦大学经济学院 金融数据


下载:Ch12.rar

  1. setwd("D:\\...\\Ch12")
  2. install.packages("fPortfolio","RQuantLib","qcc")
  3. ## 调用
  4. library(fPortfolio)
  5. library(RQuantLib)
  6. library(qcc)
  7. library(rugarch)
  8. library(fGarch)
  9. library(zoo)
  10. library(tseries)
  11. library(forecast)
  12. library(timeSeries)
  13. library(tidyverse)
  14. library(readxl)
  15. library(ggplot2)
  16. load("Ch12_data.Rdata") # 加载数据

1. 波动率期限结构

1.1 均值回转的半衰期

对于标准GARCH(1,1)模型:


时,方差的均值 ,均值回转的半衰期为:

  1. m1=garchFit(~1+garch(1,1),data=csco,trace=F)
  2. coef(m1) # 表5-1第二行:alpha1 + beta1 = 0.9770
  3. # 注:用rugarch包做,结果不同
  4. garch.sec= ugarchspec(
  5. variance.model = list(model="sGARCH",garchOrder=c(1,1)),
  6. mean.model = list(armaOrder=c(0,0),include.mean=T))
  7. garch_csco=ugarchfit(spec=garch.sec, data=csco)
  8. coef(garch_csco) # alpha1 + beta1 = 0.9927
  9. log(0.5)/log(coef(m1)[3]+coef(m1)[4]) # csco波动率冲击的半衰期
  10. m2=garchFit(~1+garch(1,1),data=cat,trace=F)
  11. log(0.5)/log(coef(m2)[3]+coef(m2)[4]) # cat波动率冲击的半衰期
  12. # 用rugarch包为cat建GARCH模型,sigma的拟合与预测
  13. garch_cat=ugarchfit(spec=garch.sec, data=cat)
  14. tdx=c(1:2515)/252+2001
  15. par(mfcol=c(2,1))
  16. plot(x=tdx,garch_cat@fit$sigma,type="l",
  17. xlab="",ylab="sigma",main="CAT波动率拟合值")
  18. fc_cat=ugarchforecast(garch_cat,n.ahead=30,data=cat) # 预测30期
  19. plot(fc_cat,which=3) # 作sigma的预测图,可以看到均值回转

1.2 波动率的期限结构

距离 时刻 期的某资产的对数收益率为:


在有效市场假设下,各期自相关系数为零,因此条件方差( 期的所有信息):

对于GARCH(1,1)模型,以 时刻为预测原点,第 期对数收益率的条件方差 等于之后各期预测值之和:

  1. par(mfcol=c(1,1))
  2. tdx=c(1:507)/252+2009
  3. H=c(1,seq(5,40,by=5)) # h=1,5,10...40
  4. t_start=which(index(csco)=="2008-12-29")
  5. # 计算时间较长,可跳过以下三行,直接加载数据作图:
  6. #source("FCmat.R") # 调用函数
  7. #FCmat_csco=FCmat(H, var=csco, t_start,t_end=length(csco))
  8. #FCmat_cat =FCmat(H, var=cat, t_start,t_end=length(cat))
  9. load("FCmat.RData")
  10. matplot(x=tdx,FCmat_csco,type="l", # P193,图5-3
  11. xlab="",ylab="",main="CSCO波动率期限结构")
  12. matplot(x=tdx,FCmat_cat,type="l", # P194,图5-4
  13. xlab="",ylab="",main="CAT波动率期限结构")

2. 期权定价和对冲

2.1 期权定价公式和Duan(1995)的方法

假定股票对数价格服从如下随机过程:


相应的离散时间形式为:

其中 。给定 ,可以通过生成随机数 模拟出 个价格序列 ,每个价格序列

在实际应用中,波动率是随时间变化的,Duan(1995)提出估计下述方程中的参数


进而可以拟合或预测随时间变化的 ,或者用:

代替收益率样本方差,进行期权价格的模拟和测算。

2.2 程序实现

  1. # 读取 Duan(1995) S&P100 数据
  2. sp100=read.csv("sp100.csv",header=T)
  3. head(sp100)
  4. rsp=diff(log(sp100$Adj.Close))
  5. rsp=c(0,rsp)
  6. tdx=(1:nrow(sp100))/252+1986
  7. plot(tdx,rsp,type="l")
  8. abline(v=tdx[which(sp100$Date=="1987/10/19")],col="red") # 黑色星期一
  9. P0=100; rf=0; K=90; # 设定期权参数
  10. source("myfun_ch12.R") # 调用程序
  1. sigma_a=0.245 # 年化收益率标准差
  2. sigma=sigma_a/sqrt(360) # 相应的日收益率标准差
  3. # 或者遵循 Duan(1995) 的方程设定,相应的年化收益率标准差也为0.245
  4. alpha0=1.5e-5; alpha1=0.19; beta1=0.72; theta=0.007; lambda=0
  5. for (fixsigma in c(1,0)){
  6. # 模拟数据
  7. P_simu=Psimu(S=1000, T=180, fixsigma)
  8. matplot(P_simu[,1:30],type="l",
  9. main=paste("价格模拟(fixsigma= ",fixsigma,")",sep=""),
  10. xlab="",ylab="")
  11. P_simu_est=Psimu(S=30, T=1000, fixsigma)
  12. # 计算期权价格
  13. max_PT_0=as.numeric(P_simu[T,]-K)*as.numeric(P_simu[T,]>K)
  14. Eurocall=exp(-rf*T)*mean(max_PT_0)
  15. Eurocall_check=EuropeanOption(type="call", underlying=100, strike=90,
  16. dividendYield = 0, riskFreeRate=0, maturity=0.5,
  17. volatility=sigma_a)$value # 用RQuantLib中命令计算
  18. print(data.frame(fixsigma,Eurocall,Eurocall_check)) # 比较欧式期权结果
  19. mnPs=apply(as.matrix(P_simu),2,mean)
  20. max_mnPs_0=as.matrix(mnPs-K)*as.numeric(mnPs>K)
  21. Asiacall=exp(-rf*T)*mean(max_mnPs_0)
  22. # AsianOption(averageType="arithmetic", type="call", underlying=100,
  23. # strike=90, div=0, riskFree=0, maturity=0.5,vol=sigma)
  24. # 亚式看涨期权计算时间较长,结果为:
  25. Asiacall_check=11.563
  26. print(data.frame(fixsigma,Asiacall,Asiacall_check)) # 比较亚式期权结果
  27. }
  1. r=rsp # 用sp100收益率
  2. for (withlambda in c(0,1)){
  3. if (withlambda==0){ # 没有lambda
  4. para0=c(alpha0=1.5e-5,alpha1=0.13,beta1=0.7,theta=0.8)*1.1
  5. }else{ # 有lambda
  6. para0=c(alpha0=1.5e-5,alpha1=0.13,beta1=0.7,theta=0.8,
  7. lambda=0.05)*1.1
  8. }
  9. report(para0)
  10. }
  11. # 或者用模拟数据估计(模拟数据没有lambda)
  12. s=floor(runif(1,min=1,max=30))
  13. r=c(0,diff(log(P_simu_est[,s])))
  14. para0=c(alpha0=1.5e-5,alpha1=0.19,beta1=0.72,theta=0.007)*1.1
  15. report(para0)

3. 随时间变化的相关系数和

3.1 随时间变化的相关系数

如果根据定义,用两种资产的收益率直接计算协方差和相关系数,则得到的 是不随时间变化的,但是在实际运用中,我们希望得到任意两种资产随时间变化的相关系数。根据统计理论,有:


进而有:

因此我们可以先为两种资产的对数收益率 以及它们的和与差 建立波动率模型,再把四个模型波动率的拟合值代入上式,得到时变的相关系数。

  1. # P198:
  2. xp=csco+cat
  3. xm=csco-cat
  4. m1=garchFit(~1+garch(1,1),data=csco,trace=F)
  5. vcsco=m1@sigma.t
  6. m2=garchFit(~1+garch(1,1),data=cat,trace=F)
  7. vcat=m2@sigma.t
  8. m3=garchFit(~1+garch(1,1),data=xp,trace=F) # 收益率之和的模型
  9. summary(m3)
  10. vxp=m3@sigma.t
  11. m4=garchFit(~1+garch(1,1),data=xm,trace=F) # 收益率之差的模型
  12. summary(m4)
  13. vxm=m4@sigma.t
  14. COV=(vxp^2-vxm^2)/4 # (5-5)式协方差
  15. COR=COV/(vcat*vcsco) # (5-6)式相关系数

此外,还可以基于指数加权移动平均方法(Exponentially Weighted Moving Average, EWMA)来计算方差和协方差。给定权重 ,样本 的EWMA预测值为:


Tsay(2015)设 qcc包中的ewma命令方程为:

其中

  1. # P200:
  2. source("EWMAvol.R")
  3. rtn=data.frame(csco,cat)
  4. M1=EWMAvol(rtn,theta=0.94)
  5. cor_ewma=M1[,3]/sqrt(M1[,1]*M1[,2])
  6. # 以均值为起点,调用qcc包中的ewma命令复制结果
  7. dacsco=c(mean(csco^2),(csco^2)[-2515])
  8. dacat=c(mean(cat^2),(cat^2)[-2515])
  9. dacscocat=c(mean(csco*cat),(csco*cat)[-2515])
  10. library(qcc)
  11. vcsco <- ewma(dacsco, lambda=0.06, plot=F)
  12. vcat <- ewma(dacat, lambda=0.06, plot=F)
  13. vcscocat <- ewma(dacscocat, lambda=0.06, plot=F)
  14. cor_ewma2=vcscocat$y/sqrt(vcsco$y*vcat$y)
  15. tdx=(1:2515)/252+2001
  16. par(mfcol=c(2,1))
  17. range(COR,cor_ewma,cor_ewma2) # 查看纵轴取值范围,设定 ylim,作图5-5
  18. plot(tdx,COR,xlab='year',ylab='cor',ylim=c(-0.35,1),type='l')
  19. title(main='(a) GARCH based')
  20. plot(tdx,cor_ewma,xlab='year',ylab='cor',ylim=c(-0.35,1),type='l')
  21. title(main='(b) EWMA with theta = 0.94')
  22. lines(tdx,cor_ewma2,col="red")

3.2 随时间变化的

回归标准CAPM模型:


可以得到资产的市场风险敏感性因子:

类似的,我们可以通过对 的波动率建模,计算可变的

  1. xp=cat+sp5
  2. xm=cat-sp5
  3. m1=garchFit(~1+garch(1,1),data=xp,trace=F) # 分别建GARCH(1,1)模型
  4. summary(m1)
  5. m2=garchFit(~1+garch(1,1),data=xm,trace=F)
  6. summary(m2)
  7. m3=garchFit(~1+garch(1,1),data=sp5,trace=F)
  8. summary(m3)
  9. vxp=m1@sigma.t # 提取标准差
  10. vxm=m2@sigma.t
  11. vsp5=m3@sigma.t
  12. beta=(vxp^2-vxm^2)/(4*vsp5^2) # 计算可变的beta
  13. par(mfcol=c(1,1))
  14. tdx=c(1:2515)/252+2001
  15. coef(lm(cat~sp5))
  16. plot(tdx,beta,xlab='year',ylab='beta',type='l') # 作图
  17. abline(h=c(1.146)) # 标记不变的beta
  18. index(cat)[which(beta==max(beta))] # 显示风险最大的一天

4. 最小方差投资组合

4.1 理论框架

种资产的协方差矩阵 ,权重 和收益率 分别为 维列向量, 为单位列向量。允许做空、目标收益率为 的最小方差投资组合问题可以描述为:


该问题的拉格朗日函数为:

将一阶条件(FOC) 和约束条件联立为线性方程组,进而求解:

如果不考虑目标收益率,问题的解简化为:

4.2 应用示例

  1. rtn=rtn5 # 5支股票
  2. colnames(rtn)=stocklist=c("BA","CAT","IBM","MSFT","PG")
  3. # 计算最小方差组合权重和波动率,绘制表5-2
  4. One=matrix(1,5,1)
  5. period=list(c(1:756),c(757:1512),c(1513:2515))
  6. Wgt=Vol=data.frame()
  7. for (p in 1:3){
  8. V=cov(rtn[period[[p]],])
  9. D=t(One) %*% solve(V) %*% One
  10. w=solve(V) %*% One / as.numeric(D)
  11. Wgt=rbind(Wgt,t(w))
  12. v=sqrt(diag(V))
  13. sd_p=sqrt(t(w)%*%V%*%w)
  14. Vol=rbind(Vol,c(v,sd_p=sd_p))
  15. }
  16. Wgt=round(t(Wgt)*100,2)
  17. rownames(Wgt)=stocklist
  18. Vol=round(t(Vol)*100,2)
  19. rownames(Vol)=c(stocklist,"Portfolio")
  20. Wgt
  21. Vol
  22. # 用fPortfolio包中的命令复制表5-2
  23. library(timeSeries)
  24. rtn=timeSeries(rtn[1:756,])
  25. Spec <- portfolioSpec()
  26. setSolver(Spec) <- "solveRshortExact"
  27. #setTargetReturn(Spec) <- 0.005 # 设定目标收益率
  28. a=efficientPortfolio(rtn, Spec)
  29. a@portfolio@portfolio$weights # 查看权重
  30. sqrt(diag(a@data@statistics$Cov)) # 资产标准差
  31. a@portfolio@portfolio$targetRisk[2] # 组合标准差
  1. rtn=rtn3 # 3只股票
  2. colnames(rtn)=c("ABT","IBM","WMT")
  3. # 下面的命令运行时间很长(约13分钟),建议跳过
  4. #timestart<-Sys.time()
  5. #source("GMVPnew.R")
  6. #M2=GMVP(rtn,start=2011)
  7. #timeend<-Sys.time()
  8. #runningtime<-timeend-timestart
  9. #print(runningtime)
  10. # 直接加载结果
  11. load("M2_Portfolio.RData")
  12. names(M2)
  13. wgt=M2$weights
  14. range(wgt)
  15. prtn=M2$returns
  16. mean(prtn)
  17. sqrt(var(prtn))
  18. Mean=apply(rtn[2012:2515,],2,mean)
  19. Mean
  20. v1=sqrt(apply(rtn[2012:2515,],2,var))
  21. print(v1)
  22. minV=sqrt(M2$minVariance)
  23. Vol=sqrt(M2$variances)
  24. range(minV,Vol)
  25. tdx=c(1:505)/505+2009
  26. par(mfcol=c(2,1))
  27. plot(tdx,wgt[1,],xlab='year',ylab='weights',type='l',ylim=c(-.75,1.5))
  28. lines(tdx,wgt[2,],lty=2,col="blue")
  29. lines(tdx,wgt[3,],lty=3,col="green")
  30. plot(tdx,Vol[,1],xlab='year',ylab='vol',type='l',ylim=c(0,0.04))
  31. lines(tdx,Vol[,2],lty=2,col="blue")
  32. lines(tdx,Vol[,3],lty=3,col="green")
  33. lines(tdx,minV,lty=4,col="red")

5. 预测

Tsay(2015)以美国1997-2010年每周的原油价格为例,介绍如何利用ARMA-GARCH模型提高对金融数据的预测力。具体步骤为:

  1. # P212
  2. da=read.table("w-petroprice.txt",header=T)
  3. head(da)
  4. price=ts(da$US,frequency=52,start=c(1997,1))
  5. dp=ts(diff(price),frequency=52,start=c(1997,2))
  6. cprice=diff(price)
  7. arima(cprice,order=c(3,0,0),
  8. seasonal=list(order=c(2,0,0),period=5)) # (5-9)式,截距项不显著
  9. arima(cprice,order=c(3,0,0),
  10. seasonal=list(order=c(2,0,0),period=5),
  11. include.mean=F) # (5-9)式,去掉截距项
  12. m2=arima(cprice,seasonal=list(order=c(2,0,0),period=5), # 季节模型
  13. include.mean=F)
  14. c_star=cprice[11:716]-coef(m2)[1]-coef(m2)[2]*cprice[1:706] # 根据季节模型进行调整
  1. m3=garchFit(~arma(3,0)+garch(1,1),data=c_star,trace=F,include.mean=F)
  2. summary(m3)
  3. par(mfcol=c(2,1))
  4. plot(m3,which=c(3,13)) # 图5-13
  5. m4=garchFit(~arma(3,0)+garch(1,1),data=c_star,trace=F,include.mean=F,cond.dist="std")
  6. summary(m4)
  7. plot(m4,which=13) # 图5-14(a)
  8. m5=garchFit(~arma(1,0)+garch(1,1),data=c_star,trace=F,include.mean=F,cond.dist="sstd")
  9. summary(m5)
  10. plot(m5,which=13) # 图5-14(b)
  11. plot(m5,which=c(3,7)) # 图5-15
  1. M3=arima(c_star,order=c(3,0,0),include.mean=F)
  2. source("backtest.R")
  3. M3F=backtest(M3,c_star,650,2,inc.mean=F) # AR(3)预测
  4. source("backtestGarch.R")
  5. M4F=backtestGarch(c_star,650,2,inc.mean=F,cdist="sstd") # AR(1)-GARCH(1,1)预测

更多应用

Ch12_paper.rar7548.3kB

添加新批注
在作者公开此批注前,只有你和作者可见。
回复批注