Python 量化金融第二版（五）_quant-research-zh-CSDN博客

原文：zh.annas-archive.org/md5/25334ce178953792df9684c784953114

译者：飞龙

协议：CC BY-NC-SA 4.0

第十二章：蒙特卡洛模拟

蒙特卡洛模拟是金融领域中非常有用的工具。例如，因为我们可以通过从对数正态分布中抽取随机数来模拟股票价格，所以著名的布莱克-斯科尔斯-梅顿期权模型可以被复制。在第九章，投资组合理论中，我们学习到，通过向投资组合中添加更多股票，企业特定风险可以被减少或消除。通过模拟，我们可以更清晰地看到多样化效应，因为我们可以反复从 5000 只股票中随机选择 50 只股票。对于资本预算，我们可以模拟几十个具有不确定未来值的变量。对于这些情况，可以通过模拟生成许多可能的未来结果、事件以及各种类型的组合。在本章中，将涵盖以下主题：

从正态分布、均匀分布和泊松分布中生成随机数
使用蒙特卡洛模拟估算π值
使用对数正态分布模拟股票价格变动
构建高效投资组合与有效前沿
通过模拟复制布莱克-斯科尔斯-梅顿期权模型
定价一些复杂期权，例如带浮动行权价格的回顾期权
自助法（有/无放回）
长期预期回报预测
效率、准蒙特卡洛模拟与 Sobol 序列

蒙特卡洛模拟的重要性

蒙特卡洛模拟，或称模拟，在金融领域发挥着非常重要的作用，具有广泛的应用。假设我们打算估算一个项目的净现值（NPV）。未来存在许多不确定性，例如借款成本、最终产品价格、原材料等。对于少数几个变量，我们仍然能够轻松处理。然而，如果面对二十多个具有不确定未来值的变量，寻找解决方案就会成为一个难题。幸运的是，这时可以应用蒙特卡洛模拟。在第十章，期权与期货中，我们了解到，布莱克-斯科尔斯-梅顿期权模型背后的逻辑是股票回报的正态性假设。正因为如此，它们的封闭形式解可以通过模拟来复制。另一个例子是从 4500 只可用股票中随机选择 50 只。与传统期权（如布莱克-斯科尔斯-梅顿模型）不同，复杂期权没有封闭形式解。幸运的是，我们可以使用模拟来为其中一些期权定价。

从标准正态分布生成随机数

正态分布在金融中起着核心作用。一个主要原因是，许多金融理论（如期权理论及其相关应用）假设股票收益服从正态分布。第二个原因是，如果我们的计量经济学模型设计得当，那么模型中的误差项应该服从零均值的正态分布。生成标准正态分布中的 n 个随机数是一个常见任务。为此，我们有以下三行代码：

import scipy as sp
x=sp.random.standard_normal(size=10)
print(x)
[-0.98350472  0.93094376 -0.81167564 -1.83015626 -0.13873015  0.33408835
  0.48867499 -0.17809823  2.1223147   0.06119195]

SciPy/NumPy 中的基本随机数是通过 numpy.random 函数中的梅森旋转算法（Mersenne Twister PRNG）生成的。numpy.random 中的分布随机数是用 Cython/Pyrex 编写的，运行速度非常快。读者不可能得到与此处相同的10个随机数。我们很快会解释如何生成相同的一组随机数。或者，我们可以使用以下代码：

>>>import scipy as sp
>>>x=sp.random.normal(size=10)

这个程序等价于以下程序：

>>>import scipy as sp 
>>>x=sp.random.normal(0,1,10)

第一个输入是均值，第二个输入是标准差，最后一个输入是随机数的个数，也就是我们期望数据集的大小。对比前两个程序，显然均值和标准差的默认设置是0和1。我们可以使用 help() 函数来查看这三个输入变量的名称。为了节省空间，这里仅显示前几行：

>>>help(sp.random.normal) 
Help on built-in function normal:
normal(...) 
normal(loc=0.0, scale=1.0, size=None)

从正态分布中抽取随机样本

正态分布的概率密度函数最早由德·莫伊夫（De Moivre）推导出来，200 年后由高斯（Gauss）和拉普拉斯（Laplace）独立推导完成，通常称为钟形曲线，因为它的特征形状；参见下图：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_01.jpg

标准正态分布的密度函数如下所示：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_17.jpg

这里，f(x) 是标准正态分布的密度函数，x 是输入值，e 是指数函数，π 是 3.1415926。以下是生成上述钟形曲线的代码：

import scipy as sp
import scipy.stats as stats
import matplotlib.pyplot as plt
x = sp.arange(-3,3,0.01)
y=stats.norm.pdf(x)
plt.plot(x,y)
plt.title("A standard normal distribution")
plt.xlabel('x')
plt.ylabel('y')
plt.show()

使用种子生成随机数

很多时候，用户希望能够反复生成相同的一组随机数。例如，当教授在讲解如何估计一组随机数的均值、标准差、偏度和峰度时，学生能够生成与教授完全相同的数值是非常有帮助的。另一个例子是，当我们在调试 Python 程序以模拟股票走势时，我们可能更希望得到相同的中间结果。对于这种情况，我们可以使用scipy.random.seed()函数，如下所示：

>>>import scipy as sp 
>>>sp.random.seed(12345) 
>>>x=sp.random.normal(0,1,20) 
>>>print x[0:5] 
[-0.20470766 0.47894334 -0.51943872 -0.5557303 1.96578057] 
>>>

这里，12345是种子。种子的值并不重要，关键是相同的种子会产生相同的随机数值。更一般的正态分布公式如下：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_18.jpg

这里，*f(x)*是正态分布的密度函数，x是输入值，e是指数函数，μ是均值，σ是标准差。

正态分布的随机数

要从正态分布中生成n个随机数，我们有以下代码：

>>>impimport scipy as sp 
>>>sp.random.seed(12345) 
>>>mean=0.05
>>>std=0.1
>>>n=50
>>>x=sp.random.normal(mean,std,n) 
>>>print(x[0:5])
[ 0.02952923 0.09789433 -0.00194387 -0.00557303 0.24657806]
>>>

这个程序与前一个程序的区别在于，均值是0.05而不是0，而标准差是0.1而不是1。

正态分布的直方图

在分析数据集属性的过程中，直方图被广泛使用。为了为从具有指定均值和标准差的正态分布中抽取的一组随机值生成直方图，我们有以下代码：

import scipy as sp 
import matplotlib.pyplot as plt 
sp.random.seed(12345) 
mean=0.1
std=0.2
n=1000
x=sp.random.normal(mean,std,n) 
plt.hist(x, 15, normed=True) 
plt.title("Histogram for random numbers drawn from a normal distribution")
plt.annotate("mean="+str(mean),xy=(0.6,1.5))
plt.annotate("std="+str(std),xy=(0.6,1.4))
plt.show()

结果图形如下所示：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_02.jpg

对数正态分布的图形展示

当股票回报率遵循正态分布时，其价格应当遵循对数正态分布。对数正态分布的定义如下：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_19.jpg

这里，*f(x;μ,σ)*是对数正态分布的密度，ln()是自然对数函数。以下代码展示了三个不同的对数正态分布，分别使用了三组参数，如(0, 0.25)、(0, 0.5)和(0, 1.0)。第一个参数是均值(μ)，第二个参数是标准差，见以下代码：

import scipy as sp
import numpy as np
import matplotlib.pyplot as plt 
from scipy import sqrt,exp,log,pi
#
x=np.linspace(0.001,3,200)
mu=0 
sigma0=[0.25,0.5,1]
color=['blue','red','green'] 
target=[(1.2,1.3),(1.7,0.4),(0.18,0.7)]
start=[(1.8,1.4),(1.9,0.6),(0.18,1.6)]
#
for i in sp.arange(len(sigma0)):
    sigma=sigma0[i]
    y=1/(x*sigma*sqrt(2*pi))*exp(-(log(x)-mu)**2/(2*sigma*sigma))
    plt.annotate('mu='+str(mu)+', sigma='+str(sigma),xy=target[i],xytext=start[i],arrowprops=dict(facecolor=color[i],shrink=0.01),) 
    plt.plot(x,y,color[i])
    plt.title('Lognormal distribution') 
    plt.xlabel('x')
    plt.ylabel('lognormal density distribution') 
#
plt.show()

这里展示了图形。显然，与正态分布的密度不同，对数正态分布的密度函数是非对称的：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_03.jpg

从均匀分布中生成随机数

当从n个可用的股票中随机选择 m 只股票时，我们可以从均匀分布中抽取一组随机数。为了从均匀分布中生成 10 个介于 1 和 100 之间的随机数，我们有以下代码。为了保证相同的数值集合，使用了seed()函数：

>>>import scipy as sp 
>>>sp.random.seed(123345) 
>>>x=sp.random.uniform(low=1,high=100,size=10)

同样，low、high 和 size 是三个输入名称。第一个指定最小值，第二个指定最大值，而 size 指定我们打算生成的随机数的数量。前五个数字如下所示：

>>>print(x[0:5])
[ 30.32749021 20.58006409 2.43703988 76.15661293 75.06929084]
>>>

下一个程序随机掷一个骰子，结果为 1、2、3、4、5 或 6：

import random
def rollDice():
    roll = random.randint(1,6)
    return roll
i =1
n=10
result=[]
random.seed(123)
while i<n:
    result.append(rollDice())
    i+=1
print(result)
[1, 1, 3, 1, 6, 1, 4, 2, 6]

在前一个程序中，应用了random.seed()函数。因此，任何读者都应得到最后一行显示的相同结果。

使用模拟估算π值

通过模拟估算π值是一个很好的练习。我们来画一个边长为 2R 的正方形。如果在正方形内部放置一个最大的圆形，那么其半径将是 R，表示为以下方程：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_20.jpg

另一方面，正方形的面积是其边长的平方：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_21.jpg

将方程 (4) 除以方程 (5)，我们得到以下结果：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_22.jpg

重新组织后，我们得到以下方程：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_23.jpg

换句话说，π的值将是4 Scircle/Square*。在运行仿真时，我们从一个均匀分布中生成n对x和y，范围为 0 到 0.5。然后我们估算一个距离，该距离是x和y的平方和的平方根，即！使用仿真估算π值。

显然，当d小于 0.5（R 值）时，它会落入圆内。我们可以想象扔一个飞镖，飞镖落入圆内时，π的值将呈现以下形式：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_24.jpg

下图展示了这些随机点在圆内和方形内的分布：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_04.jpg

估算π值的 Python 程序如下所示：

import scipy as sp 
n=100000
x=sp.random.uniform(low=0,high=1,size=n) 
y=sp.random.uniform(low=0,high=1,size=n) 
dist=sp.sqrt(x**2+y**2) 
in_circle=dist[dist<=1] 
our_pi=len(in_circle)*4./n
print ('pi=',our_pi)
print('error (%)=', (our_pi-sp.pi)/sp.pi)

每次运行前面的代码时，估算的π值会发生变化，如下所示，且其估算的准确性取决于试验次数，即n：

('pi=', 3.14168)
('error (%)=', 2.7803225891524895e-05)

从泊松分布生成随机数

为了研究私人信息的影响，Easley、Kiefer、O’Hara 和 Paperman（1996 年）设计了一个信息化交易概率（PIN）度量方法，基于买方发起交易的每日数量和卖方发起交易的数量推导而来。他们模型的基本假设是，订单到达遵循泊松分布。以下代码展示了如何从泊松分布中生成n个随机数：

import numpy as np
import scipy as sp 
import matplotlib.pyplot as plt 
x=sp.random.poisson(lam=1, size=100) 
#plt.plot(x,'o') 
a = 5\. # shape 
n = 1000 
s = np.random.power(a, n) 
count, bins, ignored = plt.hist(s, bins=30) 
x = np.linspace(0, 1, 100) 
y = a*x**(a-1.) 
normed_y = n*np.diff(bins)[0]*y 
plt.title("Poisson distribution")
plt.ylabel("y")
plt.xlabel("x")
plt.plot(x, normed_y) 
plt.show()

该图如下所示：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_05.jpg

从 n 个给定股票中随机选择 m 只股票

基于前面的程序，我们可以轻松地从 500 只可用证券中选择 20 只股票。如果我们打算研究随机选择股票数量对投资组合波动性的影响，这将是一个重要的步骤，如下所示的代码所示：

import scipy as sp 
n_stocks_available=500 
n_stocks=20 
sp.random.seed(123345) 
x=sp.random.uniform(low=1,high=n_stocks_available,size=n_stocks)
y=[] 
for i in range(n_stocks): 
    y.append(int(x[i])) 
#print y 
final=sp.unique(y) 
print(final) 
print(len(final))
[  8  31  61  99 124 148 155 172 185 205 226 275 301 334 356 360 374 379
 401 449]
20

在前面的程序中，我们从 500 个数字中选择 20 个数字。由于我们必须选择整数，可能最终会得到少于 20 个值，也就是说，一些整数在将实数转换为整数后可能会出现重复。一种解决方法是选择更多的数字，然后取前 20 个整数。另一种方法是使用randrange()和randint()函数。在下一个程序中，我们从所有可用的股票中选择n只股票。首先，我们从canisius.edu/~yany/python/yanMonthly.pkl下载数据集。假设数据集位于C:/temp/目录下：

import scipy as sp
import numpy as np
import pandas as pd
#
n_stocks=10 
x=pd.read_pickle('c:/temp/yanMonthly.pkl') 
x2=sp.unique(np.array(x.index)) 
x3=x2[x2<'ZZZZ']                        # remove all indices 
sp.random.seed(1234567) 
nonStocks=['GOLDPRICE','HML','SMB','Mkt_Rf','Rf','Russ3000E_D','US_DEBT','Russ3000E_X','US_GDP2009dollar','US_GDP2013dollar'] 
x4=list(x3) 
#
for i in range(len(nonStocks)): 
    x4.remove(nonStocks[i]) 
#
k=sp.random.uniform(low=1,high=len(x4),size=n_stocks) 
y,s=[],[] 
for i in range(n_stocks): 
    index=int(k[i]) 
    y.append(index) 
    s.append(x4[index]) 
#
final=sp.unique(y) 
print(final) 
print(s)

在前面的程序中，我们移除了非股票数据项。这些非股票项是数据项的一部分。首先，我们加载一个名为yanMonthly.pickle的数据集，其中包括 200 多只股票、黄金价格、GDP、失业率、小盘减大盘（SMB）、高估减低估（HML）、无风险利率、价格率、市场超额收益率以及罗素指数。

pandas 的一种输出格式是.pkl和.png。由于x.index会呈现每个观测值的所有索引，我们需要使用unique()函数来选择所有唯一的 ID。因为我们只考虑股票来构建我们的投资组合，所以我们必须移除所有市场指数和其他非股票证券，如 HML 和US_DEBT。因为所有股票市场指数都以插入符号（^）开头，所以我们使用小于 ZZZZ 的方式来移除它们。对于其他在 A 和 Z 之间的 ID，我们必须逐个移除它们。为此，我们使用.remove()函数，该函数适用于列表变量。最终输出如下：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_06.jpg

有/没有替换

假设我们拥有一只股票的历史数据，如价格和回报。显然，我们可以估计它们的均值、标准差和其他相关统计数据。那么，明年的预期年均值和风险是多少呢？最简单的，可能也是天真的方法是使用历史均值和标准差。更好的方法是构建年回报和风险的分布。这意味着我们必须找到一种方法，利用历史数据更有效地预测未来。在这种情况下，我们可以应用自助法（bootstrapping）方法。例如，对于一只股票，我们拥有过去 20 年的月度回报数据，即 240 个观测值。

为了估计明年 12 个月的回报，我们需要构建回报分布。首先，我们从历史回报集中随机选择 12 个回报（不放回），并估计它们的均值和标准差。我们重复这个过程 5000 次。最终输出将是我们的回报-标准差分布。基于这个分布，我们也可以估计其他属性。同样地，我们也可以在有替换的情况下进行。NumPy 中一个有用的函数是numpy.random.permutation()。假设我们有从 1 到 10 的 10 个数字（包含 1 和 10）。我们可以调用numpy.random.permutation()函数来重新洗牌，如下所示：

import numpy as np 
x=range(1,11) 
print(x) 
for i in range(5):
    y=np.random.permutation(x) 
#
print(y)

这段代码的输出如下所示：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_07.jpg

基于numpy.random.permutation()函数，我们可以定义一个函数，输入三个变量：数据、我们计划从数据中随机选择的观测值数量，以及是否选择有或没有替换的自助法，如下代码所示：

import numpy as np 
def boots_f(data,n_obs,replacement=None):
    n=len(data) 
    if (n<n_obs):
        print "n is less than n_obs" 
    else: 
        if replacement==None:
            y=np.random.permutation(data) 
            return y[0:n_obs] 
        else:
            y=[] 
    #
    for i in range(n_obs): 
        k=np.random.permutation(data) 
        y.append(k[0]) 
    return y

