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

第十二讲 波动率模型及应用（III）

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

下载：Ch12.rar

  1. setwd("D:\\...\\Ch12")
  2. install.packages("fPortfolio","RQuantLib","qcc")
  3. ## 调用
  4. library(fPortfolio)
  5. library(RQuantLib)
  6. library(qcc)
  8. library(rugarch)
  9. library(fGarch)
  10. library(zoo)
  11. library(tseries)
  12. library(forecast)
  13. library(timeSeries)
  15. library(tidyverse)
  16. library(readxl)
  17. library(ggplot2)
  19. 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
  10. log(0.5)/log(coef(m1)[3]+coef(m1)[4])                  # csco波动率冲击的半衰期
  11. m2=garchFit(~1+garch(1,1),data=cat,trace=F)
  12. log(0.5)/log(coef(m2)[3]+coef(m2)[4])                  # cat波动率冲击的半衰期
  14. # 用rugarch包为cat建GARCH模型，sigma的拟合与预测
  15. garch_cat=ugarchfit(spec=garch.sec, data=cat)  
  16. tdx=c(1:2515)/252+2001
  17. par(mfcol=c(2,1))
  18. plot(x=tdx,garch_cat@fit$sigma,type="l",
  19.      xlab="",ylab="sigma",main="CAT波动率拟合值")      
  20. fc_cat=ugarchforecast(garch_cat,n.ahead=30,data=cat)   # 预测30期
  21. 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")
  6. # 计算时间较长，可跳过以下三行，直接加载数据作图：
  7. #source("FCmat.R")                                        # 调用函数
  8. #FCmat_csco=FCmat(H, var=csco, t_start,t_end=length(csco))
  9. #FCmat_cat =FCmat(H, var=cat, t_start,t_end=length(cat))
  11. load("FCmat.RData")
  12. matplot(x=tdx,FCmat_csco,type="l",                        # P193，图5-3
  13.         xlab="",ylab="",main="CSCO波动率期限结构") 
  14. matplot(x=tdx,FCmat_cat,type="l",                         # P194，图5-4
  15.         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") # 黑色星期一
  10. P0=100; rf=0; K=90;                                      # 设定期权参数
  11. 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
  6. for (fixsigma in c(1,0)){
  7.   # 模拟数据
  8.   P_simu=Psimu(S=1000, T=180, fixsigma)
  9.   matplot(P_simu[,1:30],type="l",
  10.           main=paste("价格模拟(fixsigma= ",fixsigma,")",sep=""),
  11.           xlab="",ylab="")
  12.   P_simu_est=Psimu(S=30, T=1000, fixsigma)
  14.   # 计算期权价格
  15.   max_PT_0=as.numeric(P_simu[T,]-K)*as.numeric(P_simu[T,]>K)
  16.   Eurocall=exp(-rf*T)*mean(max_PT_0)
  17.   Eurocall_check=EuropeanOption(type="call", underlying=100, strike=90,
  18.                  dividendYield = 0, riskFreeRate=0, maturity=0.5, 
  19.                  volatility=sigma_a)$value               # 用RQuantLib中命令计算
  20.   print(data.frame(fixsigma,Eurocall,Eurocall_check))    # 比较欧式期权结果
  22.   mnPs=apply(as.matrix(P_simu),2,mean)
  23.   max_mnPs_0=as.matrix(mnPs-K)*as.numeric(mnPs>K)
  24.   Asiacall=exp(-rf*T)*mean(max_mnPs_0)
  25.   # AsianOption(averageType="arithmetic", type="call", underlying=100,
  26.   #            strike=90, div=0, riskFree=0, maturity=0.5,vol=sigma)
  27.   # 亚式看涨期权计算时间较长，结果为：
  28.   Asiacall_check=11.563
  29.   print(data.frame(fixsigma,Asiacall,Asiacall_check))    # 比较亚式期权结果
  30. }
  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. }
  12. # 或者用模拟数据估计（模拟数据没有lambda）
  13. s=floor(runif(1,min=1,max=30))
  14. r=c(0,diff(log(P_simu_est[,s])))              
  15. para0=c(alpha0=1.5e-5,alpha1=0.19,beta1=0.72,theta=0.007)*1.1 
  16. report(para0)

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

3.1 随时间变化的相关系数

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


