library(stringr)
Settle = "2020-12-31"; #当前时刻
option = "2103"
datafreerisk<-read.csv("三个月国债利率.csv",header=TRUE)
i=1
while(!as.character(as.Date(Settle)-i)%in%datafreerisk[,1]){
i=i+1;
panduan<-as.character(as.Date(Settle)-i)
}
R<-subset(datafreerisk,datafreerisk[,1]==panduan)[1,2]#3个月国债的无风险利率
R<-mean(as.numeric(as.character(R)))/100
datadetail<-read.csv("沪深300期权合约基本资料.csv",header=TRUE)
Exercise<-subset(datadetail,str_detect(datadetail[,1],option)==1)[1,12]
data<-read.csv("沪深３００期权日行情.csv",header=TRUE)
data2103_31<-merge(subset(data,data[,1]==Settle&str_detect(data[,2], option)==1),
subset(datadetail,str_detect(datadetail[,1],option)==1),
by="期权代码")
T1=as.numeric(as.Date(Exercise)-as.Date(Settle))/250
MarketStrikes = data2103_31[,29]; #读取期权敲定价数据
MarketPrices = as.numeric(as.character(data2103_31[,12])); #读取期权市场最新报价数据
dataforward<-read.csv("沪深３００指数期货行情序列.csv",header=TRUE)
CurrentForwardValue<-subset(dataforward,str_detect(dataforward[,1],Settle)==1)[1,4]; #当前远期价格
CurrentForwardValue<-as.numeric(as.character(CurrentForwardValue))
BSM=function(F0,k,T,r,sigma){
d1=(log(F0/k)+(r+0.5*sigma^2)*T)/sigma/sqrt(T)
d2=d1-sigma*sqrt(T)
call=F0*pnorm(d1)-k*exp(-r*T)*pnorm(d2)
put=k*exp(-r*T)*pnorm(-d2)-F0*pnorm(-d1)
return(c(call,put))
}
#二分法求隐含波动率implied volatility
bisection=function(P,F0,k,r,T,c){ #c=0,call;c=1,put.
value0=0#迭代初始initial value
top=100
floor=0
sigma=(floor+top)/2
i=0
while(abs(P-value0)>0.0001){
value0=BSM(F0,k,T,r,sigma)[c+1]
if ((P-value0)>0){
floor=sigma
sigma=(sigma+top)/2
}else
{
top=sigma
sigma=(sigma+floor)/2}
i=i+1
if(i>1000){break}
}
return (sigma)}
bisection_sigma=function(P,F0,k,r,T,c){
if (P==0){
return(0)
}else{
return(bisection(P,F0,k,r,T,c))}}
MarketVolatilities=c()
for(j in 1:22){
MarketVolatilities[j]=bisection_sigma(MarketPrices[j],CurrentForwardValue,MarketStrikes[j],R,T1,as.numeric(str_detect(data2103_31[j,1],"P")))
}
plot(MarketVolatilities,ylab="12月31日IO2103各行权价下的波动率",xlab="左半看涨右半看跌，行权价在同一类中依次递增")#volatility plot
ATMVolatility = MarketVolatilities[11]; #ATM implied volatility
#method for SABR model.
SABRVolitily= function(F0, k, T, alpha, beta, rho, nu)
{
z = nu/alpha*(F0*k)^((1-beta)/2)*log(F0/k)
x = log((sqrt(1-2*rho*z+z^2)+z-rho)/(1-rho))
sigma = alpha/((F0*k)^((1-beta)/2)*(1+(1-beta)^2/24*(log(F0/k))^2+(1-beta)^4/1920*(log(F0/k))^4))*z/x*(1+((1-beta)^2/24*alpha^2/(F0*k)^(1-beta)+0.25*rho*beta*nu*alpha/(F0*k)^(0.5*(1-beta))+(2-3*rho^2)*nu^2/24)*T)
return(sigma)
}
residual_square = function(parameter)
{
alpha=parameter[1]
beta=parameter[2]
rho=parameter[3]
nu=parameter[4]
residual = c()
for(i in 1:length(MarketStrikes))
{
residual[i] = SABRVolitily(CurrentForwardValue, MarketStrikes, T1, alpha, beta, rho, nu)[i]-MarketVolatilities[i]
}
result = sum(residual^2)
return(result)
}
EquivalentBlackVol = function(F0, k, T, alpha, beta, rho, nu)
{
EquivalentVol = c()
for(i in 1:length(MarketStrikes))
{
EquivalentVol = c(EquivalentVol, SABRVolitily(CurrentForwardValue, MarketStrikes[i], T1, alpha, beta, rho, nu))
}
return(EquivalentVol)
}
fit_result = optim(c(0.1,0.99,0.5,0.1),residual_square)
X = fit_result$par
Alpha1 = X[1];
Beta1 = X[2];
Rho1 = X[3];
Nu1 = X[4];
EquivalentVol=EquivalentBlackVol(CurrentForwardValue, MarketStrikes, T1, Alpha1, Beta1, Rho1, Nu1)
Option_value=c()
for(i in 1:(length(MarketStrikes)))
{
if(i<12){
Option_value = c(Option_value, BSM(CurrentForwardValue, MarketStrikes[i], R, EquivalentVol[i], T1)[1])
}else{
Option_value = c(Option_value, BSM(CurrentForwardValue, MarketStrikes[i], R, EquivalentVol[i], T1)[2])
}
}
plot(MarketStrikes[1:11], MarketVolatilities[1:11], col = 'red',xlab="12月31日IO2103不同行权价的看涨期权",ylab="波动率",ylim=c(0.1,0.25))
lines(MarketStrikes[1:11], EquivalentVol[1:11], col = 'cyan')#看涨期权的拟合
plot(MarketStrikes[12:22], MarketVolatilities[12:22], col = 'red',xlab="12月31日IO2103不同行权价的看跌期权",ylab="波动率",ylim=c(0.1,0.25))
lines(MarketStrikes[12:22], EquivalentVol[12:22], col = 'cyan')#看跌期权的拟合
plot(MarketStrikes[1:11], MarketPrices[1:11], col = 'red',xlab="12月31日IO2103不同行权价的看涨期权",ylab="期权价格",ylim=c(0,700))
lines(MarketStrikes[1:11], Option_value[1:11], col = 'cyan')
plot(MarketStrikes[12:22], MarketPrices[12:22], col = 'red',xlab="12月31日IO2103不同行权价的看跌期权",ylab="期权价格",ylim=c(0,700))
lines(MarketStrikes[12:22], Option_value[12:22], col = 'cyan')
FittingError_BS = function(sigma)
{
residual = c()
for(i in 1:(length(MarketStrikes)))
{
if(i<12){
residual[i] = BSM(CurrentForwardValue, MarketStrikes[i], R, sigma, T1)[1] - MarketPrices[i]
}else{
residual[i]=BSM(CurrentForwardValue, MarketStrikes[i], R, sigma, T1)[2] - MarketPrices[i]
}}
result = sum(residual^2)
return(result)
}
FittingError_SABR=function(value){
residual = c()
for(i in 1:(length(MarketStrikes)))
{
residual[i] = value[i] - MarketPrices[i]
}
result = sum(residual^2)
return(result)
}
BSVol = optimize(FittingError_BS, c(0, 1), tol=0.1)$minimum
Option_BS = c()
for(i in 1:length(MarketStrikes))
{
if(i<12){
Option_BS = c(Option_BS, BSM(CurrentForwardValue, MarketStrikes[i], R, BSVol, T1)[1])
}else{
Option_BS=c(Option_BS, BSM(CurrentForwardValue, MarketStrikes[i], R, BSVol, T1)[2])
}
}
plot(MarketStrikes[1:11], MarketPrices[1:11], col = 'red')
lines(MarketStrikes[1:11], Option_BS[1:11], col = 'navy')
lines(MarketStrikes[1:11], Option_value[1:11], col = 'Lightgreen')#Call option value contrary plot
plot(MarketStrikes[12:22], MarketPrices[12:22], col = 'red')
lines(MarketStrikes[12:22], Option_BS[12:22], col = 'navy')
lines(MarketStrikes[12:22], Option_value[12:22], col = 'Lightgreen')#Put option value contrary plot
FittingError1_BS = FittingError_BS(BSVol)
FittingError1_SABR=FittingError_SABR(Option_value)