4.2.1 欧式期权定价

先编写StockOption类，存放需要的参数。

""" Store common attributes of a stock option """

import math

class StockOption(object):

def __init__(self, S0, K, r, T, N, params):

self.S0 = S0

self.K = K

self.r = r

self.T = T

self.N = max(1, N) #Ensure N have at least 1 time step

self.STs = None #Declare the stock prices tree

""" Optional parameters used by derived classes """

self.pu = params.get("pu", 0) #Probability of up state

self.pd= params.get("pd", 0) #Probability of down state

self.div = params.get("div", 0) #dividend yield

self.sigma = params.get("sigma", 0) #Volatility

self.is_call = params.get("is_call", True) #Call or put

self.is_european = params.get("is_european", True) # Eu or Am

""" Compute values"""

self.dt = T/float(N) #Single time step, in years

self.df = math.exp(-(r-self.div)*self.dt) #Discount factor

Python是一种面对对象的语言，所以最重要的就是“类”的使用；本文关于面对对象和面对过程的区别不再赘述。

编写BinomialEuropeanOption类

""" Price a European option by the binomial tree model"""

from StockOption import StockOption①

import math

import numpy as np

class BinomialEuropeanOption(StockOption):②

def _setup_parameters_(self):③

""" Required calculations for the model """

self.M = self.N + 1 # Number of terminal nodes of tree

self.u = 1 + self.pu # Expected value in the up state

self.d = - self.pd # Expected value in the down state

self.qu = (math.exp((self.r-self.div)*self.dt) -

self.d) / (self.u-self.d)

self.qd = 1-self.qu

def _initialize_stock_price_tree_(self):④

# Initialize terminal price nodes to zeros

self.STs = np.zeros(self.M)

# Calculate expected stock prices for each node

for i in range(self.M):

self.STs[i] = self.S0*(self.u**(self.N-i))*(self.d**i)

def _initialize_payoffs_tree_(self):⑤

# Get payoffs when the option expires at terminal nodes

payoffs = np.maximum(

0, (self.STs-self.K) if self.is_call

else (self.K-self.STs))

return payoffs

def _traverse_tree_(self, payoffs):⑥

# Starting from the time the option expires, traverse

# backwards and calculate discounted payoffs at each node

for i in range(self.N):

payoffs = (payoffs[:-1] * self.qu +

payoffs[1:] * self.qd) * self.df

return payoffs

def __begin_tree_traversal__(self):⑦

payoffs = self._initialize_payoffs_tree_()

return self._traverse_tree_(payoffs)

def price(self):⑧

""" The pricing implementation """

self._setup_parameters_()

self._initialize_stock_price_tree_()

payoffs = self.__begin_tree_traversal__()

return payoffs[0] # Option value converges to first node

①导入先前的StockOption类和其他模块

②对StockOption类进行继承

③计算期权定价模型所需的相关参数，包括上涨的概率和上涨的幅度等。四步二叉树

为明晰作者编程的思路，按照原代码进行推演，假设进行三步二叉树，N=3。

④三步二叉树的终点有四个节点，所以要使用N+1=4。代码用np.zeros(4)生成最终结点的零矩阵 [0, 0, 0, 0]。再用range(4)生成列表[0, 1, 2, 3]，u和d分别为上涨幅度和下跌幅度。通过for循环生成最终节点的价格分布 [S0*uuu, S0*uud, S0*udd, S0*udd]

⑤对看跌期权定价，在风险中性情况下，需要持有标的物的空头和看跌期权的空头。在某节点的期权内在价值为max(0, K-S)，通过for循环将标的物价格转换为期权价格，并返回payoffs

[max(0, K-S0*uuu), max(0, K-S0*uud), max(0, K-S0*udd), max(0, K-S0*ddd)]

⑥需要将payoffs列表作为参数，乘上对应的概率。通过range控制for循环次数为三次，每次对不同部分进行乘积。payoffs[:-1]是通过切片将payoffs中除了最后一个元素的其他值选取，payoffs[1:]是除了第一个元素进行选取，最后乘上折现因子。该折现步骤比较巧妙，最后payoffs变成一个数。

⑦对⑤和⑥进行类内引用，可一步求出最终的payoffs列表

⑧对③④⑦部分进行引用，所以只要执行⑧，其他属性也将被执行。

运行原书案例：

from BinomialEuropeanOption import BinomialEuropeanOption

