精通 Python 金融第二版（二）

原文：zh.annas-archive.org/md5/8b046e39ce2c1a10ac13fd89834aaadc

译者：飞龙

协议：CC BY-NC-SA 4.0

第四章：期权定价的数值方法

衍生品是一种合同，其回报取决于某些基础资产的价值。在封闭形式衍生品定价可能复杂甚至不可能的情况下，数值程序表现出色。数值程序是使用迭代计算方法试图收敛到解的方法。其中一个基本实现是二项树。在二项树中，一个节点代表某一时间点的资产状态，与价格相关。每个节点在下一个时间步骤导致另外两个节点。同样，在三项树中，每个节点在下一个时间步骤导致另外三个节点。然而，随着树的节点数量或时间步骤的增加，消耗的计算资源也会增加。栅格定价试图通过仅在每个时间步骤存储新信息，同时在可能的情况下重复使用价值来解决这个问题。

在有限差分定价中，树的节点也可以表示为网格。网格上的终端值包括终端条件，而网格的边缘代表资产定价中的边界条件。我们将讨论有限差分方案的显式方法、隐式方法和 Crank-Nicolson 方法，以确定资产的价格。

尽管香草期权和某些奇异期权，如欧式障碍期权和回望期权，可能具有封闭形式解，但其他奇异产品，如亚洲期权，没有封闭形式解。在这些情况下，可以使用数值程序来定价期权。

在本章中，我们将涵盖以下主题：

• 使用二项树定价欧式和美式期权

• 使用 Cox-Ross-Rubinstein 二项树

• 使用 Leisen-Reimer 树定价期权

• 使用三项树定价期权

• 从树中推导希腊字母

• 使用二项和三项栅格定价期权

• 使用显式、隐式和 Crank-Nicolson 方法的有限差分

• 使用 LR 树和二分法的隐含波动率建模

期权介绍

期权是资产的衍生品，它赋予所有者在特定日期以特定价格交易基础资产的权利，称为到期日和行权价格。

认购期权赋予买方在特定日期以特定价格购买资产的权利。认购期权的卖方或写方有义务在约定日期以约定价格向买方出售基础证券，如果买方行使其权利。

认沽期权赋予买方在特定日期以特定价格出售基础资产的权利。认沽期权的卖方或写方有义务在约定日期以约定价格从买方购买基础证券，如果买方行使其权利。

最常见的期权包括欧式期权和美式期权。其他奇异期权包括百慕大期权和亚洲期权。本章主要涉及欧式期权和美式期权。欧式期权只能在到期日行使。而美式期权则可以在期权的整个生命周期内的任何时间行使。

期权定价中的二项树

在二项期权定价模型中，假设在一个时间段内，代表具有给定价格的节点的基础证券在下一个时间步骤中遍历到另外两个节点，代表上升状态和下降状态。由于期权是基础资产的衍生品，二项定价模型以离散时间为基础跟踪基础条件。二项期权定价可用于估值欧式期权、美式期权以及百慕大期权。

根节点的初始值是基础证券的现货价格S[0]，具有风险中性概率的上涨q和下跌的风险中性概率1-q，在下一个时间步骤。基于这些概率，为每个价格上涨或下跌的状态计算了证券的预期值。终端节点代表了每个预期证券价格的值，对应于上涨状态和下跌状态的每种组合。然后我们可以计算每个节点的期权价值，通过风险中性期望遍历树状图，并在从远期利率贴现后，我们可以推导出看涨期权或看跌期权的价值。

定价欧式期权

考虑一个两步二叉树。不支付股息的股票价格从 50 美元开始，在两个时间步骤中，股票可能上涨 20%，也可能下跌 20%。假设无风险利率为每年 5%，到期时间T为两年。我们想要找到行权价K为 52 美元的欧式看跌期权的价值。以下图表显示了使用二叉树定价股票和终端节点的回报：

这里，节点的计算如下：

在终端节点，行使欧式看涨期权的回报如下：

在欧式看跌期权的情况下，回报如下：

欧式看涨期权和看跌期权通常用小写字母c和p表示，而美式看涨期权和看跌期权通常用大写字母C和P表示。

从期权回报值中，我们可以向后遍历二叉树到当前时间，并在从无风险利率贴现后，我们将获得期权的现值。向后遍历树状图考虑了期权上涨状态和下跌状态的风险中性概率。

我们可以假设投资者对风险不感兴趣，并且所有资产的预期回报相等。在通过风险中性概率投资股票的情况下，持有股票的回报并考虑上涨和下跌状态的可能性将等于在下一个时间步骤中预期的连续复利无风险利率，如下所示：

投资股票的风险中性概率q可以重写如下：

这些公式与股票相关吗？期货呢？

与投资股票不同，投资者无需提前付款来持有期货合约。在风险中性意义上，持有期货合约的预期增长率为零，投资期货的风险中性概率q可以重写如下：

让我们计算前面示例中给出的股票的风险中性概率q：

In [ ]:
    import math

    r = 0.05
    T = 2
    t = T/2
    u = 1.2
    d = 0.8

    q = (math.exp(r*t)-d)/(u-d)
In [ ]:
    print('q is', q)
Out[ ]:   
    q is 0.6281777409400603

在终端节点行使欧式看跌期权的回报分别为 0 美元、4 美元和 20 美元。看跌期权的现值可以定价如下：

这给出了看跌期权价格为 4.19 美元。使用二叉树估算每个节点的欧式看跌期权的价值如下图所示：

编写 StockOption 基类

在进一步实现我们将要讨论的各种定价模型之前，让我们创建一个StockOption类，用于存储和计算股票期权的共同属性，这些属性将在本章中被重复使用：

In [ ]:
    import math

    """ 
    Stores common attributes of a stock option 
    """
    class StockOption(object):
        def __init__(
            self, S0, K, r=0.05, T=1, N=2, pu=0, pd=0, 
            div=0, sigma=0, is_put=False, is_am=False):
            """
            Initialize the stock option base class.
            Defaults to European call unless specified.

            :param S0: initial stock price
            :param K: strike price
            :param r: risk-free interest rate
            :param T: time to maturity
            :param N: number of time steps
            :param pu: probability at up state
            :param pd: probability at down state
            :param div: Dividend yield
            :param is_put: True for a put option,
                    False for a call option
            :param is_am: True for an American option,
                    False for a European option
            """
            self.S0 = S0
            self.K = K
            self.r = r
            self.T = T
            self.N = max(1, N)
            self.STs = [] # Declare the stock prices tree

            """ Optional parameters used by derived classes """
            self.pu, self.pd = pu, pd
            self.div = div
            self.sigma = sigma
            self.is_call = not is_put
            self.is_european = not is_am

        @property
        def dt(self):
            """ Single time step, in years """
            return self.T/float(self.N)

        @property
        def df(self):
            """ The discount factor """
            return math.exp(-(self.r-self.div)*self.dt)  

当前的标的价格、行权价格、无风险利率、到期时间和时间步数是定价期权的强制共同属性。时间步长dt和折现因子df的增量作为类的属性计算，并且如果需要，可以被实现类覆盖。

使用二项式树的欧式期权类

欧式期权的 Python 实现是BinomialEuropeanOption类，它继承自StockOption类的共同属性。该类中方法的实现如下：

1. BinomialEuropeanOption类的price()方法是该类所有实例的入口

2. 它调用setup_parameters()方法来设置所需的模型参数，然后调用init_stock_price_tree()方法来模拟期间内股票价格的预期值直到T

3. 最后，调用begin_tree_traversal()方法来初始化支付数组并存储折现支付值，因为它遍历二项式树回到现在的时间

4. 支付树节点作为 NumPy 数组对象返回，其中欧式期权的现值在初始节点处找到

BinomialEuropeanOption的类实现如下 Python 代码：

In [ ]:
    import math
    import numpy as np
    from decimal import Decimal

    """ 
    Price a European option by the binomial tree model 
    """
    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 = 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 init_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 init_payoffs_tree(self):
            """
            Returns the payoffs when the option 
            expires at terminal nodes
            """ 
            if self.is_call:
                return np.maximum(0, self.STs-self.K)
            else:
                return np.maximum(0, self.K-self.STs)

        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.init_payoffs_tree()
            return self.traverse_tree(payoffs)

        def price(self):
            """ Entry point of the pricing implementation """
            self.setup_parameters()
            self.init_stock_price_tree()
            payoffs = self.begin_tree_traversal()

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

让我们使用我们之前讨论的两步二项式树示例中的值来定价欧式看跌期权：

In [ ]:
    eu_option = BinomialEuropeanOption(
        50, 52, r=0.05, T=2, N=2, pu=0.2, pd=0.2, is_put=True)
In [ ]:
    print('European put option price is:', eu_option.price())
Out[ ]:    
    European put option price is: 4.1926542806038585

使用二项式期权定价模型，我们得到了欧式看跌期权的现值为 4.19 美元。

使用二项式树的美式期权类

与只能在到期时行使的欧式期权不同，美式期权可以在其寿命内的任何时候行使。

为了在 Python 中实现美式期权的定价，我们可以像BinomialEuropeanOption类一样创建一个名为BinomialTreeOption的类，该类继承自Stockoption类。setup_parameters()方法中使用的参数保持不变，只是删除了一个未使用的M参数。

美式期权中使用的方法如下：

• init_stock_price_tree：使用二维 NumPy 数组存储所有时间步的股票价格的预期回报。这些信息用于计算在每个期间行使期权时的支付值。该方法编写如下：

def init_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

• init_payoffs_tree：创建支付树作为二维 NumPy 数组，从期满时期的期权内在价值开始。该方法编写如下：

def init_payoffs_tree(self):
    if self.is_call:
        return np.maximum(0, self.STs[self.N]-self.K)
    else:
        return np.maximum(0, self.K-self.STs[self.N])

• check_early_exercise：返回在提前行使美式期权和根本不行使期权之间的最大支付值。该方法编写如下：

def check_early_exercise(self, payoffs, node):
    if self.is_call:
        return np.maximum(payoffs, self.STs[node] - self.K)
    else:
        return np.maximum(payoffs, self.K - self.STs[node])

• traverse_tree：这还包括调用check_early_exercise()方法，以检查是否在每个时间步提前行使美式期权是最优的。该方法编写如下：

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

begin_tree_traversal()和price()方法的实现保持不变。

当在类的实例化期间将is_put关键字参数设置为False或True时，BinomialTreeOption类可以定价欧式和美式期权。

以下代码是用于定价美式期权的：

In [ ]:
    am_option = BinomialTreeOption(50, 52, 
        r=0.05, T=2, N=2, pu=0.2, pd=0.2, is_put=True, is_am=True)
In [ ]:
    print('American put option price is:', am_option.price())
Out[ ]:    
    American put option price is: 5.089632474198373

美式看跌期权的价格为 5.0896 美元。由于美式期权可以在任何时候行使，而欧式期权只能在到期时行使，因此美式期权的这种灵活性在某些情况下增加了其价值。

对于不支付股息的基础资产的美式看涨期权，可能没有超过其欧式看涨期权对应的额外价值。由于时间价值的原因，今天在行权价上行使美式看涨期权的成本比在未来以相同行权价行使更高。对于实值的美式看涨期权，提前行使期权会失去对抗行权价以下不利价格波动的保护，以及其内在时间价值。没有股息支付的权利，没有动机提前行使美式看涨期权。

Cox–Ross–Rubinstein 模型

在前面的例子中，我们假设基础股价在相应的u上升状态和d下降状态分别增加 20％和减少 20％。Cox-Ross-Rubinstein（CRR）模型提出，在风险中性世界的短时间内，二项式模型与基础股票的均值和方差相匹配。基础股票的波动性，或者股票回报的标准差，如下所示：

CRR 二项式树期权定价模型的类

二项式 CRR 模型的实现与我们之前讨论的二项式树相同，唯一的区别在于u和d模型参数。

在 Python 中，让我们创建一个名为BinomialCRROption的类，并简单地继承BinomialTreeOption类。然后，我们只需要覆盖setup_parameters()方法，使用 CRR 模型中的值。

BinomialCRROption对象的实例将调用price()方法，该方法调用父类BinomialTreeOption的所有其他方法，除了被覆盖的setup_parameters()方法：

In [ ]:
    import math

    """ 
    Price an option by the binomial CRR model 
    """
    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

再次考虑两步二项式树。不支付股息的股票当前价格为 50 美元，波动率为 30％。假设无风险利率为年利率 5％，到期时间T为两年。我们想要找到 CRR 模型下的行权价为 52 美元的欧式看跌期权的价值：

