請大神講一講美式期權定價的方法？（看成同一歐式期權和early exercise premium兩部分

01-30

首先你要弄明白什麼是美式期權。別以為很簡單，很多金數人到畢業也沒懂。比如我問你，你用pde來定價，為什麼每一步要加那個不等式比較？用LSMC，為什麼要compare holding和exercising value？

美式期權因為提前行權的特性，他的定價公式不再是簡單的Q下Expectation了，而是多了一個sup。這類問題從本質上講是一個dynamic programming principle。這保證了它可以被一個dynamic programming的比較來solve出來。這解答了為什麼pde解法和lsmc都要在比較value，其實都是利用了dynamic programming。

最普遍採用的numerical method是pde和lsmc，前者好處是利於得到Greek，後者好處是可以同時得到optimal boundary，而且可以處理高維和correlated問題。除此以外，還有tree，not necessarily binomial，以及stochastic mesh，其實個人感覺pde和lsmc都算是stochastic mesh的特例。而如今，還有一種artificial neutral network的定價方法，業內也有研究。

美式期權定價有很多方法。

具體有：

1. 二叉樹。這個比較簡單直觀，但是有些greeks的計算有些問題，並不很好。

2.三叉樹。

3.BAW。是一個從歐式期權定價模型上進行調整，由PDE推出來的公式。其間做了一些簡化，適用於正利率和短期以及長期期權。

4.LSM。一個神奇的Least Square Monte Carlo模型。避免了nested MC的問題。

5. Dual Formula.選取適當的輔助變數然後進行計算，是一個比較複雜的模型

去搜longstaff的論文，用mc做美式期權定價，從最後一期行權收益一步步做回歸，每一個時間點比較。

Valuing American Options by Simulation: A Simple Least-Squares Approach
Francis A. Longstaff
UCLA
Eduardo S. Schwartz
UCLA

二叉樹就可以辦到

early exercise 必定是某一點期權行權價值大於等待到期末價值（僅限us put），那麼以此點較大的行權價值折現，肯定比用稍小的等待價值折現要大。thats where the premium comes from

近期國內商品期貨期權模擬交易開始, 不同於股票期權, 期貨期權類型為美式. 本來想申請一個個人專欄來寫點定價方面的東西, 目的也是為了鞏固所學. 但由於某些原因被知乎暫停開設個人專欄的權利(原因下面有). 美式常用定價模型有:

1, 最小二乘蒙特卡洛模擬(LSM);

2, BAW;

3, 二叉樹.

二叉樹從精度上來講是最好的, 理論本身也十分簡單明了, 而且同遞歸演算法天然匹配, 但計算速度太慢, 簡直是暴力求解了(再次表明, 暴力model不可取). 大連商品交易所採取的是BAW模型, 這裡就主要講講此模型.

BAW把美式期權價格表示為:

$American(S, T) = European(S, T) + w(S, T)$

其中S為標的價格, T為距離到期時間, w(S, T)為early exercise premium.

由於美式和歐式期權價格都滿足Black-Scholes偏微方程

$frac{sigma ^{2}S^{2}}{2} V_{SS} +bSV_{S} -rV+V_{t}=0$

其中b為標的的cost of carry, r為無風險利率.

故w亦滿足BS方程, BAW模型把w作如下表示

$w(S, T)=(1-e^{-rT})f(S, 1-e^{-rT})$

這個表示很關鍵, 因子 $1-e^{-rT}$ 的來源可以這樣理解: T時間後得到1元, 現值只值 $e^{-rT}$ 元, 假設我們要求立馬兌現這1元錢, 那麼此時有 $1-e^{-rT}$ 的early exercise premium.

在分離出 $1-e^{-rT}$ 這個因子後, w所滿足BS偏微分方程近似成了一個常微分方程. 實際上, 當期權標的的cost of carry為0時, 我們甚至可對w作分離變數處理, 令 $w(S, T)=(1-e^{-rT})f(S)$ , 這樣做的理由是當標的沒有持有成本時, early exercise premium的時間因子部分完全可由現金流貶值來刻畫! 為了說明這點我們回過頭來看BAW模型, 其把w表示為

$w(S, T) = K(T)f(S, T)$

通過取 $K(T)=1-e^{-rT}$ , BAW模型中把BS方程近似為

$left( 1 ight) S^{2}f_{SS} + NSf_{S}-frac{M}{K}f=0$

其中N為依賴cost of carry的常數, M則是與cost of carry無關的常數. 若我們假設cost of carry為0, 則N=0, 設 $w(S, T) = K(T)f(S)$ , 此時BS方程為

$left( 2 ight) S^{2}f_{SS} -M(1+frac{K_{T}}{rK})f=0$

若我們在N等於0時直接代入方程(1), 並令方程(1)與方程(2)相同, 即

$frac{1}{K} =1+frac{K_{T}}{rK}$

此時可以發現 $K=1-e^{-rT}$ 恰為這個微分方程的解.