eu_option = BinomialEuropeanOption(50, 50, 0.05, 0.5, 2, {"pu":0.2, "pd":0.2, "is_call":False})

print (eu_option.price())

>>>

4.82565175125599

4.2.4 美式期权定价

""" Price a European or American option by the binomial tree """

from StockOption import StockOption

import math

import numpy as np

class BinomialTreeOption(StockOption):

def _setup_parameters_(self):

self.u = 1 + self.pu # Expected value in the up state

self.d = 1 - self.pd # Expected value in the down state

self.qu = (math.exp((self.r-self.div)*self.dt) -

self.d)/(self.u-self.d)

self.qd = 1-self.qu

def _initialize_stock_price_tree_(self):

# Initialize a 2D tree at T=0

self.STs = [np.array([self.S0])]

# Simulate the possible stock prices path

for i in range(self.N):

prev_branches = self.STs[-1]

st = np.concatenate((prev_branches*self.u,

[prev_branches[-1]*self.d]))

self.STs.append(st) # Add nodes at each time step

def _initialize_payoffs_tree_(self):

# The payoffs when option expires

return np.maximum(

0, (self.STs[self.N]-self.K) if self.is_call

else (self.K-self.STs[self.N]))

def __check_early_exercise__(self, payoffs, node):

early_ex_payoff = \

(self.STs[node] - self.K) if self.is_call \

else (self.K - self.STs[node])

return np.maximum(payoffs, early_ex_payoff)

def _traverse_tree_(self, payoffs):

for i in reversed(range(self.N)):

# The payoffs from NOT exercising the option

payoffs = (payoffs[:-1] * self.qu +

payoffs[1:] * self.qd) * self.df

# Payoffs from exercising, for American options

if not self.is_european:

payoffs = self.__check_early_exercise__(payoffs, i)

return payoffs

def __begin_tree_traversal__(self):

payoffs = self._initialize_payoffs_tree_()

return self._traverse_tree_(payoffs)

def price(self):

self._setup_parameters_()

self._initialize_stock_price_tree_()

payoffs = self.__begin_tree_traversal__()

return payoffs[0]

相比于之前的欧式期权定价，美式期权中多了节点检查，美式期权定价的原代码也集成了欧式期权定价，运行则有：

from BinomialTreeOption import BinomialTreeOption

am_option = BinomialTreeOption(

50, 50, 0.05, 0.5, 2,

{"pu": 0.2, "pd": 0.2, "is_call": False, "is_european": False})

print(am_option.price())

>>>

5.113060082823486

4.2.5 Cox-Ross-Rubinstein模型

该模型将标的资产的回报和波动率相联系，直接继承 BinomialTreeOption，并重写属性_setup_parameters_，因为这两个模型只差一个公式，在面对对象的编程当中是不需要重新把整个代码都写一遍的。

对两种期权进行定价。

from BinomialTreeOption import BinomialTreeOption

import math

class BinomialCRROption(BinomialTreeOption):

def _setup_parameters_(self):

self.u = math.exp(self.sigma * math.sqrt(self.dt))

self.d = 1./self.u

self.qu = (math.exp((self.r-self.div)*self.dt) -

self.d)/(self.u-self.d)

self.qd = 1-self.qu

运行：

from BinomialCRROption import BinomialCRROption

eu_option = BinomialCRROption(

50, 50, 0.05, 0.5, 2,

{"sigma": 0.3, "is_call": False})

print("European put: " + str(eu_option.price()))

am_option = BinomialCRROption(

50, 50, 0.05, 0.5, 2,

{"sigma": 0.3, "is_call": False, "is_european": False})

print("American put: " +str(am_option.price()))

>>>

European put: 3.1051473412967003

American put: 3.4091814964048277

4.2.6 Leisen-Reimer模型

""" Price an option by the Leisen-Reimer tree """

from BinomialTreeOption import BinomialTreeOption

import math

class BinomialLROption(BinomialTreeOption):

def _setup_parameters_(self):

odd_N = self.N if (self.N%2 == 1) else (self.N+1)

d1 = (math.log(self.S0/self.K) +

((self.r-self.div) +

(self.sigma**2)/2.) *

self.T) / (self.sigma * math.sqrt(self.T))

d2 = (math.log(self.S0/self.K) +

((self.r-self.div) -

(self.sigma**2)/2.) *

self.T) / (self.sigma * math.sqrt(self.T))