In [ ]:
    eu_option = BinomialCRROption(
        50, 52, r=0.05, T=2, N=2, sigma=0.3, is_put=True)
In [ ]:
    print('European put:', eu_option.price())
Out[ ]:
    European put: 6.245708445206436
In [ ]:
    am_option = BinomialCRROption(50, 52, 
        r=0.05, T=2, N=2, sigma=0.3, is_put=True, is_am=True)
In [ ]:
    print('American put option price is:', am_option.price())
Out[ ]:
    American put option price is: 7.428401902704834

通过使用 CRR 两步二项式树模型，欧式看跌期权和美式看跌期权的价格分别为 6.2457 美元和 7.4284 美元。

使用 Leisen-Reimer 树

在我们之前讨论的二项式模型中，我们对上升和下降状态的概率以及由此产生的风险中性概率做出了几个假设。除了我们讨论的具有 CRR 参数的二项式模型之外，在数学金融中广泛讨论的其他形式的参数化包括 Jarrow-Rudd 参数化、Tian 参数化和 Leisen-Reimer 参数化。让我们详细看看 Leisen-Reimer 模型。

Dietmar Leisen 博士和 Matthias Reimer 提出了一个二项式树模型，旨在在步数增加时逼近 Black-Scholes 解。它被称为Leisen-Reimer（LR）树，节点不会在每个交替步骤重新组合。它使用反演公式在树遍历期间实现更好的准确性。

有关该公式的详细解释可在 1995 年 3 月的论文Binomial Models For Option Valuation - Examining And Improving Convergence中找到，网址为papers.ssrn.com/sol3/papers.cfm?abstract_id=5976。我们将使用 Peizer 和 Pratt 反演函数f的第二种方法，具有以下特征参数：

*S[0]*参数是当前股票价格，K是期权的行权价格，σ是基础股票的年化波动率，T是期权的到期时间，r是年化无风险利率，y是股息收益，Δt是每个树步之间的时间间隔。

LR 二项树期权定价模型的一个类

LR 树的 Python 实现在以下BinomialLROption类中给出。与BinomialCRROption类类似，我们只需继承BinomialTreeOption类，并用 LR 树模型的变量覆盖setup_parameters方法中的变量：

In [ ]:
    import math

    """ 
    Price an option by the Leisen-Reimer tree
    """
    class BinomialLROption(BinomialTreeOption):

        def setup_parameters(self):
            odd_N = self.N if (self.N%2 == 0) 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))

            pbar = self.pp_2_inversion(d1, odd_N)
            self.p = self.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

        def pp_2_inversion(self, z, n):
            return .5 + math.copysign(1, z)*\
                math.sqrt(.25 - .25*
                    math.exp(
                        -((z/(n+1./3.+.1/(n+1)))**2.)*(n+1./6.)
                    )
                )

使用我们之前使用的相同示例，我们可以使用 LR 树定价期权：

In [ ]:
    eu_option = BinomialLROption(
        50, 52, r=0.05, T=2, N=4, sigma=0.3, is_put=True)
In [ ]:
    print('European put:', eu_option.price())
Out[ ]:      
    European put: 5.878650106601964
In [ ]:
    am_option = BinomialLROption(50, 52, 
        r=0.05, T=2, N=4, sigma=0.3, is_put=True, is_am=True)
In [ ]:
    print('American put:', am_option.price())
Out[ ]:
    American put: 6.763641952939979

通过使用具有四个时间步长的 LR 二项树模型，欧式看跌期权的价格和美式看跌期权的价格分别为$5.87865 和$6.7636。

希腊字母免费

在我们迄今为止涵盖的二项树定价模型中，我们在每个时间点上上下遍历树来确定节点值。根据每个节点的信息，我们可以轻松地重用这些计算出的值。其中一种用途是计算希腊字母。

希腊字母衡量衍生品价格对基础资产参数变化的敏感性，例如期权，通常用希腊字母表示。在数学金融中，与希腊字母相关的常见名称包括 alpha、beta、delta、gamma、vega、theta 和 rho。

期权的两个特别有用的希腊字母是 delta 和 gamma。Delta 衡量期权价格对基础资产价格的敏感性。Gamma 衡量 delta 相对于基础价格的变化率。

如下图所示，在我们原始的两步树周围添加了额外的节点层，使其成为一个四步树，向时间向后延伸了两步。即使有额外的期末支付节点，所有节点将包含与我们原始的两步树相同的信息。我们感兴趣的期权价值现在位于树的中间，即t=0：

注意，在t=0处存在两个额外节点的信息，我们可以使用它们来计算 delta 公式，如下所示：

三角洲公式规定，期权价格在上涨和下跌状态之间的差异表示为时间t=0时各自股票价格之间的差异的单位。

相反，gamma 公式可以计算如下：

伽玛公式规定，上节点和下节点中期权价格的 delta 之间的差异与初始节点值相对于各自状态下股票价格的差异的单位进行计算。

LR 二项树的希腊字母类

为了说明在 LR 树中计算希腊字母的过程，让我们创建一个名为BinomialLRWithGreeks的新类，该类继承了BinomialLROption类，并使用我们自己的price方法的实现。

在price方法中，我们将首先调用父类的setup_parameters()方法来初始化 LR 树所需的所有变量。然而，这一次，我们还将调用new_stock_price_tree()方法，这是一个用于在原始树周围创建额外节点层的新方法。

调用begin_tree_traversal()方法执行父类中的通常 LR 树实现。返回的 NumPy 数组对象现在包含t=0处三个节点的信息，其中中间节点是期权价格。在数组的第一个和最后一个索引处是t=0处上升和下降状态的支付。

有了这些信息，price()方法计算并返回期权价格、delta 和 gamma 值：

