量化交易软件：未知概率密度函数的核密度估计

简介

随着 MQL5 的性能改善及计算机效率的稳步提升，赫兹量化平台用户能够将相当复杂和高深的数学方法用于市场分析中。这些方法可能来自经济学、计量经济学和统计学在内的多个领域，但是在任何情形下，赫兹量化使用这些方法都无法绕开概率密度函数的理念。

很多常用的分析方法基于所用模型中数据分布正态性或错误正态性的假设而开发。此外，赫兹量化通常有必要知道所用模型中各部分的分布，以评估分析结果。在两种情形下我们都有一个任务，即创建一种“工具”（在理想情况下可通用）以估计未知概率密度函数。

本文试图创建一个类，以或多或少实现用于估计未知概率密度函数的通用算法。最初的构想是在估计过程中不使用外部方法，即仅通过 MQL5 实现所有一切。但最终，最初的构想在一定程度上发生了变化。很明显，概率密度函数的目测估计任务由两个独立的部分组成。

即，估计本身的计算及其可视化（以图表或示意图显示）。计算自然已经通过 MQL5 实现，而可视化必须通过创建 HTML 页面实施和在网页浏览器中显示。该解决方案用于最终获得矢量形式的图形。

但由于计算部分和结果显示的分别实现，读者自然可使用其他任何可用的可视化方法。此外，赫兹量化希望各种库，包括图形库，将出现在赫兹量化 中（据我们所知，此项工作正在进行中）。变更建议解决方案的显示部分并不困难，因为赫兹量化提供了构建图形和示意图的进阶方法。

应初步指出的是，创建真正的通用算法用于估计序列概率密度函数经证明是一个无法实现的目标。尽管建议解决方案的专门化程度不高，但也不能称为是完全通用的。问题是，密度估计过程中的优化准则大相庭径，例如对于钟形分布，有正态分布和指数分布。

因此，如果有关于估计分布的一些初步信息，赫兹量化总是能够为每种特定情形选择最合适的解决方案。但是，尽管如此，我们将假设我们对估计密度的本质一无所知。这种方法无疑会影响估计的质量，但我们希望它能够对差异较大的密度提供估计的可能。

由于在市场数据分析的过程中我们经常要处理非平稳序列，赫兹量化最感兴趣的是较短和中等序列密度的估计。这是决定所使用估计方法的选取的关键时刻。

直方图和 P 样条曲线可成功用于包含超过一百万个值的长序列。但是，当我们尝试为包含 10 至 20 个值的序列有效地构建直方图时，出现了一些问题。因此，我们将重点放在包含约 10 到 10 000 个值的序列上。

1. 概率密度函数估计方法

目前已知有大量的受欢迎程度不一的概率密度函数估计方法。这些方法在互联网上可以很容易地搜索到，例如，通过使用诸如“概率密度估计”、“概率密度”、“密度估计”等关键词。但遗憾的是，我们没有在其中找到最佳的方法。这些方法都各有利弊。

直方图在传统上应用于密度估计 [1]。使用直方图（包括平滑直方图）可以提供高质量的概率密度估计，但这仅仅是针对长序列的处理情形。正如之前提到的，赫兹量化无法将一个短序列划分为数量庞大的组，而由 2 到 3 个柱组成的直方图无法说明这个序列的概率密度分布规律。因此，我们不得不放弃使用直方图。

另一个非常有名的估计方法是核密度估计 [2]。使用核平滑的理念在 [3] 中得到了很好的展现。因此，赫兹量化选择了该方法，尽管它有这样那样的弊端。与此方法的实现相关的一些方面将在下文中简单讨论。

还有必要提及一个非常有趣的密度估计方法，该方法使用“期望最大化”算法 [4]。这个算法允许将序列划分为具有比如正态分布规律的单独部分。在确定单独部分的参数后，有可能通过对得到的分布曲线进行求和而获得密度估计。此方法在 [5] 中进行了说明。本文未涉及此方法以及其他许多方法的实现和检验。在实践中，赫兹量化无法对各种材料包含的大量密度估计方法一一进行检验。

我们继续讨论选定实现的核密度估计方法。

2. 概率密度函数的核密度估计

概率密度函数的核密度估计基于核平滑方法。该方法的原理可在 [6]、[7] 中找到。

核平滑的基本理念相当简单。赫兹量化用户已熟悉移动平均线 (MA) 指标。该指标可以很容易地描绘成沿序列滑动的窗口，其加权值在窗口内被平滑。窗口可以是矩形、指数或具有其他形状。在核平滑过程中，赫兹量化可以很容易地看到同样的滑动窗口（例如，[3]）。但在这种情况下它是对称的。

核平滑中最常用窗口的示例可在 [8] 中找到。如果核平滑中使用了零阶回归，到达窗口（核）的序列加权值便被平滑，就像 MA 中一样。在我们处理滤波问题时，我们可以看到相同种类的窗口函数应用程序。但现在，同一过程的表示有了一点不同。当核（窗口）调用滤波器脉冲特性时，幅频特性和相频特性被使用。

这些例子表明，一个事物往往可以有不同的表现方式。自然 ，这对数学工具有所裨益。但在讨论此类问题时也可能会导致混淆。