pp_2_inversion = \

lambda z, n: \

.5 + math.copysign(1, z) * \

math.sqrt(.25 - .25 * math.exp(

-((z/(n+1./3.+.1/(n+1)))**2.)*(n+1./6.)))

pbar = pp_2_inversion(d1, odd_N)

self.p = pp_2_inversion(d2, odd_N)

self.u = 1/self.df * pbar/self.p

self.d = (1/self.df - self.p*self.u)/(1-self.p)

self.qu = self.p

self.qd = 1-self.p

运行：

from BinomialLROption import BinomialLROption

eu_option = BinomialLROption(

50, 50, 0.05, 0.5, 3,

{"sigma": 0.3, "is_call": False})

print("European put: " + str(eu_option.price()))

am_option = BinomialLROption(

50, 50, 0.05, 0.5, 3,

{"sigma": 0.3, "is_call": False, "is_european": False})

print("American put: " + str(am_option.price()))

>>>

European put: 3.5674299991832887

American put: 3.668179104133528

4.3 希腊值

用LR的二叉树模型求Delta和Gamma的值：

""" Compute option price, delta and gamma by the LR tree """

from BinomialLROption import BinomialLROption

import numpy as np

class BinomialLRWithGreeks(BinomialLROption):

def __new_stock_price_tree__(self):

"""Create additional layer of nodes to ouroriginal stock price tree"""

self.STs = [np.array([self.S0*self.u/self.d,

self.S0,

self.S0*self.d/self.u])]

for i in range(self.N):

prev_branches = self.STs[-1]

st = np.concatenate((prev_branches * self.u,

[prev_branches[-1] * self.d]))

self.STs.append(st)

def price(self):

self._setup_parameters_()

self.__new_stock_price_tree__()

payoffs = self.__begin_tree_traversal__()

""" Option value is now in the middle node at t=0"""

option_value = payoffs[int(len(payoffs)/2)]

payoff_up = payoffs[0]

payoff_down = payoffs[-1]

S_up = self.STs[0][0]

S_down = self.STs[0][-1]

dS_up = S_up - self.S0

dS_down = self.S0 - S_down

""" Get delta value """

dS = S_up - S_down

dV = payoff_up - payoff_down

delta = dV/dS

""" Get gamma value """

gamma = ((payoff_up-option_value)/dS_up -

(option_value-payoff_down)/dS_down) / \

((self.S0+S_up)/2. - (self.S0+S_down)/2.)

return option_value, delta, gamma

运行：

from BinomialLRWithGreeks import BinomialLRWithGreeks

eu_call = BinomialLRWithGreeks(

50, 50, 0.05, 0.5, 300, {"sigma": 0.3, "is_call": True})

results = eu_call.price()

print("European call values")

print("Price:%s\nDelta:%s\nGamma:%s" % results)

eu_put = BinomialLRWithGreeks(

50, 50, 0.05, 0.5, 300, {"sigma":0.3, "is_call": False,})

results = eu_put.price()

print("European put values")

print("Price:%s\nDelta:%s\nGamma:%s" % results)

>>>

European call values

Price: 4.803864657409284

Delta: 0.5888015221823649

Gamma: 0.036736782388357196

European put values

Price: 3.569360258827997

Delta: -0.41119847781760727

Gamma: 0.03673678238835338

4.4 Boyle三叉树

""" Price an option by the Boyle trinomial tree """

from BinomialTreeOption import BinomialTreeOption

import math

import numpy as np

class TrinomialTreeOption(BinomialTreeOption):

def _setup_parameters_(self):

""" Required calculations for the model """

self.u = math.exp(self.sigma*math.sqrt(2.*self.dt))

self.d = 1/self.u

self.m = 1

self.qu = ((math.exp((self.r-self.div) *

self.dt/2.) -

math.exp(-self.sigma *

math.sqrt(self.dt/2.))) /

(math.exp(self.sigma *

math.sqrt(self.dt/2.)) -

math.exp(-self.sigma *

math.sqrt(self.dt/2.))))**2

self.qd = ((math.exp(self.sigma *

math.sqrt(self.dt/2.)) -

math.exp((self.r-self.div) *

self.dt/2.)) /

(math.exp(self.sigma *

math.sqrt(self.dt/2.)) -

math.exp(-self.sigma *

math.sqrt(self.dt/2.))))**2.

self.qm = 1 - self.qu - self.qd