In [ ]:
    import numpy as np

    """ 
    Compute option price, delta and gamma by the LR tree 
    """
    class BinomialLRWithGreeks(BinomialLROption):

        def new_stock_price_tree(self):
            """
            Creates an additional layer of nodes to our
            original 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[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

            # Calculate delta value
            dS = S_up - S_down
            dV = payoff_up - payoff_down
            delta = dV/dS

            # calculate 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

使用 LR 树的相同示例，我们可以计算具有 300 个时间步的欧式看涨期权和看跌期权的期权价值和希腊值：

In [ ]:
    eu_call = BinomialLRWithGreeks(50, 52, r=0.05, T=2, N=300, sigma=0.3)
    results = eu_call.price()
In [ ]:
    print('European call values')
    print('Price: %s\nDelta: %s\nGamma: %s' % results)
Out[ ]:
    European call values
    Price: 9.69546807138366
    Delta: 0.6392477816643529
    Gamma: 0.01764795890533088

In [ ]:
    eu_put = BinomialLRWithGreeks(
        50, 52, r=0.05, T=2, N=300, sigma=0.3, is_put=True)
    results = eu_put.price()
In [ ]:
    print('European put values')
    print('Price: %s\nDelta: %s\nGamma: %s' % results)
Out[ ]:   
    European put values
    Price: 6.747013809252746
    Delta: -0.3607522183356649
    Gamma: 0.0176479589053312

从price()方法和结果中可以看出，我们成功地从修改后的二项树中获得了希腊附加信息，而没有增加计算复杂性。

期权定价中的三项树

在二项树中，每个节点导致下一个时间步中的两个其他节点。同样，在三项树中，每个节点导致下一个时间步中的三个其他节点。除了具有上升和下降状态外，三项树的中间节点表示状态不变。当扩展到两个以上的时间步时，三项树可以被视为重新组合树，其中中间节点始终保留与上一个时间步相同的值。

让我们考虑 Boyle 三项树，其中树被校准，使得上升、下降和平稳移动的概率u、d和m与风险中性概率q[u]、*q[d]和q[m]*如下：

我们可以看到  重新组合为 m =1。通过校准，无状态移动 m 以 1 的固定利率增长，而不是以无风险利率增长。变量 v 是年化股息收益，σ 是基础股票的年化波动率。

一般来说，处理更多节点时，三项树在建模较少时间步时比二项树具有更好的精度，可以节省计算速度和资源。下图说明了具有两个时间步的三项树的股价变动：

三项树期权定价模型的类

让我们创建一个TrinomialTreeOption类，继承自BinomialTreeOption类。

TrinomialTreeOption的方法如下所示：

• setup_parameters()方法实现了三项树的模型参数。该方法编写如下：

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

• init_stock_price_tree()方法设置了三项树，包括股价的平稳移动。该方法编写如下：

def init_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)

• traverse_tree()方法在打折后考虑中间节点的收益：

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

• 使用二项树的相同示例，我们得到以下结果：

In [ ]:
   eu_put = TrinomialTreeOption(
        50, 52, r=0.05, T=2, N=2, sigma=0.3, is_put=True)
In [ ]:
   print('European put:', eu_put.price())
Out[ ]:
   European put: 6.573565269142496
In [ ]:
   am_option = TrinomialTreeOption(50, 52, 
        r=0.05, T=2, N=2, sigma=0.3, is_put=True, is_am=True)
In [ ]:
   print('American put:', am_option.price())
Out[ ]:
   American put: 7.161349217272585

通过三项树模型，我们得到了欧式看跌期权和美式看跌期权的价格分别为$6.57 和$7.16。

期权定价中的栅格

在二项树中，每个节点在每个交替节点处重新组合。在三项树中，每个节点在每个其他节点处重新组合。重新组合树的这种属性也可以表示为栅格，以节省内存而无需重新计算和存储重新组合的节点。

使用二项栅格

我们将从二项 CRR 树创建一个二项栅格，因为在每个交替的上升和下降节点处，价格重新组合为相同的ud=1概率。在下图中，S[u]和S[d]与S[du] = S[ud] = *S[0]***重新组合。现在可以将树表示为单个列表：

对于N步二项式树，需要一个大小为2N +1的列表来包含关于基础股票价格的信息。对于欧式期权定价，列表的奇数节点代表到期时的期权价值。树向后遍历以获得期权价值。对于美式期权定价，随着树向后遍历，列表的两端收缩，奇数节点代表任何时间步的相关股票价格。然后可以考虑早期行权的回报。

CRR 二项式栅格期权定价模型的类

让我们通过 CRR 将二项式树定价转换为栅格。我们可以继承BinomialCRROption类（该类又继承自BinomialTreeOption类），并创建一个名为BinomialCRRLattice的新类，如下所示：

In [ ]:
    import numpy as np

    class BinomialCRRLattice(BinomialCRROption):

        def setup_parameters(self):
            super(BinomialCRRLattice, self).setup_parameters()
            self.M = 2*self.N + 1

        def init_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 init_payoffs_tree(self):
            odd_nodes = self.STs[::2]  # Take odd nodes only
            if self.is_call:
                return np.maximum(0, odd_nodes-self.K)
            else:
                return np.maximum(0, self.K-odd_nodes)

        def check_early_exercise(self, payoffs, node):
            self.STs = self.STs[1:-1]  # Shorten ends of the list
            odd_STs = self.STs[::2]  # Take odd nodes only
            if self.is_call:
                return np.maximum(payoffs, odd_STs-self.K)
            else:
                return np.maximum(payoffs, self.K-odd_STs)

以下方法被覆盖，同时保留所有其他定价函数的行为：

• setup_parameters：覆盖父类方法以初始化父类的 CRR 参数，并声明新变量M为列表大小

• init_stock_price_tree：覆盖父类方法，设置一个一维 NumPy 数组作为具有M大小的栅格

• init_payoffs_tree和check_early_exercise：覆盖父类方法，只考虑奇数节点的回报

使用我们二项式 CRR 模型示例中的相同股票信息，我们可以使用二项式栅格定价来定价欧式和美式看跌期权：

In [ ]:
    eu_option = BinomialCRRLattice(
        50, 52, r=0.05, T=2, N=2, sigma=0.3, is_put=True)
In [ ] :
    print('European put:', eu_option.price())
Out[ ]:  European put: 6.245708445206432
In [ ]:
    am_option = BinomialCRRLattice(50, 52, 
        r=0.05, T=2, N=2, sigma=0.3, is_put=True, is_am=True)
In [ ] :
    print("American put:", am_option.price())
Out[ ]:   
    American put: 7.428401902704828

通过使用 CRR 二项式树格定价模型，我们得到了欧式和美式看跌期权的价格分别为$6.2457 和$7.428。

使用三项式栅格

三项式栅格与二项式栅格的工作方式基本相同。由于每个节点在每个其他节点重新组合，而不是交替节点，因此不需要从列表中提取奇数节点。由于列表的大小与二项式栅格中的大小相同，在三项式栅格定价中没有额外的存储要求，如下图所示：

三项式栅格期权定价模型的类

在 Python 中，让我们创建一个名为TrinomialLattice的类，用于继承TrinomialTreeOption类的三项式栅格实现。

就像我们为BinomialCRRLattice类所做的那样，覆盖了setup_parameters、init_stock_price_tree、init_payoffs_tree和check_early_exercise方法，而不必考虑奇数节点的回报：

In [ ]:
    import numpy as np

    """ 
    Price an option by the trinomial lattice 
    """
    class TrinomialLattice(TrinomialTreeOption):

        def setup_parameters(self):
            super(TrinomialLattice, self).setup_parameters()
            self.M = 2*self.N + 1

        def init_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 init_payoffs_tree(self):
            if self.is_call:
                return np.maximum(0, self.STs-self.K)
            else:
                return np.maximum(0, self.K-self.STs)

        def check_early_exercise(self, payoffs, node):
            self.STs = self.STs[1:-1]  # Shorten ends of the list
            if self.is_call:
                return np.maximum(payoffs, self.STs-self.K)
            else:
                return np.maximum(payoffs, self.K-self.STs)

使用与之前相同的示例，我们可以使用三项式栅格模型定价欧式和美式期权：

In [ ]:
    eu_option = TrinomialLattice(
        50, 52, r=0.05, T=2, N=2, sigma=0.3, is_put=True)
    print('European put:', eu_option.price())
Out[ ]:
    European put: 6.573565269142496
In [ ]:
    am_option = TrinomialLattice(50, 52, 
        r=0.05, T=2, N=2, sigma=0.3, is_put=True, is_am=True)
    print('American put:', am_option.price())
Out[ ]:
    American put: 7.161349217272585

输出与从三项式树期权定价模型获得的结果一致。

期权定价中的有限差分

有限差分方案与三项式树期权定价非常相似，其中每个节点依赖于另外三个节点，即上升、下降和平移。有限差分的动机是应用 Black-Scholes偏微分方程（PDE）框架（涉及函数及其偏导数），其中价格*S(t)是f(S,t)*的函数，r是无风险利率，t是到期时间，σ是基础证券的波动率：

有限差分技术往往比栅格更快地收敛，并且很好地逼近复杂的异国期权。

通过有限差分向后工作来解决 PDE，建立大小为M乘以N的离散时间网格，以反映一段时间内的资产价格，使得S和t在网格上的每个点上取以下值：

由网格符号表示，f[i,j]=f( idS, j dt)。S[max]是一个适当大的资产价格，无法在到期时间T到达。因此dS和dt是网格中每个节点之间的间隔，分别由价格和时间递增。到期时间T的终端条件对于每个S的值是一个具有行权价K的看涨期权的max(S − K, 0)和一个具有行权价K的看跌期权的max(K − S, 0)。网格从终端条件向后遍历，遵守 PDE，同时遵守网格的边界条件，例如早期行权的支付。

边界条件是节点的两端的定义值，其中i=0和i=N对于每个时间t。边界处的值用于使用 PDE 迭代计算所有其他格点的值。

网格的可视化表示如下图所示。当i和j从网格的左上角增加时，价格S趋向于网格的右下角的S[max]（可能的最高价格）：

近似 PDE 的一些方法如下：

• 前向差分：

• 后向差分：

• 中心或对称差分：

• 二阶导数：

一旦我们设置好边界条件，现在可以使用显式、隐式或 Crank-Nicolson 方法进行迭代处理。

显式方法

用于近似*f[i,j]*的显式方法由以下方程给出：

在这里，我们可以看到第一个差分是关于t的后向差分，第二个差分是关于S的中心差分，第三个差分是关于S的二阶差分。当我们重新排列项时，我们得到以下方程：

其中：

然后：

显式方法的迭代方法可以通过以下图表进行可视化表示：

编写有限差分基类

由于我们将在 Python 中编写有限差分的显式、隐式和 Crank-Nicolson 方法，让我们编写一个基类，该基类继承了所有三种方法的共同属性和函数。

我们将创建一个名为FiniteDifferences的类，该类在__init__构造方法中接受并分配所有必需的参数。price()方法是调用特定有限差分方案实现的入口点，并将按以下顺序调用这些方法：setup_boundary_conditions()、setup_coefficients()、traverse_grid()和interpolate()。这些方法的解释如下：

• setup_boundary_conditions：设置网格结构的边界条件为 NumPy 二维数组

• setup_coefficients：设置用于遍历网格结构的必要系数

• traverse_grid：向后迭代网格结构，将计算值存储到网格的第一列

• interpolate：使用网格第一列上的最终计算值，这种方法将插值这些值以找到最接近初始股价S0的期权价格

所有这些方法都是可以由派生类实现的抽象方法。如果我们忘记实现这些方法，将抛出NotImplementedError异常类型。

基类应该包含以下强制方法：

In [ ]:
    from abc import ABC, abstractmethod
    import numpy as np

    """ 
    Base class for sharing 
    attributes and functions of FD 
    """
    class FiniteDifferences(object):

        def __init__(
            self, S0, K, r=0.05, T=1, 
            sigma=0, Smax=1, M=1, N=1, is_put=False
        ):
            self.S0 = S0
            self.K = K
            self.r = r
            self.T = T
            self.sigma = sigma
            self.Smax = Smax
            self.M, self.N = M, N
            self.is_call = not is_put

            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)

        @abstractmethod
        def setup_boundary_conditions(self):
            raise NotImplementedError('Implementation required!')

        @abstractmethod
        def setup_coefficients(self):
            raise NotImplementedError('Implementation required!')

        @abstractmethod
        def traverse_grid(self):
            """  Iterate the grid backwards in time"""
            raise NotImplementedError('Implementation required!')

        @abstractmethod
        def interpolate(self):
            """ Use piecewise linear interpolation on the initial
            grid column to get the closest price at S0.
            """
            return np.interp(
                self.S0, self.boundary_conds, self.grid[:,0])

抽象基类（ABCs）提供了定义类接口的方法。@abstractmethod()装饰器声明了子类应该实现的抽象方法。与 Java 的抽象方法不同，这些方法可能有一个实现，并且可以通过super()机制从覆盖它的类中调用。

除了这些方法，我们还需要定义dS和dt，即每单位时间内S的变化和每次迭代中T的变化。我们可以将这些定义为类属性：

@property
def dS(self):
    return self.Smax/float(self.M)

@property
def dt(self):
    return self.T/float(self.N)

最后，将price()方法添加为入口点，显示调用我们讨论的抽象方法的步骤：

def price(self):
    self.setup_boundary_conditions()
    self.setup_coefficients()
    self.traverse_grid()
    return self.interpolate()

使用有限差分的显式方法对欧式期权进行定价的类

使用显式方法的有限差分的 Python 实现如下FDExplicitEu类，它继承自FiniteDifferences类并覆盖了所需的实现方法：

In [ ]:
    import numpy as np

    """ 
    Explicit method of Finite Differences 
    """
    class FDExplicitEu(FiniteDifferences):

        def setup_boundary_conditions(self):
            if self.is_call:
                self.grid[:,-1] = np.maximum(
                    0, self.boundary_conds - self.K)
                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(
                    0, self.K-self.boundary_conds)
                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]

完成网格结构的遍历后，第一列包含t=0时刻的初始资产价格的现值。NumPy 的interp函数用于执行线性插值以近似期权价值。

除了线性插值作为插值方法的最常见选择外，还可以使用其他方法，如样条或三次插值来近似期权价值。

考虑一个欧式看跌期权的例子。标的股票价格为 50 美元，波动率为 40%。看跌期权的行权价为 50 美元，到期时间为五个月。无风险利率为 10%。

我们可以使用显式方法对该期权进行定价，Smax值为100，M值为100，N值为1000：

In [ ]:
    option = FDExplicitEu(50, 50, r=0.1, T=5./12., 
        sigma=0.4, Smax=100, M=100, N=1000, is_put=True)
    print(option.price())
Out[ ]:
    4.072882278148043

当选择其他不合适的M和N值时会发生什么？

In [ ]:
    option = FDExplicitEu(50, 50, r=0.1, T=5./12., 
        sigma=0.4, Smax=100, M=80, N=100, is_put=True)
    print(option.price())
Out[ ]:   
    -8.109445694129245e+35

显然，有限差分方案的显式方法存在不稳定性问题。

隐式方法

显式方法的不稳定问题可以通过对时间的前向差分来克服。用于近似*f[i,j]*的隐式方法由以下方程给出：

在这里，可以看到隐式和显式近似方案之间唯一的区别在于第一个差分，隐式方案中使用了对t的前向差分。当我们重新排列项时，我们得到以下表达式：

其中：

在这里：

隐式方案的迭代方法可以用以下图表进行可视化表示：

从前面的图表中，我们可以注意到需要在下一次迭代步骤中计算出j+1的值，因为网格是向后遍历的。在隐式方案中，网格可以被认为在每次迭代中代表一个线性方程组，如下所示：

通过重新排列项，我们得到以下方程：

线性方程组可以表示为Ax = B的形式，我们希望在每次迭代中解出x的值。由于矩阵A是三对角的，我们可以使用 LU 分解，其中A=LU，以加快计算速度。请记住，我们在第二章中使用 LU 分解解出了线性方程组，该章节名为《金融中的线性关系的重要性》。

使用有限差分的隐式方法对欧式期权进行定价的类

隐式方案的 Python 实现在以下FDImplicitEu类中给出。我们可以从之前讨论的FDExplicitEu类中继承显式方法的实现，并覆盖感兴趣的必要方法，即setup_coefficients和traverse_grid方法：

In [ ]:
    import numpy as np
    import scipy.linalg as linalg

    """ 
    Explicit method of Finite Differences 
    """
    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

使用与显式方案相同的示例，我们可以使用隐式方案定价欧式看跌期权：

In [ ]:
    option = FDImplicitEu(50, 50, r=0.1, T=5./12., 
        sigma=0.4, Smax=100, M=100, N=1000, is_put=True)
    print(option.price())
Out[ ]:
    4.071594188049893
In [ ]:
    option = FDImplicitEu(50, 50, r=0.1, T=5./12., 
        sigma=0.4, Smax=100, M=80, N=100, is_put=True)
    print(option.price())
Out[ ]:
    4.063684691731647

鉴于当前参数和输入数据，我们可以看到隐式方案没有稳定性问题。

Crank-Nicolson 方法

另一种避免稳定性问题的方法，如显式方法中所见，是使用 Crank-Nicolson 方法。Crank-Nicolson 方法通过使用显式和隐式方法的组合更快地收敛，取两者的平均值。这导致以下方程：

这个方程也可以重写如下：

其中：

隐式方案的迭代方法可以用以下图表形式表示：

我们可以将方程视为矩阵形式的线性方程组：

其中：

我们可以在每个迭代过程中解出矩阵M。

使用有限差分的 Crank-Nicolson 方法定价欧式期权的类

Crank-Nicolson 方法的 Python 实现在以下FDCnEu类中给出，该类继承自FDExplicitEu类，并仅覆盖setup_coefficients和traverse_grid方法：

In [ ]:
    import numpy as np
    import scipy.linalg as linalg

    """ 
    Crank-Nicolson method of Finite Differences 
    """
    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

使用与显式和隐式方法相同的示例，我们可以使用 Crank-Nicolson 方法为不同的时间点间隔定价欧式看跌期权：

In [ ]:
    option = FDCnEu(50, 50, r=0.1, T=5./12.,
        sigma=0.4, Smax=100, M=100, N=1000, is_put=True)
    print(option.price())
Out[ ]:   
    4.072238354486825
In [ ]:
    option = FDCnEu(50, 50, r=0.1, T=5./12., 
        sigma=0.4, Smax=100, M=80, N=100, is_put=True)
    print(option.price())
Out[ ]: 
    4.070145703042843

从观察到的值来看，Crank-Nicolson 方法不仅避免了我们在显式方案中看到的不稳定性问题，而且比显式和隐式方法都更快地收敛。隐式方法需要更多的迭代，或者更大的N值，才能产生接近 Crank-Nicolson 方法的值。

定价异国情调的障碍期权

有限差分在定价异国情调期权方面特别有用。期权的性质将决定边界条件的规格。

在本节中，我们将看一个使用有限差分的 Crank-Nicolson 方法定价低迷障碍期权的例子。由于其相对复杂性，通常会使用其他分析方法，如蒙特卡罗方法，而不是有限差分方案。

一种低迷的选择

让我们看一个低迷期权的例子。在期权的任何生命周期中，如果标的资产价格低于*S[barrier]*障碍价格，则认为该期权毫无价值。由于在网格中，有限差分方案代表所有可能的价格点，我们只需要考虑以下价格范围的节点：

然后我们可以设置边界条件如下：

使用有限差分的 Crank-Nicolson 方法定价低迷期权的类

让我们创建一个名为FDCnDo的类，它继承自之前讨论的FDCnEu类。我们将在构造方法中考虑障碍价格，而将FDCnEu类中的 Crank-Nicolson 实现的其余部分保持不变：

In [ ]:
    import numpy as np

    """
    Price a down-and-out option by the Crank-Nicolson
    method of finite differences.
    """
    class FDCnDo(FDCnEu):

        def __init__(
            self, S0, K, r=0.05, T=1, sigma=0, 
            Sbarrier=0, Smax=1, M=1, N=1, is_put=False
        ):
            super(FDCnDo, self).__init__(
                S0, K, r=r, T=T, sigma=sigma,
                Smax=Smax, M=M, N=N, is_put=is_put
            )
            self.barrier = Sbarrier
            self.boundary_conds = np.linspace(Sbarrier, Smax, M+1)
            self.i_values = self.boundary_conds/self.dS

        @property
        def dS(self):
            return (self.Smax-self.barrier)/float(self.M)

让我们考虑一个敲出期权的例子。标的股票价格为 50 美元，波动率为 40%。期权的行权价为 50 美元，到期时间为五个月。无风险利率为 10%。敲出价格为 40 美元。

我们可以使用 Smax 为 100，M 为 120，N 为 500 来定价看涨期权和敲出看跌期权：

In [ ]:
    option = FDCnDo(50, 50, r=0.1, T=5./12., 
        sigma=0.4, Sbarrier=40, Smax=100, M=120, N=500)
    print(option.price())
Out[ ]:   
    5.491560552934787
In [ ]:
    option = FDCnDo(50, 50, r=0.1, T=5./12., sigma=0.4, 
        Sbarrier=40, Smax=100, M=120, N=500, is_put=True)
    print(option.price())
Out[ ]:
   0.5413635028954452

敲出看涨期权和敲出看跌期权的价格分别为 5.4916 美元和 0.5414 美元。

使用有限差分定价美式期权

到目前为止，我们已经定价了欧式期权和奇异期权。由于美式期权中存在提前行权的可能性，因此定价此类期权并不那么直接。在隐式 Crank-Nicolson 方法中需要迭代过程，当前期内的提前行权收益要考虑先前期内的提前行权收益。在 Crank-Nicolson 方法中，建议使用高斯-西德尔迭代方法定价美式期权。

回顾一下，在第二章中，我们讨论了在金融中线性性的重要性，我们介绍了解决线性方程组的高斯-西德尔方法，形式为 Ax=B。在这里，矩阵 A 被分解为 A=L+U，其中 L 是下三角矩阵，U 是上三角矩阵。让我们看一个 4 x 4 矩阵 A 的例子：

然后通过迭代方式获得解决方案，如下所示：

我们可以将高斯-西德尔方法调整到我们的 Crank-Nicolson 实现中，如下所示：

这个方程满足提前行权特权方程：

使用有限差分的 Crank-Nicolson 方法定价美式期权的类

让我们创建一个名为 FDCnAm 的类，该类继承自 FDCnEu 类，后者是定价欧式期权的 Crank-Nicolson 方法的对应物。setup_coefficients 方法可以被重用，同时覆盖所有其他方法，以包括先前行权的收益，如果有的话。

__init__() 构造函数和 setup_boundary_conditions() 方法在 FDCnAm 类中给出：

In [ ]:
    import numpy as np
    import sys

    """ 
    Price an American option by the Crank-Nicolson method 
    """
    class FDCnAm(FDCnEu):

        def __init__(self, S0, K, r=0.05, T=1, 
                Smax=1, M=1, N=1, omega=1, tol=0, is_put=False):
            super(FDCnAm, self).__init__(S0, K, r=r, T=T, 
                sigma=sigma, Smax=Smax, M=M, N=N, is_put=is_put)
            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(0, 
                    self.boundary_conds[1:self.M]-self.K)
            else:
                self.payoffs = np.maximum(0, 
                    self.K-self.boundary_conds[1:self.M])

            self.past_values = self.payoffs
            self.boundary_values = self.K * np.exp(
                    -self.r*self.dt*(self.N-self.j_values))

接下来，在同一类中实现 traverse_grid() 方法：

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] = \
                self.calculate_payoff_start_boundary(
                    rhs, old_values)    

            for k in range(self.M-2)[1:]:
                new_values[k] = \
                    self.calculate_payoff(
                        k, rhs, old_values, new_values)                  

            new_values[-1] = \
                self.calculate_payoff_end_boundary(
                    rhs, old_values, new_values)

            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]))

在 while 循环的每个迭代过程中，计算收益时要考虑开始和结束边界。此外，new_values 不断根据现有和先前的值进行新的收益计算替换。

在开始边界处，索引为 0 时，通过省略 alpha 值来计算收益。在类内实现 calculate_payoff_start_boundary() 方法：

 def calculate_payoff_start_boundary(self, rhs, old_values):
    payoff = old_values[0] + \
        self.omega/(1-self.beta[1]) * \
            (rhs[0] - \
             (1-self.beta[1])*old_values[0] + \
             self.gamma[1]*old_values[1])

    return max(self.payoffs[0], payoff)       

在结束边界处，最后一个索引时，通过省略 gamma 值来计算收益。在类内实现 calculate_payoff_end_boundary() 方法：

 def calculate_payoff_end_boundary(self, rhs, old_values, new_values):
    payoff = old_values[-1] + \
        self.omega/(1-self.beta[-2]) * \
            (rhs[-1] + \
             self.alpha[-2]*new_values[-2] - \
             (1-self.beta[-2])*old_values[-1])

    return max(self.payoffs[-1], payoff)

对于不在边界的收益，通过考虑 alpha 和 gamma 值来计算收益。在类内实现 calculate_payoff() 方法：

def calculate_payoff(self, k, rhs, old_values, new_values):
    payoff = 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])

    return max(self.payoffs[k], payoff)

由于新变量 values 包含我们的终端收益值作为一维数组，因此重写父类的 interpolate 方法以考虑这一变化，使用以下代码：

def interpolate(self):
    # Use linear interpolation on final values as 1D array
    return np.interp(self.S0, self.boundary_conds, self.values)

容差参数用于高斯-西德尔方法的收敛标准。omega 变量是过度松弛参数。更高的 omega 值提供更快的收敛速度，但这也伴随着算法不收敛的可能性更高。

让我们定价一个标的资产价格为 50，波动率为 40%，行权价为 50，无风险利率为 10%，到期日为五个月的美式看涨期权和看跌期权。我们选择 Smax 值为 100，M 为 100，N 为 42，omega 参数值为 1.2，容差值为 0.001：

In [ ]:
    option = FDCnAm(50, 50, r=0.1, T=5./12., 
        sigma=0.4, Smax=100, M=100, N=42, omega=1.2, tol=0.001)
    print(option.price())
Out[ ]:
    6.108682815392217
In [ ]:
    option = FDCnAm(50, 50, r=0.1, T=5./12., sigma=0.4, Smax=100, 
        M=100, N=42, omega=1.2, tol=0.001, is_put=True)
    print(option.price())
Out[ ]:   
    4.277764229383736

使用 Crank-Nicolson 方法计算美式看涨和看跌期权的价格分别为 6.109 美元和 4.2778 美元。

将所有内容整合在一起-隐含波动率建模

到目前为止，我们学到的期权定价方法中，有一些参数被假定为常数：利率、行权价、股息和波动率。在这里，感兴趣的参数是波动率。在定量研究中，波动率比率被用来预测价格趋势。

要得出隐含波动率，我们需要参考第三章金融中的非线性，在那里我们讨论了非线性函数的根查找方法。在我们的下一个示例中，我们将使用数值程序的二分法来创建一个隐含波动率曲线。

AAPL 美式看跌期权的隐含波动率

让我们考虑股票苹果（AAPL）的期权数据，这些数据是在 2014 年 10 月 3 日收集的。以下表格提供了这些细节。期权到期日为 2014 年 12 月 20 日。所列价格为买入价和卖出价的中间价：

行权价	看涨期权价格	看跌期权价格
75	30	0.16
80	24.55	0.32
85	20.1	0.6
90	15.37	1.22
92.5	10.7	1.77
95	8.9	2.54
97.5	6.95	3.55
100	5.4	4.8
105	4.1	7.75
110	2.18	11.8
115	1.05	15.96
120	0.5	20.75
125	0.26	25.8

AAPL 的最后交易价格为 99.62 美元，利率为 2.48%，股息率为 1.82%。美式期权在 78 天后到期。

利用这些信息，让我们创建一个名为ImpliedVolatilityModel的新类，它在构造函数中接受股票期权的参数。如果需要，导入我们在本章前面部分介绍的用于 LR 二项树的BinomialLROption类。还需要导入我们在第三章金融中的非线性中介绍的bisection函数。

option_valuation()方法接受K行权价和sigma波动率值，计算期权的价值。在这个例子中，我们使用BinomialLROption定价方法。

get_implied_volatilities()方法接受一个行权价和期权价格的列表，通过bisection方法计算每个价格的隐含波动率。因此，这两个列表的长度必须相同。

ImpliedVolatilityModel类的 Python 代码如下所示：

In [ ]:
    """
    Get implied volatilities from a Leisen-Reimer binomial
    tree 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(strikes)):
                # 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

