R codes

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)