Python金融分析（七）：常用数学工具

发表于 分类于 阅读次数： Disqus：

引言

每年的9月份，杨振宁会为清华的大一新生上一堂课，回顾他波澜壮阔的一生，回味物理学尤其是高能物理学的黄金岁月，然后告诫大家，个人的奋斗必须要与历史的进程结合起来。

事实的确如此，在经历了20世纪50~70年代的黄金时代以后，物理学家纷纷投身到金融行业，即“Rocket Scientists on Wall Street”，强大的智力支持确实让现代金融在数理基础上实现了腾飞。

接下来，我们将介绍一些比较常用的金融数学方法，以及这些方法的Python实现：

  • 近似
  • 优化
  • 积分
  • 符号计算

近似

回归

当涉及到函数逼近时，回归是一种非常有效的工具，不仅适用于一维函数，对于高维情况同样有效。回归的本质就是辨识，也就是用一个模型来近似实际的系统，如果用优化的形式来表达，就是极小化系统观测值与模型近似值之间的误差。

max⁡a1,…,aD1I∑i=1I(yi−∑d=1Dad⋅bd(xi))2(1)\max \limits_{a_1,\ldots,a_D}\frac{1}{I}\sum_{i=1}^{I}\left(y_i-\sum\limits_{d=1}^{D}a_d\cdot b_d(x_i)\right)^2 \tag{1}

既然是模型，那么模型的结构也就至关重要。这里模型的结构可以认为由两部分:

  • b(x)b(x)的形式，也就是基函数的类型；
  • DD的数目，也就是模型的阶次。

模型的结果与系统实际的结构一致当然是最理想的，但是很多时候，我们并不知道系统究竟具有怎样的结构，往往我们也不需要知道。但是，为了达到近似的目的，选择一些具有代表性的模型还是非常有必要的。

多项式基函数

正如我们在信号处理里面选择正余弦函数作为基函数一样，多项式回归采用幂函数系作为基函数。在Python中，可以用Numpy中的 polyfit()和polyval()估计参数并进行拟合。下面我们通过例子来说明：

1
2
3
4
5
6
7
8
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib as mpl

plt.style.use('ggplot')
mpl.rcParams['figure.figsize']=(15, 9)
%matplotlib inline
1
2
3
4
5
6
7
def create_plot(x, y, colors, linestyles,labels, axlabels):
    plt.figure(figsize=(10, 6))
    for i in range(len(x)):
        plt.plot(x[i], y[i], color=colors[i], linestyle=linestyles[i],label=labels[i])
        plt.xlabel(axlabels[0])
        plt.ylabel(axlabels[1])
    plt.legend(loc=0)

我们以f(x)=sin(x)+0.5xf(x)=sin(x)+0.5x作为例子：

1
2
def f(x):
    return np.sin(x) + 0.5 * x

确定使用多项式进行回归以后，我们还需要确定多项式的阶次，可以发现选用的阶次不同，拟合效果也显著不同：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
models = []
x = np.linspace(-2 * np.pi, 2 * np.pi, 50)
y = [f(x)]
legend = ['f(x)']
# cm = mpl.cm.Set1.colors
color_cycle = list(plt.rcParams['axes.prop_cycle'])
cm = [color_cycle[0]['color']]
ls = ['-']
for i in range(1,6):
    models.append(np.polyfit(x, f(x), deg=i, full=True))
    y.append(np.polyval(models[i-1][0],x))
    legend.append(f'order={i}')
    cm.append(color_cycle[i]['color'])
    ls.append('--')
create_plot([x]*6, y, cm[:6], ls, legend, ['x', 'f(x)'])

output_12_0

显然，模型的阶次越高，模型的参数越多，拟合的结果就越好。但是，实际的系统并没有这么多的参数，从这个意义上来讲，只用一次多项式去逼近的时候，反而是贴近真实的情况（反映了长期的趋势）。这也是系统辨识有意思的地方。

自定义基底

在numpy中除了使用幂函数基底去做回归以外，我们还可以自定义一组基，然后用np.linalg.lstsq()方法进行参数估计。由于回归十分常见，也就不多做介绍。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# 定义由幂函数和正弦函数组成的基
def func0(x):
    return x-x+1
def func1(x):
    return x
def func2(x):
    return x**2
def func3(x):
    return np.sin(x)