导入我们在上一章讨论过的bisection函数：

In [ ]:
    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

利用这个模型，让我们使用这组特定数据找出美式看跌期权的隐含波动率：

In [ ]:
    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]
In [ ]:
    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)

隐含波动率值现在以list对象的形式存储在impvols_put变量中。让我们将这些值绘制成隐含波动率曲线：

In [ ]:
    %matplotlib inline
    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()

这将给我们提供波动率微笑，如下图所示。在这里，我们建立了一个包含 77 个步骤的 LR 树，每一步代表一天：

当然，每天定价一个期权可能并不理想，因为市场变化是以毫秒为单位的。我们使用了二分法来解决隐含波动率，这是由二项树隐含的，而不是直接从市场价格观察到的实现波动率值。

我们是否应该将这条曲线与多项式曲线拟合，以确定潜在的套利机会？或者推断曲线，以从远离实值和虚值期权的隐含波动率中获得更多见解？好吧，这些问题是供像你这样的期权交易员去发现的！

总结

在本章中，我们研究了衍生品定价中的一些数值程序，最常见的是期权。其中一种程序是使用树，二叉树是最简单的结构来建模资产信息，其中一个节点在每个时间步长延伸到另外两个节点，分别代表上升状态和下降状态。在三叉树中，每个节点在每个时间步长延伸到另外三个节点，分别代表上升状态、下降状态和无移动状态。随着树向上遍历，基础资产在每个节点处被计算和表示。然后期权采用这棵树的结构，并从期末回溯并向根部遍历，收敛到当前折现期权价格。除了二叉树和三叉树，树还可以采用 CRR、Jarrow-Rudd、Tian 或 LR 参数的形式。

通过在我们的树周围添加另一层节点，我们引入了额外的信息，从中我们可以推导出希腊字母，如 delta 和 gamma，而不会产生额外的计算成本。

晶格被引入是为了节省存储成本，相比二叉树和三叉树。在晶格定价中，只保存具有新信息的节点一次，并在以后需要不改变信息的节点上重复使用。

我们还讨论了期权定价中的有限差分方案，包括期末和边界条件。从期末条件开始，网格使用显式方法、隐式方法和 Crank-Nicolson 方法向后遍历时间。除了定价欧式和美式期权，有限差分定价方案还可以用于定价异国期权，我们看了一个定价下触及障碍期权的例子。

通过引入在第三章中学到的二分根查找方法，以及本章中的二叉 LR 树模型，我们使用美式期权的市场价格来创建隐含波动率曲线以进行进一步研究。

在下一章中，我们将研究利率和衍生品建模。

第五章：建模利率和衍生品

利率影响各个层面的经济活动。包括美联储（俗称为联邦储备系统）在内的中央银行将利率作为一种政策工具来影响经济活动。利率衍生品受到投资者的欢迎，他们需要定制的现金流需求或对利率变动的特定观点。