在之前的程序中指定的约束条件是，给定的观察次数应大于我们计划选择的随机回报次数。这对于无替换的自助法是成立的。对于有替换的自助法，我们可以放宽这个约束；请参考相关习题。

年度回报分布

估算年化回报分布并将其表示为图形是一个很好的应用。为了使我们的练习更有意义，我们下载了微软的每日价格数据。然后，我们估算了其每日回报，并将其转换为年回报。基于这些年回报，我们通过应用带替换的自助法进行 5,000 次模拟，从而生成其分布，如下代码所示：

import numpy as np 
import scipy as sp
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.finance import quotes_historical_yahoo_ochl as getData 
# Step 1: input area
ticker='MSFT'          # input value 1 
begdate=(1926,1,1)      # input value 2 
enddate=(2013,12,31)    # input value 3 
n_simulation=5000       # input value 4
# Step 2: retrieve price data and estimate log returns
x=getData(ticker,begdate,enddate,asobject=True)
logret = sp.log(x.aclose[1:]/x.aclose[:-1])
# Step 3: estimate annual returns 
date=[]
d0=x.date
for i in range(0,sp.size(logret)): 
    date.append(d0[i].strftime("%Y"))
y=pd.DataFrame(logret,date,columns=['logret']) 
ret_annual=sp.exp(y.groupby(y.index).sum())-1 
ret_annual.columns=['ret_annual']
n_obs=len(ret_annual)
# Step 4: estimate distribution with replacement 
sp.random.seed(123577) 
final=sp.zeros(n_obs,dtype=float)
for i in range(0,n_obs):
    x=sp.random.uniform(low=0,high=n_obs,size=n_obs) 
    y=[]
    for j in range(n_obs): 
        y.append(int(x[j]))
        z=np.array(ret_annual)[y] 
    final[i]=sp.mean(z)
# step 5: graph
plt.title('Mean return distribution: number of simulations ='+str(n_simulation))
plt.xlabel('Mean return')
plt.ylabel('Frequency')
mean_annual=round(np.mean(np.array(ret_annual)),4) 
plt.figtext(0.63,0.8,'mean annual='+str(mean_annual)) 
plt.hist(final, 50, normed=True)
plt.show()

相应的图表如下所示：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_08.jpg

股票价格波动模拟

我们在前面的章节中提到，在金融领域，回报假定服从正态分布，而价格则服从对数正态分布。股票在时间 t+1 的价格是时间 t 股票价格、均值、标准差和时间间隔的函数，如下公式所示：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_25.jpg

在此公式中，St + 1 是时间 t+1 的股票价格，ˆ μ 是预期股票回报，t_ 是时间间隔（T t n_= ），T 是时间（以年为单位），n 是步数，ε 是均值为零的分布项，σ 是标的股票的波动率。通过简单的操作，方程（4）可以推导出以下我们将在程序中使用的方程：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_26.jpg

在一个风险中性世界中，投资者不要求承担风险的补偿。换句话说，在这样的世界里，任何证券（投资）的预期回报率都是无风险利率。因此，在风险中性世界中，前面的方程式变为以下方程式：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_27.jpg

如果你想了解更多关于风险中性概率的内容，请参考*《期权、期货及其他衍生品》第七版，约翰·赫尔，皮尔森，2009 年*。模拟股票价格运动（路径）的 Python 代码如下：

import scipy as sp 
import matplotlib.pyplot as plt
# input area
stock_price_today = 9.15 # stock price at time zero 
T =1\.                    # maturity date (in years) 
n_steps=100\.             # number of steps 
mu =0.15                 # expected annual return 
sigma = 0.2              # annualized volatility
sp.random.seed(12345)    # fixed our seed 
n_simulation = 5         # number of simulations 
dt =T/n_steps 
#
S = sp.zeros([n_steps], dtype=float) 
x = range(0, int(n_steps), 1) 
for j in range(0, n_simulation): 
    S[0]= stock_price_today 
    for i in x[:-1]: 
        e=sp.random.normal() 
        S[i+1]=S[i]+S[i]*(mu-0.5*pow(sigma,2))*dt+sigma*S[i]*sp.sqrt(dt)*e; 
    plt.plot(x, S)
#
plt.figtext(0.2,0.8,'S0='+str(S[0])+',mu='+str(mu)+',sigma='+str(sigma)) 
plt.figtext(0.2,0.76,'T='+str(T)+', steps='+str(int(n_steps))) 
plt.title('Stock price (number of simulations = %d ' % n_simulation +')') 
plt.xlabel('Total number of steps ='+str(int(n_steps))) 
plt.ylabel('stock price') 
plt.show()

为了使我们的图表更具可读性，我们故意只选择了五次模拟。由于应用了scipy.random.seed()函数，你可以通过运行之前的代码来复制以下图表。图表如下所示：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_09.jpg

期权到期日股票价格的图形展示

到目前为止，我们已经讨论了期权实际上是路径无关的，这意味着期权价格取决于终值。因此，在定价此类期权之前，我们需要知道终期股票价格。为了扩展之前的程序，我们有以下代码来估算一组给定值的终期股票价格：S0（初始股票价格），n_simulation（终期价格的数量），T（到期日，按年计算），n_steps（步骤数），mu（预期年股票回报率），sigma（波动率）：

import scipy as sp 
import matplotlib.pyplot as plt
from scipy import zeros, sqrt, shape 
#input area
S0 = 9.15               # stock price at time zero 
T =1\.                   # years
n_steps=100\.            # number of steps 
mu =0.15                # expected annual return 
sigma = 0.2             # volatility (annual) 
sp.random.seed(12345)   # fix those random numbers 
n_simulation = 1000     # number of simulation 
dt =T/n_steps 
#
S = zeros([n_simulation], dtype=float) 
x = range(0, int(n_steps), 1) 
for j in range(0, n_simulation): 
    tt=S0 
    for i in x[:-1]: 
        e=sp.random.normal() 
        tt+=tt*(mu-0.5*pow(sigma,2))*dt+sigma*tt*sqrt(dt)*e; 
        S[j]=tt 
#
plt.title('Histogram of terminal price') 
plt.ylabel('Number of frequencies') 
plt.xlabel('Terminal price') 
plt.figtext(0.5,0.8,'S0='+str(S0)+',mu='+str(mu)+',sigma='+str(sigma)) 
plt.figtext(0.5,0.76,'T='+str(T)+', steps='+str(int(n_steps))) 
plt.figtext(0.5,0.72,'Number of terminal prices='+str(int(n_simulation))) 
plt.hist(S) 
plt.show()

我们模拟的终期价格的直方图如下所示：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_10.jpg

正如我们在第九章中提到的，投资组合理论，为了生成两个相关的随机时间序列，需要进行两个步骤：生成两个零相关的随机时间序列x1和x2；然后应用以下公式：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_33.jpg

这里，ρ是这两个时间序列之间的预定相关性。现在，y1和y2与预定的相关性相关。以下 Python 程序将实现上述方法：

import scipy as sp
sp.random.seed(123)
n=1000
rho=0.3
x1=sp.random.normal(size=n)
x2=sp.random.normal(size=n)
y1=x1
y2=rho*x1+sp.sqrt(1-rho**2)*x2
print(sp.corrcoef(y1,y2))
[[ 1\.          0.28505213]
 [ 0.28505213  1\.        ]]

使用模拟复制布莱克-斯科尔斯-默顿看涨期权

在知道终期价格后，如果给定行权价格，我们可以估算看涨期权的支付。使用无风险利率作为贴现率计算的这些贴现支付的均值将是我们的看涨期权价格。以下代码帮助我们估算看涨期权价格：

import scipy as sp 
from scipy import zeros, sqrt, shape 
#
S0 = 40\.              # stock price at time zero 
X= 40\.                # exercise price 
T =0.5                # years 
r =0.05               # risk-free rate 
sigma = 0.2           # annualized volatility 
n_steps=100          # number of steps 
#
sp.random.seed(12345) # fix those random numbers 
n_simulation = 5000   # number of simulation 
dt =T/n_steps 
call = sp.zeros([n_simulation], dtype=float) 
x = range(0, int(n_steps), 1) 
for j in range(0, n_simulation): 
    sT=S0 
    for i in x[:-1]: 
        e=sp.random.normal() 
        sT*=sp.exp((r-0.5*sigma*sigma)*dt+sigma*e*sqrt(dt)) 
        call[j]=max(sT-X,0) 
#
call_price=sp.mean(call)*sp.exp(-r*T) 
print('call price = ', round(call_price,3))

估算的看涨期权价格为$2.748。相同的逻辑适用于定价看跌期权。

异型期权 #1 – 使用蒙特卡洛模拟来定价平均值

到目前为止，我们在第九章中已经讨论了欧式期权和美式期权，投资组合理论。布莱克-斯科尔斯-默顿期权模型，也叫做普通期权。其特征之一是路径无关。另一方面，异型期权更复杂，因为它们可能有多个触发因素与其支付的确定相关。例如，一家炼油厂担心未来三个月内原油价格的波动。它们计划对潜在的原油价格跳跃进行对冲。公司可以购买看涨期权。然而，由于公司每天消耗大量的原油，自然更关心的是平均价格，而不仅仅是普通看涨期权所依赖的终期价格。对于这种情况，平均期权会更有效。平均期权是一种亚洲期权。对于平均期权，其支付是由一段预设时间内的基础价格的平均值决定的。平均值有两种类型：算术平均和几何平均。亚洲看涨期权（平均价格）的支付函数如下所示：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_29.jpg

这里给出了亚洲看跌期权（平均价格）的支付函数：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_30.jpg

亚洲期权是另类期权的一种基本形式。亚洲期权的另一个优势是，与欧洲和美式普通期权相比，它们的成本更低，因为平均值的波动要比终值的波动小得多。下面的 Python 程序是针对具有算术平均价格的亚洲期权的：

import scipy as sp
s0=40\.                 # today stock price 
x=40\.                  # exercise price 
T=0.5                  # maturity in years 
r=0.05                 # risk-free rate 
sigma=0.2              # volatility (annualized) 
sp.random.seed(123)    # fix a seed here 
n_simulation=100       # number of simulations 
n_steps=100\.           # number of steps
#	
dt=T/n_steps 
call=sp.zeros([n_simulation], dtype=float) 
for j in range(0, n_simulation): 
    sT=s0 
    total=0 
    for i in range(0,int(n_steps)): 
         e=sp.random.normal()
         sT*=sp.exp((r-0.5*sigma*sigma)*dt+sigma*e*sp.sqrt(dt)) 
         total+=sT 
         price_average=total/n_steps 
    call[j]=max(price_average-x,0) 
#
call_price=sp.mean(call)*sp.exp(-r*T) 
print('call price based on average price = ', round(call_price,3))
('call price based on average price = ', 1.699)

根据上述结果，这个平均价格看涨期权的权利金为$1.70。

另类期权 #2 – 使用蒙特卡罗模拟法定价障碍期权

与黑-肖尔斯-莫顿期权模型中的看涨和看跌期权不同，后者是路径无关的，障碍期权则是路径相关的。障碍期权在许多方面与普通期权类似，但存在触发条件。敲入期权的生命周期从无价值开始，除非标的股票达到预定的敲入障碍。而相反，敲出障碍期权从一开始就有效，只有在价格突破敲出障碍时才会变得无效。此外，如果障碍期权到期时没有激活，它可能变得一文不值，或者可能支付一部分权利金作为现金返还。四种障碍期权如下所示：

向上敲出期权：在这种障碍期权中，价格从较低的障碍水平开始。如果它达到障碍，则被“敲出”。
跌破障碍期权：在这种障碍期权中，价格从一个较高的障碍开始。如果价格达到障碍，则被“敲出”。
向上敲入期权：在这种障碍期权中，价格从较低的障碍开始，必须达到该障碍才能激活。
向下敲入期权：在这种障碍期权中，价格从较高的障碍开始，必须达到该障碍才能激活。

接下来的 Python 程序用于一个向上敲出障碍期权，类型为欧洲看涨期权：

import scipy as sp 
from scipy import log,exp,sqrt,stats 
#
def bsCall(S,X,T,r,sigma):
    d1=(log(S/X)+(r+sigma*sigma/2.)*T)/(sigma*sqrt(T)) 
    d2 = d1-sigma*sqrt(T)
    return S*stats.norm.cdf(d1)-X*exp(-r*T)*stats.norm.cdf(d2)
#
def up_and_out_call(s0,x,T,r,sigma,n_simulation,barrier):
    n_steps=100\. 
    dt=T/n_steps 
    total=0 
    for j in sp.arange(0, n_simulation): 
        sT=s0 
        out=False
        for i in range(0,int(n_steps)): 
            e=sp.random.normal() 
            sT*=sp.exp((r-0.5*sigma*sigma)*dt+sigma*e*sp.sqrt(dt)) 
            if sT>barrier: 
               out=True 
        if out==False: 
            total+=bsCall(s0,x,T,r,sigma) 
    return total/n_simulation

基本设计是我们模拟股票价格波动n次，例如模拟 100 次。对于每次模拟，我们有 100 个步骤。每当股票价格达到障碍时，支付将为零。否则，支付将为一种普通的欧洲看涨期权。最终价值将是所有未被敲出的看涨期权价格的总和，再除以模拟次数，如下所示的代码：

s0=40\.              # today stock price 
x=40\.               # exercise price 
barrier=42          # barrier level 
T=0.5               # maturity in years 
r=0.05              # risk-free rate 
sigma=0.2           # volatility (annualized) 
n_simulation=100    # number of simulations 
sp.random.seed(12)  # fix a seed
#
result=up_and_out_call(s0,x,T,r,sigma,n_simulation,barrier) 
print('up-and-out-call = ', round(result,3))
('up-and-out-call = ', 0.937)

根据上述结果，我们知道这个向上敲出看涨期权的价格为$0.94。

喜欢使用模拟法计算 VaR 的两种方法

在上一章，第十一章，风险价值中，我们学习了可以使用两种方法来估算个股或投资组合的 VaR：这取决于正态性假设和基于历史回报排序的方法。蒙特卡罗模拟法能够将这两种方法结合起来，见下述代码：

import numpy as np
import numpy as np
import scipy as sp
import pandas as pd
from scipy.stats import norm
#
position=1e6              # portfolio value
std=0.2                   # volatility
mean=0.08                 # mean return
confidence=0.99           # confidence level
nSimulations=50000        # number of simulations
# Method I
z=norm.ppf(1-confidence)
VaR=position*(mean+z*std)
print("Holding=",position, "VaR=", round(VaR,2), "tomorrow")
#
# Method II: Monte Carlo simulaiton 
sp.random.seed(12345) 
ret2=sp.random.normal(mean,std,nSimulations) 
ret3=np.sort(ret2) 
m=int(nSimulations*(1-confidence))
VaR2=position*(ret3[m])
print("Holding=",position, "VaR2=", round(VaR2,2), "tomorrow")
('Holding=', 1000000.0, 'VaR=', -385270.0, 'tomorrow')
('Holding=', 1000000.0, 'VaR2=', -386113.0, 'tomorrow')

蒙特卡洛模拟的结果为$386,113，而基于公式计算的结果为$385,270（假设投资组合当前价值为 100 万美元）。

使用蒙特卡洛模拟进行资本预算

正如我们在本章开头提到的，当变量的数量有许多不同的值时，我们可以使用蒙特卡洛模拟来进行资本预算。我们的目标是通过对所有未来自由现金流进行折现，来估算给定预算的净现值（NPV）：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_31.jpg

这里，NPV 是某一提案的净现值，FCF0 将是零时点的自由现金流，FCFt 将是第 I 年末的自由现金流，R 是折现率。计算第 t 年末自由现金流的公式如下：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_32.jpg

这里，FCTt 是第 t 年的自由现金流，Dt 是第 t 年的折旧，CaptExt 是第 t 年的净资本支出，NWC 是净营运资金，即流动资产减去流动负债，Δ表示变化。让我们看一个简单的例子。假设公司购买一项长期等效资产，总成本为 50 万美元，使用年限为五年：

项目	0	1	2	3	4	5
价格	0	28	28	28	28	28
单位	0	100000	100000	100000	100000	100000
销售	0	2800000	2800000	2800000	2800000	2800000
销售成本	0	840000	840000	840000	840000	840000
其他成本	0	100000	100000	100000	100000	100000
销售、一般管理和行政费用	15000	15000	15000	15000	15000	15000
研发费用	20000
折旧		1000000	1000000	1000000	1000000	1000000
息税前利润（EBIT）	-35000	845000	845000	845000	845000	845000
税率 35%	-12250	295750	295750	295750	295750	295750
净收入（NI）	-47250	1140750	1140750	1140750	1140750	1140750
加上折旧	-47250	2140750	2140750	2140750	2140750	2140750

表 12.1 每年的现金流

我们有以下等效代码：