basis_func = [func0, func1, func2, func3]
variables = np.asarray([f(x) for f in basis_func])

# 线性回归
reg = np.linalg.lstsq(variables.T, f(x), rcond=None)[0]
reg.round(4)
ry = np.dot(reg, variables)
create_plot([x, x], [f(x), ry], ['b', 'r'], ['-', '--'],
            ['f(x)', 'regression'], ['x', 'f(x)'])

output_16_0

插值

插值与拟合有很多相似的地方，一个显著的区别是，插值的曲线必定经过样本点，而拟合只需要接近即可，从某种程度上来说，我们可以认为插值是加了过样本点约束的拟合。

插值的方法也有很多种，线性插值、三次样条插值、B样条插值等等……其中金融中经常用的是样条插值。样条插值，是指用“样条”的特殊分段多项式进行插值的形式。由于样条插值可以使用低阶多项式样条实现较小的插值误差，这样就避免了使用高阶多项式所出现的“龙格现象”。线性插值实际上就是一次样条插值，而最常用的当属三次样条插值。三次样条插值之所以常用，主要是因为这样的道德插值曲线，具有二阶可导的性质。正可谓是：“三生万物”。

在Python中，scipy包提供了插值模块interpolate，其中最常用到的两个方法是splrep()和splev()，前者用来估计插值函数的参数，后者用来计算插值的点。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
import scipy.interpolate as spi

xd = np.linspace(-2*np.pi, 2 * np.pi, 6)
models = []
y = [f(x)]
legend = ['f(x)']
# cm = mpl.cm.Set1.colors
color_cycle = list(plt.rcParams['axes.prop_cycle'])
cm = [color_cycle[0]['color']]
ls = ['-']
for i in range(1,4):
    models.append(spi.splrep(xd, f(xd), k=i))
    y.append(spi.splev(x, models[i-1]))
    legend.append(f'interploation by {i}')
    cm.append(color_cycle[i]['color'])
    ls.append('--')
create_plot([x]*4, y, cm[:4], ls, legend, ['x', 'f(x)'])

output_19_0

可以看到，三次样条已经取得了比较不错的结果。需要注意的是，在可以应用样条插值的情况下，与最小二乘回归方法相比，可以预期更好的近似结果。 但是，请记住，需要排序（和“非噪声”）数据，并且该方法仅限于低维问题。 它在计算上也要求更高，因此在某些用例中可能比回归更长（更长）

此外，我试验了一下，如果调整控制点的个数，在本例里面，当点数小于6的时候会出现比较有意思的现象。反过来，我又想到了如果是做趋势分析的话，对于分段有没有什么启发意义？（每一段可以用一个三次样条的系数进行表征）

优化

在Python中，scipy包的optimize模块提供了优化所需要的方法。根据是否能够找到全局最优解，优化方法可以分为全局优化和局部优化方法，实际上全局优化是非常困难的事情，局部优化则相对容易得多。

optimize

针对下面的例子，分别用全局和局部方法来寻优：

1
2
3
4
5
6
7
8
def fm(p):
    x, y = p
    return np.sin(x) + 0.05 * x ** 2 + np.sin(y) + 0.05 * y ** 2

x = np.linspace(-10, 10, 100)
y = np.linspace(-10, 10, 100)
X, Y = np.meshgrid(x,y)
Z = fm((X, Y))
1
2
3
4
5
6
7
8
9
10
11
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(12, 6))
ax = Axes3D(fig)
surf = ax.plot_surface(X, Y, Z, rstride=2, cstride=2,
                                cmap='coolwarm', linewidth=0.5,
                                antialiased=True)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('f(x, y)')
fig.colorbar(surf, shrink=0.5, aspect=5)
plt.show()

output_24_0

全局优化

在scipy.optimize中，调用brute()函数可以通过暴力搜索的方式寻找全局最优解：

1
2
3
4
import scipy.optimize as sco
opt1 = sco.brute(fm, ((-10, 10.1, 0.02), (-10, 10.1, 0.02)), finish=None)
print(opt1)
print(fm(opt1))
[-1.42 -1.42]
-1.7756635257034394

局部优化

在全局搜索的过程中，由于暴力搜索计算复杂度较高，所以优化的尺度也比较粗，为了获得更加精确的结果，通过局部优化进一步提高精度。在scipy.optimize中，调用fmin()，将最小化函数和起始参数值作为输入。 可选参数值是输入参数容差和功能值容差，以及最大迭代次数和函数调用次数。