言歸正傳! 從常微分方程中得到f的形式解後, 需利用邊界條件來確定唯一解, 實際上, BAW模型中該解唯一依賴的未知量是美式期權行權時刻標的資產價格 $S^{ast }$ ,我們知道美式的行權時間為停時, 即應該在美式期權價值第一次等於內在價值時行權, 也就是American(S, T)與max(S - StrikePrice, 0)相切的點執行. 自然, 我們就可選擇牛頓迭代法來數值計算這個點. 至此, 整個BAW模型就講完了.

最後Talk is cheap. Show me the code. 我用C++來實現這個演算法, 當作一個小小的練習. 建立兩個頭文件, 一個是關於定價和計算波動率用到的數學工具, 一個是定價模型本身. C++知識點主要用了三個(言必稱三):

1, lambda表達式和template, 此處運用只為形式上簡潔一些;

2, namespace, 為的是防止命名衝突;

3, function object, 不使用函數而使用函數對象一是為了代碼表現得更像C++而不是C, 更重要的是BAW模型有一些臨時變數的計算很繁瑣, 並且有共同變數, 故可用函數對象直接封裝起來.

頭文件norm_distribution.h

#pragma once


#ifndef LIUTING_NORM_DISTTRIBUTION_H

#define LIUTING_NORM_DISTTRIBUTION_H
#include &

#include &
namespace AmurTiger

{

	const double PI = 3.141592653589793238462643;
	auto NormDensity = [](const double d) { return (1.0 / sqrt(2.0 * PI)) * exp(-0.5 * d * d); };

double NormDist(const double d) { if (d &> 6.0) { return 1.0; }; if (d &< -6.0) { return 0.0; }; double b1 = 0.31938153; double b2 = -0.356563782; double b3 = 1.781477937; double b4 = -1.821255978; double b5 = 1.330274429; double p = 0.2316419; double c2 = 0.3989423; double a = fabs(d); double t = 1.0 / (1.0 + a*p); double b = c2 * exp((-d) * (d / 2.0)); double n = ((((b5 * t + b4) * t + b3) * t + b2) * t + b1) * t; n = 1.0 - b * n; if (d &< 0.0) n = 1.0 - n; return n; }; template& T CalcEXP(const std::vector& v) { size_t sz = v.size(); if (sz == 0) return 0; T s = 0; for (auto i : v) s += i; return s / sz; } template& T CalcVAR(const std::vector& v) { size_t sz = v.size(); if (sz &<= 1) return 0; T e = CalcEXP(v); T s = 0; for (auto i : v) s += (i - e) * (i - e); return s / (sz - 1); } } #endif // !LIUTING_NORM_DISTTRIBUTION_H

頭文件option.h

#pragma once


#ifndef LIUTING_OPTION_H

#define LIUTING_OPTION_H
#include &

#include "norm_distribution.h"
namespace AmurTiger

{

	enum class CALLPUT

	{

		CALL,

		PUT

	};

	enum class UNDERLYING {

		COMMODITY,

		FUTURE

	};

	enum class OPTIONTYPE {

		EUROPEAN,

		AMERICAN

	};
	struct OPTION {

		UNDERLYING underlying;

		OPTIONTYPE option_type;

		CALLPUT call_put;

		double underlying_price;//標的資產現價

		double strike_price;//執行價

		double sigma;//波動率

		double t;//距離到期時間

		double r;//無風險利率

		double b;//cost of carry, 若無分紅等, 持有成本一般等於無風險利率r, 期貨持有成本為0
		OPTION() = delete;

		OPTION(UNDERLYING u, OPTIONTYPE ot, CALLPUT cp, double up, double sp, double sig, double x, double y, double z) :

			underlying(u), option_type(ot), call_put(cp), underlying_price(up), strike_price(sp), sigma(sig), t(x), r(y), b(z) {}

		OPTION(const OPTION op) :

			underlying(op.underlying), option_type(op.option_type), call_put(op.call_put),

			underlying_price(op.underlying_price), strike_price(op.strike_price), sigma(op.sigma),

			t(op.t), r(op.r), b(op.b) {}

	};
	class BLACKSCHOLES {//只針對歐式期權

	public:

		double operator()(const OPTION op)