import scipy as sp
nYear=5                 # number of years
costEquipment=5e6       # 5 million 
n=nYear+1               # add year zero
price=28                # price of the product
units=100000            # estimate number of units sold 
otherCost=100000        # other costs
sellingCost=1500        # selling and administration cost 
R_and_D=200000          # Research and development
costRawMaterials=0.3    # percentage cost of raw materials
R=0.15                  # discount rate
tax=0.38                # corporate tax rate
#
sales=sp.ones(n)*price*units
sales[0]=0              # sales for 1st year is zero
cost1=costRawMaterials*sales
cost2=sp.ones(n)*otherCost
cost3=sp.ones(n)*sellingCost
cost4=sp.zeros(n)
cost4[0]=costEquipment
RD=sp.zeros(n)
RD[0]=R_and_D                     # assume R&D at time zero
D=sp.ones(n)*costEquipment/nYear  # straight line depreciation 
D[0]=0                            # no depreciation at time 0
EBIT=sales-cost1-cost2-cost3-cost4-RD-D
NI=EBIT*(1-tax)
FCF=NI+D                         # add back depreciation
npvProject=sp.npv(R,FCF)         # estimate NPV
print("NPV of project=",round(npvProject,0))
('NPV of project=', 1849477.0)

该项目的净现值（NPV）为$1,848,477。由于它是正值，如果我们的标准是基于 NPV 规则，那么我们应当接受该提案。现在，让我们加入一些不确定性。假设我们有三个不确定因素：价格、预计销售的产品单位数以及折现率，请参见以下代码：

import scipy as sp
import matplotlib.pyplot as plt
nYear=5                 # number of years
costEquipment=5e6       # 5 million 
n=nYear+1               # add year zero
otherCost=100000        # other costs
sellingCost=1500        # selling and administration cost 
R_and_D=200000          # Research and development
costRawMaterials=0.3    # percentage cost of raw materials
tax=0.38                # corporate tax rate
thousand=1e3            # unit of thousand 
million=1e6             # unit of million 
#
# three uncertainties: price, unit and discount rate
nSimulation=100         # number of simulation
lowPrice=10             # low price
highPrice=30            # high price
lowUnit=50*thousand     # low units expected to sell 
highUnit=200*thousand   # high units expected to sell 
lowRate=0.15            # lower discount rate
highRate=0.25           # high discount rate 
#
n2=nSimulation
sp.random.seed(123)
price0=sp.random.uniform(low=lowPrice,high=highPrice,size=n2)
units0=sp.random.uniform(low=lowUnit,high=highUnit,size=n2)
R0=sp.random.uniform(lowRate,highRate,size=n2)
#
npv=[]
for i in sp.arange(nSimulation):
    units=sp.ones(n)*units0[i]
    price=price0[i]
    R=R0[i]
    sales=units*price
    sales[0]=0              # sales for 1st year is zero
    cost1=costRawMaterials*sales
    cost2=sp.ones(n)*otherCost
    cost3=sp.ones(n)*sellingCost
    cost4=sp.zeros(n)
    cost4[0]=costEquipment
    RD=sp.zeros(n)
    RD[0]=R_and_D                     # assume R&D at time zero
    D=sp.ones(n)*costEquipment/nYear  # straight line depreciation 
    D[0]=0                            # no depreciation at time 0
    EBIT=sales-cost1-cost2-cost3-cost4-RD-D
    NI=EBIT*(1-tax)
    FCF=NI+D                          # add back depreciation
    npvProject=sp.npv(R,FCF)/million  # estimate NPV
    npv.append(npvProject)
print("mean NPV of project=",round(sp.mean(npv),0))
print("min  NPV of project=",round(min(npv),0))
print("max  NPV of project=",round(max(npv),0))
plt.title("NPV of the project: 3 uncertainties")
plt.xlabel("NPV (in million)")
plt.hist(npv, 50, range=[-3, 6], facecolor='blue', align='mid')
plt.show()

NPV 分布的直方图如下所示：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_11.jpg

Python SimPy 模块

SimPy 是一个基于标准 Python 的基于过程的离散事件仿真框架。它的事件调度器基于 Python 的生成器，也可以用于异步网络通信或实现多智能体系统（包括模拟和真实通信）。SimPy 中的进程是简单的 Python 生成器函数，用于建模活跃的组件，如顾客、车辆或代理。SimPy 还提供了各种类型的共享资源，用于建模有限容量的拥堵点（如服务器、结账柜台和隧道）。从 3.1 版本开始，它还将提供监控功能，帮助收集有关资源和进程的统计数据：

import simpy
def clock(env, name, tick):
     while True:
         print(name, env.now)
         yield env.timeout(tick)
#
env = simpy.Environment()
env.process(clock(env, 'fast', 0.5))
env.process(clock(env, 'slow', 1))
env.run(until=2)
('fast', 0)
('slow', 0)
('fast', 0.5)
('slow', 1)
('fast', 1.0)
('fast', 1.5)

两种社会政策的比较——基本收入与基本工作

这个例子借鉴自 Stucchhio（2013）。在过去几十年的发展中，各国的财富持续积累，尤其是在发达国家中尤为明显。支持公平的基本论点之一是每个公民都应拥有基本的生活标准。基于这一论点，许多国家为其公民提供了大量福利，例如普及医疗、免费教育等。一项政策建议是基本收入，即每个公民每年都会获得一笔无附加条件的基本收入。例如，如果我们假设基本时薪为 7.50 美元，每周工作 40 小时，每年工作 50 周，那么基本收入应该为 15,000 美元。Zhong（2017）报告称，印度正在考虑通过普遍基本收入计划来对抗贫困。显而易见的优势是行政成本会相当低。此外，腐败侵占政府为穷人发放的资金的可能性也较小。2017 年，芬兰启动了一个试点项目，加拿大和荷兰的地方当局也宣布了相关实验。2016 年，瑞士选民拒绝了一项最低收入提案。

一种替代方案是所谓的基本工作，在这种工作中，政府保证为无法找到体面工作的人提供一份低薪工作。这些方法各有优缺点。基于一组假设，如时薪、每周工作小时数、每年工作周数、人口、劳动力等，Stucchhio（2013）比较了这两种提案的成本和收益。存在若干不确定性；请参阅下表中的列表：

政策	命令	描述
基本收入	`unitAdmCost = norm(250,75)`	每个人的行政费用
`binom(nNonWorkers,tiny).rvs()`	从二项分布中随机生成一个数字
`nonWorkerMultiplier = uniform(-0.10, 0.15).rvs()`	非劳动者的乘数
基本工作	`unitAdmCost4disabled= norm(500,150).rvs()`	每个残疾成年人所需的行政费用
`unitAdmCost4worker = norm(5000, 1500).rvs()`	每个工人的行政费用
`nonWorkerMultiplier = uniform(-0.20, 0.25).rvs()`	非工作者的乘数
`hourlyProductivity = uniform(0.0,hourlyPay).rvs()`	每小时生产力

表 12.2：两种提案的成本与收益

程序使用三种分布：正态分布、均匀分布和二项分布。uniform(a,b).rvs()命令生成一个在a和b之间均匀分布的随机数。norm(mean,std).rvs()命令生成一个来自具有指定均值和标准差的正态分布的随机数。binom(n,k).rvs()命令生成一个来自二项分布的随机数，输入值为n和k：

import scipy as sp
import scipy.stats as stats
sp.random.seed(123)
u=stats.uniform(-1,1).rvs()
n=stats.norm(500,150).rvs()
b=stats.binom(10000,0.1).rvs()
x='random number from a '
print(x+"uniform distribution ",u)
print(x+" normal distribution ",n)
print(x+" binomial distribution ",b)
('random number from a uniform distribution ', -0.30353081440213836)
('random number from a  normal distribution ', 357.18541897080166)
('random number from a  binomial distribution', 1003)

Stucchhio 的 Python 程序，经过少许修改，显示如下：

from pylab import *
from scipy.stats import *
#input area
million=1e6                        # unit of million 
billion=1e9                        # unit of billion 
trillion=1e12                      # unit of trillion 
tiny=1e-7                          # a small number 
hourlyPay = 7.5                    # hourly wage
workingHoursPerWeek=40             # working hour per week                                
workingWeeksPerYear=50             # working weeks per year
nAdult           = 227*million     # number of adult
laborForce       = 154*million     # labor force
disabledAdults   =  21*million     # disability 
nSimulations     = 1024*32         # number of simulations 
#
basicIncome = hourlyPay*workingHoursPerWeek*workingWeeksPerYear
# define a few function
def geniusEffect(nNonWorkers):
    nGenious = binom(nNonWorkers,tiny).rvs()
    return nGenious* billion
#
def costBasicIncome():
    salaryCost= nAdult * basicIncome
    unitAdmCost = norm(250,75)
    nonWorkerMultiplier = uniform(-0.10, 0.15).rvs()
    nonWorker0=nAdult-laborForce-disabledAdults
    nNonWorker = nonWorker0*(1+nonWorkerMultiplier)
    marginalWorkerHourlyProductivity = norm(10,1)
    admCost = nAdult * unitAdmCost.rvs()
    unitBenefitNonWorker=40*52*marginalWorkerHourlyProductivity.rvs()
    benefitNonWorkers = 1 * (nNonWorker*unitBenefitNonWorker)
    geniusBenefit=geniusEffect(nNonWorker)
    totalCost=salaryCost + admCost - benefitNonWorkers-geniusBenefit
    return totalCost
#
def costBasicJob():
    unitAdmCost4disabled= norm(500,150).rvs()
    unitAdmCost4worker = norm(5000, 1500).rvs()
    nonWorkerMultiplier = uniform(-0.20, 0.25).rvs()
    hourlyProductivity = uniform(0.0, hourlyPay).rvs()
    cost4disabled=disabledAdults * (basicIncome + unitAdmCost4disabled)
    nBasicWorkers=((nAdult-disabledAdults-laborForce)*(1+nonWorkerMultiplier))
    annualCost=workingHoursPerWeek*workingWeeksPerYear*hourlyProductivity
    cost4workers=nBasicWorkers * (basicIncome+unitAdmCost4worker-annualCost)
    return cost4disabled + cost4workers
#
N = nSimulations
costBI = zeros(shape=(N,),dtype=float)
costBJ = zeros(shape=(N,),dtype=float)
for k in range(N):
    costBI[k] = costBasicIncome()
    costBJ[k] = costBasicJob()
#
def myPlot(data,myTitle,key):
    subplot(key)
    width = 4e12
    height=50*N/1024
    title(myTitle)
    #xlabel("Cost (Trillion = 1e12)")
    hist(data, bins=50)
    axis([0,width,0,height])
#
myPlot(costBI,"Basic Income",211)
myPlot(costBJ,"Basic Job",212)
show()

根据这里展示的图表，他得出结论，基本工作提案的成本低于基本收入提案。为了节省篇幅，我们不再详细阐述该程序。更多详细解释及相关假设，请参阅 Stucchhio（2013）发布的博客：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_12.jpg

通过模拟找到基于两只股票的有效前沿

以下程序旨在基于已知均值、标准差和相关性的两只股票生成有效前沿。我们只有六个输入值：两个均值，两个标准差，相关性（ρ）和模拟次数。为了生成相关的* y1 和 y2 时间序列，我们首先生成不相关的 x1 和 x2 *序列。然后，我们应用以下公式：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_33.jpg

另一个重要的问题是如何构造一个目标函数来最小化。我们的目标函数是投资组合的标准差，外加一个惩罚项，该惩罚项定义为与目标投资组合均值的绝对偏差的缩放。

换句话说，我们最小化投资组合的风险以及投资组合回报与目标回报之间的偏差，代码如下所示：

import numpy as np 
import scipy as sp 
import pandas as pd
import matplotlib.pyplot as plt
from datetime import datetime as dt 
from scipy.optimize import minimize
#
# Step 1: input area
mean_0=(0.15,0.25)   # mean returns for 2 stocks
std_0= (0.10,0.20)   # standard deviations for 2 stocks 
corr_=0.2       # correlation between 2 stocks
nSimulations=1000    # number of simulations 
#
# Step 2: Generate two uncorrelated time series 
n_stock=len(mean_0)
n=nSimulations
sp.random.seed(12345) # to get the same random numbers 
x1=sp.random.normal(loc=mean_0[0],scale=std_0[0],size=n) 
x2=sp.random.normal(loc=mean_0[1],scale=std_0[1],size=n) 
if(any(x1)<=-1.0 or any(x2)<=-1.0):
    print ('Error: return is <=-100%')
#
# Step 3: Generate two correlated time series 
index_=pd.date_range(start=dt(2001,1,1),periods=n,freq='d') 
y1=pd.DataFrame(x1,index=index_) 
y2=pd.DataFrame(corr_*x1+sp.sqrt(1-corr_**2)*x2,index=index_)
#
# step 4: generate a return matrix called R 
R0=pd.merge(y1,y2,left_index=True,right_index=True) 
R=np.array(R0)
#
# Step 5: define a few functions 
def objFunction(W, R, target_ret):
    stock_mean=np.mean(R,axis=0) 
    port_mean=np.dot(W,stock_mean)            # portfolio mean
    cov=np.cov(R.T)                           # var-covar matrix 
    port_var=np.dot(np.dot(W,cov),W.T)        # portfolio variance 
    penalty = 2000*abs(port_mean-target_ret)  # penalty 4 deviation
    return np.sqrt(port_var) + penalty        # objective function
#
# Step 6: estimate optimal portfolio for a given return 
out_mean,out_std,out_weight=[],[],[] 
stockMean=np.mean(R,axis=0)
#
for r in np.linspace(np.min(stockMean),np.max(stockMean),num=100): 
    W = sp.ones([n_stock])/n_stock             # start equal w
    b_ = [(0,1) for i in range(n_stock)]       # bounds
    c_ = ({'type':'eq', 'fun': lambda W: sum(W)-1\. })# constraint 
    result=minimize(objFunction,W,(R,r),method='SLSQP',constraints=c_,bounds=b_)
    if not result.success:                     # handle error 
        raise BaseException(result.message)
    out_mean.append(round(r,4))                # decimal places
    std_=round(np.std(np.sum(R*result.x,axis=1)),6) 
    out_std.append(std_) 
    out_weight.append(result.x)
#
# Step 7: plot the efficient frontier
plt.title('Simulation for an Efficient Frontier from given 2 stocks') 
plt.xlabel('Standard Deviation of the 2-stock Portfolio (Risk)') 
plt.ylabel('Return of the 2-stock portfolio')
plt.figtext(0.2,0.80,' mean = '+str(stockMean)) 
plt.figtext(0.2,0.75,' std  ='+str(std_0)) 
plt.figtext(0.2,0.70,' correlation ='+str(corr_))
plt.plot(np.array(std_0),np.array(stockMean),'o',markersize=8) 
plt.plot(out_std,out_mean,'--',linewidth=3)
plt.show()

输出如下所示：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_13.jpg

构建带有 n 只股票的有效前沿

当股票数量* n 增加时，每对股票之间的相关性会显著增加。对于 n 只股票，我们有 n(n-1)/2 个相关性。例如，如果* n 是 10，则我们有 45 个相关性。由于这个原因，手动输入这些值并不是一个好主意。相反，我们通过从多个均匀分布中抽取随机数来生成均值、标准差和相关性。为了产生相关的回报，首先我们生成 n *个不相关的股票回报时间序列，然后应用 Cholesky 分解，具体如下：

import numpy as np
import scipy as sp
import pandas as pd
import matplotlib.pyplot as plt
from datetime import datetime as dt
from scipy.optimize import minimize
#
# Step 1: input area
nStocks=20
sp.random.seed(1234)                        # produce the same random numbers 
n_corr=nStocks*(nStocks-1)/2                # number of correlation 
corr_0=sp.random.uniform(0.05,0.25,n_corr)  # generate correlations 
mean_0=sp.random.uniform(-0.1,0.25,nStocks) # means
std_0=sp.random.uniform(0.05,0.35,nStocks)  # standard deviation 
nSimulations=1000                           # number of simulations 
#
# Step 2: produce correlation matrix: Cholesky decomposition
corr_=sp.zeros((nStocks,nStocks))
for i in range(nStocks):
    for j in range(nStocks):
        if i==j:
            corr_[i,j]=1
        else:
            corr_[i,j]=corr_0[i+j]
U=np.linalg.cholesky(corr_)
#
# Step 3: Generate two uncorrelated time series 
R0=np.zeros((nSimulations,nStocks))
for i in range(nSimulations):
    for j in range(nStocks):
        R0[i,j]=sp.random.normal(loc=mean_0[j],scale=std_0[j],size=1)
if(R0.any()<=-1.0):
    print ('Error: return is <=-100%')