利率衍生品交易员面临的一个关键挑战是为这些产品制定一个良好而稳健的定价程序。这涉及理解单个利率变动的复杂行为。已经提出了几种利率模型用于金融研究。金融学中研究的一些常见模型是 Vasicek、CIR 和 Hull-White 模型。这些利率模型涉及对短期利率的建模，并依赖于因素（或不确定性来源），其中大多数只使用一个因素。已经提出了双因素和多因素利率模型。

在本章中，我们将涵盖以下主题：

• 理解收益曲线

• 估值零息债券

• 引导收益曲线

• 计算远期利率

• 计算债券的到期收益率和价格

• 使用 Python 计算债券久期和凸性

• 短期利率建模

• Vasicek 短期利率模型

• 债券期权的类型

• 定价可赎回债券期权

固定收益证券

公司和政府发行固定收益证券是为了筹集资金。这些债务的所有者借钱，并期望在债务到期时收回本金。希望借钱的发行人可能在债务寿命期间的预定时间发行固定金额的利息支付。

债务证券持有人，如美国国债、票据和债券，面临发行人违约的风险。联邦政府和市政府被认为面临最低的违约风险，因为他们可以轻松提高税收并创造更多货币来偿还未偿还的债务。

大多数债券每半年支付固定金额的利息，而有些债券每季度或每年支付。这些利息支付也被称为票息。它们以债券的面值或票面金额的百分比来报价，以年为单位。

例如，一张面值为 10,000 美元的 5 年期国库券，票面利率为 5%，每年支付 500 美元的票息，或者每 6 个月支付 250 美元的票息，直到到期日。如果利率下降，新的国库券只支付 3%的票息，那么新债券的购买者每年只会收到 300 美元的票息，而 5%债券的现有持有人将继续每年收到 500 美元的票息。由于债券的特性影响其价格，它们与当前利率水平呈反向关系：随着利率的上升，债券的价值下降。随着利率的下降，债券价格上升。

收益曲线

在正常的收益曲线环境中，长期利率高于短期利率。投资者希望在较长时间内借出资金时获得更高的回报，因为他们面临更高的违约风险。正常或正收益曲线被认为是向上倾斜的，如下图所示：

在某些经济条件下，收益曲线可能会倒挂。长期利率低于短期利率。当货币供应紧缩时会出现这种情况。投资者愿意放弃长期收益，以保护他们的短期财富。在通货膨胀高涨的时期，通货膨胀率超过票息利率时，可能会出现负利率。投资者愿意在短期内支付费用，只是为了保护他们的长期财富。倒挂的收益曲线被认为是向下倾斜的，如下图所示：

估值零息债券

零息债券是一种除到期日外不支付任何定期利息的债券，到期时偿还本金或面值。零息债券也被称为纯贴现债券。

零息债券的价值可以如下计算：

在这里，y是债券的年复利收益率或利率，t是债券到期的剩余时间。

让我们来看一个面值 100 美元的五年期零息债券的例子。年收益率为 5%，年复利。价格可以如下计算：

一个简单的 Python 零息债券计算器可以用来说明这个例子：

In [ ]:
    def zero_coupon_bond(par, y, t):
        """
        Price a zero coupon bond.

        :param par: face value of the bond.
        :param y: annual yield or rate of the bond.
        :param t: time to maturity, in years.
        """
        return par/(1+y)**t

使用上述例子，我们得到以下结果：

In [ ]:
    print(zero_coupon_bond(100, 0.05, 5))
Out[ ]:
    78.35261664684589

在上述例子中，我们假设投资者能够以 5%的年利率在 5 年内年复利投资 78.35 美元。

现在我们有了一个零息债券计算器，我们可以用它通过抽靴法收益率曲线来确定零息率，如下一节所述。

即期和零息率

随着复利频率的增加（比如，从年复利到日复利），货币的未来价值达到了一个指数极限。也就是说，今天的 100 美元在以连续复利率R投资T时间后将达到未来价值 100e^(RT)。如果我们对一个在未来时间T支付 100 美元的证券进行贴现，使用连续复利贴现率R，其在时间零的价值为。这个利率被称为即期利率。

即期利率代表了当前的利率，用于几种到期日，如果我们现在想要借款或出借资金。零息率代表了零息债券的内部收益率。

通过推导具有不同到期日的债券的即期利率，我们可以使用零息债券通过抽靴过程构建当前收益率曲线。

抽靴法收益率曲线

短期即期利率可以直接从各种短期证券（如零息债券、国库券、票据和欧元存款）中推导出来。然而，长期即期利率通常是通过一个抽靴过程从长期债券的价格中推导出来的，考虑到与票息支付日期相对应的到期日的即期利率。在获得短期和长期即期利率之后，就可以构建收益率曲线。

抽靴法收益率曲线的一个例子

让我们通过一个例子来说明抽靴法收益率曲线。以下表格显示了具有不同到期日和价格的债券列表：

债券面值（美元）	到期年限（年）	年息（美元）	债券现金价格（美元）
100	0.25	0	97.50
100	0.50	0	94.90
100	1.00	0	90.00
100	1.50	8	96.00
100	2.00	12	101.60

今天以 97.50 美元购买三个月的零息债券的投资者将获得 2.50 美元的利息。三个月的即期利率可以如下计算：

因此，3 个月的零息率是 10.127%，采用连续复利。零息债券的即期利率如下表所示：

到期年限（年）	即期利率（百分比）
0.25	10.127
0.50	10.469
1.00	10.536

使用这些即期利率，我们现在可以如下定价 1.5 年期债券：

y的值可以通过重新排列方程轻松求解，如下所示：

有了 1.5 年期债券的即期利率为 10.681%，我们可以使用它来定价 2 年期债券，年息为 6 美元，半年付息，如下所示：

重新排列方程并解出y，我们得到 2 年期债券的即期利率为 10.808。

通过这种迭代过程，按照到期日递增的顺序计算每个债券的即期利率，并在下一次迭代中使用它，我们得到了一个可以用来构建收益率曲线的不同到期日的即期利率列表。

编写收益率曲线引导类

编写 Python 代码引导收益率曲线并生成图表输出的步骤如下：

1. 创建一个名为BootstrapYieldCurve的类，它将在 Python 代码中实现收益率曲线的引导：

import math

class BootstrapYieldCurve(object):    

    def __init__(self):
        self.zero_rates = dict()
        self.instruments = dict()

1. 在构造函数中，声明了两个zero_rates和instruments字典变量，并将被几种方法使用，如下所示：

• 添加一个名为add_instrument()的方法，该方法将债券信息的元组附加到以到期时间为索引的instruments字典中。此方法的编写如下：

def add_instrument(self, par, T, coup, price, compounding_freq=2):
    self.instruments[T] = (par, coup, price, compounding_freq)

• 添加一个名为get_maturities()的方法，它简单地按升序返回一个可用到期日列表。此方法的编写如下：

def get_maturities(self):
    """ 
    :return: a list of maturities of added instruments 
    """
    return sorted(self.instruments.keys())

• 添加一个名为get_zero_rates()的方法，该方法对收益率曲线进行引导，计算沿该收益率曲线的即期利率，并按到期日升序返回零息率的列表。该方法的编写如下：

def get_zero_rates(self):
    """ 
    Returns a list of spot rates on the yield curve.
    """
    self.bootstrap_zero_coupons()    
    self.get_bond_spot_rates()
    return [self.zero_rates[T] for T in self.get_maturities()]

• 添加一个名为bootstrap_zero_coupons()的方法，该方法计算给定零息债券的即期利率，并将其添加到以到期日为索引的zero_rates字典中。此方法的编写如下：

def bootstrap_zero_coupons(self):
    """ 
    Bootstrap the yield curve with zero coupon instruments first.
    """
    for (T, instrument) in self.instruments.items():
        (par, coup, price, freq) = instrument
        if coup == 0:
            spot_rate = self.zero_coupon_spot_rate(par, price, T)
            self.zero_rates[T] = spot_rate  

• 添加一个名为zero_coupon_spot_rate()的方法，该方法计算零息债券的即期利率。此方法由bootstrap_zero_coupons()调用，并编写如下：

def zero_coupon_spot_rate(self, par, price, T):
    """ 
    :return: the zero coupon spot rate with continuous compounding.
    """
    spot_rate = math.log(par/price)/T
    return spot_rate

• 添加一个名为get_bond_spot_rates()的方法，它计算非零息债券的即期利率，并将其添加到以到期日为索引的zero_rates字典中。此方法的编写如下：

def get_bond_spot_rates(self):
    """ 
    Get spot rates implied by bonds, using short-term instruments.
    """
    for T in self.get_maturities():
        instrument = self.instruments[T]
        (par, coup, price, freq) = instrument
        if coup != 0:
            spot_rate = self.calculate_bond_spot_rate(T, instrument)
            self.zero_rates[T] = spot_rate

• 添加一个名为calculate_bond_spot_rate()的方法，该方法由get_bond_spot_rates()调用，用于计算特定到期期间的即期利率。此方法的编写如下：

def calculate_bond_spot_rate(self, T, instrument):
    try:
        (par, coup, price, freq) = instrument
        periods = T*freq
        value = price
        per_coupon = coup/freq
        for i in range(int(periods)-1):
            t = (i+1)/float(freq)
            spot_rate = self.zero_rates[t]
            discounted_coupon = per_coupon*math.exp(-spot_rate*t)
            value -= discounted_coupon

        last_period = int(periods)/float(freq)        
        spot_rate = -math.log(value/(par+per_coupon))/last_period
        return spot_rate
    except:
        print("Error: spot rate not found for T=", t)

1. 实例化BootstrapYieldCurve类，并从前表中添加每个债券的信息：

In [ ]:
    yield_curve = BootstrapYieldCurve()
    yield_curve.add_instrument(100, 0.25, 0., 97.5)
    yield_curve.add_instrument(100, 0.5, 0., 94.9)
    yield_curve.add_instrument(100, 1.0, 0., 90.)
    yield_curve.add_instrument(100, 1.5, 8, 96., 2)
    yield_curve.add_instrument(100, 2., 12, 101.6, 2)
In [ ]:
    y = yield_curve.get_zero_rates()
    x = yield_curve.get_maturities()

1. 在类实例中调用get_zero_rates()方法会返回一个即期利率列表，顺序与分别存储在x和y变量中的到期日相同。发出以下 Python 代码以在图表上绘制x和y：

In [ ]:
    %pylab inline

    fig = plt.figure(figsize=(12, 8))
    plot(x, y)
    title("Zero Curve") 
    ylabel("Zero Rate (%)")
    xlabel("Maturity in Years");

1. 我们得到以下的收益率曲线：

在正常的收益率曲线环境中，随着到期日的增加，利率也会增加，我们得到一个上升的收益率曲线。

远期利率

计划在以后投资的投资者可能会好奇知道未来的利率会是什么样子，这是由当今的利率期限结构所暗示的。例如，您可能会问，“一年后的一年期即期利率是多少？”为了回答这个问题，您可以使用以下公式计算*T[1]和T[2]*之间的期间的远期利率：

在这里，*r[1]和r[2]分别是T[1]和T[2]*时期的连续复利年利率。

以下的ForwardRates类帮助我们从即期利率列表生成远期利率列表：

class ForwardRates(object):

    def __init__(self):
        self.forward_rates = []
        self.spot_rates = dict()

    def add_spot_rate(self, T, spot_rate):
        self.spot_rates[T] = spot_rate

    def get_forward_rates(self):
        """
        Returns a list of forward rates
        starting from the second time period.
        """
        periods = sorted(self.spot_rates.keys())
        for T2, T1 in zip(periods, periods[1:]):
            forward_rate = self.calculate_forward_rate(T1, T2)
            self.forward_rates.append(forward_rate)

        return self.forward_rates

    def calculate_forward_rate(self, T1, T2):
        R1 = self.spot_rates[T1]
        R2 = self.spot_rates[T2]
        forward_rate = (R2*T2-R1*T1)/(T2-T1)
        return forward_rate        

使用从前述收益率曲线派生的即期利率，我们得到以下结果：

In [ ]:
    fr = ForwardRates()
    fr.add_spot_rate(0.25, 10.127)
    fr.add_spot_rate(0.50, 10.469)
    fr.add_spot_rate(1.00, 10.536)
    fr.add_spot_rate(1.50, 10.681)
    fr.add_spot_rate(2.00, 10.808)
In [ ]:
    print(fr.get_forward_rates())
Out[ ]:
    [10.810999999999998, 10.603, 10.971, 11.189]

调用ForwardRates类的get_forward_rates()方法返回一个从下一个时间段开始的远期利率列表。

计算到期收益率