def _initialize_stock_price_tree_(self):

""" Initialize a 2D tree at t=0 """

self.STs = [np.array([self.S0])]

for i in range(self.N):

prev_nodes = self.STs[-1]

self.ST = np.concatenate(

(prev_nodes*self.u, [prev_nodes[-1]*self.m,

prev_nodes[-1]*self.d]))

self.STs.append(self.ST)

def _traverse_tree_(self, payoffs):

""" Traverse the tree backwards """

for i in reversed(range(self.N)):

payoffs = (payoffs[:-2] * self.qu +

payoffs[1:-1] * self.qm +

payoffs[2:] * self.qd) * self.df

if not self.is_european:

payoffs = self.__check_early_exercise__(payoffs,

i)

return payoffs

构造原理和二叉树完全一样，运行可得：

from TrinomialTreeOption import TrinomialTreeOption

print("European put:", TrinomialTreeOption(

50, 50, 0.05, 0.5, 2, {"sigma": 0.3, "is_call": False}).price())

print("American put:", TrinomialTreeOption(

50, 50, 0.05, 0.5, 2, {"sigma": 0.3, "is_call": False, "is_european": False}).price())

>>>

European put: 3.330905491759248

American put: 3.48241453902267

4.5.1 二叉树Lattice方法

这种方法在效率上会更高。

""" Price an option by the binomial CRR lattice """

from BinomialCRROption import BinomialCRROption

import numpy as np

class BinomialCRRLattice(BinomialCRROption):

def _setup_parameters_(self):

super(BinomialCRRLattice, self)._setup_parameters_()

self.M = 2*self.N + 1

def _initialize_stock_price_tree_(self):

self.STs = np.zeros(self.M)

self.STs[0] = self.S0 * self.u**self.N

for i in range(self.M)[1:]:

self.STs[i] = self.STs[i-1]*self.d

def _initialize_payoffs_tree_(self):

odd_nodes = self.STs[::2]

return np.maximum(

0, (odd_nodes - self.K) if self.is_call

else(self.K - odd_nodes))

def __check_early_exercise__(self, payoffs, node):

self.STs = self.STs[1:-1] # Shorten the ends of the list

odd_STs = self.STs[::2]

early_ex_payoffs = \

(odd_STs-self.K) if self.is_call \

else (self.K-odd_STs)

payoffs = np.maximum(payoffs, early_ex_payoffs)

return payoffs

运行：

from BinomialCRRLattice import BinomialCRRLattice

eu_option = BinomialCRRLattice(

50, 50, 0.05, 0.5, 2,

{"sigma": 0.3, "is_call": False})

print("European put:%s" % eu_option.price())

am_option = BinomialCRRLattice(

50, 50, 0.05, 0.5, 2,

{"sigma": 0.3, "is_call": False, "is_european": False})

print("American put:%s" % am_option.price())

>>>

European put: 3.1051473412967057

American put: 3.409181496404833

4.5.1 三叉树Lattice方法

""" Price an option by the trinomial lattice """

from TrinomialTreeOption import TrinomialTreeOption

import numpy as np

class TrinomialLattice(TrinomialTreeOption):

def _setup_parameters_(self):

super(TrinomialLattice, self)._setup_parameters_()

self.M = 2*self.N+1

def _initialize_stock_price_tree_(self):

self.STs = np.zeros(self.M)

self.STs[0] = self.S0 * self.u**self.N

for i in range(self.M)[1:]:

self.STs[i] = self.STs[i-1]*self.d

def _initialize_payoffs_tree_(self):

return np.maximum(

0, (self.STs-self.K) if self.is_call

else(self.K-self.STs))

def __check_early_exercise__(self, payoffs, node):

self.STs = self.STs[1:-1] # Shorten the ends of the list

early_ex_payoffs = \

(self.STs-self.K) if self.is_call \

else(self.K-self.STs)

payoffs = np.maximum(payoffs, early_ex_payoffs)

return payoffs

运行：

from TrinomialLattice import TrinomialLattice

eu_option = TrinomialLattice(

50, 50, 0.05, 0.5, 2,

{"sigma": 0.3, "is_call":False})

print("European put: " + str(eu_option.price()))

am_option = TrinomialLattice(

50, 50, 0.05, 0.5, 2,

{"sigma": 0.3, "is_call": False, "is_european": False})

print ("American put: " + str(am_option.price()))

>>>

European put: 3.330905491759248