#
# Step 4: generate correlated return matrix: Cholesky     
R=np.dot(R0,U)
R=np.array(R)
#
# Step 5: define a few functions
def objFunction(W, R, target_ret): 
    stock_mean=np.mean(R,axis=0)  
    port_mean=np.dot(W,stock_mean)           # portfolio mean
    cov=np.cov(R.T)                          # var-covar matrix
    port_var=np.dot(np.dot(W,cov),W.T)       # portfolio variance
    penalty = 2000*abs(port_mean-target_ret) # penalty 4 deviation 
    return np.sqrt(port_var) + penalty       # objective function 
#
# Step 6: estimate optimal portfolo for a given return 
out_mean,out_std,out_weight=[],[],[] 
stockMean=np.mean(R,axis=0)    
#
for r in np.linspace(np.min(stockMean), np.max(stockMean), num=100):
    W = sp.ones([nStocks])/nStocks             # starting:equal w 
    b_ = [(0,1) for i in range(nStocks)]       # bounds
    c_ = ({'type':'eq', 'fun': lambda W: sum(W)-1\. })# constraint
    result=minimize(objFunction,W,(R,r),method='SLSQP',constraints=c_, bounds=b_)    
    if not result.success:                    # handle error
        raise BaseException(result.message) 
    out_mean.append(round(r,4))               # a few decimal places
    std_=round(np.std(np.sum(R*result.x,axis=1)),6)
    out_std.append(std_)
    out_weight.append(result.x) 
#
# Step 7: plot the efficient frontier
plt.title('Simulation for an Efficient Frontier: '+str(nStocks)+' stocks')
plt.xlabel('Standard Deviation of the Porfolio')
plt.ylabel('Return of the2-stock portfolio')
plt.plot(out_std,out_mean,'--',linewidth=3)
plt.show()

图表如下所示：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_14.jpg

当 n 是一个很大的数字时，模拟一个 n 股票组合是困难的。原因在于生成方差-协方差矩阵非常耗时，请看这里的协方差（相关性）数量：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_34.jpg

假设我们的投资组合中有 500 只股票。那么我们需要估算 124,750 对相关性。为了简化这个计算，我们可以应用资本资产定价模型（CAPM），见下式：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_35.jpg

这里 R[i,t] 是股票 i 在时间 t 的收益，α[i] 和 β[i] 分别是股票 i 的截距和斜率，R[M,t] 是市场指数在时间 t 的收益，e[i], t 是时间 t 的误差项。由于单只股票的总风险包含两个成分：系统性风险和公司特定风险。因此，股票 i 的方差与市场指数的关系如下所示：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_36.jpg

股票 i 和 j 之间的协方差如下所示：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_37.jpg

因此，我们可以将估算从 124,750 减少到仅 1,000。首先估算 500 个 β 值。然后我们应用之前的公式来估算协方差。类似地，估算股票 i 和 j 之间的相关性的公式如下：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_38.jpg

长期收益预测

许多研究人员和实践者认为，如果基于过去收益的算术平均数来进行长期收益预测，会导致过高的估计；而如果基于几何平均数，则会导致过低的估计。Jacquier、Kane 和 Marcus（2003）建议使用以下加权方案，利用 80 年的历史收益来预测未来 25 年的收益：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_39.jpg

以下程序反映了上述方程：

import numpy as np
import pandas as pd
from matplotlib.finance import quotes_historical_yahoo_ochl as getData 
#
# input area
ticker='IBM'           # input value 1 
begdate=(1926,1,1)     # input value 2 
enddate=(2013,12,31)   # input value 3 
n_forecast=25          # input value 4
#
def geomean_ret(returns): 
    product = 1
    for ret in returns: 
        product *= (1+ret)
    return product ** (1.0/len(returns))-1
#
x=getData(ticker,begdate,enddate,asobject=True, adjusted=True)
logret = np.log(x.aclose[1:]/x.aclose[:-1]) 
date=[]
d0=x.date
for i in range(0,np.size(logret)):
    date.append(d0[i].strftime("%Y"))
#
y=pd.DataFrame(logret,date,columns=['logret'],dtype=float)
ret_annual=np.exp(y.groupby(y.index).sum())-1 
ret_annual.columns=['ret_annual']
n_history=len(ret_annual) 
a_mean=np.mean(np.array(ret_annual))
g_mean=geomean_ret(np.array(ret_annual))
w=n_forecast/n_history
future_ret=w*g_mean+(1-w)*a_mean
print('Arithmetric mean=',round(a_mean,3), 'Geomean=',round(g_mean,3),'forecast=',future_ret)

输出结果如下所示：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_15.jpg

效率、拟蒙特卡洛法和 Sobol 序列

当使用蒙特卡洛模拟来解决各种金融相关问题时，会生成一组随机数。当精度要求非常高时，我们必须生成大量的随机数。例如，在定价期权时，我们使用非常小的间隔或大量的步长来提高解决方案的精度。因此，蒙特卡洛模拟的效率在计算时间和成本方面至关重要。如果需要定价几千个期权，这一点尤其重要。一种提高效率的方法是应用更好的算法，即优化我们的代码。另一种方法是使用一些分布更均匀的特殊类型的随机数，这被称为准蒙特卡洛模拟。一个典型的例子是所谓的 Sobol 序列。Sobol 序列属于低差异序列，满足随机数的属性，但分布更加均匀：

import numpy as np
import matplotlib.pyplot as plt
np.random.seed(12345)
n=200
a = np.random.uniform(size=(n*2))
plt.scatter(a[:n], a[n:])
plt.show()

附录 A – 数据案例#8 - 蒙特卡洛模拟与二十一点

二十一点是一种双人游戏，包括发牌员和玩家。这里，我们假设你是玩家。

规则#1：2 到 10 的牌按面值计分，而杰克、皇后和国王的点数为 10，A 的点数为 1 或 11（由玩家选择）。

术语：

二十一点：一张 A 牌加任何一张值 10 点的牌
输：玩家的赌注被发牌员收走
胜：玩家赢得与赌注相等的金额
二十一点（自然牌）：玩家赢得 1.5 倍赌注
平局：玩家保持赌注，既不赢也不输
步骤 1：发牌员抽两张牌，一张面朝上，而玩家抽两张牌（面朝上）
步骤 2：玩家可以抽第三张牌
胜或负：如果你的牌的点数少于 21 且大于发牌员的点数，你获胜。查看 www.pagat.com/banking/ blackjack.html

参考文献

请参考以下文章：

布鲁诺，朱塞佩，蒙特卡洛模拟在定价欧洲和美式篮子 期权中的应用，意大利银行，www.r-project.org/conferences/ useR-2010/abstracts/Bruno.pdf
Easley, D., Kiefer, N.M., O’Hara, M., Paperman, J.B., 1996，流动性、信息和不常交易的股票，《金融学杂志》51 卷，1405–1436，www.defaultrisk.com/pa_liqty_03.htm
Jacquier, Eric, Alex Kane, and Alan J. Marcus, 2003, 几何平均数与算术平均数：一种重新考虑, www2.bc.edu/alan-marcus/papers/FAJ_2003.pdf
Stucchio, Chris, 2013, 使用 Python 和蒙特卡罗模拟模型构建基本收入, www.chrisstucchio.com/blog/2013/basic_income_vs_basic_job.html
Zhong, Raymond, 2017, 印度考虑通过普遍基本收入来对抗贫困*, 《华尔街日报》*, blogs.wsj.com/indiarealtime/2017/01/31/india-considers-fighting-poverty-witha-universal-basic-income/

练习

finance.yahoo.com，下载一些公司（如 IBM、WMT 和 C（花旗集团））过去五年的价格数据。测试它们的每日回报是否符合正态分布。
编写一个 Python 程序，使用scipy.permutation()函数，从过去五年的数据中随机选择 12 个月的回报，不放回。为了测试该程序，你可以使用花旗集团的数据，时间范围为 2012 年 1 月 2 日到 2016 年 12 月 31 日，数据来源于 Yahoo! Finance。
编写一个 Python 程序，使用给定的n个回报进行自助法抽样。每次选择m个回报，其中m>n。
为了将均匀分布的随机数转换为正态分布，我们有以下公式：https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_12_40.jpg

根据公式，生成 5,000 个正态分布的随机数；估算它们的均值、标准差，并进行测试。
假设当前股票价格为$10.25，过去五年的均值为$9.35，标准差为 4.24。编写一个 Python 程序，生成 1,000 个未来的股票价格。
下载过去 10 年内 10 只股票的价格数据。构建一个等权重投资组合，并对该组合的每日回报进行 Shapiro-Wilk 检验：

公司名称股票代码戴尔公司 DELL
国际商业机器公司 IBM 通用电气 GE
微软 MSFT 谷歌 GOOG
Family Dollar Stores FDO Apple AAPL
Wal-Mart Stores WMT eBay EBAY
麦当劳 MCD
前往 Yahoo! Finance 查找今天的 IBM 股票价格，然后下载其历史价格信息，估算过去五年的均值和标准差。生成未来一年每日价格的预测。
对于 20 个股票代码，下载并保存它们的每日价格为 20 个不同的 CSV 文件。编写一个 Python 程序，随机选择五只股票，估算它们的等权重投资组合回报和风险。
重复前一个程序，但将其保存为一个文件，而不是 20 个单独的 CSV 文件。

提示

生成一个额外的变量，称为 ticker。
班级里有 30 名学生。编写一个程序，从中随机选择七个学生。
测试提取ffMonthly.pkl、ffDaily.pkl或ffMonthly.csv、ffDaily.csv的时间差，并进行一些测试。
通常我们观察到投资组合的波动率与投资组合中股票数量之间存在负相关关系。编写一个程序，显示投资组合方差与其中股票数量之间的关系。
从标记为 1 到 10 的 10 个球中，抽取 1、2、3 和 4 的概率是多少？使用两种方法：a. 使用公式。b. 编写一个程序生成五个随机数。
编写一个程序生成 176 百万组 Mega Millions 游戏的组合。赢得（1、2、3、4、5）和（1）的机会是多少？
对于 Powerball 游戏，我们从 59 个标有 1 到 59 的白球中选择 5 个白球，并从 39 个标有 1 到 39 的红球中选择 1 个红球。编写一个程序随机选择这六个球。
从 20 只股票中选择 7 只，选择前七只股票的概率是多少？使用模拟来证明你的结果。

公司名称	股票代码	戴尔公司	DELL
国际商业机器公司	IBM	通用电气	GE
微软	MSFT	谷歌	GOOG
Family Dollar Stores	FDO	Apple	AAPL
Wal-Mart Stores	WMT	eBay	EBAY
麦当劳	MCD

总结

在本章中，我们讨论了几种类型的分布：正态分布、标准正态分布、对数正态分布和泊松分布。由于假设股票遵循对数正态分布且回报遵循正态分布是期权理论的基石，因此蒙特卡洛模拟被用于定价欧式期权。在某些情况下，亚洲期权可能在对冲方面更为有效。奇异期权比标准期权更为复杂，因为前者没有封闭解，而后者可以通过 Black-Scholes-Merton 期权模型定价。定价这些奇异期权的一种方法是使用蒙特卡洛模拟。我们还讨论了定价亚洲期权和回溯期权的 Python 程序。

第十三章：信用风险分析

信用风险分析的目标是衡量潜在违约支付承诺金额的概率。信用评级反映了一个公司或债券的信用状况。公司的评级不同于其债券的评级，因为后者取决于债券的到期日以及某些特征，比如是否具有可赎回或可出售的选项。在第五章《债券与股票估值》中，我们学到了到期收益率（YTM），或简称收益率，它与信用质量相关。信用质量越低，要求的回报率越高，即收益率越高。在本章中，我们将讨论与信用风险相关的许多基本概念，如信用评级、信用利差、一年期信用评级迁移矩阵、违约概率、违约损失率、回收率和 KMV 模型。具体来说，将涵盖以下主题：

穆迪（Moody’s）、标准普尔（Standard and Poor’s）和惠誉（Fitch）的信用评级
信用利差、一年期和五年期迁移矩阵
利率的期限结构
未来利率的模拟
阿尔特曼 Z 评分预测企业破产
KMV 模型估算总资产及其波动性
违约概率和违约距离
信用违约掉期

信用风险分析简介

在本章中，我们将讨论与信用风险相关的基本概念，如信用评级、信用利差、一年期和五年期评级迁移矩阵、违约概率、回收率和违约损失率。信用利差，即债券收益率与基准收益率（无风险利率）之间的差异，反映了其信用风险或违约风险。例如，要估算一只 AA 评级债券在两年后的票息支付现值，折现率（收益率）将是无风险收益率（国债收益率）加上相应的利差。在分析公司或债券的信用状况时，有很多工具可供使用。第一个工具是信用评级，由信用评级机构提供，如穆迪（Moody’s）或标普（Standard and Poor’s）。其显著的优点是潜在用户可以花费更少的时间和精力评估公司或债券的信用风险。明显的缺点是，信用评级对于大多数用户而言是一个“黑匣子”。换句话说，用户无法复制一个信用评级。因此，想要揭示这种简单字母评级系统背后的逻辑（例如 AA 或 A1）是相当困难的。还有其他方式来评估公司（债券）的信用状况，例如利差，这是一个容易获取的数据。最量化的模型之一是所谓的 KMV 模型，它应用了我们在第十章中学到的期权理论来评估公司的信用风险。

信用评级

目前，美国有三大信用评级机构：穆迪、标准普尔和惠誉。它们的网站分别是 www.moodys.com/，www.standardandpoors.com/en_US/web/guest/home 和 www.fitchratings.com/site/home。尽管它们的评级符号（字母）不同，但将一个评级机构的字母评级转换为另一个评级机构的评级是非常容易的。根据以下链接 www.quadcapital.com/Rating%20Agency%20Credit%20Ratings.pdf，生成了一个名为 creditRatigs3.pkl 的数据集，可以从作者的网站下载，网址是 http://canisius.edu/~yany/python/creditRatings3.pkl。假设它位于 C:/temp/ 目录下。

以下代码显示其内容：

import pandas as pd
x=pd.read_pickle("c:/temp/creditRatings3.pkl")
print(x)
       Moody's S&P Fitch  NAIC  InvestmentGrade
0      Aaa   AAA   AAA     1                1
1      Aa1   AA+   AA+     1                1
2      Aa2    AA    AA     1                1
3      Aa3   AA-   AA-     1                1
4       A1    A+    A+     1                1
5       A2     A     A     1                1
6       A3    A-    A-     1                1
7     Baa1  BBB+  BBB+     2                1
8     Baa2   BBB   BBB     2                1
9     Baa3  BBB-  BBB-     2                1
10     Ba1   BB+   BB+     3                0
11     Ba2    BB    BB     3                0
12     Ba3   BB-   BB-     3                0
13      B1    B+    B+     3                0
14      B2     B     B     3                0
15      B3    B-    B-     3                0

第一列是行号，没有特定含义。接下来的三列分别是穆迪、标准普尔 和惠誉的信用评级。NAIC 代表 国家保险委员会协会。任何评级等于或高于 BBB 的被归类为投资级别，参见最后一列（变量），其值为 1 或 0。许多共同基金和养老金基金只能投资评级为投资级别的债券。

当一家公司今年获得 Aaa 评级时，它明年保持相同信用评级的概率是多少？根据以下表格，它明年保持 Aaa 评级的概率为 89%，穆迪（2007 年）。另一方面，降级的概率是 3%，即从 Aaa 降级为 Aa1。对于 B1 评级的债券，保持相同评级的概率是 65%。同时，它有 12% 的概率升级，9% 的概率降级。B1 评级债券的违约概率是 3%，请参见以下图表的最后一列，给出了这一年的信用评级迁移矩阵：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_13_01.jpg

一年期信用评级迁移矩阵

注意以下缩写：

WR 表示穆迪已撤回其评级
DEF 表示违约概率

同样，Aaa 评级公司在明年变成 Aa2 评级的概率是 3%。主对角线上的数值（从西北到东南）是明年保持相同评级的概率。主对角线下方的数值（左下三角和底部三角）是降级的概率，而主对角线上方的数值（右上三角）是升级的概率。最后一列提供了各种评级的违约概率。例如，Ba2 评级的债券违约概率为 1%，而 Caa3 评级的债券违约概率为 25%。可以使用名为 migration1year.pkl 的 Python 数据集，代码如下所示。数据集可以通过以下网址获取：http://canisius.edu/~yany/python/migration1year.pkl：

import pandas as pd
x=pd.read_pickle("c:/temp/migration1year.pkl")
print(x.head(1))
print(x.tail(1))
    Aaa   Aa1   Aa2  Aa3   A1   A2   A3  Baa1  Baa2  Baa3 ...   Ba3   B1  \
Aaa  0.89  0.03  0.03  0.0  0.0  0.0  0.0   0.0   0.0   0.0 ...   0.0  0.0   
      B2   B3  Caa1  Caa2  Caa3  Ca-C    WR  DEF  
Aaa  0.0  0.0   0.0   0.0   0.0   0.0  0.05  0.0  
[1 rows x 22 columns]
      Aaa  Aa1  Aa2  Aa3   A1   A2   A3  Baa1  Baa2  Baa3 ...   Ba3   B1   B2\
