Python金融分析(七):常用数学工具

引言

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

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

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

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

近似

回归

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

maxa1,,aD1Ii=1I(yid=1Dadbd(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.png

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

自定义基底

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.png

插值

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

插值的方法也有很多种,线性插值、三次样条插值、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.png

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

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

优化

在Python中,scipy包的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.png

全局优化

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,出于风险分散和边际递减的角度,其效用函数为各地收益的平方根之和,综上可以得到优化问题为:

maxa,bE(u(w1))=pw1u+(1p)w1dw1u=15a+5bw1d=5a+12bsubject to10010a+10ba,b0(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.png

数值积分

对于上面的函数,由于非常简单,我们可以直接求出积分的函数,然后计算精确值,但是也可以利用数值分析的方法求解数值解,数值积分的方法也有多种:辛普森积分、梯形公式、龙贝格积分等,这些方法在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.25x2cos(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中的实现,例如回归会在多因子模型中经常用到,而数值积分则经常在金融衍生品的定价中用到。实际用到的会更加复杂,但是万变不离其宗。

张da统帅 wechat
扫码订阅我的微信公众号