American put: 3.48241453902267

4.6.1 显式有限差分法(PDE)

编写FiniteDifferences类，用于分配参数。

""" Shared attributes and functions of FD """

import numpy as np

class FiniteDifferences(object):

def __init__(self, S0, K, r, T, sigma, Smax, M, N,

is_call=True):

self.S0 = S0

self.K = K

self.r = r

self.T = T

self.sigma = sigma

self.Smax = Smax

self.M, self.N = int(M), int(N) # Ensure M&N are integers

self.is_call = is_call

self.dS = Smax / float(self.M)

self.dt = T / float(self.N)

self.i_values = np.arange(self.M)

self.j_values = np.arange(self.N)

self.grid = np.zeros(shape=(self.M+1, self.N+1))

self.boundary_conds = np.linspace(0, Smax, self.M+1)

def _setup_boundary_conditions_(self):

pass

def _setup_coefficients_(self):

pass

def _traverse_grid_(self):

""" Iterate the grid backwards in time """

pass

def _interpolate_(self):

"""Use piecewise linear interpolation on the initialgrid column to get the closest price at S0."""

return np.interp(self.S0,

self.boundary_conds,

self.grid[:, 0])

def price(self):

self._setup_boundary_conditions_()

self._setup_coefficients_()

self._traverse_grid_()

return self._interpolate_()

编写显式PDE类：

""" Explicit method of Finite Differences """

import numpy as np

from FiniteDifferences import FiniteDifferences

class FDExplicitEu(FiniteDifferences):

def _setup_boundary_conditions_(self):

if self.is_call:

self.grid[:, -1] = np.maximum(

self.boundary_conds - self.K, 0)

self.grid[-1, :-1] = (self.Smax - self.K) * \

np.exp(-self.r *

self.dt *

(self.N-self.j_values))

else:

self.grid[:, -1] = \

np.maximum(self.K-self.boundary_conds, 0)

self.grid[0, :-1] = (self.K - self.Smax) * \

np.exp(-self.r *

self.dt *

(self.N-self.j_values))

def _setup_coefficients_(self):

self.a = 0.5*self.dt*((self.sigma**2) *

(self.i_values**2) -

self.r*self.i_values)

self.b = 1 - self.dt*((self.sigma**2) *

(self.i_values**2) +

self.r)

self.c = 0.5*self.dt*((self.sigma**2) *

(self.i_values**2) +

self.r*self.i_values)

def _traverse_grid_(self):

for j in reversed(self.j_values):

for i in range(self.M)[2:]:

self.grid[i,j] = self.a[i]*self.grid[i-1,j+1] +\

self.b[i]*self.grid[i,j+1] + \

self.c[i]*self.grid[i+1,j+1]

运行：

from FDExplicitEu import FDExplicitEu

option = FDExplicitEu(50, 50, 0.1, 5./12., 0.4, 100, 100, 1000, False)

print("European_put: ", option.price())

option = FDExplicitEu(50, 50, 0.1, 5./12., 0.4, 100, 100, 100, False)

print("European_put: ", option.price())

>>>

European_put: 4.072882278148043

European_put: -1.6291077072251005e+53

如果参数选择不正确会导致显式方法的不稳定。

4.6.2 隐式有限差分法(PDE)

"""Price a European option by the implicit methodof finite differences."""

import numpy as np

import scipy.linalg as linalg

from FDExplicitEu import FDExplicitEu

class FDImplicitEu(FDExplicitEu):

def _setup_coefficients_(self):

self.a = 0.5*(self.r*self.dt*self.i_values -

(self.sigma**2)*self.dt*(self.i_values**2))

self.b = 1 + \

(self.sigma**2)*self.dt*(self.i_values**2) + \

self.r*self.dt

self.c = -0.5*(self.r * self.dt*self.i_values +

(self.sigma**2)*self.dt*(self.i_values**2))

self.coeffs = np.diag(self.a[2:self.M], -1) + \

np.diag(self.b[1:self.M]) + \

np.diag(self.c[1:self.M-1], 1)

def _traverse_grid_(self):

""" Solve using linear systems of equations """

P, L, U = linalg.lu(self.coeffs)

aux = np.zeros(self.M-1)

for j in reversed(range(self.N)):

aux[0] = np.dot(-self.a[1], self.grid[0, j])

x1 = linalg.solve(L, self.grid[1:self.M, j+1]+aux)