Ca-C  0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0   0.0   0.0 ...   0.0  0.0  0.0   
       B3  Caa1  Caa2  Caa3  Ca-C    WR  DEF  
Ca-C  0.0  0.01  0.01  0.01  0.35  0.13  0.2  
[1 rows x 22 columns]

以下表格展示了穆迪的 5 年过渡（迁移）矩阵。请注意DEF列（用于违约概率）：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_13_02.jpg

穆迪的平均 5 年评级过渡矩阵（1920-1992）

来源：穆迪（2007 年）。

注意以下缩写：

WR 表示穆迪已撤销其评级
DEF 表示违约概率

一个名为migration5year.pkl的数据集已生成。该数据集可以在canisius.edu/~yany/python/migration5year.pkl下载。以下代码将打印其首尾行：

import pandas as pd
x=pd.read_pickle("c:/temp/migration5year.pkl")
print(x.head(1))
print(x.tail(1))
    Aaa   Aa1  Aa2   Aa3    A1    A2   A3  Baa1  Baa2  Baa3 ...   Ba3   B1  \
Aaa  0.56  0.07  0.1  0.03  0.01  0.01  0.0   0.0   0.0   0.0 ...   0.0  0.0   
      B2   B3  Caa1  Caa2  Caa3  Ca-C   WR  DEF  
Aaa  0.0  0.0   0.0   0.0   0.0   0.0  0.2  0.0  
[1 rows x 22 columns]
      Aaa  Aa1  Aa2  Aa3   A1   A2   A3  Baa1  Baa2  Baa3  ...   Ba3   B1  \
Ca-C  0.0  0.0  0.0  0.0  0.0  0.0  0.0   0.0   0.0   0.0  ...   0.0  0.0   
        B2    B3  Caa1  Caa2  Caa3  Ca-C    WR   DEF  
Ca-C  0.02  0.02  0.01  0.01  0.01  0.04  0.43  0.46

评级与违约呈负相关。评级越高，违约概率越低。以下是累积历史违约率（以百分比表示）：

	违约率（%）
	穆迪
评级类别	市政
`Aaa`/`AAA`	0.00
`Aa`/`AA`	0.06
`A`/`A`	0.03
`Baa`/`BBB`	0.13
`Ba`/`BB`	2.65
`B`/`B`	11.86
`Caa-C`/`CCC-C`	16.58

平均值
投资级	0.07
非投资级	4.29
所有	0.10

表 13.3 信用评级与 DP（违约概率）之间的关系

数据来自于网站monevator.com/bond-default-rating-probability/。

例如，对于穆迪的Aaa级相关公司债券，其违约概率为 0.52%。标准普尔的对应违约概率为 0.60%。给定违约后的回收率是一个重要概念。债务的状态（资历）对回收率有重大影响。根据 Altman 和 Kishore（1997 年）的研究，我们得到了以下表格：

	回收率（占面值的百分比）
高级担保债务	58%
高级-非担保债务	48%
高级-次级债务	35%
次级债务	32%
折扣债券和零息债券	21%

表 13.4 基于资历的回收率

担保债务是指由资产担保支付的债务。高级债务和次级债务是指优先级结构。另一方面，不同的行业有不同的回收率，这是由于它们各自的行业特征，如固定长期资产和无形资产的比例：

行业	平均回收率	观察数量
公用事业	70.5%	56
化工、石油、橡胶和塑料制品	62.7%	35
机械、仪器及相关产品	48.7%	36
服务业 - 商业和个人	46.2%	14
食品及相关产品	45.3%	18
批发和零售贸易	44.0%	12
多元化制造业	42.3%	20
赌场、酒店和娱乐	40.2%	21
建筑材料、金属和加工制品	38.8%	68
运输和运输设备	38.4%	52
通信、广播、电影制作	37.1%	65
印刷和出版	NA	NA
金融机构	35.7%	66
建筑和房地产	35.3%	35
一般商品商店	33.2%	89
矿业和石油钻探	33.0%	45
纺织和服装产品	31.7%	31
木材、纸张和皮革制品	29.8%	11
住宿、医院和护理设施	26.5%	22
总计	41.0%	696

表 13.5 基于行业的回收率

请参见关于回收率的文章： Recovery Rates PDF。

前面的表格是根据回收率从高到低排序的。对于印刷和出版行业，根据原始数据没有相关数据。违约损失（LGD）等于 1 减去回收率：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_13_09.jpg

在这里，我们通过一个假设性的例子来解释默认概率和回收率的使用，以计算债券的价格。假设一只一年期债券的面值为 100 美元，票息率为 6%，到期收益率（YTM）为 7%。我们有以下四种情况：

情况 #1：无违约。今天的价格将是其折现后的未来现金流，(6+100)/(1+0.07)。
情况 #2：确定违约并且无法回收任何金额。对于这种情况，其价格将为零。
情况 #3：如果发生违约，我们将无法收回任何金额。
情况 #4：如果发生违约，我们将收到某些金额。

以下表格总结了前面四种情况：

#	条件	违约概率回收率	今天的价格
1	无违约	P=0, 回收率（NA）	$99.07
2	100% 违约/无回收	P=100%，回收率=0	0
3	如果违约，无法回收任何金额	O<P<100%，回收率=0	$99.07 *(1-P)
4	如果违约，回收部分金额	O<P<100%，回收率>0	$99.07 1-P(1- ![信用评级)]

表 13.6 不同违约概率和回收率的四种情况

债券的价格是其所有预期未来现金流现值的总和：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_13_10.jpg

如果P是违约概率，我们有以下预期未来现金流：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_13_11.jpg

折现所有未来现金流可以得出债券的价格：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_13_12.jpg

假设根据穆迪评级，信用评级为 A。根据表 13.3，其违约率为 1.29%。进一步假设它是一家公用事业公司。因此，根据表 13.5，违约时的回收率为 70.5%。该债券的面值为 100 美元，要求回报率（YTM）为 5%。根据前述公式，若没有违约，一年期债券的价格将为 95.24 美元，即 100/(1+0.05)。我们这只债券的卖出价格，将在 1.29%的违约概率下为 94.88 美元，即95.24(1-0.0129*(1-0.705))*。

信用利差

信用利差（违约风险溢价）反映了其违约风险。例如，要估算一个 AA 评级债券在两年后的票息支付现值，折现率（收益率）将是无风险利率加上相应的利差。对于给定的信用评级，可以通过使用历史数据来找到其信用利差。以下是一个典型的表格，显示了信用风险溢价（利差）与信用评级之间的关系，请参见下表：

我们感谢亚当·阿莫多尔教授在其网站上提供的数据集，people.stern.nyu.edu/adamodar/pc/datasets/：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_13_03.jpg

基于信用评级的信用利差

除了前面表格中的最后一行外，利差的单位是基点，即百分之一的百分之一。例如，对于一个 A-（A 减）评级的五年期债券，其利差为 83.6 个基点。由于无风险利率为 1.582%（五年期国债利率），该债券的到期收益率（YTM）将是 2.418%，即 0.01582+83.6/100/100。根据前面的表格，我们生成了一个名为bondSpread2014.p的 Python 数据集，数据集可以在作者的网站上找到，canisius.edu/~yany/python/creditSpread2014.pkl：

import pandas as pd
x=pd.read_pickle("c:/temp/creditSpread2014.pkl")
print(x.head())
print(x.tail())
  Rating     1     2     3     5     7    10     30
0  Aaa/AAA   5.0   8.0  12.0  18.0  28.0  42.0   65.0
1  Aa1/AA+  11.2  20.0  27.0  36.6  45.2  56.8   81.8
2   Aa2/AA  16.4  32.8  42.6  54.8  62.8  71.2   97.8
3  Aa3/AA-  21.6  38.6  48.6  59.8  67.4  75.2   99.2
4    A1/A+  26.2  44.0  54.2  64.6  71.4  78.4  100.2
               Rating        1        2        3        5        7       10  \
13              B1/B+  383.600  409.600  431.400  455.600  477.600  500.800   
14               B2/B  455.800  481.600  505.200  531.000  555.400  581.400   
15              B3/B-  527.800  553.800  579.400  606.400  633.600  661.800   
16           Caa/CCC+  600.000  626.000  653.000  682.000  712.000  743.000   
17  US Treasury Yield    0.132    0.344    0.682    1.582    2.284    2.892

经过仔细研究前面的表格后，我们会发现两个单调趋势。首先，利差是信用质量的递减函数。信用评级越低，其利差越高。其次，对于相同的信用评级，其利差每年都会增加。例如，对于 AAA 评级的债券，一年期的利差为 5 个基点，而五年期的利差为 18 个基点。

AAA 评级债券的收益率，Altman Z 分数

从前面的章节中，我们已经了解到，债券的收益率与同到期国债收益率之间的差额是违约风险溢价。为了获取AAA和AA债券的收益率，我们使用以下代码。穆迪的Aaa公司债券收益率可以在 https://fred.stlouisfed.org/series/AAA 下载。数据集可以在canisius.edu/~yany/python/moodyAAAyield.p下载。请注意，.p格式的.png文件适用于.pickle格式：

import pandas as pd
x=pd.read_pickle("c:/temp/moodyAAAyield.p")
print(x.head())
print(x.tail())

输出如下所示：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_13_04.jpg

请注意，名为moodyAAAyield.p的数据集中的第二列值是年化的。因此，如果我们想估计 1919 年 1 月的月度收益率（回报率），该收益率应为 0.4458333%，即 0.0535/12。

艾尔特曼 Z-score 广泛应用于金融领域的信用分析，用于预测公司破产的可能性。该分数是基于公司资产负债表和损益表的五个比率的加权平均值。对于上市公司，艾尔特曼（1968）提供了以下公式：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_13_13.jpg

在此，X1、X2、X3、X4 和 X5 的定义列在下表中：

变量	定义
X1	息税前利润/总资产
X2	净销售额/总资产
X3	股本市场价值/总负债
X4	营运资金/总资产
X5	留存收益/总资产

表 13.8 Z-scores 估算中变量的定义

基于 Z-score 的范围，我们可以将上市公司分为以下四类。Eidlenan（1995）发现，Z-score 正确预测了 72%的破产事件，这些破产事件发生在事件发生前两年：

Z-score 范围	描述
> 3.0	安全
2.7 至 2.99	警戒状态
1.8 至 2.7	两年内破产的高概率
< 1.80	财务困境的概率非常高

艾尔特曼 Z-score 属于信用评分（方法）类别。另一方面，更先进的模型，例如 KMV 模型，是基于现代金融理论，如期权理论。

使用 KMV 模型估算总资产的市场价值及其波动性

KMV 代表Kealhofer、McQuown和Vasicek，他们创办了一家公司，专注于衡量违约风险。KMV 方法是通过使用公司资产负债表信息和股市信息来估计公司违约概率的最重要方法之一。本节的目标是展示如何估算总资产（A）的市场价值及其相应的波动性（σA）。该结果将在本章后续使用。基本思路是将公司的股本视为一个看涨期权，其债务的账面价值视为行权价。让我们看一个最简单的例子。对于一家公司，如果其债务为$70，股本为$30，则总资产为$100，见下表：

100	70
30

假设总资产跳升至$110，债务保持不变。现在，股本的价值增加到$40。另一方面，如果资产下降到$90，股本将被评估为$20。由于股东是剩余索偿人，他们的价值满足以下表达式：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_13_14.jpg

这里，E 是股权的价值，A 是总资产，D 是总债务水平。回顾欧式看涨期权，我们有以下支付函数：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_13_15.jpg

这里，ST 是到期日的终端股票价格，T 是到期日，K 是执行价格，max() 是最大值函数。前两个公式之间的相似性表明，我们可以将股权视为以债务水平为执行价格的看涨期权。通过适当的符号表示，我们将得到以下公司的股权公式。KMV 模型在此定义：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_13_16.jpg

另一方面，以下是股权与总资产波动率之间的关系。在下面的公式中，我们有：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_13_27.jpg https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_13_17.jpg

由于 d1 和 d2 是通过前面公式定义的，我们有两个关于两个未知数 (A 和 https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_13_28.jpg) 的方程；请参见以下公式。因此，我们可以使用试错法或联立方程法来求解这两个未知数。最终，我们要解以下两个关于 A 和 https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_13_29.jpg 的方程：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_13_18.jpg

我们应该注意前面公式中估算的 A（总资产市场价值），因为它不同于资产市场价值与债务账面价值的总和。

以下 Python 程序用于估算给定 E（股权）、D（债务）、T（到期日）、r（无风险利率）和股权波动率（sigmaE）下的总资产（A）及其波动率（sigmA）。该程序的基本逻辑是，我们输入大量的 (A, sigmaE) 配对数据，然后根据前面的公式估算 E 和 sigmaE。由于 E 和 sigmaE 都是已知的，我们可以估算出差异，diff4E=estimatedE – knownE 和 diff4sigmaE = estimatedSigmaE – knownSigmaE。最小化这两个绝对差值和的 (A, sigmaE) 配对即为我们的解：

import scipy as sp
import pandas as pd
import scipy.stats as stats
from scipy import log,sqrt,exp
# input area 
D=30\.            # debt
E=70\.            # equity 
T=1\.             # maturity 
r=0.07           # risk-free
sigmaE=0.4       # volatility of equity 
#
# define a function to siplify notations later 
def N(x):
    return stats.norm.cdf(x)
#
def KMV_f(E,D,T,r,sigmaE):
    n=10000
    m=2000
    diffOld=1e6     # a very big number
    for i in sp.arange(1,10):
        for j in sp.arange(1,m):
            A=E+D/2+i*D/n
            sigmaA=0.05+j*(1.0-0.001)/m
            d1 = (log(A/D)+(r+sigmaA*sigmaA/2.)*T)/(sigmaA*sqrt(T))
            d2 = d1-sigmaA*sqrt(T)
            diff4A= (A*N(d1)-D*exp(-r*T)*N(d2)-E)/A  # scale by assets
            diff4sigmaE= A/E*N(d1)*sigmaA-sigmaE     # a small number already
            diffNew=abs(diff4A)+abs(diff4sigmaE)
            if diffNew<diffOld:
               diffOld=diffNew
               output=(round(A,2),round(sigmaA,4),round(diffNew,5))
    return output
#
print("KMV=", KMV_f(D,E,T,r,sigmaE))
print("KMV=", KMV_f(D=65e3,E=110e3,T=1,r=0.01,sigmaE=0.2))

输出如下所示：

print("KMV=", KMV_f(D,E,T,r,sigmaE))

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_13_05.jpg

请注意结果，因为债务的账面价值与股本的市场价值总和为 175,000，而我们估算的结果为 142,559。由于公司的股本是看涨期权，我们可以使用 Black-Scholes-Merton 模型来再次验证我们的结果。

利率期限结构

在第五章《债券与股票估值》中，我们讨论了利率期限结构的概念。利率期限结构被定义为无风险利率与时间之间的关系。无风险利率通常被定义为无违约的国债利率。从许多来源，我们可以获取当前的利率期限结构。例如，在 2017 年 2 月 27 日，我们可以从 http://finance.yahoo.com/bonds 获取以下信息：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_13_06.jpg

绘制的利率期限结构可能更吸引眼球；请参见以下代码：

import matplotlib.pyplot as plt
time=[3./12.,6./12.,2.,3.,5.,10.,30.]
rate=[0.45,0.61,1.12,1.37,1.78,2.29,2.93]
plt.title("Term Structure of Interest Rate ")
plt.xlabel("Time (in years) ")
plt.ylabel("Risk-free rate (%)")
plt.plot(time,rate)
plt.show()

违约距离

违约距离 (DD) 由以下公式定义；这里 A 是总资产的市场价值，https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_13_22.jpg

就 违约点 而言，目前没有理论指导如何选择理想的违约点。然而，我们可以将所有短期债务加上长期债务的一半作为我们的违约点。在得到资产的市场价值及其波动性后，我们可以使用前述公式来估算违约距离。A 和 https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_13_31.jpg 来自 公式（10） 的输出。另一方面，如果违约点等于 E，我们将得到以下公式：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_13_23.jpg

根据 Black-Scholes-Merton 看涨期权模型，DD 与 DP（违约概率） 之间的关系如下：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_13_24.jpg

信用违约掉期