到期收益（YTM）衡量了债券隐含的利率，考虑了所有未来票息支付和本金的现值。假设债券持有人可以以 YTM 利率投资收到的票息，直到债券到期；根据风险中性预期，收到的支付应与债券支付的价格相同。

让我们来看一个例子，一个 5.75％的债券，将在 1.5 年内到期，面值为 100。债券价格为 95.0428 美元，票息每半年支付一次。定价方程可以陈述如下：

！

这里：

• c是每个时间段支付的票面金额

• T是以年为单位的支付时间段

• n是票息支付频率

• y是我们感兴趣的 YTM 解决方案

解决 YTM 通常是一个复杂的过程，大多数债券 YTM 计算器使用牛顿法作为迭代过程。

债券 YTM 计算器由以下bond_ytm()函数说明：

import scipy.optimize as optimize

def bond_ytm(price, par, T, coup, freq=2, guess=0.05):
    freq = float(freq)
    periods = T*2
    coupon = coup/100.*par
    dt = [(i+1)/freq for i in range(int(periods))]
    ytm_func = lambda y: \
        sum([coupon/freq/(1+y/freq)**(freq*t) for t in dt]) +\
        par/(1+y/freq)**(freq*T) - price

    return optimize.newton(ytm_func, guess)

请记住，我们在第三章中介绍了牛顿法和其他非线性函数根求解器的使用，金融中的非线性。对于这个 YTM 计算器函数，我们使用了scipy.optimize包来解决 YTM。

使用债券示例的参数，我们得到以下结果：

In [ ] :
    ytm = bond_ytm(95.0428, 100, 1.5, 5.75, 2)
In [ ]:
    print(ytm)
Out[ ]:
    0.09369155345239522

债券的 YTM 为 9.369％。现在我们有一个债券 YTM 计算器，可以帮助我们比较债券的预期回报与其他证券的回报。

计算债券价格

当 YTM 已知时，我们可以以与使用定价方程相同的方式获得债券价格。这是由bond_price()函数实现的：

In [ ]:
    def bond_price(par, T, ytm, coup, freq=2):
        freq = float(freq)
        periods = T*2
        coupon = coup/100.*par
        dt = [(i+1)/freq for i in range(int(periods))]
        price = sum([coupon/freq/(1+ytm/freq)**(freq*t) for t in dt]) + \
            par/(1+ytm/freq)**(freq*T)
        return price

插入先前示例中的相同值，我们得到以下结果：

In [ ]:
    price = bond_price(100, 1.5, ytm, 5.75, 2)
    print(price)
Out[ ]:   
    95.04279999999997

这给我们了先前示例中讨论的相同原始债券价格，计算到期收益。使用bond_ytm()和bond_price()函数，我们可以将这些应用于债券定价的进一步用途，例如查找债券的修改后持续时间和凸性。债券的这两个特征对于债券交易员来说非常重要，可以帮助他们制定各种交易策略并对冲风险。

债券持续时间

持续时间是债券价格对收益变化的敏感度度量。一些持续时间度量是有效持续时间，麦考利持续时间和修改后的持续时间。我们将讨论的持续时间类型是修改后的持续时间，它衡量了债券价格相对于收益变化的百分比变化（通常为 1％或 100 个基点（bps））。

债券的持续时间越长，对收益变化的敏感度就越高。相反，债券的持续时间越短，对收益变化的敏感度就越低。

债券的修改后持续时间可以被认为是价格和收益之间关系的第一导数：

！

这里：

• • dY *是给定的收益变化

• P^−是债券因* dY *减少而导致的价格

• *P^+是债券因 dY *增加而导致的价格

• *P[0]*是债券的初始价格

应该注意，持续时间描述了Y的小变化对价格-收益关系的线性关系。由于收益曲线不是线性的，使用较大的dY不能很好地近似持续时间度量。

修改后的持续时间计算器的实现在以下bond_mod_duration()函数中给出。它使用了本章前面讨论的bond_ytm()函数，计算到期收益，来确定具有给定初始值的债券的收益。此外，它使用bond_price()函数来确定具有给定收益变化的债券的价格：

In [ ]:
    def bond_mod_duration(price, par, T, coup, freq, dy=0.01):
        ytm = bond_ytm(price, par, T, coup, freq)

        ytm_minus = ytm - dy    
        price_minus = bond_price(par, T, ytm_minus, coup, freq)

        ytm_plus = ytm + dy
        price_plus = bond_price(par, T, ytm_plus, coup, freq)

        mduration = (price_minus-price_plus)/(2*price*dy)
        return mduration

我们可以找出之前讨论的 5.75％债券的修改后持续时间，计算到期收益，它将在 1.5 年内到期，面值为 100，债券价格为 95.0428：

In [ ]:
    mod_duration = bond_mod_duration(95.0428, 100, 1.5, 5.75, 2)
In [ ]:
    print(mod_duration)
Out[ ]:
    1.3921935426561034

债券的修正久期为 1.392 年。

债券凸性

凸性是债券久期对收益率变化的敏感度度量。将凸性视为价格和收益率之间关系的二阶导数：

债券交易员使用凸性作为风险管理工具，以衡量其投资组合中的市场风险。相对于债券久期和收益率相同的低凸性投资组合，高凸性投资组合受利率波动的影响较小。因此，其他条件相同的情况下，高凸性债券比低凸性债券更昂贵。

债券凸性的实现如下：

In [ ]:
    def bond_convexity(price, par, T, coup, freq, dy=0.01):
        ytm = bond_ytm(price, par, T, coup, freq)

        ytm_minus = ytm - dy    
        price_minus = bond_price(par, T, ytm_minus, coup, freq)

        ytm_plus = ytm + dy
        price_plus = bond_price(par, T, ytm_plus, coup, freq)

        convexity = (price_minus + price_plus - 2*price)/(price*dy**2)
        return convexity

现在我们可以找到之前讨论的 5.75%债券的凸性，它将在 1.5 年后到期，票面价值为 100，债券价格为 95.0428：

In [ ]:
    convexity = bond_convexity(95.0428, 100, 1.5, 5.75, 2)
In [ ]:
    print(convexity)
Out[ ] :    
    2.633959390331875

债券的凸性为 2.63。对于两个具有相同票面价值、票息和到期日的债券，它们的凸性可能不同，这取决于它们在收益率曲线上的位置。相对于收益率的变化，高凸性债券的价格变化更大。

短期利率建模

在短期利率建模中，短期利率*r(t)*是特定时间的即期利率。它被描述为收益率曲线上无限短的时间内的连续复利化年化利率。短期利率在利率模型中采用随机变量的形式，其中利率可能在每个时间点上以微小的变化。短期利率模型试图模拟利率随时间的演变，并希望描述特定时期的经济状况。

短期利率模型经常用于评估利率衍生品。债券、信用工具、抵押贷款和贷款产品对利率变化敏感。短期利率模型被用作利率组成部分，结合定价实现，如数值方法，以帮助定价这些衍生品。

利率建模被认为是一个相当复杂的话题，因为利率受到多种因素的影响，如经济状态、政治决策、政府干预以及供求法则。已经提出了许多利率模型，以解释利率的各种特征。

在本节中，我们将研究金融研究中使用的一些最流行的一因素短期利率模型，即瓦西切克、考克斯-英格索尔-罗斯、伦德尔曼和巴特尔、布伦南和施瓦茨模型。使用 Python，我们将执行一条路径模拟，以获得对利率路径过程的一般概述。金融学中常讨论的其他模型包括何厉、赫尔-怀特和布莱克-卡拉辛基。

瓦西切克模型

在一因素瓦西切克模型中，短期利率被建模为单一随机因素：

在这里，K、θ和σ是常数，σ是瞬时标准差。W(t)是随机维纳过程。瓦西切克遵循奥恩斯坦-乌伦贝克过程，模型围绕均值θ回归，K是均值回归速度。因此，利率可能变为负值，这在大多数正常经济条件下是不希望的特性。

为了帮助我们理解这个模型，以下代码生成了一组利率：

In [ ]:
    import math
    import numpy as np

    def vasicek(r0, K, theta, sigma, T=1., N=10, seed=777):    
        np.random.seed(seed)
        dt = T/float(N)    
        rates = [r0]
        for i in range(N):
            dr = K*(theta-rates[-1])*dt + \
                sigma*math.sqrt(dt)*np.random.normal()
            rates.append(rates[-1]+dr)

        return range(N+1), rates

vasicek()函数返回瓦西切克模型的一组时间段和利率。它接受一些输入参数：r0是t=0时的初始利率；K、theta和sigma是常数；T是以年为单位的期间；N是建模过程的间隔数；seed是 NumPy 标准正态随机数生成器的初始化值。

假设当前利率接近零，为 0.5%，长期均值水平theta为0.15，瞬时波动率sigma为 5%。我们将使用T值为10和N值为200来模拟不同均值回归速度K的利率，使用值为0.002，0.02和0.2：

In [ ]:
    %pylab inline

    fig = plt.figure(figsize=(12, 8))

    for K in [0.002, 0.02, 0.2]:
        x, y = vasicek(0.005, K, 0.15, 0.05, T=10, N=200)
        plot(x,y, label='K=%s'%K)
        pylab.legend(loc='upper left');

    pylab.legend(loc='upper left')
    pylab.xlabel('Vasicek model');

运行前述命令后，我们得到以下图表：

在这个例子中，我们只运行了一个模拟，以查看 Vasicek 模型的利率是什么样子。请注意，利率在某个时候变为负值。当均值回归速度K较高时，该过程更快地达到其长期水平 0.15。

Cox-Ingersoll-Ross 模型

Cox-Ingersoll-Ross（CIR）模型是一个一因素模型，旨在解决 Vasicek 模型中发现的负利率。该过程如下：

术语随着短期利率的增加而增加标准差。现在vasicek()函数可以重写为 Python 中的 CIR 模型：

In [ ]:
    import math
    import numpy as np

    def CIR(r0, K, theta, sigma, T=1.,N=10,seed=777):        
        np.random.seed(seed)
        dt = T/float(N)    
        rates = [r0]
        for i in range(N):
            dr = K*(theta-rates[-1])*dt + \
                sigma*math.sqrt(rates[-1])*\
                math.sqrt(dt)*np.random.normal()
            rates.append(rates[-1] + dr)

        return range(N+1), rates

使用Vasicek 模型部分中给出的相同示例，假设当前利率为 0.5%，theta为0.15，sigma为0.05。我们将使用T值为10和N值为200来模拟不同均值回归速度K的利率，使用值为0.002，0.02和0.2：

In [ ] :
    %pylab inline

    fig = plt.figure(figsize=(12, 8))

    for K in [0.002, 0.02, 0.2]:
        x, y = CIR(0.005, K, 0.15, 0.05, T=10, N=200)
        plot(x,y, label='K=%s'%K)

    pylab.legend(loc='upper left')
    pylab.xlabel('CRR model');

以下是前述命令的输出：

请注意，CIR 利率模型没有负利率值。

Rendleman 和 Bartter 模型

在 Rendleman 和 Bartter 模型中，短期利率过程如下：

这里，瞬时漂移是θr(t)，瞬时标准差为σr(t)。Rendleman 和 Bartter 模型可以被视为几何布朗运动，类似于对数正态分布的股价随机过程。该模型缺乏均值回归属性。均值回归是利率似乎被拉回到长期平均水平的现象。

以下 Python 代码模拟了 Rendleman 和 Bartter 的利率过程：

In [ ]:
    import math
    import numpy as np

    def rendleman_bartter(r0, theta, sigma, T=1.,N=10,seed=777):        
        np.random.seed(seed)
        dt = T/float(N)    
        rates = [r0]
        for i in range(N):
            dr = theta*rates[-1]*dt + \
                sigma*rates[-1]*math.sqrt(dt)*np.random.normal()
            rates.append(rates[-1] + dr)

        return range(N+1), rates

我们将继续使用前几节中的示例并比较模型。

假设当前利率为 0.5%，sigma为0.05。我们将使用T值为10和N值为200来模拟不同瞬时漂移theta的利率，使用值为0.01，0.05和0.1：

In [ ]:
    %pylab inline

    fig = plt.figure(figsize=(12, 8))

    for theta in [0.01, 0.05, 0.1]:
        x, y = rendleman_bartter(0.005, theta, 0.05, T=10, N=200)
        plot(x,y, label='theta=%s'%theta)

    pylab.legend(loc='upper left')
    pylab.xlabel('Rendleman and Bartter model');

以下图表是前述命令的输出：

总的来说，该模型缺乏均值回归属性，并向长期平均水平增长。

Brennan 和 Schwartz 模型