x2 = linalg.solve(U, x1)

self.grid[1:self.M, j] = x2

运行：

from FDImplicitEu import FDImplicitEu

option = FDImplicitEu(50, 50, 0.1, 5./12., 0.4, 100, 100, 100, False)

print(option.price())

option = FDImplicitEu(50, 50, 0.1, 5./12., 0.4, 100, 100, 1000, False)

print(option.price())

>>>

4.065801939431454

4.071594188049893

4.6.3 Crank-Nicolson方法

""" Crank-Nicolson method of Finite Differences """

import numpy as np

import scipy.linalg as linalg

from FDExplicitEu import FDExplicitEu

class FDCnEu(FDExplicitEu):

def _setup_coefficients_(self):

self.alpha = 0.25*self.dt*(

(self.sigma**2)*(self.i_values**2) -

self.r*self.i_values)

self.beta = -self.dt*0.5*(

(self.sigma**2)*(self.i_values**2) +

self.r)

self.gamma = 0.25*self.dt*(

(self.sigma**2)*(self.i_values**2) +

self.r*self.i_values)

self.M1 = -np.diag(self.alpha[2:self.M], -1) + \

np.diag(1-self.beta[1:self.M]) - \

np.diag(self.gamma[1:self.M-1], 1)

self.M2 = np.diag(self.alpha[2:self.M], -1) + \

np.diag(1+self.beta[1:self.M]) + \

np.diag(self.gamma[1:self.M-1], 1)

def _traverse_grid_(self):

""" Solve using linear systems of equations """

P, L, U = linalg.lu(self.M1)

for j in reversed(range(self.N)):

x1 = linalg.solve(L,

np.dot(self.M2,

self.grid[1:self.M, j+1]))

x2 = linalg.solve(U, x1)

self.grid[1:self.M, j] = x2

运行：

from FDCnEu import FDCnEu

option = FDCnEu(50, 50, 0.1, 5./12., 0.4, 100, 100, 100, False)

print(option.price())

option = FDCnEu(50, 50, 0.1, 5./12., 0.4, 100, 100, 1000, False)

print(option.price())

>>>

4.072254507998114

4.072238354486825

4.6.4 障碍期权定价之下降出局期权

"""Price a down-and-out option by the Crank-Nicolsonmethod of finite differences."""

import numpy as np

from FDCnEu import FDCnEu

class FDCnDo(FDCnEu):

def __init__(self, S0, K, r, T, sigma, Sbarrier, Smax, M, N,

is_call=True):

super(FDCnDo, self).__init__(

S0, K, r, T, sigma, Smax, M, N, is_call)

self.dS = (Smax-Sbarrier)/float(self.M)

self.boundary_conds = np.linspace(Sbarrier,

Smax,

self.M+1)

self.i_values = self.boundary_conds/self.dS

运行：

from FDCnDo import FDCnDo

option = FDCnDo(50, 50, 0.1, 5./12., 0.4, 40, 100, 120, 500)

print(option.price())

option = FDCnDo(50, 50, 0.1, 5./12., 0.4, 40, 100, 120, 500, False)

print(option.price())

>>>

5.491560552934787

0.5413635028954449

5.6.5 美式期权定价的有限差分

""" Price an American option by the Crank-Nicolson method """

import numpy as np

import sys

from FDCnEu import FDCnEu

class FDCnAm(FDCnEu):

def __init__(self, S0, K, r, T, sigma, Smax, M, N, omega, tol,

is_call=True):

super(FDCnAm, self).__init__(

S0, K, r, T, sigma, Smax, M, N, is_call)

self.omega = omega

self.tol = tol

self.i_values = np.arange(self.M+1)

self.j_values = np.arange(self.N+1)

def _setup_boundary_conditions_(self):

if self.is_call:

self.payoffs = np.maximum(

self.boundary_conds[1:self.M]-self.K, 0)

else:

self.payoffs = np.maximum(

self.K-self.boundary_conds[1:self.M], 0)

self.past_values = self.payoffs

self.boundary_values = self.K * \

np.exp(-self.r *

self.dt *

(self.N-self.j_values))

def _traverse_grid_(self):

""" Solve using linear systems of equations """

aux = np.zeros(self.M-1)

new_values = np.zeros(self.M-1)

for j in reversed(range(self.N)):

aux[0] = self.alpha[1]*(self.boundary_values[j] +

self.boundary_values[j+1])

