最小二乗モンテカルロ( Least Squares Monte Carlo）法でアメリカンオプションの評価
Longstaff and Schwartzが開発した最小二乗モンテカルロ法( Least Squares Monte Carlo)でアメリカンオプションのプライスを出してみたというお話。
numpy使えば楽に書けるかなと思ったけど意外とめんどかった。
答えはWolfram Alphaのとチェック。だいたい同じであることを確認（Wolfram|Alpha: Computational Intelligence）
計算条件は
| option | vanilla |
| option name | American |
| option type | call |
| strike price | Au$50 (Australian dollars) |
| time to expiration | 1 year |
| underlying price | Au$50 (Australian dollars) |
| volatility per period | 20% |
| dividend yield | 0% |
| risk-free interest rate | 1% |
としておいた。なぜかデフォルトの通貨がオージードルになる・・・
んで、Wolfram Alphaが弾き出した値（リスク指標付）は
| value | delta | gamma | vega | theta | rho |
| Au$4.24 | 0.560 | 0.039 | 19.850 | -2.217 | 23.802 |
一方、私が出したプライスは
4.2009
となってだいたいよろしい感じか。
以下ソース、なんかおかしい箇所があったら教えてください。
# -*- coding: utf-8 -*- import numpy #株価パスの生成 def create_path(current_price, risk_free, volatility, maturity, num_of_step, num_of_path, seed): #時間刻み dt = maturity / num_of_step #乱数生成器 rand = numpy.random.mtrand.RandomState(seed) #確率過程のドリフト drift = (risk_free - 0.5 * volatility**2.0) * dt #標準正規乱数の生成 norm_rand = rand.standard_normal(num_of_step * num_of_path).reshape((num_of_path, num_of_step)) #１ステップだけ時間進める演算子的なもの＆その累積 operateor_one_step = numpy.exp(drift + volatility * numpy.sqrt(dt) * norm_rand) operateor_cumulated = numpy.cumprod(operateor_one_step, axis = 1) return current_price * operateor_cumulated #回帰用の説明変数の設定。ハルの本に合わせてある def create_regressors(x): return numpy.vstack([numpy.ones(len(x)), x, x**2]).T #オプションを行使すべきかそのまま持つべきかを最小二乗モンテカルロで推定 def execute_or_not(cash_flow_execute, cash_flow_continue, cash_flow_continue_future, stock_price, discount_factor): #回帰係数を決定＆継続価値の予測値の算出。回帰には０以上の値の所だけ使う index_positive_cash_flow = cash_flow_execute > 0 regressors = create_regressors(stock_price) regressand = cash_flow_continue_future * discount_factor coefficient = numpy.linalg.lstsq(regressors[index_positive_cash_flow], regressand[index_positive_cash_flow])[0] cash_flow_continue_regressed = numpy.dot(regressors, coefficient.T) #行使した方が所持し続けるより儲かるなら行使するフラグ return index_positive_cash_flow & (cash_flow_execute > cash_flow_continue_regressed) def price_american_call_option(current_price, risk_free, volatility, strike_price, maturity, num_of_step, num_of_path, seed): #株価パスの生成 stock_prices = create_path(current_price, risk_free, volatility, maturity, num_of_step, num_of_path, seed) #生成する行列のサイズ size = (num_of_path, num_of_step) #時間刻み dt = maturity / num_of_step #オプションを実行した場合のキャッシュフロー cash_flow_execute = numpy.zeros(size) #オプションを持ったままにした場合のキャッシュフロー cash_flow_continue = numpy.zeros(size) #オプションを行使すべきか否かのフラグ flag_execute = numpy.zeros(size) #行使時のCF・継続価値のCF・行使すべきかフラグの満期時点での初期化 cash_flow_execute[:, -1] = numpy.maximum(stock_prices[:, -1] - strike_price, 0) cash_flow_continue[:, -1] = numpy.maximum(stock_prices[:, -1] - strike_price, 0) flag_execute[:, -1] = 1 #満期時点-1時間単位からスタートして時刻dtまでループ for col in xrange(num_of_step - 2, 0, -1) : stock_price = stock_prices[:, col] cash_flow_execute[:, col] = numpy.maximum(stock_price - strike_price, 0) is_executed = execute_or_not(cash_flow_execute[:, col], cash_flow_continue[:, col], cash_flow_continue[:, col + 1], stock_price, numpy.exp(-risk_free * dt)) #各パスに対して #-行使する場合は行使のCFを代入＆行使フラグの位置を現在の位置に修正 #-行使しない場合は継続価値（1期間割引き） #を代入 for row in xrange(0, num_of_path): if(is_executed[row]): cash_flow_continue[row, col] = cash_flow_execute[row, col] flag_execute[row, col] = 1 flag_execute[row, (col + 1):] = 0 else: cash_flow_continue[row, col] = cash_flow_continue[row, col + 1] * numpy.exp(-risk_free * dt) #ディスカウントファクター（全時間分）算出 discount_factor = numpy.cumprod(numpy.exp(-risk_free * dt) * numpy.ones(size), axis = 1) #各パスの行使価値（payoff)の現在価値（ディスカウント）を行使フラグに応じた時点で算出。平均値を計算して価格とする return numpy.sum(numpy.maximum(stock_prices - strike_price, 0) * discount_factor * flag_execute) / num_of_path if __name__ == '__main__': #現時点の株価 current_price = 50.0 #リスク・フリーレート risk_free = 0.01 #インプライド・ボラティリティ volatility = 0.2 #ストライク strike_price = 50.0 #オプションの満期 maturity = 1.0 #ステップ数（＝満期を何ステップに分割してシミュレーションするのか） num_of_step = 250 #モンテカルロパス数 num_of_path = 5000 #乱数のシード seed = 100 #アメリカンオプションの計算結果 print price_american_call_option(current_price, risk_free, volatility, strike_price, maturity, num_of_step, num_of_path, seed)