Brennan 和 Schwartz 模型是一个双因素模型，其中短期利率向长期利率作为均值回归，也遵循随机过程。短期利率过程如下：

可以看出，Brennan 和 Schwartz 模型是几何布朗运动的另一种形式。

我们的 Python 代码现在可以这样实现：

In [ ]:
    import math
    import numpy as np

    def brennan_schwartz(r0, K, theta, sigma, T=1., N=10, seed=777):    
        np.random.seed(seed)
        dt = T/float(N)    
        rates = [r0]
        for i in range(N):
            dr = K*(theta-rates[-1])*dt + \
                sigma*rates[-1]*math.sqrt(dt)*np.random.normal()
            rates.append(rates[-1] + dr)

        return range(N+1), rates

假设当前利率保持在 0.5%，长期均值水平theta为 0.006。sigma为0.05。我们将使用T值为10和N值为200来模拟不同均值回归速度K的利率，使用值为0.2，0.02和0.002：

In [ ]:
    %pylab inline

    fig = plt.figure(figsize=(12, 8))

    for K in [0.2, 0.02, 0.002]:
        x, y = brennan_schwartz(0.005, K, 0.006, 0.05, T=10, N=200)
        plot(x,y, label='K=%s'%K)

    pylab.legend(loc='upper left')
    pylab.xlabel('Brennan and Schwartz model');

运行前述命令后，我们将得到以下输出：

当 k 为 0.2 时，均值回归速度最快，达到长期均值 0.006。

债券期权

当债券发行人，如公司，发行债券时，他们面临的风险之一是利率风险。利率下降时，债券价格上升。现有债券持有人会发现他们的债券更有价值，而债券发行人则处于不利地位，因为他们将发行高于市场利率的利息支付。相反，当利率上升时，债券发行人处于有利地位，因为他们能够继续按债券合同规定的低利息支付。

为了利用利率变化，债券发行人可以在债券中嵌入期权。这使得发行人有权利，但没有义务，在特定时间段内以预定价格买入或卖出发行的债券。美式债券期权允许发行人在债券的任何时间内行使期权的权利。欧式债券期权允许发行人在特定日期行使期权的权利。行使日期的确切方式因债券期权而异。一些发行人可能选择在债券在市场上流通超过一年后行使债券期权的权利。一些发行人可能选择在几个特定日期中的一个上行使债券期权的权利。无论债券的行使日期如何，您可以按以下方式定价嵌入式期权的债券：

债券价格=没有期权的债券价格-嵌入式期权的价格

没有期权的债券定价相当简单：未来日期收到的债券的现值，包括所有票息支付。必须对未来的理论利率进行一些假设，以便将票息支付再投资。这样的假设可能是短期利率模型所暗示的利率变动，我们在前一节中介绍了短期利率建模。另一个假设可能是在二项或三项树中的利率变动。为简单起见，在债券定价研究中，我们将定价零息债券，这些债券在债券的寿命期间不发行票息。

要定价期权，必须确定可行的行使日期。从债券的未来价值开始，将债券价格与期权的行使价格进行比较，并使用数值程序（如二项树）回溯到现在的时间。这种价格比较是在债券期权可能行使的时间点进行的。根据无套利理论，考虑到行使时债券的现值超额值，我们得到了期权的价格。为简单起见，在本章后面的债券定价研究中，定价可赎回债券期权，我们将把零息债券的嵌入式期权视为美式期权。

可赎回债券

在利率较高的经济条件下，债券发行人可能面临利率下降的风险，并不得不继续发行高于市场利率的利息支付。因此，他们可能选择发行可赎回债券。可赎回债券包含一项嵌入式协议，约定在约定日期赎回债券。现有债券持有人被认为已经向债券发行人出售了一项认购期权。

如果利率下降，公司有权在特定价格回购债券的期间行使该选择权，他们可能会选择这样做。公司随后可以以较低的利率发行新债券。这也意味着公司能够以更高的债券价格形式筹集更多资本。

可回售债券

与可赎回债券不同，可赎回债券的持有人有权利，但没有义务，在一定期限内以约定价格将债券卖回给发行人。可赎回债券的持有人被认为是从债券发行人购买了一个认沽期权。当利率上升时，现有债券的价值变得更不值钱，可赎回债券持有人更有动力行使以更高行使价格出售债券的权利。由于可赎回债券对买方更有利而对发行人不利，它们通常比可赎回债券更少见。可赎回债券的变体可以在贷款和存款工具的形式中找到。向金融机构存入固定利率存款的客户在指定日期收到利息支付。他们有权随时提取存款。因此，固定利率存款工具可以被视为带有内嵌美式认沽期权的债券。

希望从银行借款的投资者签订贷款协议，在协议的有效期内进行利息支付，直到债务连同本金和约定利息全部偿还。银行可以被视为在债券上购买了一个看跌期权。在某些情况下，银行可能行使赎回贷款协议全部价值的权利。

因此，可赎回债券的价格可以如下所示：

可赎回债券的价格=无期权债券价格+认沽期权价格

可转换债券

公司发行的可转换债券包含一个内嵌期权，允许持有人将债券转换为一定数量的普通股。债券转换为股票的数量由转换比率定义，该比率被确定为使股票的美元金额与债券价值相同。

可转换债券与可赎回债券有相似之处。它们允许债券持有人在约定时间以约定的转换比率行使债券，换取相等数量的股票。可转换债券通常发行的票面利率低于不可转换债券，以补偿行使权利的附加价值。

当可转换债券持有人行使其股票权利时，公司的债务减少。另一方面，随着流通股数量的增加，公司的股票变得更加稀释，公司的股价预计会下跌。

随着公司股价的上涨，可转换债券价格往往会上涨。相反，随着公司股价的下跌，可转换债券价格往往会下跌。

优先股

优先股是具有债券特性的股票。优先股股东在普通股股东之前对股利支付具有优先权，通常作为其面值的固定百分比进行协商。虽然不能保证股利支付，但所有股利都优先支付给优先股股东而不是普通股股东。在某些优先股协议中，未按约支付的股利可能会累积直到以后支付。这些优先股被称为累积。

优先股的价格通常与其普通股同步变动。它们可能具有与普通股股东相关的表决权。在破产情况下，优先股在清算时对其票面价值具有第一顺位权。

定价可赎回债券期权

在本节中，我们将研究定价可赎回债券。我们假设要定价的债券是一种带有内嵌欧式认购期权的零息支付债券。可赎回债券的价格可以如下所示：

可赎回债券的价格=无期权债券价格-认购期权价格

通过 Vasicek 模型定价零息债券

在时间t和当前利率r下，面值为 1 的零息债券的价值定义如下：

由于利率r总是在变化，我们将零息债券重写如下：

现在，利率r是一个随机过程，考虑从时间t到T的债券价格，其中T是零息债券的到期时间。

为了对利率r进行建模，我们可以使用短期利率模型作为随机过程之一。为此，我们将使用 Vasicek 模型来对短期利率过程进行建模。

对于对数正态分布变量X的期望值如下：

对对数正态分布变量X的矩：

我们得到了对数正态分布变量的期望值，我们将在零息债券的利率过程中使用。

记住 Vasicek 短期利率过程模型：

然后，*r(t)*可以如下导出：

利用特征方程和 Vasicek 模型的利率变动，我们可以重写零息债券价格的期望值：

这里：

零息债券价格的 Python 实现在exact_zcb函数中给出：

In [ ]:
    import numpy as np
    import math

    def exact_zcb(theta, kappa, sigma, tau, r0=0.):
        B = (1 - np.exp(-kappa*tau)) / kappa
        A = np.exp((theta-(sigma**2)/(2*(kappa**2)))*(B-tau) - \
                   (sigma**2)/(4*kappa)*(B**2))
        return A * np.exp(-r0*B)

例如，我们有兴趣找出多种到期日的零息债券价格。我们使用 Vasicek 短期利率过程，theta值为0.5，kappa值为0.02，sigma值为0.03，初始利率r0为0.015进行建模。

将这些值代入exact_zcb函数，我们得到了从 0 到 25 年的时间段内以 0.5 年为间隔的零息债券价格，并绘制出图表：

In [ ]:    
    Ts = np.r_[0.0:25.5:0.5]
    zcbs = [exact_zcb(0.5, 0.02, 0.03, t, 0.015) for t in Ts]
In [ ]:
    %pylab inline

    fig = plt.figure(figsize=(12, 8))
    plt.title("Zero Coupon Bond (ZCB) Values by Time")
    plt.plot(Ts, zcbs, label='ZCB')
    plt.ylabel("Value ($)")
    plt.xlabel("Time in years")
    plt.legend()
    plt.grid(True)
    plt.show()

以下图表是上述命令的输出：

提前行使权的价值

可赎回债券的发行人可以按合同规定的约定价格赎回债券。为了定价这种债券，可以将折现的提前行使价值定义如下：

这里，k是行权价格与票面价值的价格比率，r是行权价格的利率。

然后，提前行使期权的 Python 实现可以写成如下形式：

In [ ]:
    import math

    def exercise_value(K, R, t):
        return K*math.exp(-R*t)

在上述示例中，我们有兴趣对行权比率为 0.95 且初始利率为 1.5%的认购期权进行定价。然后，我们可以将这些值作为时间函数绘制，并将它们叠加到零息债券价格的图表上，以更好地呈现零息债券价格与可赎回债券价格之间的关系：

In [ ]:
    Ts = np.r_[0.0:25.5:0.5]
    Ks = [exercise_value(0.95, 0.015, t) for t in Ts]
    zcbs = [exact_zcb(0.5, 0.02, 0.03, t, 0.015) for t in Ts]
In [ ]:
    import matplotlib.pyplot as plt

    fig = plt.figure(figsize=(12, 8))
    plt.title("Zero Coupon Bond (ZCB) and Strike (K) Values by Time")
    plt.plot(Ts, zcbs, label='ZCB')
    plt.plot(Ts, Ks, label='K', linestyle="--", marker=".")
    plt.ylabel("Value ($)")
    plt.xlabel("Time in years")
    plt.legend()
    plt.grid(True)
    plt.show()

以下是上述命令的输出：

从上述图表中，我们可以近似计算可赎回零息债券的价格。由于债券发行人拥有行权，可赎回零息债券的价格可以如下所述：

这种可赎回债券价格是一种近似值，考虑到当前的利率水平。下一步将是通过进行一种政策迭代来处理提前行使，这是一种用于确定最佳提前行使价值及其对其他节点的影响的循环，并检查它们是否到期提前行使。在实践中，这样的迭代只会发生一次。

通过有限差分的政策迭代

到目前为止，我们已经在我们的短期利率过程中使用了 Vasicek 模型来模拟零息债券。我们可以通过有限差分进行政策迭代，以检查提前行使条件及其对其他节点的影响。我们将使用有限差分的隐式方法进行数值定价程序，如第四章*，期权定价的数值程序*中所讨论的那样。

让我们创建一个名为VasicekCZCB的类，该类将包含用于实现 Vasicek 模型定价可赎回零息债券的所有方法。该类及其构造函数定义如下：

import math
import numpy as np
import scipy.stats as st

class VasicekCZCB:

    def __init__(self):
        self.norminv = st.distributions.norm.ppf
        self.norm = st.distributions.norm.cdf    

在构造函数中，norminv和normv变量对于所有需要计算 SciPy 的逆正态累积分布函数和正态累积分布函数的方法都是可用的。

有了这个基类，让我们讨论所需的方法并将它们添加到我们的类中：

• 添加vasicek_czcb_values()方法作为开始定价过程的入口点。r0变量是时间t=0的短期利率；R是债券价格的零利率；ratio是债券的面值每单位的行权价格；T是到期时间；sigma是短期利率r的波动率；kappa是均值回归率；theta是短期利率过程的均值；M是有限差分方案中的步数；prob是vasicek_limits方法后用于确定短期利率的正态分布曲线上的概率；max_policy_iter是用于找到提前行使节点的最大政策迭代次数；grid_struct_const是确定calculate_N()方法中的N的dt移动的最大阈值；rs是短期利率过程遵循的利率列表。

该方法返回一列均匀间隔的短期利率和一列期权价格，写法如下：