1
sco.fmin(fm, opt1, xtol=0.001, ftol=0.001, maxiter=100, maxfun=100)
Optimization terminated successfully.
         Current function value: -1.775726
         Iterations: 15
         Function evaluations: 31

array([-1.42737105, -1.42763835])

但是，如果直接使用局部优化算法的话，则可能会陷入到局部极小当中：

1
sco.fmin(fm, (1.5, 1.5), maxiter=100)
Optimization terminated successfully.
         Current function value: -0.879950
         Iterations: 51
         Function evaluations: 98

array([ 4.27110118, -1.4275218 ])

带约束的优化

我们前面提到的都是无约束的优化，如果考虑资源的稀缺性，也就变成了带约束的优化问题，也就是运筹的由来。

考虑投资中的期望效用最大化问题，现有两个地方u和c，两种投资产品a和b，其价格和收益如下表：

项目 u地 u地 价格
a 15 5 10
b 5 12 10

假设投资者目前共有资金100，出于风险分散和边际递减的角度，其效用函数为各地收益的平方根之和，综上可以得到优化问题为：

max⁡a,bE(u(w1))=pw1u+(1−p)w1dw1u=15a+5bw1d=5a+12bsubject to100≥10a+10ba,b≥0(2)\begin{aligned} \max \limits_{a,b}E(u(w_1)) &= p\sqrt{w_{1u}}+(1-p)\sqrt{w_{1d}}\\ w_{1u} &= 15a+5b\\ w_{1d} &= 5a+12b\\ \text{subject to}\\ 100 &\geq10a+10b\\ a,b&\geq0\\ \end{aligned} \tag{2}

Python代码具体实现如下：

1
2
3
4
5
6
7
8
9
10
11
import math
def Eu(p):
    a, b = p
    return -(0.5 * math.sqrt(a*15 + b*5)+
            0.5 * math.sqrt(a*5+ b*12))
constraints = ({'type': 'ineq',
                'fun': lambda p: 100 - p[0] * 10 - p[1]*10})
bounds = ((0,1000), (0,1000))
result = sco.minimize(Eu, [10,10], method='SLSQP',
                     bounds=bounds, constraints=constraints)
result
     fun: -9.700883611651133
     jac: array([-0.48507607, -0.48491502])
 message: 'Optimization terminated successfully.'
    nfev: 25
     nit: 6
    njev: 6
  status: 0
 success: True
       x: array([8.02543856, 1.97456144])

积分

我们用最开始回归的例子去做积分，首先看一下积分的区域：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
from matplotlib.patches import Polygon

x = np.linspace(0, 10)
y = f(x)
a = 0.5  
b = 9.5  
Ix = np.linspace(a, b)  
Iy = f(Ix)

fig, ax = plt.subplots(figsize=(10, 6))
plt.plot(x, y, 'b', linewidth=2)
plt.ylim(bottom=0)
Ix = np.linspace(a, b)
Iy = f(Ix)
verts = [(a, 0)] + list(zip(Ix, Iy)) + [(b, 0)]
poly = Polygon(verts, facecolor='0.7', edgecolor='0.5')
ax.add_patch(poly)
plt.text(0.7 * (a + b), 1.8, r"$\int_a^b f(x)dx$",
      horizontalalignment='center', fontsize=20)
plt.figtext(0.9, 0.075, '$x$')
plt.figtext(0.075, 0.9, '$f(x)$')
ax.set_xticks((a, b))
ax.set_xticklabels(('$a$', '$b$'))
ax.set_yticks([f(a), f(b)]);

output_38_0

数值积分

对于上面的函数，由于非常简单，我们可以直接求出积分的函数，然后计算精确值，但是也可以利用数值分析的方法求解数值解，数值积分的方法也有多种：辛普森积分、梯形公式、龙贝格积分等，这些方法在Python的scipy.integrate中都有实现：

