python barrier option pricing_《Python金融数据分析》书内代码实战与讲解（二）金融衍生物定价...
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()波動率微笑
總結
以上是生活随笔為你收集整理的python barrier option pricing_《Python金融数据分析》书内代码实战与讲解（二）金融衍生物定价...的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： java guava 使用_Java8-
- 下一篇： java自带工具_深入了解Java JD