def vasicek_czcb_values(self, r0, R, ratio, T, sigma, kappa, theta,
                        M, prob=1e-6, max_policy_iter=10, 
                        grid_struct_const=0.25, rs=None):
    (r_min, dr, N, dtau) = \
        self.vasicek_params(r0, M, sigma, kappa, theta,
                            T, prob, grid_struct_const, rs)
    r = np.r_[0:N]*dr + r_min
    v_mplus1 = np.ones(N)

    for i in range(1, M+1):
        K = self.exercise_call_price(R, ratio, i*dtau)
        eex = np.ones(N)*K
        (subdiagonal, diagonal, superdiagonal) = \
            self.vasicek_diagonals(
                sigma, kappa, theta, r_min, dr, N, dtau)
        (v_mplus1, iterations) = \
            self.iterate(subdiagonal, diagonal, superdiagonal,
                         v_mplus1, eex, max_policy_iter)
    return r, v_mplus1

• 添加vasicek_params()方法来计算 Vasicek 模型的隐式方案参数。它返回一个元组r_min，dr，N和dt。如果未向rs提供值，则r_min到r_max的值将由vasicek_limits()方法自动生成，作为prob的正态分布函数的函数。该方法的写法如下：

def vasicek_params(self, r0, M, sigma, kappa, theta, T,
                  prob, grid_struct_const=0.25, rs=None):
    if rs is not None:
        (r_min, r_max) = (rs[0], rs[-1])
    else:
        (r_min, r_max) = self.vasicek_limits(
            r0, sigma, kappa, theta, T, prob)      

    dt = T/float(M)
    N = self.calculate_N(grid_struct_const, dt, sigma, r_max, r_min)
    dr = (r_max-r_min)/(N-1)

    return (r_min, dr, N, dt)

• 添加calculate_N()方法，该方法由vasicek_params()方法使用，用于计算网格大小参数N。该方法的写法如下：

def calculate_N(self, max_structure_const, dt, sigma, r_max, r_min):
    N = 0
    while True:
        N += 1
        grid_structure_interval = \
            dt*(sigma**2)/(((r_max-r_min)/float(N))**2)
        if grid_structure_interval > max_structure_const:
            break
    return N

• 添加vasicek_limits()方法来计算 Vasicek 利率过程的最小值和最大值，通过正态分布过程。Vasicek 模型下短期利率过程r(t)的期望值如下：

方差定义如下：

该方法返回一个元组，其中包括由正态分布过程的概率定义的最小和最大利率水平，写法如下：

def vasicek_limits(self, r0, sigma, kappa, theta, T, prob=1e-6):
    er = theta+(r0-theta)*math.exp(-kappa*T)
    variance = (sigma**2)*T if kappa==0 else \
                (sigma**2)/(2*kappa)*(1-math.exp(-2*kappa*T))
    stdev = math.sqrt(variance)
    r_min = self.norminv(prob, er, stdev)
    r_max = self.norminv(1-prob, er, stdev)
    return (r_min, r_max)

• 添加vasicek_diagonals()方法，该方法返回有限差分隐式方案的对角线，其中：

边界条件是使用诺伊曼边界条件实现的。该方法的写法如下：

def vasicek_diagonals(self, sigma, kappa, theta, r_min,
                      dr, N, dtau):
    rn = np.r_[0:N]*dr + r_min
    subdiagonals = kappa*(theta-rn)*dtau/(2*dr) - \
                    0.5*(sigma**2)*dtau/(dr**2)
    diagonals = 1 + rn*dtau + sigma**2*dtau/(dr**2)
    superdiagonals = -kappa*(theta-rn)*dtau/(2*dr) - \
                    0.5*(sigma**2)*dtau/(dr**2)

    # Implement boundary conditions.
    if N > 0:
        v_subd0 = subdiagonals[0]
        superdiagonals[0] = superdiagonals[0]-subdiagonals[0]
        diagonals[0] += 2*v_subd0
        subdiagonals[0] = 0

    if N > 1:
        v_superd_last = superdiagonals[-1]
        superdiagonals[-1] = superdiagonals[-1] - subdiagonals[-1]
        diagonals[-1] += 2*v_superd_last
        superdiagonals[-1] = 0

    return (subdiagonals, diagonals, superdiagonals)

诺伊曼边界条件指定了给定常规或偏微分方程的边界。更多信息可以在mathworld.wolfram.com/NeumannBoundaryConditions.html找到。

• 添加check_exercise()方法，返回一个布尔值列表，指示建议从提前行使中获得最佳回报的索引。该方法的写法如下：

def check_exercise(self, V, eex):
    return V > eex

• 添加exercise_call_price()方法，该方法返回折现值的行权价比率，写法如下：

def exercise_call_price(self, R, ratio, tau):
    K = ratio*np.exp(-R*tau)
    return K

• 添加vasicek_policy_diagonals()方法，该方法被政策迭代过程调用，用于更新一个迭代的子对角线、对角线和超对角线。在进行早期行权的索引中，子对角线和超对角线的值将被设置为 0，对角线上的剩余值。该方法返回新的子对角线、对角线和超对角线值的逗号分隔值。该方法的写法如下：

 def vasicek_policy_diagonals(self, subdiagonal, diagonal, \
                             superdiagonal, v_old, v_new, eex):
    has_early_exercise = self.check_exercise(v_new, eex)
    subdiagonal[has_early_exercise] = 0
    superdiagonal[has_early_exercise] = 0
    policy = v_old/eex
    policy_values = policy[has_early_exercise]
    diagonal[has_early_exercise] = policy_values
    return (subdiagonal, diagonal, superdiagonal)

• 添加iterate()方法，该方法通过执行政策迭代来实现有限差分的隐式方案，其中每个周期都涉及解决三对角方程组，调用vasicek_policy_diagonals()方法来更新三个对角线，并在没有更多早期行权机会时返回可赎回零息债券价格。它还返回执行的政策迭代次数。该方法的写法如下：

def iterate(self, subdiagonal, diagonal, superdiagonal,
            v_old, eex, max_policy_iter=10):
    v_mplus1 = v_old
    v_m = v_old
    change = np.zeros(len(v_old))
    prev_changes = np.zeros(len(v_old))

    iterations = 0
    while iterations <= max_policy_iter:
        iterations += 1

        v_mplus1 = self.tridiagonal_solve(
                subdiagonal, diagonal, superdiagonal, v_old)
        subdiagonal, diagonal, superdiagonal = \
            self.vasicek_policy_diagonals(
                subdiagonal, diagonal, superdiagonal, 
                v_old, v_mplus1, eex)

        is_eex = self.check_exercise(v_mplus1, eex)
        change[is_eex] = 1

        if iterations > 1:
            change[v_mplus1 != v_m] = 1

        is_no_more_eex = False if True in is_eex else True
        if is_no_more_eex:
            break

        v_mplus1[is_eex] = eex[is_eex]
        changes = (change == prev_changes)

        is_no_further_changes = all((x == 1) for x in changes)
        if is_no_further_changes:
            break

        prev_changes = change
        v_m = v_mplus1

    return v_mplus1, iterations-1

• 添加tridiagonal_solve()方法，该方法实现了 Thomas 算法来解决三对角方程组。方程组可以写成如下形式：

这个方程可以用矩阵形式表示：

这里，a是子对角线的列表，b是对角线的列表，c是矩阵的超对角线。

Thomas 算法是一个矩阵算法，用于使用简化的高斯消元法解决三对角方程组。更多信息可以在faculty.washington.edu/finlayso/ebook/algebraic/advanced/LUtri.htm找到。

tridiagonal_solve()方法的写法如下：

def tridiagonal_solve(self, a, b, c, d):
    nf = len(a)  # Number of equations
    ac, bc, cc, dc = map(np.array, (a, b, c, d))  # Copy the array
    for it in range(1, nf):
        mc = ac[it]/bc[it-1]
        bc[it] = bc[it] - mc*cc[it-1] 
        dc[it] = dc[it] - mc*dc[it-1]

    xc = ac
    xc[-1] = dc[-1]/bc[-1]

    for il in range(nf-2, -1, -1):
        xc[il] = (dc[il]-cc[il]*xc[il+1])/bc[il]

    del bc, cc, dc  # Delete variables from memory

    return xc

有了这些定义的方法，我们现在可以运行我们的代码，并使用 Vasicek 模型定价可赎回零息债券。

假设我们使用以下参数运行此模型：r0为0.05，R为0.05，ratio为0.95，sigma为0.03，kappa为0.15，theta为0.05，prob为1e-6，M为250，max_policy_iter为10，grid_struc_interval为0.25，我们对 0%到 2%之间的利率感兴趣。

以下 Python 代码演示了这个模型的 1 年、5 年、7 年、10 年和 20 年到期的情况：

In [ ]:
    r0 = 0.05
    R = 0.05
    ratio = 0.95
    sigma = 0.03
    kappa = 0.15
    theta = 0.05
    prob = 1e-6
    M = 250
    max_policy_iter=10
    grid_struct_interval = 0.25
    rs = np.r_[0.0:2.0:0.1]
In [ ]:
    vasicek = VasicekCZCB()
    r, vals = vasicek.vasicek_czcb_values(
        r0, R, ratio, 1., sigma, kappa, theta, 
        M, prob, max_policy_iter, grid_struct_interval, rs)
In [ ]:
    %pylab inline

    fig = plt.figure(figsize=(12, 8))
    plt.title("Callable Zero Coupon Bond Values by r")
    plt.plot(r, vals, label='1 yr')

    for T in [5., 7., 10., 20.]:
        r, vals = vasicek.vasicek_czcb_values(
            r0, R, ratio, T, sigma, kappa, theta, 
            M, prob, max_policy_iter, grid_struct_interval, rs)
        plt.plot(r, vals, label=str(T)+' yr', linestyle="--", marker=".")

    plt.ylabel("Value ($)")
    plt.xlabel("r")
    plt.legend()
    plt.grid(True)
    plt.show()

运行上述命令后，您应该得到以下输出：

我们得到了各种到期日和各种利率下可赎回零息债券的理论价值。

可赎回债券定价的其他考虑

在定价可赎回零息债券时，我们使用 Vasicek 利率过程来模拟利率的变动，借助正态分布过程。在Vasicek 模型部分，我们演示了 Vasicek 模型可以产生负利率，这在大多数经济周期中可能不太实际。定量分析师通常在衍生品定价中使用多个模型以获得现实的结果。CIR 和 Hull-White 模型是金融研究中常讨论的模型之一。这些模型的限制在于它们只涉及一个因素，或者说只有一个不确定性来源。

我们还研究了有限差分的隐式方案，用于早期行权的政策迭代。另一种考虑方法是有限差分的 Crank-Nicolson 方法。其他方法包括蒙特卡洛模拟来校准这个模型。

最后，我们得到了一份短期利率和可赎回债券价格的最终清单。为了推断特定短期利率的可赎回债券的公平价值，需要对债券价格清单进行插值。通常使用线性插值方法。其他考虑的插值方法包括三次和样条插值方法。

总结

在本章中，我们专注于使用 Python 进行利率和相关衍生品定价。大多数债券，如美国国债，每半年支付固定利息，而其他债券可能每季度或每年支付。债券的一个特点是它们的价格与当前利率水平密切相关，但是呈现出相反的关系。正常或正斜率的收益曲线，即长期利率高于短期利率，被称为向上倾斜。在某些经济条件下，收益曲线可能会倒挂，被称为向下倾斜。

零息债券是一种在其存续期内不支付利息的债券，只有在到期时偿还本金或面值时才支付。我们用 Python 实现了一个简单的零息债券计算器。

收益曲线可以通过零息债券、国债、票据和欧元存款的短期零点利率推导出来，使用引导过程。我们使用 Python 使用大量债券信息来绘制收益曲线，并从收益曲线中推导出远期利率、到期收益率和债券价格。

对债券交易员来说，两个重要的指标是久期和凸性。久期是债券价格对收益率变化的敏感度度量。凸性是债券久期对收益率变化的敏感度度量。我们在 Python 中实现了使用修正久期模型和凸性计算器进行计算。

短期利率模型经常用于评估利率衍生品。利率建模是一个相当复杂的话题，因为它们受到诸多因素的影响，如经济状态、政治决策、政府干预以及供求法则。已经提出了许多利率模型来解释利率的各种特征。我们讨论的一些利率模型包括 Vasicek、CIR 和 Rendleman 和 Bartter 模型。

债券发行人可能在债券中嵌入期权，以使他们有权利（但非义务）在规定的时间内以预定价格购买或出售发行的债券。可赎回债券的价格可以被视为不带期权的债券价格与嵌入式认购期权价格之间的价格差异。我们使用 Python 来通过有限差分的隐式方法来定价可赎回的零息债券，应用了 Vasicek 模型。然而，这种方法只是量化分析师在债券期权建模中使用的众多方法之一。

在下一章中，我们将讨论时间序列数据的统计分析。