债权人可以购买一种所谓的信用违约掉期（CDS）来在违约发生时进行保护。CDS 的买方向卖方支付一系列款项，作为交换，如果贷款发生违约，买方将获得赔偿。我们来看一个简单的例子。一个基金刚刚购买了 1 亿美元的公司债券，债券的到期时间为 15 年。如果发行公司没有发生违约，养老基金将每年享受利息支付，并在到期时拿回 1 亿美元本金。为了保护他们的投资，基金与一家金融机构签订了一份 15 年的 CDS 合约。根据债券发行公司的信用状况，约定的利差为 80 个基点，按年支付。这意味着每年，养老基金（CDS 买方）将在未来 10 年里支付给金融机构（CDS 卖方）80,000 美元。如果发生信用事件，CDS 卖方将根据其损失向 CDS 买方进行赔偿。如果合同规定的是实物结算，CDS 买方可以将债券以 1 亿美元的价格卖给 CDS 卖方。如果合同规定的是现金结算，CDS 卖方将支付*Max（$100m-X,0）*给 CDS 买方，其中 X 是债券的市场价值。如果债券的市场价值为 7000 万美元，则 CDS 卖方将向 CDS 买方支付 3000 万美元。在上述案例中，利差或费用与发行公司违约的概率密切相关。违约概率越高，CDS 的利差越高。下表显示了这种关系：

CDS	P	CDS	P	CDS	P	CDS	P	CDS	P	CDS	P	CDS	P
0	0.0%	100	7.8%	200	13.9%	300	19.6%	500	30.2%	500	30.2%	1000	54.1%
5	0.6%	105	8.1%	205	14.2%	310	20.2%	510	30.7%	525	31.4%	1025	55.2%
10	1.1%	110	8.4%	210	14.5%	320	20.7%	520	31.2%	550	32.7%	1050	56.4%
15	1.6%	115	8.7%	215	14.8%	330	21.2%	530	31.7%	575	33.9%	1075	57.5%
20	2.0%	120	9.1%	220	15.1%	340	21.8%	540	32.2%	600	35.2%	1100	58.6%
25	2.4%	125	9.4%	225	15.4%	350	22.3%	550	32.7%	625	36.4%	1125	59.7%
30	2.8%	130	9.7%	230	15.7%	360	22.9%	560	33.2%	650	37.6%	1150	60.9%
35	3.2%	135	10.0%	235	16.0%	370	23.4%	570	33.7%	675	38.8%	1175	62.0%
40	3.6%	140	10.3%	240	16.2%	380	23.9%	580	34.2%	700	40.0%	1200	63.1%
45	4.0%	145	10.6%	245	16.5%	390	24.5%	590	34.7%	725	41.2%	1225	64.2%
50	4.3%	150	10.9%	250	16.8%	400	25.0%	600	35.2%	750	42.4%	1250	65.3%
55	4.7%	155	11.2%	255	17.1%	410	25.5%	610	35.7%	775	43.6%	1275	66.4%
60	5.0%	160	11.5%	260	17.4%	420	26.0%	620	36.1%	800	44.8%	1300	67.5%
65	5.4%	165	11.8%	265	17.7%	430	26.6%	630	36.6%	825	46.0%	1325	68.6%
70	5.7%	170	12.1%	270	17.9%	440	27.1%	640	37.1%	850	47.2%	1350	69.7%
75	6.1%	175	12.4%	275	18.2%	450	27.6%	650	37.6%	875	48.3%	1375	70.7%
80	6.4%	180	12.7%	280	18.5%	460	28.1%	660	38.1%	900	49.5%	1400	71.8%
85	6.8%	185	13.0%	285	18.8%	470	28.6%	670	38.6%	925	50.6%	1425	72.9%
90	7.1%	190	13.3%	290	19.1%	480	29.1%	680	39.1%	950	51.8%	1450	74.0%
95	7.4%	195	13.6%	295	19.3%	490	29.6%	690	39.6%	975	52.9%	1475	75.1%
100	7.8%	200	13.9%	300	19.6%	500	30.2%	700	40.0%	1000	54.1%	1500	76.1%

表 13.9：违约概率和信用违约掉期。

违约概率估计的五年累计违约概率（P）

和五年期信用违约掉期（5Y CDS）

附录 A – 数据案例#8 - 使用 Z 分数预测破产

Altman 的 Z 分数用于预测公司破产的可能性。该分数是基于公司资产负债表和利润表的五个比率的加权平均值。对于上市公司，Altman（1968）提供了以下公式：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_13_25.jpg

在这里，X1、X2、X3、X4 和 X5 的定义在下面的表格中给出：

变量	定义
X1	息税前利润/总资产
X2	净销售额/总资产
X3	市场价值的股本/总负债
X4	营运资金/总资产
X5	留存收益/总资产

根据 Z 分数的范围，我们可以将 20 家上市公司分为以下四类。Eidlenan（1995）发现，Z 分数能够准确预测事件发生前两年的 72%的破产情况：

Z 分数范围	描述
> 3.0	安全
2.7 到 2.99	处于警戒状态
1.8 到 2.7	两年内破产的可能性较大
< 1.80	财务困境的可能性非常高

参考文献

Altman, Edward I, 1968，财务比率、判别分析与公司破产预测，金融学杂志，189–209，onlinelibrary.wiley.com/doi/10.1111/j.1540-6261.1968.tb00843.x/abstract
Altman, E.I., Kishore, V., 1997。高收益债务市场的违约与回报，1991-1996，纽约大学 Salomon 中心特别报告
Altman, Edward I., 2000，公司财务困境预测，pages.stern.nyu.edu/~ealtman/Zscores.pdf
Eidleman, Gregory J.,1995,Z-Scores – A Guide to Failure Prediction, CPA Journal Online，www.easycalculation.com/statistics/altman-z-score.php
Fitch，www.fitchratings.com/site/home。
KMV 模型，github.com/ghlingjun/kmv-model
穆迪网站，www.moodys.com/
穆迪，2007，介绍穆迪信用迁移模型，www.moodysanalytics.com/~/media/Brochures/Credit-Research-Risk-Measurement/Quantative-Insight/Credit-Transition-Model/Introductory-Article-Credit-Transition-Model.pdf
标准普尔，www.standardandpoors.com/en_US/web/guest/home

练习

美国有多少家信用评级机构？哪些是主要的机构？
风险的定义有几种？信用风险和市场风险有什么区别？
如何估计一个公司的总风险和市场风险？相关的数学公式是什么？
如何估计一个公司的信用风险？相关的数学公式是什么？
为什么债券的信用风险可能与其公司的信用评级不同？
如果一切条件相等，哪种债券的风险更高，是长期债券还是短期债券？
信用利差的定义是什么？为什么它有用？
利率期限结构的用途是什么？
对于 Altman 的 Z-score，X1、X2、X3、X4和X5的定义是什么？解释为什么 Z-score 越高，破产的概率越低：https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_13_26.jpg
确定 Z-score 的一个问题并找到解决方法。
一年期迁移（转移）矩阵是什么？
信用评级与违约概率之间有什么关系？
使用债券现值的概念，如何通过折现预计的未来现金流来推导方程式（1）？
信用迁移矩阵中（从西北到东南）的主对角线上的值是什么？
沃尔玛计划发行 5000 万美元（总面值）的公司债，每个债券的面值为 1000 美元。这些债券将在 10 年后到期，票面利率为 8%，每年支付一次利息。沃尔玛今天能筹集多少资金？如果沃尔玛将其信用评级提高一个等级，该公司能筹集多少额外资金？

下表展示了评级、违约风险（利差）与时间之间的关系。编写一个 Python 程序来插值缺失的利差，如从第 11 年到第 29 年的 S。Python 数据集可以从canisius.edu/~yany/python/creditSpread2014.p下载：

import matplotlib.pyplot as plt
import pandas as pd
x=pd.read_pickle("c:/temp/creditSpread2014.p")
print(x.head())
    Rating     1     2     3     5     7    10     30
0  Aaa/AAA   5.0   8.0  12.0  18.0  28.0  42.0   65.0
1  Aa1/AA+  11.2  20.0  27.0  36.6  45.2  56.8   81.8
2   Aa2/AA  16.4  32.8  42.6  54.8  62.8  71.2   97.8
3  Aa3/AA-  21.6  38.6  48.6  59.8  67.4  75.2   99.2
4    A1/A+  26.2  44.0  54.2  64.6  71.4  78.4  100.2

总结

在本章中，我们从信用风险分析的基础知识开始，涵盖信用评级、信用利差、1 年期评级迁移矩阵、违约概率（PD）、违约损失（LGD）、利率期限结构、Altman Z 分数、KMV 模型、违约概率、违约距离和信用违约掉期。在第十章，期权与期货，讨论了一些基本的香草期权，如 Black-Scholes-Merton 期权及其相关应用。此外，在第十二章，蒙特卡洛模拟中，解释了两种特殊期权。

在下一章中，我们将讨论更多的特殊期权，因为它们对于缓解许多类型的金融风险非常有用。

第十四章：另类期权

在第十章，期权与期货中，我们讨论了著名的布莱克-斯科尔斯-梅顿期权模型以及涉及各种类型期权、期货和基础证券的各种交易策略。布莱克-斯科尔斯-梅顿的闭式解法适用于只能在到期日行使的欧式期权。美式期权可以在到期日之前或当天行使。通常，这些类型的期权被称为香草期权。另一方面，也存在各种类型的另类期权，它们具有各种特性，使其比常见的香草期权更为复杂。

例如，如果期权买方可以在到期日前的多个时刻行使权利，那么它就是一个百慕大期权。在第十二章，蒙特卡洛模拟中，讨论了两种类型的另类期权。许多另类期权（衍生品）可能有多个触发条件与其支付有关。另类期权还可能包括非标准的基础证券或工具，专为特定客户或特定市场开发。另类期权通常是场外交易（OTC）。

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

欧式、美式和百慕大期权
简单选择期权
喊叫期权、彩虹期权和双向期权
平均价格期权
障碍期权——上涨入场期权和上涨退出期权
障碍期权——下跌入场期权和下跌退出期权

欧式、美式和百慕大期权

在第十章，期权与期货中，我们已经学习了对于欧式期权，期权买方只能在到期日行使权利，而对于美式期权买方，他们可以在到期日之前的任何时候行使权利。因此，美式期权的价值通常高于其对应的欧式期权。百慕大期权可以在几个预定的日期内行使一次或多次。因此，百慕大期权的价格应介于欧式期权和美式期权之间，前提是它们具有相同的特征，例如相同的到期日和相同的行权价格，见以下两个关于看涨期权的不等式：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_14_01.jpg

这里是一个关于百慕大期权的例子。假设一家公司发行了一张 10 年期债券。七年后，公司可以在接下来的三年中的每年年底，选择赎回即退还债券。这种可赎回特性最终形成了一个嵌入式百慕大期权，其行使日期为第 8、9 和 10 年的 12 月。

首先，我们来看一下使用二项模型的美式看涨期权的 Python 程序：

def binomialCallAmerican(s,x,T,r,sigma,n=100):
    from math import exp,sqrt
    import numpy as np
    deltaT = T /n
    u = exp(sigma * sqrt(deltaT)) 
    d = 1.0 / u
    a = exp(r * deltaT)
    p = (a - d) / (u - d)
    v = [[0.0 for j in np.arange(i + 1)] for i in np.arange(n + 1)] 
    for j in np.arange(n+1):
        v[n][j] = max(s * u**j * d**(n - j) - x, 0.0) 
    for i in np.arange(n-1, -1, -1):
        for j in np.arange(i + 1):
            v1=exp(-r*deltaT)*(p*v[i+1][j+1]+(1.0-p)*v[i+1][j]) 
            v2=max(v[i][j]-x,0)         # early exercise 
            v[i][j]=max(v1,v2)
    return v[0][0]
#
s=40\.        # stock price today 
x=40\.        # exercise price
T=6./12      # maturity date ii years
tao=1/12     # when to choose
r=0.05       # risk-free rate
sigma=0.2    # volatility 
n=1000       # number of steps
#
price=binomialCallAmerican(s,x,T,r,sigma,n)
print("American call =", price)
('American call =', 2.7549263174936502)

这个美式看涨期权的价格是$2.75。为了修改之前的程序以仅满足少数几个行权价格，以下两行是关键：

            v2=max(v[i][j]-x,0)         # early exercise 
            v[i][j]=max(v1,v2)

这是一个适用于伯穆达看涨期权的 Python 程序。关键的不同点是名为 T2 的变量，它包含伯穆达期权可以行权的日期：

def callBermudan(s,x,T,r,sigma,T2,n=100):
    from math import exp,sqrt
    import numpy as np
    n2=len(T2)
    deltaT = T /n
    u = exp(sigma * sqrt(deltaT)) 
    d = 1.0 / u
    a = exp(r * deltaT)
    p = (a - d) / (u - d)
    v =[[0.0 for j in np.arange(i + 1)] for i in np.arange(n + 1)] 
    for j in np.arange(n+1):
        v[n][j] = max(s * u**j * d**(n - j) - x, 0.0) 
    for i in np.arange(n-1, -1, -1):
        for j in np.arange(i + 1):
            v1=exp(-r*deltaT)*(p*v[i+1][j+1]+(1.0-p)*v[i+1][j])
            for k in np.arange(n2):
                if abs(j*deltaT-T2[k])<0.01:
                    v2=max(v[i][j]-x,0)  # potential early exercise 
                else: 
                    v2=0
            v[i][j]=max(v1,v2)
    return v[0][0]
#
s=40\.                 # stock price today 
x=40\.                 # exercise price
T=6./12               # maturity date ii years
r=0.05                # risk-free rate
sigma=0.2             # volatility 
n=1000                # number of steps
T2=(3./12.,4./12.)    # dates for possible early exercise 
#
price=callBermudan(s,x,T,r,sigma,T2,n)
print("Bermudan call =", price)
('Bermudan call =', 2.7549263174936502)

选择权选项

对于选择权期权，它允许期权买方在期权到期之前的某个预定时间点选择是欧洲看涨期权还是欧洲看跌期权。对于一个简单的选择权期权，标的的看涨和看跌期权具有相同的到期日和行权价格。我们来看两个极端案例。期权买方必须在今天做出决策，即当他们进行此类购买时。这个选择权期权的价格应当是看涨和看跌期权的最大值，因为期权买方没有更多的信息。第二个极端情况是投资者可以在到期日做出决策。由于看涨和看跌期权具有相同的行权价格，如果看涨期权处于价内，看跌期权应该是价外的，反之亦然。因此，选择权期权的价格应当是看涨期权和看跌期权的总和。这相当于购买一个看涨期权和一个看跌期权，行权价格和到期日相同。在第十章《期权与期货》中，我们知道这样的交易策略被称为跨式期权（Straddle）。通过这种交易策略，我们押注标的资产会偏离当前的价格，但我们不确定其方向。

首先，我们来看一个简单的选择权期权的定价公式，假设看涨和看跌期权有相同的到期日和行权价格，并假设到期之前没有股息。一个简单的选择权期权的定价公式如下：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_14_04.jpg

这里，Pchooer 是选择权的价格或溢价，call (T) 是一个到期日为 T 的欧洲看涨期权，put(τ) 将很快定义。对于第一个 call (T) 期权，我们有以下定价公式：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_14_05.jpg

这里，call (T) 是看涨期权的溢价，S 是今天的价格，K 是行权价格，T 是到期日（年），σ 是波动率，N() 是累积标准正态分布。实际上，这与 Black-Scholes-Merton 看涨期权模型完全相同。put (τ) 的公式如下：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_14_06.jpg

同样，put(τ) 是看跌期权的溢价，τ 是选择权买方可以做出决策的时刻。为了使 d1 和 d2 与前面公式中的两个值有所区分，使用了 https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_14_26.jpg 和 https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_14_27.jpg 来代替 d1 和 d2。请注意，前面的公式与 Black-Scholes-Merton 看跌期权模型不同，因为我们有 T 和 τ 而不仅仅是 T。现在，我们来看一个极端情况：期权买方可以在到期日做出决策，即 τ=T。从前面的公式可以明显看出，选择权的价格将是这两种期权的总和：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_14_07.jpg

以下 Python 程序是针对选择期权的。为了节省空间，我们可以将看涨期权和看跌期权合并，见下文 Python 代码。为此，我们有两个时间变量输入，分别是T和tao：

from scipy import log,exp,sqrt,stats 
def callAndPut(S,X,T,r,sigma,tao,type='C'):
    d1=(log(S/X)+r*T+0.5*sigma*sigma*tao)/(sigma*sqrt(tao)) 
    d2 = d1-sigma*sqrt(tao)
    if type.upper()=='C':
        c=S*stats.norm.cdf(d1)-X*exp(-r*T)*stats.norm.cdf(d2)
        return c
    else:
        p=X*exp(-r*T)*stats.norm.cdf(-d2)-S*stats.norm.cdf(-d1)
        return p
#
def chooserOption(S,X,T,r,sigma,tao):
    call_T=callAndPut(S,X,T,r,sigma,T)
    put_tao=callAndPut(S,X,T,r,sigma,tao,type='P')
    return call_T- put_tao
#
s=40\.        # stock price today 
x=40\.        # exercise price
T=6./12      # maturity date ii years
tao=1./12\.   # when to choose
r=0.05       # risk-free rate
sigma=0.2    # volatility 
#
price=chooserOption(s,x,T,r,sigma,tao)
print("price of a chooser option=",price)
('price of a chooser option=', 2.2555170735574421)