1
2
3
4
5
6
7
8
9
import scipy.integrate as sci
xi = np.linspace(0.5,9.5,100)
print(f'Simpon\'s rule:            {sci.simps(f(xi), xi)}')
print(f'Trapezodial rule:         {sci.trapz(f(xi),xi)}')
print(f'Fixed Guassian quadrature:{sci.fixed_quad(f,a,b)[0]}')
print(f'Adpative quadrature:      {sci.quad(f,a,b)[0]}')
print(f'Romberg integration:      {sci.romberg(f,a,b)}')
print('-------------------------------------')
print(f'The accurate result       {-np.cos(9.5)+np.cos(0.5)+(9.5**2-0.5**2)/4}')
Simpon's rule:            24.37474011548407
Trapezodial rule:         24.373463386819818
Fixed Guassian quadrature:24.366995967084602
Adpative quadrature:      24.374754718086752
Romberg integration:      24.374754718086713
-------------------------------------
The accurate result       24.374754718086752

符号计算（Sympy）

许多科学计算的软件都具有符号计算的功能，如Mathematica、Matlab等，Python中可以使用SymPy进行符号计算。

基本用法

Sympy为数学表达式提供三种渲染方法：

  • 基于LaTeX
  • 基于Unicode
  • 基于ASCII

如果是在Jupyter Notebook或者Markdown中使用，推荐采用LaTex渲染。

1
2
3
4
5
6
7
import sympy as sy
sy.init_session()
# Define symbols to work with
x = sy.Symbol('x')
y = sy.Symbol('y')
f = sy.Integral(x ** 2 + sy.sqrt(x))
f
IPython console for SymPy 1.3 (Python 3.7.1-64-bit) (ground types: python)

These commands were executed:
>>> from __future__ import division
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)
>>> init_printing()

Documentation can be found at http://docs.sympy.org/1.3/

∫(x+x2) dx\int \left(\sqrt{x} + x^{2}\right)\, dx

解方程

Sympy中的solve()方法可以用来解方程，并且支持复数：

1
sy.solve(x**2+y**0.5-3, x)

[−−y0.5+3.0,−y0.5+3.0]\left [ - \sqrt{- y^{0.5} + 3.0}, \quad \sqrt{- y^{0.5} + 3.0}\right ] 

1
sy.solve(sy.sin(x)-1)

[π2]\left [ \frac{\pi}{2}\right ]

微积分

前面我们用数值积分的方法求解了积分，我们也可以通过符号计算直接求解：

1
2
int_func = sy.integrate(sy.sin(x)+0.5*x, x)
int_func

0.25x2−cos⁡(x)0.25 x^{2} - \cos{\left (x \right )}

上面是不定积分的结果，如果我们计算[a,b]区间的定积分：

1
2
3
4
a = sy.Symbol('a')
b = sy.Symbol('b')
int_func_limits = sy.integrate(sy.sin(x) + 0.5 * x, (x, a, b))
int_func_limits

−0.25a2+0.25b2+cos⁡(a)−cos⁡(b)- 0.25 a^{2} + 0.25 b^{2} + \cos{\left (a \right )} - \cos{\left (b \right )}

那么，令a=0.5,b=9.5a=0.5,b=9.5可得结果：

1
int_func_limits.subs({a:0.5, b:9.5}).evalf()

24.374754718086824.3747547180868

类似地，我们可以进行微分的符号计算，并且可以求解导数为0的点：

1
2
3
4
5
6
f = (sy.sin(x) + 0.05 * x ** 2 +
    sy.sin(y) + 0.05 * y ** 2)
# Derive partial derivative
del_x = sy.diff(f,x)
del_y = sy.diff(f,y)
del_x, del_y

(0.1x+cos⁡(x),0.1y+cos⁡(y))\left ( 0.1 x + \cos{\left (x \right )}, \quad 0.1 y + \cos{\left (y \right )}\right ) 

1
2
3
4
5
6
7
8
9
xo = sy.nsolve(del_x, 1.5)
yo = sy.nsolve(del_y, 1.5)
print(xo,yo)
print(f.subs({x:xo, y:yo}).evalf())
print('----------------------')
xo = sy.nsolve(del_x, -1.5)
yo = sy.nsolve(del_y, -1.5)
print(xo,yo)
print(f.subs({x:xo, y:yo}).evalf())
1.74632928225285 1.74632928225285
2.27423381055640
----------------------
-1.42755177876459 -1.42755177876459
-1.77572565314742

小结

本文简单介绍金融中经常用到的一些数学工具在Python中的实现，例如回归会在多因子模型中经常用到，而数值积分则经常在金融衍生品的定价中用到。实际用到的会更加复杂，但是万变不离其宗。