共计 6894 个字符，预计需要花费 18 分钟才能阅读完成。
金融市场的工夫序列数据是出了名的芜杂，并且很难解决。这也是为什么人们都对金融数学畛域如此乏味的局部起因!
咱们能够用来更好地了解趋势 (或帮忙模式识别 / 预测算法) 的一种办法是工夫序列平滑。以下传统的办法:
挪动平均线——简略、容易、无效(但会给工夫序列数据一个“滞后”的观测)，Savitzky-Golay 过滤器——无效但更简单，它蕴含了有一些直观的超参数
还有一个不太传统的办法是解热方程，它有更直观的超参数，而且也十分快! 在本文中，咱们将思考一个略微简单一些的方程，但它具备保留边缘的成果。
这个方程叫做 Perona-Malik PDE (偏微分方程)，它的平滑成果能够在上面的动图中看到:
上图是该放弃边缘平滑办法在用于于特斯拉 (TSLA) 在 2022 年的收盘价的成果。题目中的“t=x”对应于咱们平滑级数的工夫(以非维度单位)。
如果你对下面的成果感兴趣，那么本文将解释以下内容:
- Perona-Malik PDE(偏微分方程)，以及为什么要应用它
- 如何求解偏微分方程。
- 和热方程的比拟
Perona-Malik PDE
上面是将要解决的方程公式:
Perona-Malik PDE。式中 u 是咱们要平滑的工夫序列，α 是管制边保的参数(α 越小对应的边保越多)。
看着有点简单，咱们持续解释。当咱们试图解决这个 PDE 的原始模式时，它会导致一些问题; 所以咱们要思考一种批改后的模式:
基本上，函数 g 的外部进行了一次高斯函数卷积(也就是说，它变得更平滑了)。这被称为正则化，咱们只有晓得它是可解的就能够了
这个一个可怕的等式比下面更简单了，然而这咱们没有多个空间维度，咱们在平滑的是一个工夫序列，所以它只有一个维度。所以须要解的方程是:
这样看着就简略多了，至多所有的字母咱们都意识了😁，然而如果想理解细节，能够通过扩大和简化失去下面的公式进行推导，尽管不倡议这么做，然而如果你喜爱的话那随便。
咱们刚提到解决的工夫序列是一维的，然而为什么偏微分方程是二维的?
这个偏微分方程是依据工夫来求解的。从实质上讲工夫上的每一步都使数据进一步平滑。所以 t 越大，工夫序列越平滑，这意味着空间变量 x 示意工夫序列中的“工夫”，前面的求解会具体解释。
为什么要用这个方程呢?
热方程的问题是它不能很好地保留边。保留这些边缘来捕获价格的大幅疾速稳定可能是可取的，但要去除任何小但高频的噪声。
这种办法比热方程更难，因为 Perona-Malik PDE 是非线性的(不像热方程是线性的)。一般来说，非线性方程不像线性方程那么容易求解。
如何求解这个偏微分方程
咱们将应用一种称为无限差分（finite differences）的办法。它是一种求偏微分（或常微分）方程和方程组定解问题的数值解的办法。你能够将其视为每次咱们在下图中遇到交加时找到解决方案：
随着工夫和空间决裂成离散距离的图示。这里空间中的离散区间是从 [0, 1] 开始的，工夫上的离散区间是从 t=0 到 t=sk，其中 s 是咱们获取的区间。线的交点是咱们找到偏微分方程解的地位。
在解决数字之前，咱们须要用数学方法来定义整个问题。因为方程在空间上是二阶的，在工夫上是一阶的，所以须要两个边界条件和一个初始条件:
咱们将求解以平滑工夫序列的方程组（这个方程看起来比代码简单得多！），咱们的终点是股票价格工夫序列，并且起点总是具备雷同的价格。
那么咱们如何从数值上开始求解呢? 咱们最后的办法是用这些导数的无限差分近似，Perona-Malik PDE 中导数的近似值，这些导数的推导超出了本文的范畴，所以就不具体写了。
下面公式中，h 和 k 别离是空间和工夫离散点之间的间隔。这里能够应用雷同的公式计算 c 的导数。
在咱们的问题中，空间中的离散点被离开一天，所以 h = 1。倡议应用 k < 0.1 的值来放弃无限差分公式的一些准确性。因为如果 k 值太大，则无限差分计划会变得不稳固。
咱们当初须要将这些近似值放入偏微分方程中……这会让公式看起来更加简单😢，这是我的计算代数软件给出的后果:
这就是 Perona-Malik PDE 的离散模式，越来越简单了。有没有更好的办法呢？
咱们能够偷懒并应用微分矩阵。因为工夫序列是一组离散点，所以能够应用矩阵向量乘积进行微分。
如果以前从未据说过这个，那么这里有一些代码展现了如何应用矩阵向量积计算多项式的简略导数：
import numpy as np
import plotly.io as pio
import plotly.graph_objects as go
pio.renderers.default='browser'
if __name__ == '__main__':
n = 100 # The number of discrete points
a = -4 # The interval starting point
b = 4 # The interval ending point
# The discrete points and spacing
x = np.linspace(a, b, n)
h = (b-a)/n
# Function to differentiate, and its analytical derivative for checking
f = 10*x**3 + x**2
fx = 30*x**2 + 2*x
# Form the differentiation matrix
Dx = (np.diag(np.ones(n-1), 1)
- np.diag(np.ones(n-1), -1)
)/(2*h)
# Calculate the numerical derivative on the interior only
Ux = Dx[1:-1, :]@f
# Plot the figure
fig = go.Figure()
fig.add_trace(go.Scatter(x = x, y = f, name = 'Function'))
fig.add_trace(go.Scatter(x = x, y = fx, name = 'Analyitcal Derivative'))
fig.add_trace(
go.Scatter(x = x[1:-1],
y = Ux,
name='Finite Difference Derivative',
line={'dash': 'dash'}
)
)
fig.update_layout(yaxis = {'title': 'f(x), f_x(x)'},
xaxis = {'title': 'x'},
title = 'Finite Difference Derivative Test',
legend = {'x': 0, 'y': 1.075, 'orientation': 'h'},
width = 700,
height = 700,
)
fig.show()以上代码应该产生以下输入:
应用矩阵向量积能够对简略多项式求导。它实质上是一阶导数的无限差分迫近
已转化为矩阵向量乘积，应用上面的代码
Dx = (np.diag(np.ones(n-1), 1) # u_{r+1, s} terms
- np.diag(np.ones(n-1), -1) # u_{r-1, s} terms
)/(2*h)# Calculate the numerical derivative on the interior only
Ux = Dx[1:-1, :]@f通过应用这种无限差分公式，咱们只能找到外部点的导数。比方在域的第一个点 (x = r = 0) 有近似值：
尽管这是没有意义的，因为须要的计算点在域之外。然而这对咱们来说不是一个问题——因为咱们只解外部点的偏微分方程，而这些解在端点处是固定的。
咱们应用一个简略的小零碎的离散方程(比方有 5 个离散点)，下面的解释可能会清晰得多。
还有最初一个问题卷积是如何执行的? 最简略的办法是应用 scipy. nimage.gaussian_filter，然而这里我抉择的办法是解热方程，通过一点数学技巧，能够证实高斯卷积能够解决热方程（https://dsp.stackexchange.com…）。换句话说，咱们要解
这能够用离散模式示意为
高斯滤波中的标准差 (σ) 与咱们通过 σ²(τ) = 2τ 求解上述方程的“工夫”量无关，所以，要解的工夫越长，标准差越大，工夫序列就越平滑。
为什么要用这种形式进行卷积? 偏微分方程到卷积的连贯十分简洁! 并且因为能够将偏微分方程求解逻辑硬编码为循环，所以将其包装在 @numba.jit 装璜器中，进步了计算效率。
Python 实现
当初咱们曾经在数学方面做了艰辛的工作，编码就变得十分间接了!
卷积的实现如下所示:
import numpy as np
def convolve_PDE(U: np.array,
sigma: float = 1,
k: float = 0.05) -> np.array:
'''
Perform Gaussian convolution by solving the heat equation with Neumann
boundary conditions
Parameters
----------
U : np.array
The array to perform convolution on.
sigma : float, optional
The standard deviation of the guassian convolution.
k : float, optional
The step-size for the finite difference scheme (keep < 0.1 for accuracy)
Returns
-------
U : np.array
The convolved function
'''
t = 0
t_end = sigma**2/2
while t < t_end:
# Implementing the nuemann boundary conditions
U[0] = 2*k*U[1] + (1-2*k)*U[0]
U[-1] = 2*k*U[-2] + (1-2*k)*U[-1]
# Scheme on the interior nodes
U[1:-1] = k*(U[2:] + U[:-2]) + (1-2*k)*U[1:-1]
t += k
return U代码看起来就很简略了，对吧?Perona-Malik 求解器的实现也不简单：
def get_diff_mat(n: int) -> Tuple[np.array, np.array]:
'''
Get the first and second differentiation matrices, which are truncated
to find the derivative on the interior points.
Parameters
----------
n : int
The total number of discrete points
Returns
-------
[Dx, Dxx] : np.array
The first and second differentiation matrices, respecitvely.
'''
Dx = (np.diag(np.ones(n-1), 1)
- np.diag(np.ones(n-1), -1)
)/2
Dxx = (np.diag(np.ones(n-1), 1)
- 2*np.diag(np.ones(n), 0)
+ np.diag(np.ones(n-1), -1)
)
# Truncate the matrices so that we only determine the derivative on the
# interior points (i.e. we don't calculate the derivative on the boundary)
return Dx[1:-1, :], Dxx[1:-1, :]
def perona_malik_smooth(p: np.array,
alpha: float = 10.0,
k: float = 0.05,
t_end: float = 5.0) -> np.array:
'''
Solve the Gaussian convolved Perona-Malik PDE using a basic finite
difference scheme.
Parameters
----------
p : np.array
The price array to smoothen.
alpha : float, optional
A parameter to control how much the PDE resembles the heat equation,
the perona malik PDE -> heat equation as alpha -> infinity
k : float, optional
The step size in time (keep < 0.1 for accuracy)
t_end : float, optional
When to termininate the algorithm the larger the t_end, the smoother
the series
Returns
-------
U : np.array
The Perona-Malik smoothened time series
'''
Dx, Dxx = get_diff_mat(p.shape[0])
U = deepcopy(p)
t = 0
while t < t_end:
# Find the convolution of U with the guassian, this ensures that the
# PDE problem is well posed
C = convolve_PDE(deepcopy(U))
# Determine the derivatives by using matrix multiplication
Cx = [email protected]
Cxx = [email protected]
Ux = [email protected]
Uxx = [email protected]
# Find the spatial component of the PDE
PDE_space = (alpha*Uxx/(alpha + Cx**2)
- 2*alpha*Ux*Cx*Cxx/(alpha + Cx**2)**2
)
# Solve the PDE for the next time-step
U = np.hstack((np.array([p[0]]),
U[1:-1] + k*PDE_space,
np.array([p[-1]]),
))
t += k
return U因为应用了微分矩阵，所以基本上能够依照咱们查看方程的格局中写出偏微分方程的离散模式。这使得调试比硬编码无限差分方程更简略。
然而这会不会引入数据透露？
如果平滑一个大的工夫序列，而后将该序列宰割成更小的局部，那么相对会有数据透露。所以最好的办法是先切碎工夫序列，而后平滑每个较小的序列。这样基本不会有数据泄露！
与热方程的比拟
到目前为止，还没有说 Perona-Malik PDE 参数 α。如果你取 α 十分大，趋于无穷，就能够将 Perona-Malik PDE 化简为热方程。
对于大的 α，基本上有一个扩散主导的机制，其中边缘保留是无限的。咱们最终会失去这个方程组：
这里一维的热方程，以及问题的适当的初始 / 边界条件。用咱们的微分矩阵法能够很好很容易地写他的代码:
def heat_smooth(p: np.array,
k: float = 0.05,
t_end: float = 5.0) -> np.array:
'''
Solve the heat equation using a basic finite difference approach.
Parameters
----------
p : np.array
The price array to smoothen.
k : float, optional
The step size in time (keep < 0.1 for accuracy)
t_end : float, optional
When to termininate the algorithm the larger the t_end, the smoother
the series
Returns
-------
U : np.array
The heat equation smoothened time series
'''
# Obtain the differentiation matrices
_, Dxx = get_diff_mat(p.shape[0])
U = deepcopy(p)
t = 0
while t < t_end:
U = np.hstack((np.array([p[0]]),
U[1:-1] + k*[email protected],
np.array([p[-1]]),
))
t += k
return U运行 Perona-Malik PDE 其中设 α = 1,000,000, 和热方程, 根本失去了雷同的后果。
最初让咱们看看一些比照后果的动图
Perona-Malik PDE 与热方程的比拟。参数 σ = 1，α = 20, k=0.05。
上图是比拟 Perona-Malik、热方程和指数挪动均匀办法对 MSFT 股价在 2022 年期间的工夫序列数据进行平滑解决。
总结
总的来说，Perona-Malik 办法更好一些。尽管他的数学求解要简单的多，但它的确对数据产生了十分好的后果。就集体而言，倡议在开发过程中同时思考 Perona Malik 和热方程办法，看看哪种办法能够为咱们解决的问题提供更好的后果。
本文的残缺代码：https://avoid.overfit.cn/post…
作者：Danny Groves