オプションの価格を機械学習アプローチで求めてみる。

みなさん、ごきげんよう。次世代システム研究室のK.N.です。今日は、オプションの価格を通常とは少し違う機械学習アプローチで求めてみたいと思います。参考文献はIgorの論文になります。

そもそもオプションって？

詳しい説明は証券会社や取引所のHPでわかりやすく解説されているますので、ここでは概要だけ述べます。まず、オプションとはデリバティブとよばれる金融派生商品の内の一つです。派生というくらいなので、何か他の商品の将来の価値（価格）に紐付いてます。そのデリバティブの対象資産を原資産と呼びます。例えば、原資産には金利や為替、株価、経済指数、コモディティ、不動産価格が含まれます。また、原資産の候補は金融資産だけではなく、未来に渡って不確実な変数なら何でも候補になり得るため、天候の結果に依存した天候デリバティブとよばれるものも存在します。

加えて、原資産に対しての依存方法によって、商品性が別れます。代表的なデリバティブは他にも先物/先渡やスワップ取引があります。つまり、デリバティブは原資産ｘ商品性によって特定することができます。

そして、オプションとは「ある原資産を将来の特定の時点で、あらかじめ設定した行使価格で買う（または、売る）ことのできる権利」のことを言います。まず、ポイントとなるのは権利であって、義務ではないということです。つまり、権利を行使することが損失にならば、放棄されます。また、用語の意味として、将来の権利行使時点を「満期」、行使価格を「ストライク」と呼びます。

以下、簡単な例を使って、状況を説明しましょう。

ここでは、Kさんが報酬としてA社の株式で支払われたと仮定しましょう。そして、1年間は売却制限があるとします。株価が上がり続ければ良いですが、下がれば当然損失です。また、最終的に株価が今より高くなっていたとしても、その間の株価の上げ下げに一喜一憂するのも好ましくありません。そこで、ある程度の対価を払ってでも、株高の恩恵を受けつつ下落による損失を抑えることをオプションを買うことで実現させます。ここでは、「現時点で１株100円として、１年後にも100円で売ることができるプットオプション」を金融機関から購入します。ちなみに、売る権利のことを「プット」、買う権利のことを「コール」と呼びます。

例えば、１年後に１株80円になった場合、オプションが無かったら100-80=20円の損失ですが、オプション権利を行使することでA株を100円で売却できるます。つまり、オプションにより20円の利益を得たことを意味します。

一方、1年後に120円になった場合を考えます。もしオプションを行使したとすると、120円の価値のものを100円で売ることになってしまいます。そのような場合、オプション権利を行使する必要はありません。すなわち、オプションから利益はゼロですが、持っているA株は1株120円で売れるため、株高による利益を享受することができます。つまり、株価上昇/下落のいずれのシナリオにおいても損失は発生しませんでした。このように、資産価格の変動、つまり市場リスクを軽減する行為をヘッジと呼びます。

それでは、このプットオプションの価値はどう見積もれば良いでしょうか？少なくともタダでないことは明らかでしょう。株価がどのように動いても、オプションによる損失は発生しません。つまり、もしオプションが無料だとすると、元手がゼロで、損失する確率がゼロで確実に利益を得る（このような状況を裁定と呼ぶ）ことになるので、あり得ません。かといって、いくらでも高くてよいわけでもないでしょう。結局のところ、将来の原資産価格の変動の仕方、つまり未来の株価の確率分布によって決まると考えて良さそうです。

