@fanxy
2020-05-17T19:13:24.000000Z
字数 11251
阅读 15632
樊潇彦
复旦大学经济学院
金融数据
setwd("D:\\..\\Ch13")
install.packages("sm","quantreg","evir")
## 调用
library(sm)
library(quantreg)
library(evir)
library(stats)
library(fPortfolio)
library(RQuantLib)
library(qcc)
library(rugarch)
library(fGarch)
library(zoo)
library(tseries)
library(forecast)
library(timeSeries)
library(tidyverse)
library(readxl)
library(ggplot2)
根据巴塞尔协议,金融风险被分为三类:市场风险、信用风险和操作风险。一个与基本金融理论一致的市场风险测度 应满足下面四个条件( 为随机变量,):
1. 次可加性:
2. 单调性:
3. 正齐性:
4. 变换不变性:
一个金融头寸在时刻 的价值记为 ,在接下来的 期的损失
mn=0.01; sd=2; p=0.05; q=1-p; nu=5 # 参数设定
set.seed(1)
rn=rnorm(1000) # 模拟数据
rt=rt(1000,nu) # 可以用 plot(density(rn)) 命令画密度线
y = c(rn,rt)
g = factor(rep(c("r_n","r_t"), each=1000)) # 用因子变量标记两组
library(sm)
sm.density.compare(y, g) # 合并做图,进行比较
legend("topright", levels(g), fill=2+(0:nlevels(g)))
title(main="正态分布和学生分布密度曲线")
z=qnorm(q); t=qt(q,nu)
VaR_n=mn+z*sd # (7-2)正态分布VaR
VaR_t=mn+t*sd # (7-3)学生分布VaR
data.frame(VaR_n, VaR_t)
定义 ES 为损失 时的期望值,也称为条件或尾部风险值(CVaR,TVaR)。
fz=dnorm(z); ft=dt(t,nu)
ES_n=mn+fz/p*sd # (7-6)正态分布ES
ES_t=mn+ft/p*sd*(nu+t^2)/(nu-1) # (7-7)学生分布ES
data.frame(ES_n,ES_t)
定义第 期的损失率:
# P266:IBM股票
da=read.table("d-ibm-0110.txt",header=T)
head(da)
ibm=log(da[,2]+1)
source("RMfit.R")
mm=RMfit(ibm)
# 用rugarch包复制结果
library(rugarch)
igarch=ugarchspec(variance.model = list(model="iGARCH"),
mean.model = list(armaOrder=c(0,0),include.mean=F), # 均值方程截距为0
fixed.pars=list(omega=0)) # 方差方程截距为0
igarch.fit=ugarchfit(spec=igarch, data=ibm)
igarch.fcst = ugarchforecast(igarch.fit, n.ahead=1)
ibm_sd=as.numeric(igarch.fcst@forecast$sigmaFor) # 预测sigma_T(1)
names(igarch.fit@fit)
c(alpha=1-as.numeric(igarch.fit@fit$coef[2]),
sd=as.numeric(igarch.fit@fit$se.coef),
ibm_sd=ibm_sd) # 回归结果
ibm_mn=0; p=c(0.05,0.01,0.001); q=1-p
z=qnorm(q)
ibm_VaR_n=ibm_mn+z*ibm_sd # VaR
fz=dnorm(z)
ibm_ES_n=ibm_mn+fz/p*ibm_sd # ES
data.frame(p, ibm_VaR_n, ibm_ES_n)
# P267:汇率
da1=read.table("d-useu0111.txt",header=T)
head(da1)
rt=diff(log(da1[,4]))
m2=RMfit(rt)
为损失率 建立ARIMA(p,q)-GARCH(m,n)模型:
# P271:IBM股票
xt=-ibm # 计算多头损失率
library(fGarch)
m1=garchFit(~garch(1,1),data=xt,trace=F)
m1 # P270估计结果
m1pre1=as.numeric(predict(m1,1)) # 提前一步预测
source("RMeasure.R")
m11=RMeasure(mu=m1pre1[1],sigma=m1pre1[2])
m2=garchFit(~garch(1,1),data=xt,trace=F,cond.dist="std") # 学生分布
m2 # P271估计结果
m2pre1=as.numeric(predict(m2,1))
m22=RMeasure(mu=m2pre1[1],sigma=m2pre1[2],cond.dist="std",
df=as.numeric(m2@fit$coef[5])) # 与P271结果不符
# 用rugarch包复制结果
sgarch_n=ugarchspec(variance.model = list(model="sGARCH"),
mean.model = list(armaOrder=c(0,0),include.mean=T),
distribution.model = "norm") # 正态分布
sgarch_n.fit=ugarchfit(spec=sgarch_n, data=xt)
sgarch_n.fcst = ugarchforecast(sgarch_n.fit, n.ahead=1)
xt_nsd=as.numeric(sgarch_n.fcst@forecast$sigmaFor)
xt_nmn=as.numeric(sgarch_n.fit@fit$coef[1])
p=c(0.05,0.01); q=1-p
z=qnorm(q)
xt_VaR_n=xt_nmn+z*xt_nsd
fz=dnorm(z)
xt_ES_n=xt_nmn+fz/p*xt_nsd
data.frame(p, xt_VaR_n, xt_ES_n) # P270结果
sgarch_t=ugarchspec(variance.model = list(model="sGARCH"),
mean.model = list(armaOrder=c(0,0),include.mean=T),
distribution.model = "std") # 学生分布
sgarch_t.fit=ugarchfit(spec=sgarch_t, data=xt)
sgarch_t.fcst = ugarchforecast(sgarch_t.fit, n.ahead=1)
xt_tsd=as.numeric(sgarch_t.fcst@forecast$sigmaFor)
xt_tmn=as.numeric(sgarch_t.fit@fit$coef[1])
nu=as.numeric(sgarch_t.fit@fit$coef[5])
t=qt(q,nu)
xt_VaR_t=xt_tmn+ t*xt_tsd
fz=dnorm(z); ft=dt(t,nu)
xt_ES_t=xt_tmn+ft/p*xt_tsd*(nu+t^2)/(nu-1)
data.frame(p, xt_VaR_t, xt_ES_t) # P271结果
# P274
M1=predict(m1,15) # 正态分布
pmean=sum(M1$meanForecast)
pvar=sum((M1$meanError)^2)
pstd=sqrt(pvar)
M11=RMeasure(pmean,pstd)
# P275
vol=volatility(m2) # 学生分布
a1=c(1.922*10^(-6),0.06448); b1=0.9286; mu=-4.113*10^(-4)
ini=c(ibm[2515],vol[2515])
set.seed(1) # 如果设不同的随机种子,结果与教材差别很大
source("SimGarcht.R") # 模拟数据
mm=SimGarcht(h=15,mu=mu,alpha=a1,b1=b1,df=5.751,ini=ini,nter=30000)
rr=mm$rtn
mean(rr)
q=quantile(rr,c(0.95,0.99)) # VaR
idx=c(1:30000)[rr>q[1]]
mean(rr[idx]) # ES_0.95
idx=c(1:30000)[rr>q[2]]
mean(rr[idx]) # ES_0.99
# P277:样本分位数
quantile(xt,0.95)
sxt=sort(xt)
0.95*2515
es=sum(sxt[2390:2515])/(2515-2389)
es
# P279:分位数回归
dd=read.table("d-ibm-rq.txt",header=T)
head(dd)
dd[,3]=dd[,3]/100
library(quantreg)
mm=rq(nibm~vol+vix,tau=0.95,data=dd) # q=0.95
summary(mm)
names(mm)
fit=mm$fitted.values
tdx=c(2:2515)/252+2001
plot(tdx,dd$nibm,type='l',xlab='year',ylab='neg-log-rtn') # 图7-7
lines(tdx,fit,col='red')
mm=rq(xt~vol+vix,tau=0.99,data=dd) # q=0.99,与P280结果不同
summary(mm)
极值理论(Extreme Value Theory, EVT)讨论 时服从独立同分布的随机变量的极大值 的统计性质。
根据 Fisher–Tippett–Gnedenko 定理,当 时,如果存在均值和标准差数列 和 ,使正规化的极大值随机变量 有渐进分布:
则称 服从广义极值分布(Generalized Extreme Value, GEV),相应的累积分布函数为:
极值分布的三个重要参数:形状参数 (shape),位置参数 (location) 和尺度参数 (scale)。
evir
包中的gev
命令就用这一方法。
# P285表7-1第一列
source("Hill.R") # compile R script
Hill(ibm,110); Hill(xt,110)
# P286表7-2,子区间长度为21个交易日
library(evir)
gev(ibm,block=21)$par.ests
gev(xt,block=21)$par.ests
# P287图7-11
par(mfcol=c(1,2))
m1=gev(xt,block=21)
plot(m1)
# P291:收益率水平
m1=gev(xt,block=21) # n=21
rl.21.12=rlevel.gev(m1,k.block=12) # g=12
rl.21.12 # L=5.43% in (4.65%,6.76%)
# P293:图7-12
par(mfcol=c(2,1))
qplot(xt,threshold=0.01,pch='*',cex=0.8,main="Loss variable of daily IBM log returns")
meplot(ibm)
title(main="Daily IBM log returns")
# P294:表7-3
m1=pot(xt,threshold=0.01) # pot命令
m1$par.ests
m2=pot(xt,threshold=0.012)
m2$par.ests
m3=pot(xt,threshold=0.008)
m3$par.ests
riskmeasures(m1,0.95) # VaR和ES
riskmeasures(m2,0.95)
riskmeasures(m3,0.95)
# P296
m1gpd=gpd(xt,threshold=0.01) # gpd命令
m1gpd
names(m1gpd)
par(mfcol=c(2,2))
plot(m1gpd) # 图7-13
riskmeasures(m1gpd,0.95) # VaR和ES
m1=exindex(xt,10) # 图7-10
m1[which(m1[,2]==21),] # 子区间为21天,theta=0.725
1.966-(1.029/.251)*(1-(-21*.72*log(.99))^(-.251)) # VaR