rhs = np.dot(self.M2, self.past_values) + aux

old_values = np.copy(self.past_values)

error = sys.float_info.max

while self.tol < error:

new_values[0] = \

max(self.payoffs[0],

old_values[0] +

self.omega/(1-self.beta[1]) *

(rhs[0] -

(1-self.beta[1])*old_values[0] +

(self.gamma[1]*old_values[1])))

for k in range(self.M-2)[1:]:

new_values[k] = \

max(self.payoffs[k],

old_values[k] +

self.omega/(1-self.beta[k+1]) *

(rhs[k] +

self.alpha[k+1]*new_values[k-1] -

(1-self.beta[k+1])*old_values[k] +

self.gamma[k+1]*old_values[k+1]))

new_values[-1] = \

max(self.payoffs[-1],

old_values[-1] +

self.omega/(1-self.beta[-2]) *

(rhs[-1] +

self.alpha[-2]*new_values[-2] -

(1-self.beta[-2])*old_values[-1]))

error = np.linalg.norm(new_values - old_values)

old_values = np.copy(new_values)

self.past_values = np.copy(new_values)

self.values = np.concatenate(([self.boundary_values[0]],

new_values,

[0]))

def _interpolate_(self):

# Use linear interpolation on final values as 1D array

return np.interp(self.S0,

self.boundary_conds,

self.values)

运行：

from FDCnAm import FDCnAm

option = FDCnAm(50, 50, 0.1, 5./12., 0.4, 100, 100, 42, 1.2, 0.001)

print(option.price())

option = FDCnAm(50, 50, 0.1, 5./12., 0.4, 100, 100, 42, 1.2, 0.001, False)

print(option.price())

>>>6.108682815392218

4.2777642293837355

4.7 通过二叉树计算隐含波动率

先编写二分法：

def bisection(f, a, b, tol=0.1, maxiter=10):

""":param f: The function to solve:param a: The x-axis value where f(a)<0:param b: The x-axis value where f(b)>0:param tol: The precision of the solution:param maxiter: Maximum number of iterations:return: The x-axis value of the root,number of iterations used"""

c = (a+b)*0.5 # Declare c as the midpoint ab

n = 1 # Start with 1 iteration

while n <= maxiter:

c = (a+b)*0.5

if f(c) == 0 or abs(a-b)*0.5 < tol:

# Root is found or is very close

return c, n

n += 1

if f(c) < 0:

a = c

else:

b = c

return c, n

再通过LR模型求解隐含波动率：

from BinomialLROption import BinomialLROption

from bisection import bisection

"""Get implied volatilities from a Leisen-Reimer binomialtree using the bisection method as the numerical procedure."""

class ImpliedVolatilityModel(object):

def __init__(self, S0, r=0.05, T=1, div=0,

N=1, is_put=False):

self.S0 = S0

self.r = r

self.T = T

self.div = div

self.N = N

self.is_put = is_put

def option_valuation(self, K, sigma):

""" Use the binomial Leisen-Reimer tree """

lr_option = BinomialLROption(

self.S0, K, r=self.r, T=self.T, N=self.N,

sigma=sigma, div=self.div, is_put=self.is_put

)

return lr_option.price()

def get_implied_volatilities(self, Ks, opt_prices):

impvols = []

for i in range(len(Ks)):

# Bind f(sigma) for use by the bisection method

f = lambda sigma: \

self.option_valuation(Ks[i], sigma)-\

opt_prices[i]

impv = bisection(f, 0.01, 0.99, 0.0001, 100)[0]

impvols.append(impv)

return impvols

运行：

from ImpliedVolatilityModel import ImpliedVolatilityModel

strikes = [75, 80, 85, 90, 92.5, 95, 97.5,

100, 105, 110, 115, 120, 125]

put_prices = [0.16, 0.32, 0.6, 1.22, 1.77, 2.54, 3.55,

4.8, 7.75, 11.8, 15.96, 20.75, 25.81]

model = ImpliedVolatilityModel(

99.62, r=0.0248, T=78/365., div=0.0182, N=77, is_put=True)

impvols_put = model.get_implied_volatilities(strikes, put_prices)

import matplotlib.pyplot as plt

plt.plot(strikes, impvols_put)

plt.xlabel('Strike Prices')

plt.ylabel('Implied Volatilities')

plt.title('AAPL Put Implied Volatilities expiring in 78 days')

plt.show()波动率微笑