SVモデルでのオプションのプライシング

ボラティリティを予測して儲けましょうというお話。

本当は資産の上がり下がり（リターン）を見極めて運用したいのだけれども、将来の上がり下がりの予想は価格に織り込まれている（つまり上がりそうなものは既に高い価格になってしまっている）からなかなかそうもいかない。
そこで資産価格の変動のしやすさ（ボラティリティ）に注目する。実はデリバティブ（たとえばオプション）を使えば、ボラティリティを正しく予想することで儲けることができる。たとえばコールオプションとプットオプションを買い建てるロングスストラドル・ストラングル戦略(【参考】http://lounge.monex.co.jp/advance/fop/strategy-long-straddle.html)をとれば、上がり下がりに関わらず大きく価格が変動しさえすれば儲けることができる。

オプションのプライシングにはブラックショールズモデルが一般的に使われているが、ボラティリティに関して非常に単純な仮定を置いている。そこで、確率的ボラティリティ変動モデルという少し込み入ったモデルでボラティリティを予測して、オプションが買いか売りかを判断しましょうというのが本題。

従来のブラックショールズモデルとの違いは

ブラックショールズモデルはボラティリティを定数と仮定
確率的ボラティリティ変動モデルはボラティリティがランダムウォークすると仮定

$y_t=log(x_t/x_{t-1}) \\y_t =exp(h_t/2)\epsilon_t \\h_{t+1}=\mu+\phi(h_t-\mu)+\sigma_{\eta}\eta_t \\\eta_t,\epsilon_t \sim_{i.i.d}N(0,1)$

ここから少し技術的になるが、プライシングは

SVモデルの潜在パラメーター推定自体は時間がかかる（かつ学習期間によって急激に変わることは想定されない）ので、先行研究から借用し([http://www.imf.org/external/pubs/ft/wp/2003/wp03125.pdf:title=http://www.imf.org/external/pubs/ft/wp/2003/wp03125.pdf])
状態変数（すなわち現在のボラティリティ）の推定のみを粒子フィルターで行う

【参考】

お前の逐モン、GETだぜ！ from tera monagi

将来の資産価格の変動(および満期日でのオプションの行使結果)については1000個の粒子でボラティリティの変動を考慮してシミュレーションし、リスク中立で単なる期待割引現在価値で測定した

原資産(日経225)の推移

ボラティリティの推移

現在ボラティリティは小さく、長期的には本来の水準(より高い水準)に戻ると予想される。

※logスケールなので、リターンの標準偏差に変換するには、sd=exp(h/2)とする必要。直近100営業日のヒストリカルからのフィルタリング結果と、次の100営業日の予測値。真ん中のラインが期待値で、上下がそれぞれ95%信用区間

オプションの理論価格

sv=確率的ボラティリティ変動モデル,bs=ブラックショールズ
実際の価格は大阪証券取引所の指数日報(http://www.ose.or.jp/market/daily_report/pdf/20140704/n09.pdf)で見ることができる。
行使日が7月の場合ほとんど差がないが、8月の場合かなり差がある。ボラティリティの上昇予想を反映して、svモデルではより高い理論価格をつけている。

2014-07-11　がSQ(清算日)のコールオプションの理論価格

## | strike| sv| bs|
## |------:|---------:|---------:|
## | 11250| 4195.7330| 4190.9821|
## | 11500| 3945.8186| 3941.0677|
## | 11750| 3695.9042| 3691.1533|
## | 12000| 3445.9898| 3441.2389|
## | 12250| 3196.0754| 3191.3245|
## | 12500| 2946.1610| 2941.4101|
## | 12750| 2696.2466| 2691.4957|
## | 13000| 2446.3322| 2441.5813|
## | 13250| 2196.4178| 2191.6669|
## | 13500| 1946.5034| 1941.7525|
## | 13750| 1696.5890| 1691.8381|
## | 14000| 1446.6746| 1441.9240|
## | 14250| 1196.7602| 1192.0200|
## | 14500| 947.0345| 942.3026|
## | 14750| 701.2854| 694.4499|
## | 15000| 467.3743| 456.9516|
## | 15250| 269.0365| 252.3625|
## | 15500| 129.4894| 109.1894|
## | 15750| 52.0619| 34.8709|
## | 16000| 17.3579| 7.8874|
## | 16250| 4.3944| 1.2344|
## | 16500| 1.2284| 0.1324|
## | 16750| 0.4325| 0.0097|
## | 17000| 0.1636| 0.0005|
## | 17250| 0.0000| 0.0000|
## | 17500| 0.0000| 0.0000|
## | 17750| 0.0000| 0.0000|
## | 18000| 0.0000| 0.0000|
## | 18250| 0.0000| 0.0000|
## | 18500| 0.0000| 0.0000|
## | 18750| 0.0000| 0.0000|
## | 19000| 0.0000| 0.0000|
## | 19250| 0.0000| 0.0000|

2014-08-08 の分のコールオプション

## | strike| sv| bs|
## |------:|---------:|---------:|
## | 11250| 4233.3824| 4206.3772|
## | 11500| 3983.8101| 3956.8049|
## | 11750| 3734.2378| 3707.2326|
## | 12000| 3484.6655| 3457.6604|
## | 12250| 3235.0932| 3208.0882|
## | 12500| 2985.5209| 2958.5169|
## | 12750| 2736.1552| 2708.9507|
## | 13000| 2487.3576| 2459.4089|
## | 13250| 2239.4128| 2209.9656|
## | 13500| 1992.7853| 1960.8567|
## | 13750| 1750.0035| 1712.7134|
## | 14000| 1513.0153| 1466.9587|
## | 14250| 1282.6914| 1226.3058|
## | 14500| 1062.5552| 995.1290|
## | 14750| 859.2617| 779.3620|
## | 15000| 673.7090| 585.6741|
## | 15250| 512.9096| 420.0272|
## | 15500| 377.2256| 286.1252|
## | 15750| 269.1247| 184.4295|
## | 16000| 187.4256| 112.1631|
## | 16250| 127.6027| 64.2334|
## | 16500| 85.9461| 34.5990|
## | 16750| 56.4683| 17.5210|
## | 17000| 36.6471| 8.3425|
## | 17250| 22.4926| 3.7369|
## | 17500| 13.5932| 1.5762|
## | 17750| 7.4098| 0.6267|
## | 18000| 3.3847| 0.2353|
## | 18250| 1.3087| 0.0835|
## | 18500| 0.6687| 0.0281|
## | 18750| 0.3579| 0.0089|
## | 19000| 0.1084| 0.0027|
## | 19250| 0.0000| 0.0008|

ということで、ボラティリティが本来の水準に戻るべく上昇していくにも関わらず、市場は織り込んでいないようなので、SVモデルによると「コールオプションの買い」という結論になる。

以下コード。

### Price of European Option

### Under Stochastic Volatility Model Using Particle Filter
## Model

# y_t=log(x_t/x_{t-1})

# y_t~N(0,exp(h_t/2))*epsilon

# h_{t+1}=mu+phi(h_t-mu)+sig.eta*epsilon.eta

# epsilon,epsilon.eta~i.i.d N(0,1)
library(quantmod)

library(reshape2)

library(scales)

library(ggplot2)
# number of particle

Nx=1000
# sv parameters

# source: http://www.imf.org/external/pubs/ft/wp/2003/wp03125.pdf

sv_params=list(

  mu= -8.55,

  phi=0.944,

  sig.eta=0.178

  )
# simulator
state_sim=function(particles,mu,phi,sig.eta,Nx){

  mu+(particles-mu)*phi+sig.eta*rnorm(Nx)

}
measurement_eq=function(particles,obs,Nx){

  volatility=exp(particles/2)

  dnorm(obs/volatility)/volatility

}
pf_sim=function(dat,sv_params,Nx,pred_t){

  

  #sv params

  mu=sv_params$mu

  phi=sv_params$phi

  sig.eta=sv_params$sig.eta

  particles=replicate(Nx,mu)

  

  volatility.hist=xts(matrix(0,ncol=3,nrow=nrow(dat)),index(dat))
  for(i in 1:nrow(dat)){

    new_particles=state_sim(particles,mu,phi,sig.eta,Nx)

    obs=dat[i]

    particles=sample(new_particles,prob=measurement_eq(new_particles,obs,Nx),replace=T)

    volatility.hist[i,]=c(mean(particles),quantile(particles,c(0.05,0.95)))

  }

  

  result=list()

  result$volatility.hist=volatility.hist

  

  pred_p=index(dat)[nrow(dat)]+(1:(pred_t*3))

  pred_p=pred_p[weekdays(pred_p)%in%c("月曜日","火曜日","水曜日","木曜日","金曜日")][1:pred_t]
  

  result$volatility.pred=xts(matrix(0,ncol=3,nrow=pred_t),pred_p)

  result$return.pred=xts(matrix(0,ncol=Nx,nrow=pred_t),pred_p)

  

  for(i in 1:nrow(result$return.pred)){

    particles=state_sim(particles,mu,phi,sig.eta,Nx)

    mhl=c(mean(particles),quantile(particles,c(0.05,0.95)))

    result$volatility.pred[i,]=mhl

    result$return.pred[i,]=exp(particles/2)*rnorm(Nx)

  }

  result$return_cum.pred=cumsum(result$return.pred)

  

  result

}
# simulate option

sim_option=function(underlying_name,pred_t=100,CP="C",r=0.01,prev_t=100){

  getSymbols(underlying_name)

  underlying_dat=get(gsub("\\^","",underlying_name))

  underlying_return=dailyReturn(underlying_dat,type="log")  

  result=pf_sim(underlying_return,sv_params,Nx,pred_t)

  
  prev_p=index(underlying_return)

  prev_p=prev_p[(length(prev_p)-prev_t+1):length(prev_p)]

  pred_p=index(result$volatility.pred)

  s=as.numeric(tail(Cl(underlying_dat),1))

  

  # plot underlying

  #underlying.plot=data.frame(list(t=index(underlying_dat),ret=as.numeric(Cl(underlying_dat))))

  #result$underlying.plot=ggplot(underlying.plot,aes(x=as.Date(t),y=ret))+geom_line()+ scale_x_date(labels = date_format("%m/%d"))

  

  result$underlying.plot=underlying_dat[prev_p,]
  

  # plot volatility

  volatility.plot=data.frame(rbind(result$volatility.hist[prev_p],result$volatility.pred[pred_p]))

  volatility.plot$t=rownames(volatility.plot)

  volatility.plot <- melt(volatility.plot,id="t",measure=c("X1","X2","X3"))

  result$volatility.plot=ggplot(volatility.plot,aes(x=as.Date(t),y=value,group=variable,color=variable))+geom_line()+ scale_x_date(labels = date_format("%m/%d"))+ geom_vline(xintercept = as.numeric(as.Date(prev_p[length(prev_p)])),linetype=4)

  

  # date of strike

  d=as.data.frame(list(index=pred_p,month=as.factor(months(pred_p)),d=as.numeric(pred_p),dofw=weekdays(pred_p)))

  selected=d[d$dofw=="金曜日" & as.numeric(date_format("%d")(d$index))%in%8:14,"index"]

  

  result$theoretical_price.table=list()

  result$option.plot=list()

  

  for(i in 1:length(selected)){

    result_index=as.character(selected[i])

    num_of_days=length(pred_p[pred_p<=selected[i]])

    

    # define pricing function in standardized form

    pricing=function(rn_value,r,K,term){

      if(CP=="C"){

        sapply(rn_value,function(x){max(x-K,0)/(1+r)^term})

      }else{

        sapply(rn_value,function(x){max(K-x,0)/(1+r)^term})

      }

    }  

    BS_pricing=function(S,K,q,r,sigma,term){

      d1=(log(S/K)+(r-q+.5*sigma^2)*term)/(sigma*term^.5)

      d2=d1-sigma*term^.5

      exp(-1*q*term)*S*pnorm(d1)-exp(-1*r*term)*K*pnorm(d2)

    }

    

    f=function(x){num_of_days;rn_value=exp(result$return_cum.pred[selected[i],]);

                  sapply(x/s,function(y){s*mean(pricing(rn_value,r/365,y,num_of_days))})}

    fbs=function(x){num_of_days;volatility=exp(result$volatility.hist[nrow(result$volatility.hist),1]/2);

                    sapply(x/s,function(y){s*mean(BS_pricing(1,y,0,r/365,volatility,num_of_days))})}

    

    strikes=  floor(s/250)*250+(-16:16)*250

    

    

    result$option.plotresult_index=ggplot(data.frame(list(x=s*c(0.75,1.25))),aes(x=x))+

      stat_function(fun=f,n=30,aes(color="sv"))+

      stat_function(fun=fbs,n=30,aes(color="bs"))+

      geom_vline(xintercept = s)

    

    result$theoretical_price.tableresult_index=data.frame(list(strike=strikes,sv=f(strikes),bs=fbs(strikes)))

  }

  result$result_indices=as.character(selected)

  result

}
test=function(){

  sim_option("^N225")

}

レポート用

Daily Report `r format(Sys.time(), "%Y %b%d%a %H:%M:%S")`

========================================================
```{r ,echo=FALSE,warning=FALSE}

source("./pf.R")

library("xtable")

library(quantmod)

result=sim_option("^N225",pred_t=100,r=0.025,prev_t=100)
#chartSeries(result$underlying.plot)

#chartSeries(tail(result$underlying.plot,10))
kable(tail(underlying.plot,10))


plot(result$volatility.plot)
for(i in result$result_indices){

  print(i)

  plot(result$option.ploti)

  kable(result$theoretical_price.tablei)

}

```{r}

メールで配信

library("EasyHTMLReport")
easyHtmlReport(rmd.file="report.Rmd",

               from="hogehoge@piyo.com",

               to="piyopiyo@hoge.com",subject="DailyFinancialReport")