尽管核密度估计与前面提到的核平滑采用相同的原则，其算法稍有不同。

接下来我们看看定义密度估计的表达式。

​

添加图片注释，不超过 140 字（可选）

其中

• x - 长度为 n 的序列；

• K - 对称核；

• h - 平滑参数范围。

在接下来的密度估计中，我们仅使用高斯核：

​

添加图片注释，不超过 140 字（可选）

根据上述表达式，点 X 的密度通过对点 X 和序列之间的差值定义的数量的核值求和计算得出。此外，用于密度计算的 X 点可能与序列本身的值不一致。

下面是实现核密度估计算法的基本步骤。

1. 估计平均值和输入序列标准偏差。

2. 正态化输入序列。从它的每个值中减去先前获得的平均值，再除以标准偏差值。在经过这样的正态化后，原始序列将具有 0 平均值和等于 1 的标准偏差。这种正态化不是计算密度所必需的，但它可以使结果图表一致，因为就 X 刻度上的任何序列而言，都是用标准偏差单位来表示的值。

3. 找出正态化序列的最大值和最小值。

4. 创建 2 个数组，数组的大小要和在所得图表上显示的预计点数相对应。例如，如果图表使用 200 个点构建，则每个数组必须相应地包含 200 个值。

5. 保留一个创建的数组用于存储结果。另一个数组用于形成为其执行密度估计的点的值。为此，赫兹量化必须在之前准备的最大值和最小值之间形成 200（在本例中）个等距的值，并将它们保存在准备好的数组中。

6. 使用上文提到的表达式，我们应在 200 个（在本例中）检验点中执行密度估计，将结果保存在步骤 4 准备的数组中。

下面是这个算法的软件实现。

//+------------------------------------------------------------------+ //| CDens.mqh | //| 2012, victorg | //| https://www.mql5.com | //+------------------------------------------------------------------+ #property copyright "2012, victorg" #property link "https://www.mql5.com" #include <Object.mqh> //+------------------------------------------------------------------+ //| Class Kernel Density Estimation | //+------------------------------------------------------------------+ class CDens:public CObject { public: double X[]; // Data int N; // Input data length (N >= 8) double T[]; // Test points for pdf estimating double Y[]; // Estimated density (pdf) int Np; // Number of test points (Npoint>=10, default 200) double Mean; // Mean (average) double Var; // Variance double StDev; // Standard deviation double H; // Bandwidth public: void CDens(void); int Density(double &x[],double hh); void NTpoints(int n); private: void kdens(double h); }; //+------------------------------------------------------------------+ //| Constructor | //+------------------------------------------------------------------+ void CDens::CDens(void) { NTpoints(200); // Default number of test points } //+------------------------------------------------------------------+ //| Setting number of test points | //+------------------------------------------------------------------+ void CDens::NTpoints(int n) { if(n<10)n=10; Np=n; // Number of test points ArrayResize(T,Np); // Array for test points ArrayResize(Y,Np); // Array for result (pdf) } //+------------------------------------------------------------------+ //| Density | //+------------------------------------------------------------------+ int CDens::Density(double &x[],double hh) { int i; double a,b,min,max,h; N=ArraySize(x); // Input data length if(N<8) // If N is too small { Print(__FUNCTION__+": Error! Not enough data length!"); return(-1); } ArrayResize(X,N); // Array for input data ArrayCopy(X,x); // Copy input data ArraySort(X); Mean=0; for(i=0;i<N;i++)Mean=Mean+(X[i]-Mean)/(i+1.0); // Mean (average) Var=0; for(i=0;i<N;i++) { a=X[i]-Mean; X[i]=a; Var+=a*a; } Var/=N; // Variance if(Var<1.e-250) // Variance is too small { Print(__FUNCTION__+": Error! The variance is too small or zero!"); return(-1); } StDev=MathSqrt(Var); // Standard deviation for(i=0;i<N;i++)X[i]=X[i]/StDev; // Data normalization (mean=0,stdev=1) min=X[ArrayMinimum(X)]; max=X[ArrayMaximum(X)]; b=(max-min)/(Np-1.0); for(i=0;i<Np;i++)T[i]=min+b*(double)i; // Create test points //-------------------------------- Bandwidth selection h=hh; if(h<0.001)h=0.001; H=h; //-------------------------------- Density estimation kdens(h); return(0); } //+------------------------------------------------------------------+ //| Gaussian kernel density estimation | //+------------------------------------------------------------------+ void CDens::kdens(double h) { int i,j; double a,b,c; c=MathSqrt(M_PI+M_PI)*N*h; for(i=0;i<Np;i++) { a=0; for(j=0;j<N;j++) { b=(T[i]-X[j])/h; a+=MathExp(-b*b*0.5); } Y[i]=a/c; // pdf } } //--------------------------------------------------------------------

NTpoints() 方法允许设置所需的等距检验点数，将为这些点执行密度估计。此方法必须在调用 Density() 方法之前调用。当调用 Density() 方法时，指向包含输入数据和范围值（平滑参数）的数组的链接作为该方法的参数传递至该方法。