		{

			double _S = 1, _t = sqrt(op.t), d1 = 0, d2 = 0;
			if (op.underlying == UNDERLYING::COMMODITY) _S = op.underlying_price;

			else _S = op.underlying_price * exp(-op.b * op.t);
			d1 = (log(_S / op.strike_price) + (op.b + 0.5 * op.sigma * op.sigma) * op.t) / (_t * op.sigma);

			d2 = d1 - _t * op.sigma;
			if (op.call_put == CALLPUT::CALL)

				return _S * exp((op.b - op.r) * op.t) * NormDist(d1) - op.strike_price * exp(-op.r * op.t) * NormDist(d2);

			else

				return -_S * exp((op.b - op.r) * op.t) * NormDist(-d1) + op.strike_price * exp(-op.r * op.t) * NormDist(-d2);

		}

	};

class BARONEADESIWHALEY {//歐美皆可 public: double operator()(const OPTION op) { if (op.option_type == OPTIONTYPE::EUROPEAN) return bs(op); if (op.call_put == CALLPUT::CALL op.b &>= op.r) return bs(op); else { double first_passage_price = FirstPassagePrice(op); if (op.call_put == CALLPUT::CALL) { if (op.underlying_price &< first_passage_price) return bs(op) + _A(op, first_passage_price) * pow(op.underlying_price / first_passage_price, _q(op)); else return op.underlying_price - op.strike_price; } else { if (op.underlying_price &> first_passage_price) return bs(op) + _A(op, first_passage_price) * pow(op.underlying_price / first_passage_price, _q(op)); else return op.strike_price - op.underlying_price; } } } private: BLACKSCHOLES bs; const double ACCURACY = 1.0e-6; const int MAXITERNUM = 1000; double FirstPassagePrice(const OPTION op) { double fpp = _seed(op); for (int i = 0; i &< MAXITERNUM; ++i) { fpp = FirstPassagePriceIterator(op, fpp); if (abs(_lhs(op, fpp) - _rhs(op, fpp)) / op.strike_price &< ACCURACY) return fpp; } return fpp; } double FirstPassagePriceIterator(const OPTION op, double p) { double temp = _b(op, p); if (op.call_put == CALLPUT::CALL) return (op.strike_price + _rhs(op, p) - temp * p) / (1 - temp); else return (op.strike_price - _rhs(op, p) + temp * p) / (1 + temp); } double _d(const OPTION op, double x) { return (log(x / op.strike_price) + (op.b + 0.5 * op.sigma * op.sigma) * op.t) / (op.sigma * sqrt(op.t)); } double _q(const OPTION op) { double _m = 2 * op.r / op.sigma / op.sigma, _n = 2 * op.b / op.sigma / op.sigma, _k = 1 - exp(-op.r * op.t); if (op.call_put == CALLPUT::CALL) return (-_n + 1 + sqrt((_n - 1) * (_n - 1) + 4 * _m / _k)) / 2; else return (-_n + 1 - sqrt((_n - 1) * (_n - 1) + 4 * _m / _k)) / 2; } double _q_inf(const OPTION op) { double _m = 2 * op.r / op.sigma / op.sigma, _n = 2 * op.b / op.sigma / op.sigma, _k = 1 - exp(-op.r * op.t); if (op.call_put == CALLPUT::CALL) return (-_n + 1 + sqrt((_n - 1) * (_n - 1) + 4 * _m)) / 2; else return (-_n + 1 - sqrt((_n - 1) * (_n - 1) + 4 * _m)) / 2; } double _lhs(const OPTION op, double x) { if (op.call_put == CALLPUT::CALL) return x - op.strike_price; else return op.strike_price - x; } double _rhs(const OPTION op, double x) { OPTION _op(op); _op.underlying_price = x; if (op.call_put == CALLPUT::CALL) return bs(_op) + (1 - exp((op.b - op.r) * op.t) * NormDist(_d(op, x))) * x / _q(op); else return bs(_op) - (1 - exp((op.b - op.r) * op.t) * NormDist(-_d(op, x))) * x / _q(op); } double _b(const OPTION op, double x) { double temp = 0; if (op.call_put == CALLPUT::CALL) temp = exp((op.b - op.r) * op.t) * NormDist(_d(op, x)) * (1 - 1 / _q(op)) + (1 - exp((op.b - op.r) * op.t) * NormDensity(_d(op, x)) / op.sigma / sqrt(op.t)) / _q(op); else temp = -exp((op.b - op.r) * op.t) * NormDist(_d(op, x)) * (1 - 1 / _q(op)) - (1 + exp((op.b - op.r) * op.t) * NormDensity(_d(op, x)) / op.sigma / sqrt(op.t)) / _q(op); return temp; } double _seed(const OPTION op) { double s_inf = op.strike_price / (1 - 1 / _q_inf(op)), h = 0; if (op.call_put == CALLPUT::CALL) { h = -(op.b * op.t + 2 * op.sigma * sqrt(op.t)) * (op.strike_price / (s_inf - op.strike_price)); return op.strike_price + (s_inf - op.strike_price) * (1 - exp(h)); } else { h = (op.b * op.t - 2 * op.sigma * sqrt(op.t)) * (op.strike_price / (op.strike_price - s_inf)); return s_inf + (op.strike_price - s_inf) * exp(h); } } double _A(const OPTION op, double x) { if (op.call_put == CALLPUT::CALL) return x / _q(op) * (1 - exp((op.b - op.r) * op.t) * NormDist(_d(op, x))); else return -x / _q(op) * (1 - exp((op.b - op.r) * op.t) * NormDist(-_d(op, x))); } }; } #endif // !LIUTING_OPTION_H