这个选择权的价格是$2.26。

呼叫期权

呼叫期权是一种标准的欧洲期权，不同之处在于，期权买方可以在到期日前对期权卖方呼叫，将最小支付设置为Sτ-X，其中Sτ是买方呼叫时的股票价格，X是行权价。行权价的设定可以与现货价格的特定关系相关，例如设定为现货价的 3%或 5%（高于或低于）。Python 代码如下：

def shoutCall(s,x,T,r,sigma,shout,n=100):
    from math import exp,sqrt
    import numpy as np
    deltaT = T /n
    u = exp(sigma * sqrt(deltaT)) 
    d = 1.0 / u
    a = exp(r * deltaT)
    p = (a - d) / (u - d)
    v =[[0.0 for j in np.arange(i + 1)] for i in np.arange(n + 1)] 
    for j in np.arange(n+1):
        v[n][j] = max(s * u**j * d**(n - j) - x, 0.0) 
    for i in np.arange(n-1, -1, -1):
        for j in np.arange(i + 1):
            v1=exp(-r*deltaT)*(p*v[i+1][j+1]+(1.0-p)*v[i+1][j]) 
            v2=max(v[i][j]-shout,0)   # shout  
            v[i][j]=max(v1,v2)
    return v[0][0]
#
s=40\.              # stock price today 
x=40\.              # exercise price
T=6./12            # maturity date ii years
tao=1/12           # when to choose
r=0.05             # risk-free rate
sigma=0.2          # volatility 
n=1000             # number of steps
shout=(1+0.03)*s   # shout out level 
#
price=shoutCall(s,x,T,r,sigma,shout,n)
print("Shout call =", price)

二元期权

二元期权，或资产无条件期权，是一种期权类型，其支付结构为：若期权到期时处于实值，则支付固定金额的补偿；若期权到期时处于虚值，则支付零。由于这个特性，我们可以应用蒙特卡洛模拟来找到解决方案。Python 代码如下：

import random
import scipy as sp
#
def terminalStockPrice(S, T,r,sigma):
    tao=random.gauss(0,1.0)
    terminalPrice=S * sp.exp((r - 0.5 * sigma**2)*T+sigma*sp.sqrt(T)*tao)
    return terminalPrice
#
def binaryCallPayoff(x, sT,payoff):
    if sT >= x:
        return payoff
    else:
        return 0.0
# input area 
S = 40.0            # asset price
x = 40.0            # exercise price 
T = 0.5             # maturity in years 
r = 0.01            # risk-free rate 
sigma = 0.2         # vol of 20%
fixedPayoff = 10.0  # payoff 
nSimulations =10000 # number of simulations 
#
payoffs=0.0
for i in xrange(nSimulations):
    sT = terminalStockPrice(S, T,r,sigma) 
    payoffs += binaryCallPayoff(x, sT,fixedPayoff)
#
price = sp.exp(-r * T) * (payoffs / float(nSimulations))
print('Binary options call= %.8f' % price)

请注意，由于前述程序未固定种子，因此每次运行时，用户应获得不同的结果。

彩虹期权

许多金融问题可以总结为或与多项资产的最大值或最小值相关。我们来看一个简单的例子：基于两个资产最大值或最小值的期权。这种类型的期权称为彩虹期权。由于涉及两项资产，我们需要熟悉所谓的二元正态分布。以下代码展示了其图形。原始代码可以在scipython.com/blog/visualizing-the-bivariate-gaussian-distribution/网站找到：

import numpy as np
from matplotlib import cm
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
#
# input area
n   = 60                      # number of intervals
x   = np.linspace(-3, 3, n)   # x dimension
y   = np.linspace(-3, 4, n)   # y dimension 
x,y = np.meshgrid(x, y)       # grid 
#
# Mean vector and covariance matrix
mu = np.array([0., 1.])
cov= np.array([[ 1\. , -0.5], [-0.5,  1.5]])
#
# combine x and y into a single 3-dimensional array
pos = np.empty(x.shape + (2,))
pos[:, :, 0] = x
pos[:, :, 1] = y
#
def multiNormal(pos, mu, cov):
    n = mu.shape[0]
    Sigma_det = np.linalg.det(cov)
    Sigma_inv = np.linalg.inv(cov)
    n2 = np.sqrt((2*np.pi)**n * Sigma_det)
    fac=np.einsum('...k,kl,...l->...', pos-mu, Sigma_inv, pos-mu)
    return np.exp(-fac/2)/n2
#
z    = multiNormal(pos, mu, cov)
fig  = plt.figure()
ax   = fig.gca(projection='3d')
ax.plot_surface(x, y, z, rstride=3, cstride=3,linewidth=1, antialiased=True,cmap=cm.viridis)
cset = ax.contourf(x, y, z, zdir='z', offset=-0.15, cmap=cm.viridis)
ax.set_zlim(-0.15,0.2)
ax.set_zticks(np.linspace(0,0.2,5))
ax.view_init(27, -21)
plt.title("Bivariate normal distribtuion")
plt.ylabel("y values ")
plt.xlabel("x values")
plt.show()

图形如下：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_14_02.jpg

假设这两项资产的收益率遵循具有相关系数ρ的二元正态分布。为了简化我们的估算，我们假设在到期日前没有分红。对于基于两个资产最小值的看涨期权，其支付为：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_14_08.jpg

这里，https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_14_24.jpg是股票 1（2）的终值股价，T 是到期日（以年为单位）。基于两项资产最小值的看涨期权定价公式如下：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_14_09.jpg

这里，S1（S2）是股票 1（2）的当前股价，N2(a,b,ρ)是具有上限 a 和 b、资产间相关系数ρ的二元正态分布，K是行权价。d11、d12、d21、d22、ρ1和ρ2的参数在此定义：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_14_10.jpg

首先，我们应学习这里描述的二元累积正态分布N2_f(d1,d2,rho)：

def N2_f(d1,d2,rho):
    """cumulative bivariate standard normal distribution 
       d1: the first value
       d2: the second value
       rho: correlation

       Example1:
               print(N2_f(0,0,1.)) => 0.5
       Example2:
               print(N2_f(0,0,0)  => 0.25
     """
    import statsmodels.sandbox.distributions.extras as extras
    muStandardNormal=0.0    # mean of a standard normal distribution 
    varStandardNormal=1.0   # variance of standard normal distribution 
    upper=([d1,d2])         # upper bound for two values
    v=varStandardNormal     # simplify our notations
    mu=muStandardNormal     # simplify our notations
    covM=([v,rho],[rho,v])
    return extras.mvnormcdf(upper,mu,covM)
#

让我们看一些特殊情况。从单变量标准正态分布中，我们知道当输入值为0时，我们期望累积标准正态分布为0.5，因为基础正态分布是对称的。当两个时间序列完全正相关时，累积标准正态分布也应该为0.5，请参见前面的结果。另一方面，如果两个时间序列不相关，当输入都为零时，它们的累积标准正态分布期望为重叠，即，0.5 * 0.5 = 0.25。这一点通过调用之前的N2_f()函数得到验证。对于这一类期权，相关的 Python 程序如下：

from math import exp,sqrt,log
import statsmodels.sandbox.distributions.extras as extras
#
def dOne(s,k,r,sigma,T):
    #print(s,k,r,sigma,T)
    a=log(s/k)+(r-0.5*sigma**2)*T
    b=(sigma*sqrt(T))
    return a/b
#
def sigmaA_f(sigma1,sigma2,rho):
    return sqrt(sigma1**2-2*rho*sigma1*sigma2+sigma2**2)
#
def dTwo(d1,sigma,T):
    return d1+sigma*sqrt(T)
#
def rhoTwo(sigma1,sigma2,sigmaA,rho):
    return (rho*sigma2-sigma1)/sigmaA
#
def N2_f(d1,d2,rho):
    import statsmodels.sandbox.distributions.extras as extras
    muStandardNormal=0.0    # mean of a standard normal distribution 
    varStandardNormal=1.0   # variance of standard normal distribution 
    upper=([d1,d2])         # upper bound for two values
    v=varStandardNormal     # simplify our notations
    mu=muStandardNormal     # simplify our notations
    covM=([v,rho],[rho,v])
    return extras.mvnormcdf(upper,mu,covM)
#
def dOneTwo(s1,s2,sigmaA,T):
    a=log(s2/s1)-0.5*sigmaA**2*T
    b=sigmaA*sqrt(T)
    return a/b
#
def rainbowCallOnMinimum(s1,s2,k,T,r,sigma1,sigma2,rho):
    d1=dOne(s1,k,r,sigma1,T)
    d2=dOne(s2,k,r,sigma2,T)
    d11=dTwo(d1,sigma1,T)
    d22=dTwo(d2,sigma2,T)
    sigmaA=sigmaA_f(sigma1,sigma2,rho)
    rho1=rhoTwo(sigma1,sigma2,sigmaA,rho)
    rho2=rhoTwo(sigma2,sigma1,sigmaA,rho)
    d12=dOneTwo(s1,s2,sigmaA,T)
    d21=dOneTwo(s2,s1,sigmaA,T)
    #
    part1=s1*N2_f(d11,d12,rho1)
    part2=s2*N2_f(d21,d22,rho2)
    part3=k*exp(-r*T)*N2_f(d1,d2,rho)
    return part1 + part2 - part3
#
s1=100.
s2=95.
k=102.0
T=8./12.
r=0.08
rho=0.75
sigma1=0.15
sigma2=0.20
price=rainbowCallOnMinimum(s1,s2,k,T,r,sigma1,sigma2,rho)
print("price of call based on the minimum of 2 assets=",price)
('price of call based on the minimum of 2 assets=', 3.747423936156629)

另一种定价各种类型彩虹期权的方法是使用蒙特卡洛模拟。如我们在第十二章中提到的，蒙特卡洛模拟，我们可以生成两个相关的随机数时间序列。这个过程分为两步：生成两个零相关的随机时间序列x1和x2；然后应用以下公式：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_14_11.jpg

在这里，ρ是这两个时间序列之间的预定相关性。现在，y1和y2具有预定的相关性。以下 Python 程序将实现上述方法：

import scipy as sp
sp.random.seed(123)
n=1000
rho=0.3
x1=sp.random.normal(size=n)
x2=sp.random.normal(size=n)
y1=x1
y2=rho*x1+sp.sqrt(1-rho**2)*x2
print(sp.corrcoef(y1,y2))
[[ 1\.          0.28505213]
 [ 0.28505213  1\.        ]]

接下来，我们应用我们在第十二章中所学的同样技术，蒙特卡洛模拟，来定价两个资产最小值的彩虹期权看涨：

import scipy as sp 
from scipy import zeros, sqrt, shape 
#
sp.random.seed(123)  # fix our random numbers
s1=100\.              # stock price 1 
s2=95\.               # stock price 2
k=102.0              # exercise price
T=8./12\.             # maturity in years
r=0.08               # risk-free rate
rho=0.75             # correlation between 2
sigma1=0.15          # volatility for stock 1
sigma2=0.20          # volatility for stock 1
nSteps=100\.          # number of steps 
nSimulation=1000     # number of simulations 
#
# step 1: generate correlated random number
dt =T/nSteps 
call = sp.zeros([nSimulation], dtype=float) 
x = range(0, int(nSteps), 1) 
#
# step 2: call call prices 
for j in range(0, nSimulation): 
    x1=sp.random.normal(size=nSimulation)
    x2=sp.random.normal(size=nSimulation)
    y1=x1
    y2=rho*x1+sp.sqrt(1-rho**2)*x2
    sT1=s1
    sT2=s2 
    for i in x[:-1]: 
        e1=y1[i]
        e2=y2[i]
        sT1*=sp.exp((r-0.5*sigma1**2)*dt+sigma1*e1*sqrt(dt)) 
        sT2*=sp.exp((r-0.5*sigma2**2)*dt+sigma2*e2*sqrt(dt)) 
        minOf2=min(sT1,sT2)
        call[j]=max(minOf2-k,0) 
#
# Step 3: summation and discount back 
call=sp.mean(call)*sp.exp(-r*T) 
print('Rainbow call on minimum of 2 assets = ', round(call,3))
('Rainbow call on minimum of 2 assets = ', 4.127)

如果我们增加更多的资产，获得封闭式解将变得更加困难。这里我们展示如何使用蒙特卡洛模拟来定价基于最大终端股价的彩虹看涨期权。基本逻辑相当简单：生成三个终端股价，然后应用以下公式计算看涨期权收益：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_14_12.jpg

最终价格将是折现收益的平均值。关键是如何生成一组三个相关的随机数。这里，我们应用了著名的乔尔斯基分解。假设我们有一个相关矩阵叫做C，一个使得https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_14_25.jpg的乔尔斯基分解矩阵L。进一步假设无关的回报矩阵叫做U。现在，相关回报矩阵R = UL。以下是对应的 Python 代码：

import numpy as np
# input area
nSimulation=5000              # number of simulations
c=np.array([[1.0, 0.5, 0.3],  # correlation matrix
            [0.5, 1.0, 0.4],
            [0.3, 0.4, 1.0]])
np.random.seed(123)           # fix random numbers 
#
# generate uncorrelated random numbers
x=np.random.normal(size=3*nSimulation)
U=np.reshape(x,(nSimulation,3))
#
# Cholesky decomposition 
L=np.linalg.cholesky(c)
# generate correlated random numbers
r=np.dot(U,L)
#check the correlation matrix
print(np.corrcoef(r.T))
[[ 1\.          0.51826188  0.2760649 ]
 [ 0.51826188  1\.          0.35452286]
 [ 0.2760649   0.35452286  1\.        ]]

定价平均期权

在第十二章，蒙特卡罗模拟中，我们讨论了两种奇异期权。为了方便起见，我们也将它们包括在这一章中。因此，读者可能会看到一些重复的内容。欧洲期权和美式期权是路径无关的期权。这意味着期权的支付仅取决于终端股票价格和执行价格。路径相关期权的一个相关问题是到期日的市场操控。另一个问题是，一些投资者或对冲者可能更关注平均价格而非终端价格。

例如，一家炼油厂担心其主要原料——石油的价格波动，尤其是在接下来的三个月内。他们计划对冲原油价格的潜在跳动。该公司可能会购买看涨期权。然而，由于公司每天消耗大量的原油，它显然更关心的是平均价格，而不是仅仅依赖于普通看涨期权的终端价格。在这种情况下，平均期权会更加有效。平均期权是一种亚洲期权。对于平均期权，其支付取决于某一预定时间段内的基础资产价格的平均值。平均数有两种类型：算术平均数和几何平均数。以下给出了亚洲看涨期权（平均价格）的支付函数：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_14_13.jpg

这里给出了亚洲看跌期权（平均价格）的支付函数：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_14_14.jpg

亚洲期权是奇异期权的一种基本形式。亚洲期权的另一个优点是，相比于欧洲和美式的普通期权，它们的成本更低，因为平均价格的波动会比终端价格的小得多。以下 Python 程序适用于一个具有算术平均价格的亚洲期权：

import scipy as sp 
s0=30\.                 # today stock price 
x=32\.                  # exercise price 
T=3.0/12\.              # maturity in years 
r=0.025                # risk-free rate 
sigma=0.18             # volatility (annualized) 
sp.random.seed(123)    # fix a seed here 
n_simulation=1000      # number of simulations 
n_steps=500\.           # number of steps
#
dt=T/n_steps 
call=sp.zeros([n_simulation], dtype=float) 
for j in range(0, n_simulation): 
    sT=s0 
    total=0 
    for i in range(0,int(n_steps)): 
         e=sp.random.normal()
         sT*=sp.exp((r-0.5*sigma*sigma)*dt+sigma*e*sp.sqrt(dt)) 
         total+=sT 
         price_average=total/n_steps 
    call[j]=max(price_average-x,0) 
#
call_price=sp.mean(call)*sp.exp(-r*T) 
print('call price based on average price = ', round(call_price,3))
('call price based on average price = ', 0.12)

定价障碍期权

不同于 Black-Scholes-Merton 期权模型中的看涨和看跌期权，这些期权是路径无关的，障碍期权是路径相关的。障碍期权与普通期权在许多方面相似，唯一的区别是存在一个触发条件。敲入期权在基础资产触及预定的敲入障碍之前是无效的。相反，敲出障碍期权在其生命周期一开始是有效的，只有当价格突破敲出障碍价时，它才变得无效。此外，如果障碍期权到期时未激活，它可能一文不值，或者会按溢价的一定比例支付现金回扣。四种类型的障碍期权如下：

上涨敲出：在这种障碍期权中，价格从障碍下方开始。如果价格触及障碍，则会被敲出。
下跌敲出：在这种障碍期权中，价格从障碍上方开始。如果价格触及障碍，则会被敲出。
敲入上涨：在这种障碍期权中，价格从障碍下方开始，必须触及障碍才能激活。
下行入场：在这种障碍期权中，价格始于一个障碍点，并且必须达到该障碍点才能被激活。