进而有：

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

  1. # P198：
  2. xp=csco+cat
  3. xm=csco-cat  
  5. m1=garchFit(~1+garch(1,1),data=csco,trace=F) 
  6. vcsco=m1@sigma.t 
  7. m2=garchFit(~1+garch(1,1),data=cat,trace=F) 
  8. vcat=m2@sigma.t  
  9. m3=garchFit(~1+garch(1,1),data=xp,trace=F)          # 收益率之和的模型
  10. summary(m3)                       
  11. vxp=m3@sigma.t   
  12. m4=garchFit(~1+garch(1,1),data=xm,trace=F)          # 收益率之差的模型
  13. summary(m4)
  14. vxm=m4@sigma.t   
  15. COV=(vxp^2-vxm^2)/4                                 # (5-5)式协方差
  16. 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])
  7. # 以均值为起点，调用qcc包中的ewma命令复制结果
  8. dacsco=c(mean(csco^2),(csco^2)[-2515])
  9. dacat=c(mean(cat^2),(cat^2)[-2515])
  10. dacscocat=c(mean(csco*cat),(csco*cat)[-2515])
  11. library(qcc)
  12. vcsco <- ewma(dacsco, lambda=0.06, plot=F)
  13. vcat <- ewma(dacat, lambda=0.06, plot=F)
  14. vcscocat <- ewma(dacscocat, lambda=0.06, plot=F)
  15. cor_ewma2=vcscocat$y/sqrt(vcsco$y*vcat$y)
  17. tdx=(1:2515)/252+2001
  18. par(mfcol=c(2,1))
  19. range(COR,cor_ewma,cor_ewma2)        # 查看纵轴取值范围，设定 ylim，作图5-5
  20. plot(tdx,COR,xlab='year',ylab='cor',ylim=c(-0.35,1),type='l')
  21. title(main='(a) GARCH based')
  22. plot(tdx,cor_ewma,xlab='year',ylab='cor',ylim=c(-0.35,1),type='l')
  23. title(main='(b) EWMA with theta = 0.94') 
  24. 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
  14. par(mfcol=c(1,1))
  15. tdx=c(1:2515)/252+2001
  16. coef(lm(cat~sp5))
  17. plot(tdx,beta,xlab='year',ylab='beta',type='l')            # 作图
  18. abline(h=c(1.146))                                         # 标记不变的beta
  19. 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
  23. # 用fPortfolio包中的命令复制表5-2
  24. library(timeSeries)
  25. rtn=timeSeries(rtn[1:756,])
  26. Spec <- portfolioSpec()
  27. setSolver(Spec) <- "solveRshortExact"
  28. #setTargetReturn(Spec) <- 0.005                  # 设定目标收益率
  29. a=efficientPortfolio(rtn, Spec)
  30. a@portfolio@portfolio$weights                    # 查看权重
  31. sqrt(diag(a@data@statistics$Cov))                # 资产标准差
  32. a@portfolio@portfolio$targetRisk[2]              # 组合标准差
  1. rtn=rtn3                                         # 3只股票
  2. colnames(rtn)=c("ABT","IBM","WMT")
  4. # 下面的命令运行时间很长（约13分钟），建议跳过
  5. #timestart<-Sys.time()
  6. #source("GMVPnew.R")
  7. #M2=GMVP(rtn,start=2011)  
  8. #timeend<-Sys.time()
  9. #runningtime<-timeend-timestart
  10. #print(runningtime) 
  12. # 直接加载结果
  13. load("M2_Portfolio.RData")
  14. names(M2)
  15. wgt=M2$weights
  16. range(wgt)
  17. prtn=M2$returns
  18. mean(prtn)
  19. sqrt(var(prtn))
  20. Mean=apply(rtn[2012:2515,],2,mean)
  21. Mean
  22. v1=sqrt(apply(rtn[2012:2515,],2,var))
  23. print(v1)
  24. minV=sqrt(M2$minVariance)
  25. Vol=sqrt(M2$variances)
  26. range(minV,Vol)
  27. tdx=c(1:505)/505+2009
  28. par(mfcol=c(2,1))
  29. plot(tdx,wgt[1,],xlab='year',ylab='weights',type='l',ylim=c(-.75,1.5))
  30. lines(tdx,wgt[2,],lty=2,col="blue")
  31. lines(tdx,wgt[3,],lty=3,col="green") 
  32. plot(tdx,Vol[,1],xlab='year',ylab='vol',type='l',ylim=c(0,0.04))
  33. lines(tdx,Vol[,2],lty=2,col="blue")
  34. lines(tdx,Vol[,3],lty=3,col="green")
  35. 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)式，去掉截距项
  13. m2=arima(cprice,seasonal=list(order=c(2,0,0),period=5),       # 季节模型
  14.       include.mean=F)
  15. 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

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