実際には、オプション価格を簡単に算出する大変有名な公式が存在します。考案者の名前を取ってブラック・ショールズ(・マートン）(以下、BS）の価格式と呼ばれます（一般的にはブラック・ショールズ方程式が有名ですが、実務的には、パラメーターを入れればオプション価値が出てくる価格式の方が便利です）

上式は、コール・オプションの式です。プット・オプションの価値Pは以下の関係式から簡単に決まります。

N()は標準正規分布の分布関数です。Stは原資産価格を表します。パラメータは、Tが満期、Kがストライク（行使価格）、rが無リスク金利、Sが現時点tでの原資産価格で、σが原資産のボラティリティ(時間単位の収益率の標準偏差）です。ここで、T, Kはオプションの条件で、r, S,σはマーケットから決まる値です。r, Sは現時点の値をそのまま使えるので良いですが、ボラティリティは未来の原資産の分布に関する量なので現時点で求める値には精度に限度があります。ボラティリティが分からないとオプション価格が決まらないという意味では、オプション価格をボラティリティというパラメーターで表しただけにも見えます。実際、マーケットで値が決まるオプション価格からBS公式を逆算して求めたボラティリティをインプライド・ボラティリティ(IV)と呼び、オプション価値の高低を判断するにはこちらのほうが都合が良いことが多いです。ですが、オプションの価格を少ないパラメータで簡単に計算することができ、過去の収益率の分散から計算したボラティリティを使っても、割と近い値になるという意味でBSの公式は画期的でした。

BSの価格式は、原資産の収益率が正規分布に従うという仮定の下で導出されています。次式のように書き表されます。これをBSモデルと呼びます。

dztは標準ブラウン運動の変動項です。μは原資産の期待収益率です。上式のBS方程式にμが含まれていないことに注意してください。オプション価値は原資産の期待収益率に依存しません。

しかしながら、原資産の収益率が正規分布に従うというのは計算都合上の仮定の話で、実際の原資産の収益率の分布は正規分布に対してずれがあります。横軸にストライクを、縦軸にIVをプロットすると、低いストライクと高いストライクの領域においてIVが割高になる現象が確認できます。これを、その形状からボラティリティ・スマイルと呼びます。オプションの条件を表すストライクによって、原資産の分布を司るパラメータが変化するのは変な話ですが、これは原資産の分布を無理やりBSモデルで記述することに由来します。より実際の原資産の振る舞いを再現できるような複雑なスマイル・モデルが様々提案されています。ちなみに、上のBS価格式は、行使機会が満期時のみのヨーロピアンオプションに対する式です。満期までの任意の時点で行使可能なオプションをアメリカンオプションと呼びます。

プライシング

以上前置きが長くなりましたが、ここから、実際にオプション価格を計算してみたいと思います。ここでは、BSモデルを仮定せず、多数発生させたモンテカルロを使って、オプションの価格を再現するような複製ポートフォリオを最適化で求める方法を考えます。詳細ははじめに挙げた論文を参考にしていただくとして、ここでは概要だけを述べます。（ちなみに、論文では強化学習(Q学習）の枠組みでオプションプライシングの方法を示してましたが、単なる教師あり学習（最適化）だけで済む話のように思えました。ここは、私もまだしっくりきていないというか、理解できていない箇所ですので、整理がついたら追記するかもしれません）

まず、複製ポートフォリオを次式で表します。Stは原資産価値、ut、Btは原資産、預金（無リスク資産）の保有量です。

満期時には、支払われるペイオフがオプション価値そのものになり、このとき所持している原資産は精算するとします。これは次式のように表されます。

そして、次式は、リバランス時の条件式です。リバランス時は外部からのお金の出し入れは無いとします。

これらの関係式から、次の関係が求まります。

上式より、もしポジション保有量{ut}がわかっていたら、原資産のパス(Stの列）に対して、満期から時間T→T-1→T-2と割り戻すことで、ヘッジポートフォリオ価値が計算できます。そして、ポジション保有量utは、モンテカルロパスごとに計算したヘッジポートフォリオの分散が最小となるように決定します。

上式はヘッジポートフォリオの変化量を原資産でヘッジした２乗誤差を表しています。また、utに対する２次式になっているので、解を解析的に計算できます。

utは原資産Stの関数になっているのですが、ここでStを別な確率変数Xtで表現し直します（標準化）

そして、関数utを以下のように、Xtを引数とするatで表されるとします。

これから、atの関数を決定するのですが、atはある基底関数上の１次関数で表現されるとします。

上式を使って、ポートフォリオの分散をモンテカルロパスに対してあらわに計算すると、下の式になります。

これはφに対して２次式となっているので、結果として、最小解は次の連立方程式を解く問題に帰着されます。

以上、モンテカルロパスを発生させ、後ろから逐次的に連立方程式を繰り返し、計算することでヘッジポジションの量、即ちオプションデルタが求まられるので、結果としてオプションの価値を得ることを示しました。このアプローチで特筆すべきは原資産過程の分布に対しては、何ら仮定をおいていないということです。つまり、モンテカルロサンプルだけではなく、ヒストリカルなデータに対しても適用できる可能性があります。その意味ではモデルフリーなプライシングであるので（基底関数に何を選ぶかいう恣意性はあるが）、幅広く応用できる可能性がありそうです。

とはいっても、まずここでは、簡単な検証として、BSモデルで発生させたモンテカルロパスに対して上記の機械学習アプローチでオプション価格を求め、それがBS価格式と整合性があるのか検証してみたいと思います。

数値検証

まず、BSモデルでモンテカルロパスの発生させます。

基本条件

np.random.seed(100)
S0 = 100      # initial stock price
mu = 0.05     # drift
sigma = 0.15  # volatility
r = 0.03     # risk-free rate
M = 1         # maturity

K = S0 #strike
N_STEPS = 250        # number of time steps (business dates per year)
N_MC = 10000  # number of paths

パス生成

def calc_BS_paths(s0, mu, sigma, maturity, N_step, N_path):
    dt = maturity / N_step
    paths= np.ones((N_path, N_step+1))*s0
    paths[:, 1:] = np.exp(np.random.randn(N_path, N_step) * sigma * dt**0.5  
                        + (mu - 0.5*sigma**2)* dt)
                                                                                                              
    paths= np.cumprod(paths, axis=1)
    return np.linspace(0, maturity, N_step+1), paths
ts, paths = calc_BS_paths(S0, mu, sigma, M, N_STEPS, N_MC)
for i in range(20): plt.plot(ts, paths[i,:])

BSモデルで計算したデルタで複製ポートフォリオ作成を試みます。

class BlackScholesEuropean(object):

    @staticmethod
    def calc_price(s, k, r, T, sigma, isCall):
        T = max(1e-8, T)
        vt = sigma * np.sqrt(T)
        gamma = np.exp(-r * T)
        d1 = (np.log(s/k) + (r + 0.5 * sigma**2) * T) / vt
        d2 =d1 - vt
        if isCall: price = s * norm.cdf(d1) - gamma*k * norm.cdf(d2)
        else: price = gamma*k*norm.cdf(-d2) - s * norm.cdf(-d1) 
        return price
    
    @staticmethod
    def calc_delta(s, k, r, T, sigma, isCall):
        T = max(1e-8, T)
        vt = sigma * np.sqrt(T)
        gamma = np.exp(-r * T)
        d1 = (np.log(s/k) + (r + 0.5 * sigma**2) * T) / vt
        if isCall: delta = norm.cdf(d1)
        else: delta = -norm.cdf(-d1)
        return delta

def calc_hedged_portfolio(tsteps, paths, tstep_delta, deltas, spot_opt_price):
    delta = deltas[0]
    deposit = spot_opt_price - delta*paths[0]
    hedged_pfs= np.zeros(len(paths))
    hedged_pfs[0] = spot_opt_price 
    delta_idx = np.searchsorted(tstep_delta, tsteps)
    for i in range(1, len(paths)):
        dt =tsteps[i] - tsteps[i-1]
        pf = delta*paths[i] + deposit * np.exp(r*dt)
        delta = deltas[delta_idx[i]-1]
        deposit = pf - delta*paths[i]
        hedged_pfs [i] = pf
    return hedged_pfs


idx = 32 # one path is chosen
deltas = [BlackScholesEuropean.calc_delta(s, K, r, M-ts[i], sigma, isCall=False) for i, s in enumerate(paths[idx])]
bs_opts = [BlackScholesEuropean.calc_price(s, K, r, M-ts[i], sigma, isCall=False) for i, s in enumerate(paths[idx])]
spot_opt_price = BlackScholesEuropean.calc_price(S0, K, r, M, sigma, isCall=False)
plt.plot(ts, bs_opts, label='BS')
for rehedge_inv in [5, 50]:
    pfs = calc_hedged_portfolio(ts, paths[idx], ts[::rehedge_inv], deltas[::rehedge_inv], spot_opt_price)
    plt.plot(ts, pfs, label='HEDGE({})'.format(rehedge_inv))
plt.legend()

 

上図は、発生させたBSモデルのあるモンテカルロパス（株価）に対して、ATMプットオプションの満期までの価格変化を表しています。青線はBS価格式で求められるプライス、オレンジ線、緑線はBSモデルのデルタから計算した複製ポートフォリオ価値になります。オレンジは5日ごと、緑は50日ごとリバランス（デルタの更新）をしています。緑線のほうが青線に対するずれが大きくヘッジ誤差が大きいことが分かります。ですが、概ねオプションの複製に成功していると言えるでしょう。

続いて、本題に入ります。今回紹介した方法でポジション量を計算し、オプションを複製できるか確認してみたいと思います。

class OptimizeHedgePortfolioFromMCPaths(object):
    
    def __init__(self, tsteps, paths, r, mu, sigma, reg_param=1e-3, order_basis=4, num_basis=12):
        
        self.reg_param = reg_param
        self.paths = paths
        self.tsteps = tsteps
        self.num_basis = num_basis
        self.num_path, self.num_step = self.paths.shape
        
        dts = tsteps[1:] - tsteps[:-1]
        self.dpaths = self.paths[:, 1:] - np.exp(r * dts) * self.paths[:, :-1]
        self.dpaths_hat = self.dpaths - np.mean(self.dpaths, axis=0)

        # state variable
        X = - (mu - 1/2 * sigma**2) * self.tsteps + np.log(paths)  
        X_min = np.min(X)
        X_max = np.max(X)

        # Spline basis of order p on knots k
        tau = np.linspace(X_min, X_max, self.num_basis)
        k = bspline.splinelab.aptknt(tau, order_basis)
        self.basis = bspline.Bspline(k, order_basis)  
        self.data_mat_t = np.array([[ self.basis(el) for el in X[:,i]] for i in range(self.num_step)]) #num_step x num_mc x num_basis
       
        self.gammas = np.exp(-r*dts)
         
    def _calc_Amat(self, t):
        pX = self.dpaths_hat[:,t,np.newaxis] * self.data_mat_t[t,:,:]
        A_mat = np.matmul(pX.T, pX) + self.reg_param * np.eye(self.num_basis) 
        return A_mat
        
    def _calc_Bvec(self, t,  pi_hat):
        
        # coef = 1.0/(2 * self.gammas[t]  * risk_lambda)
        coef = 0.0
        tmp = pi_hat  * self.dpaths_hat[:,t] + coef * self.dpaths[:, t]
        X_mat = self.data_mat_t[t, :, :]  
        B_vec = np.matmul(X_mat.T, tmp)
        return B_vec
           
    def calc_rollback(self, payoff_func):
        self.pi = np.zeros((self.num_path, self.num_step))
        self.opt_hedges, self.deposits = np.zeros((self.num_path, self.num_step)), np.zeros((self.num_path, self.num_step))
        self.phi = np.zeros((self.num_step, self.num_basis))
        self.deposits[:, -1] = np.array([payoff_func(term_price) for term_price in self.paths[:,-1]])
        
        for t in range(self.num_step-2, -1, -1):
            self.pi[:, t+1] = self.deposits[:, t+1] + self.opt_hedges[:, t+1] * self.paths[:, t+1]
            pi_hat = self.pi[:, t+1] - np.mean(self.pi[:, t+1])
            Amat = self._calc_Amat(t)
            Bvec = self._calc_Bvec(t, pi_hat)
            self.phi[t,:] = np.matmul(np.linalg.inv(Amat), Bvec)
            self.opt_hedges[:,t] = np.matmul(self.data_mat_t[t,:,:], self.phi[t,:])
            self.deposits[:, t] = self.gammas[t]*(self.deposits[:, t+1] + 
                                                  (self.opt_hedges[:,t+1]-self.opt_hedges[:,t])*self.paths[:,t+1])
            #self.pi[:,t] = self.gammas[t] * (self.pi[:,t+1] - self.opt_hedges[:,t] * self.dpaths[:,t])
        self.pi[:, 0] = self.deposits[:, 0] + self.opt_hedges[:, 0] * self.paths[:, 0]
           
    def calc_pos(self, new_paths):
        # time grid must be identical for that used in the training
        new_x= - (mu - 1/2 * sigma**2) * self.tsteps + np.log(new_paths)
        psi_x = np.array([[self.basis(el) for el in new_x[:, i]] for i in range(new_paths.shape[1])]) #num_step x num_mc x num_basis
        pos = np.matmul(psi_x, self.phi[:,:,np.newaxis])
        pos = np.squeeze(pos).T
        return pos

ts_train, paths_train = calc_BS_paths(S0, mu, sigma, M, 10, 10000) # path for training
obj = OptimizeHedgePortfolioFromMCPaths(ts_train, paths_train, r, mu, sigma) 
def terminal_payoff(x): return max(K - x, 0) # put payoff
obj.calc_rollback(terminal_payoff) #calc optimized position 

mc_price_mean, mc_price_std = np.mean(obj.pi[:,0]), np.std(obj.pi[:,0])
bs_price = BlackScholesEuropean.calc_price(S0, K, r, M, sigma, isCall=False)
print("Replicating portfolio value: MEAN: {}, STD: {}".format(mc_price_mean, mc_price_std))
print("BS Model option value {}".format(bs_price))

上の実行結果は以下のようになります。
Replicating portfolio value: MEAN: 4.499, STD: 1.738
BS Model option value 4.530
割と再現できているようです。

次に、デルタの値を比べてみましょう。trainingの場合とは異なるpathで評価します。

ts_test = ts[::25] # time step is adjusted to match that used in the training.
paths_test = paths[:,::25]
pos_test = obj.calc_pos(paths_test)

idx=32
bs_pos = [BlackScholesEuropean.calc_delta(s, K, r, M-ts_test[i], sigma, isCall=False) for i, s in enumerate(paths_test[idx])]
plt.plot(ts_train[:-1], pos_test[idx, :-1], 'o-', label="MC")
plt.plot(ts_train[:-1], bs_pos[:-1], 'o-', label="BS")
plt.legend()

デルタ値も概ね一致している様子が分かります。先と同様に複製ポートフォリオの評価をしてみましょう。

pfs_mc =  calc_hedged_portfolio(ts, paths[idx], ts_test, pos_test[idx], spot_opt_price)
pfs_bs =  calc_hedged_portfolio(ts, paths[idx], ts_test, bs_pos, spot_opt_price)

plt.plot(ts, bs_opts, label="BS(MODEL)")
plt.plot(ts, pfs_bs, label="BS(HEDGE)")
plt.plot(ts, pfs_mc, label="MC(HEDGE)")
plt.legend()

オレンジ線がBSデルタから、緑線がMCから学習したデルタでの複製ポートフォリオですが、概ね一致している様子が分かります。
実は、これはちょっと出来過ぎです。下図の違うパターン（idx=20)では、少し複製誤差が見られます。

最後に、100本のパスにおける両者(BSデルタ、MCデルタ）のヘッジ誤差を確認したいと思います。

bs_pos = [[BlackScholesEuropean.calc_delta(s, K, r, M-ts_test[i], sigma, False) for i, s in enumerate(paths_test[ipath])] for ipath in range(100)]
pfs_bs =  np.array([calc_hedged_portfolio(ts, paths[ipath], ts_test, bs_pos[ipath], spot_opt_price) for ipath in range(100)])
pfs_mc =  np.array([calc_hedged_portfolio(ts, paths[ipath], ts_test, pos_test[ipath], spot_opt_price) for ipath in range(100)])
bs_prices = np.array([[BlackScholesEuropean.calc_price(s, K, r, M-ts[i], sigma, False) for i, s in enumerate(paths[ipath])] for ipath in range(100)])
mse_bs =  np.mean((pfs_bs - bs_prices)**2, axis=1)
mse_mc =  np.mean((pfs_mc - bs_prices)**2, axis=1)
print("Mean Squared Error:  BSDelta  {:.3f}, MCDelta {:.3f}".format(np.mean(mse_bs), np.mean(mse_mc)))

上コードを実行すると、下記のようになります。概ねBSデルタと同様の精度です。

Mean Squared Error: BSDelta 0.941, MCDelta 1.052

まとめ

以上、通常のBSモデルとは異なるアプローチでオプション のプライシングを試みました。今回の数値検証では、BSモデルで生成したモンテカルロパスで学習させるとBSデルタを十分に再現できることが分かりました。BSモデルのパラメータは、ボラティリティ一つであるの簡潔なモデルであるのに対して、今回の機械学習アプローチは、リバランス回数 x 基底関数の数個のパラメータを必要とする分モデルの複雑度が増すため、表現力があるのは当然といえます。しかし、この方法では資産の分布に対する仮定を置いていないため、より複雑なモデルから生成されたサンプルパス、ないしはヒストリカルデータを直接使って、同様のアプローチで評価することができるため、広く応用できる可能性があります。今後は、そのようなデータで本当に応用可能か検討したいと思います。

最後に

次世代システム研究室では、ビッグデータ解析プラットホームの設計・開発を行うアーキテクトとデータサイエンティストを募集しています。興味を持って頂ける方がいらっしゃいましたら、ぜひ 募集職種一覧からご応募をお願いします。