以下 Python 程序适用于一个带有欧洲看涨期权的上行出场障碍期权：

import scipy as sp 
from scipy import log,exp,sqrt,stats 
#
def bsCall(S,X,T,r,sigma):
    d1=(log(S/X)+(r+sigma*sigma/2.)*T)/(sigma*sqrt(T)) 
    d2 = d1-sigma*sqrt(T)
    return S*stats.norm.cdf(d1)-X*exp(-r*T)*stats.norm.cdf(d2)
#
def up_and_out_call(s0,x,T,r,sigma,n_simulation,barrier):
    n_steps=100\. 
    dt=T/n_steps 
    total=0 
    for j in sp.arange(0, n_simulation): 
        sT=s0 
        out=False
        for i in range(0,int(n_steps)): 
            e=sp.random.normal() 
            sT*=sp.exp((r-0.5*sigma*sigma)*dt+sigma*e*sp.sqrt(dt)) 
            if sT>barrier: 
               out=True 
        if out==False: 
            total+=bsCall(s0,x,T,r,sigma) 
    return total/n_simulation 
#

基本设计是模拟股票运动 n 次，例如 100 次。对于每次模拟，我们有 100 步。每当股票价格达到障碍点时，收益为零。否则，收益将是一个普通的欧洲看涨期权。最终的值将是所有未被打击的看涨期权价格的总和，再除以模拟次数，如下代码所示：

s0=30\.              # today stock price 
x=30\.               # exercise price 
barrier=32          # barrier level 
T=6./12\.            # maturity in years 
r=0.05              # risk-free rate 
sigma=0.2           # volatility (annualized) 
n_simulation=100    # number of simulations 
sp.random.seed(12)  # fix a seed
#
result=up_and_out_call(s0,x,T,r,sigma,n_simulation,barrier) 
print('up-and-out-call = ', round(result,3))
('up-and-out-call = ', 0.93)

下述是下行入场看跌期权的 Python 代码：

def down_and_in_put(s0,x,T,r,sigma,n_simulation,barrier): 
    n_steps=100.
    dt=T/n_steps 
    total=0
    for j in range(0, n_simulation): 
        sT=s0
        in_=False
        for i in range(0,int(n_steps)): 
            e=sp.random.normal()
            sT*=sp.exp((r-0.5*sigma*sigma)*dt+sigma*e*sp.sqrt(dt)) 
            if sT<barrier:
                in_=True
            #print 'sT=',sT
            #print 'j=',j ,'out=',out if in_==True:
            total+=p4f.bs_put(s0,x,T,r,sigma) 
    return total/n_simulation
#

进出障碍平价

如果我们购买一个上行出场的欧洲看涨期权和一个上行入场的欧洲看涨期权，则应满足以下平价关系：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_14_15.jpg

逻辑非常简单——如果股票价格触及障碍点，则第一个看涨期权失效，第二个看涨期权被激活。如果股票价格从未触及障碍点，第一个看涨期权保持有效，而第二个则永远不会被激活。无论哪种情况，其中一个期权是有效的。以下 Python 程序演示了这种情境：

def upCall(s,x,T,r,sigma,nSimulation,barrier):
    import scipy as sp
    import p4f 
    n_steps=100
    dt=T/n_steps 
    inTotal=0 
    outTotal=0
    for j in range(0, nSimulation): 
        sT=s
        inStatus=False 
        outStatus=True
        for i in range(0,int(n_steps)):
            e=sp.random.normal()
            sT*=sp.exp((r-0.5*sigma*sigma)*dt+sigma*e*sp.sqrt(dt)) 
            if sT>barrier:
                outStatus=False 
                inStatus=True
        if outStatus==True:
            outTotal+=p4f.bs_call(s,x,T,r,sigma) 
        else:
            inTotal+=p4f.bs_call(s,x,T,r,sigma) 
    return outTotal/nSimulation, inTotal/nSimulation
#

我们输入一组值来测试一个上行出场看涨期权与一个上行入场看涨期权的总和是否与一个普通看涨期权相同：

import p4f
s=40\.                 # today stock price 
x=40\.                 # exercise price 
barrier=42.0          # barrier level 
T=0.5                 # maturity in years 
r=0.05                # risk-free rate 
sigma=0.2             # volatility (annualized) 
nSimulation=500       # number of simulations 
#
upOutCall,upInCall=upCall(s,x,T,r,sigma,nSimulation,barrier) 
print 'upOutCall=', round(upOutCall,2),'upInCall=',round(upInCall,2) 
print 'Black-Scholes call', round(p4f.bs_call(s,x,T,r,sigma),2)

上行出场与上行入场平价图

使用蒙特卡罗模拟来呈现这种平价是一个不错的主意。以下代码旨在实现这一点。为了让我们的模拟更清晰，我们故意选择了仅进行五次模拟：

import p4f
import scipy as sp
import matplotlib.pyplot as plt
#
s =9.25              # stock price at time zero
x =9.10              # exercise price
barrier=10.5         # barrier
T =0.5               # maturity date (in years)
n_steps=30           # number of steps
r =0.05              # expected annual return
sigma = 0.2          # volatility (annualized) 
sp.random.seed(125)  # seed()
n_simulation = 5     # number of simulations 
#
dt =T/n_steps
S = sp.zeros([n_steps], dtype=float) 
time_= range(0, int(n_steps), 1) 
c=p4f.bs_call(s,x,T,r,sigma) 
sp.random.seed(124)
outTotal, inTotal= 0.,0\. 
n_out,n_in=0,0

for j in range(0, n_simulation):
    S[0]= s
    inStatus=False
    outStatus=True
    for i in time_[:-1]:
        e=sp.random.normal()
        S[i+1]=S[i]*sp.exp((r-0.5*pow(sigma,2))*dt+sigma*sp.sqrt(dt)*e) 
        if S[i+1]>barrier:
            outStatus=False 
            inStatus=True
    plt.plot(time_, S) 
    if outStatus==True:
        outTotal+=c;n_out+=1 
    else:
        inTotal+=c;n_in+=1 
        S=sp.zeros(int(n_steps))+barrier 
        plt.plot(time_,S,'.-') 
        upOutCall=round(outTotal/n_simulation,3) 
        upInCall=round(inTotal/n_simulation,3) 
        plt.figtext(0.15,0.8,'S='+str(s)+',X='+str(x))
        plt.figtext(0.15,0.76,'T='+str(T)+',r='+str(r)+',sigma=='+str(sigma)) 
        plt.figtext(0.15,0.6,'barrier='+str(barrier))
        plt.figtext(0.40,0.86, 'call price  ='+str(round(c,3)))
        plt.figtext(0.40,0.83,'up_and_out_call ='+str(upOutCall)+' (='+str(n_out)+'/'+str(n_simulation)+'*'+str(round(c,3))+')') 
        plt.figtext(0.40,0.80,'up_and_in_call ='+str(upInCall)+' (='+str(n_in)+'/'+ str(n_simulation)+'*'+str(round(c,3))+')')
#
plt.title('Up-and-out and up-and-in parity (# of simulations = %d ' % n_simulation +')')
plt.xlabel('Total number of steps ='+str(int(n_steps))) 
plt.ylabel('stock price')
plt.show()

对应的图形如下所示。请注意，在前面的程序中，由于使用了种子，若使用相同的种子，不同的用户应该得到相同的图形：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_14_03.jpg

使用浮动行权价定价回望期权

回望期权取决于标的资产所经历的路径（历史）。因此，它们也被称为路径依赖的奇异期权。其中一种被称为浮动行权价。当行使价格为期权生命周期内达到的最低价格时，给定的看涨期权的收益函数如下所示：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_14_16.jpg

这是该回望期权的 Python 代码：

plt.show()
def lookback_min_price_as_strike(s,T,r,sigma,n_simulation): 
    n_steps=100
    dt=T/n_steps
    total=0
    for j in range(n_simulation): 
        min_price=100000\.  # a very big number 
        sT=s
        for i in range(int(n_steps)): 
            e=sp.random.normal()
            sT*=sp.exp((r-0.5*sigma*sigma)*dt+sigma*e*sp.sqrt(dt)) 
            if sT<min_price:
                min_price=sT
                #print 'j=',j,'i=',i,'total=',total 
                total+=p4f.bs_call(s,min_price,T,r,sigma)
    return total/n_simulation

记住，前面的函数需要两个模块。因此，我们必须在调用函数之前导入这些模块，如下代码所示：

import scipy as sp
import p4f
s=40\.             # today stock price
T=0.5               # maturity in years
r=0.05              # risk-free rate
sigma=0.2           # volatility (annualized)
n_simulation=1000   # number of simulations
result=lookback_min_price_as_strike(s,T,r,sigma,n_simulation)
print('lookback min price as strike = ', round(result,3))

一次运行的结果如下所示：

('lookback min price as strike = ', 53.31)t(

附录 A – 数据案例 7 – 对冲原油

假设一家炼油厂每天使用原油。因此，他们必须面对主要原材料——原油价格的不确定性风险。保护底线和确保生产顺利之间存在权衡；公司研究所有可能的结果，例如是否对冲油价，或者完全不对冲。假设总的年原油消耗量为 2000 万加仑。再次强调，公司每天都必须运营。比较以下几种策略，并指出其优缺点：

不进行对冲
使用期货
使用期权
使用外汇期权

存在几种策略，例如美国式期权；请参见下表中的规格。以下表格展示了部分原油期权合约的规格：

合约单位	在交易所交易的轻质甜油看跌（看涨）期权代表着在交易所交易的轻质甜油期货中持有空头（多头）仓位的期权。
最小价格波动	每桶$0.01。
价格报价	每桶以美元和美分计。
产品代码	CME Globex: LO, CME ClearPort: LO, 清算: LO。
上市合约	当前年份和接下来五个日历年份的每月合约，以及三个额外年份的 6 月和 12 月合约。新日历年的剩余每月合约将在当前年 12 月合约交易终止后加入。
交易终止	交易在相关期货合约的交易终止前三个工作日终止。
交易风格	美国式。
结算方式	可交割。
标的	轻质甜油期货。

表 1：原油期权合约的一些规格

如果我们使用期货进行对冲，我们有以下公式：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_14_17.jpg

N 是期货合约的数量，VA 是我们的投资组合价值（我们想要对冲的金额），β 是基于我们材料和标的工具的回归斜率（注意如果我们的材料与标的对冲工具相同，那么 β 等于 1），VF 是一个期货合约的价值：

来源: www.cmegroup.com/trading/energy/crude-oil/light-sweet-crude_contractSpecs_options.html?gclid=CLjWq92Yx9ICFU1MDQodP5EDLg&gclsrc=aw.ds
数据来源: 原油价格：西德克萨斯中质油（WTI）- 俄克拉荷马州库欣（DCOILWTICO），fred.stlouisfed.org/series/DCOILWTICO/downloaddata

一个相关的数据集叫做crudeOilPriceDaily.pkl。这里显示的是该数据集的前几条和最后几条观测值。该数据集可以在canisius.edu/~yany/python/crudeOilPriceDaily.pkl下载：

import scipy as sp
import pandas as pd
x=pd.read_pickle("c:/temp/cruideOilPriceDaily.pkl")
print(x.head())
print(x.tail())
            PRICE
1986-01-02  25.56
1986-01-03  26.00
1986-01-06  26.53
1986-01-07  25.85
1986-01-08  25.87
            PRICE
2017-02-28  54.00
2017-03-01  53.82
2017-03-02  52.63
2017-03-03  53.33
2017-03-06  53.19

参考文献

Clewlow，Les 和 Chris Strickland，1997 年，《异型期权，技术前沿》，Thomaston Business Press
Kurtverstegen，《模拟：模拟无相关和相关的随机变量》， kurtverstegen.wordpress.com/2013/12/07/simulation/
张，彼得，1998 年，《异型期权》，世界科学出版社，第 2 版

练习

异型期权的定义是什么？
为什么有人说可赎回债券相当于普通债券加上一种 Bermudan 期权（发行公司是这个 Bermudan 期权的买方，而债券买方是卖方）？
编写一个 Python 程序，根据算术平均数定价一个亚洲平均价格看跌期权。
编写一个 Python 程序，根据几何平均数定价一个亚洲平均价格看跌期权。
编写一个 Python 程序，定价一个向上障碍期权（up-and-in call）。
编写一个 Python 程序，定价一个向下障碍期权（down-and-out put）。
编写一个 Python 程序，展示向下障碍期权和向下敲入期权的平价关系。
编写一个 Python 程序，使用 SciPy 中的permutation()函数，从过去五年的数据中随机选择 12 个月的回报，并且不放回。测试程序时，可以使用花旗银行的数据，时间段为 2009 年 1 月 1 日到 2014 年 12 月 31 日，数据来自 Yahoo Finance。
编写一个 Python 程序，使用 n 个给定的回报进行自助法。每次选择 m 个回报，其中 m>n。
在本章中，我们学习了一个简单的 chooser 期权的定价公式：https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_14_18.jpg

这里，T 是到期日（以年为单位），τ是期权决定选择看涨还是看跌的时间。是否可以使用以下公式？

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_14_19.jpg
当股票支付连续复利的股息，股息收益率为δ时，我们有以下的 Chooser 期权定价公式：https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_14_20.jpg

其中Pchooser是 chooser 期权的价格或溢价，*call (T)*是到期时间为 T 的欧洲看涨期权，*put(τ)将在后续定义。对于第一个call (T)*期权，我们有以下定价公式：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_14_21.jpg

其中*call(T)*是看涨期权价格或期权溢价，S是今天的价格，K是执行价格，T 是期权到期时间（以年为单位），σ是波动率，*N()*是累积分布函数。实际上，这正是 Black-Scholes-Merton 看涨期权模型。看跌期权(τ)的公式如下：

https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_14_22.jpg

编写一个相关的 Python 程序。
如果今天两只股票的价格分别是 $40 和 $55，这两只股票的收益率标准差分别为 0.1 和 0.2。它们的相关系数为 0.45。基于这两只股票终端股价的最大值，彩虹看涨期权的价格是多少？行权价为 $60，期限为六个月，无风险利率为 4.5%。

解释单变量累积标准正态分布与双变量累积标准正态分布之间的异同。对于单变量累积标准正态分布 N_f() 和双变量累积标准正态分布 N2_f()，我们有以下代码：

def N_f(x):
    from scipy import stats
    return stats.norm.cdf(x)
#
def N2_f(x,y,rho):
    import statsmodels.sandbox.distributions.extras as extras
    muStandardNormal=0.0    # mean of a standard normal distribution 
    varStandardNormal=1.0   # variance of standard normal distribution 
    upper=([x,y])           # upper bound for two values
    v=varStandardNormal     # simplify our notations
    mu=muStandardNormal     # simplify our notations
    covM=([v,rho],[rho,v])
return extras.mvnormcdf(upper,mu,covM)

编写一个 Python 程序来定价基于两只相关资产终端价格最大值的看涨期权：https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/py-fin-2e/img/B06175_14_23.jpg

注意

S1、S2、d1、d2、d11、d12、d21、d22 和 N2() 函数的定义在本章中有说明。
基于蒙特卡洛模拟，编写一个 Python 程序来定价两只相关资产的最小终端价格的卖出期权。
在本章中，涉及美式期权和百慕大期权的两个程序，其输入集为 s=40、x=40、T=6./12、r=0.05、sigma=0.2、n=1000、T2=(3./12.,4./1)；对于潜在的提前行使日期，提供相同的结果。为什么？
编写一个 Python 程序来定价百慕大卖出期权。
编写一个 Python 程序来定价基于五个资产最小终端价格的彩虹看涨期权。

总结

我们在第十章中讨论的期权，期权与期货 通常称为普通期权，它们有封闭形式的解，即 Black-Scholes-Merton 期权模型。除了这些普通期权外，还有许多异域期权存在。在本章中，我们讨论了几种类型的异域期权，如百慕大期权、简单选择期权、喊叫期权和二元期权、平均价格期权、上进期权、上出期权、下进期权和下出期权。对于欧式看涨期权，期权买方只能在到期日行使权利，而对于美式期权买方，他们可以在到期日前或到期日当天随时行使权利。百慕大期权可以在到期前的几个时间点行使。

在下一章，我们将讨论各种波动性度量方法，例如我们常用的标准差，下偏标准差（LPSD）。将收益的标准差作为风险度量依据了一个关键假设，即股票收益服从正态分布。基于这一点，我们引入了几种正态性检验。此外，我们还通过图示展示了波动性聚集现象——高波动性通常会跟随一个高波动期，而低波动性则通常会跟随一个低波动期。为了解决这一现象，Angel（1982）提出了自回归条件异方差（ARCH）过程，Bollerslev（1986）则提出了广义自回归条件异方差（GARCH）过程，这些是 ARCH 过程的扩展。它们的图示及相关的 Python 程序将在下一章中讨论。