在 MQL5 中实现广义赫斯特指数和方差比检验

MetaTrader 5 — 交易系统 | 12 九月 2024, 09:18

概述

在计算赫斯特指数一文中，我们了解了分形分析的概念以及如何将其应用于金融市场。在那篇文章中，作者介绍了估计赫斯特指数的重标定范围法 (R/S)。在本文中，我们采用一种不同的方法，展示广义赫斯特指数（GHE）的应用，以对序列的性质进行分类。我们将重点利用 GHE 来识别表现出均值回复趋势的外汇交易品种，希望能利用这种行为。

首先，我们简要讨论一下 GHE 的基础知识，以及它与最初的赫斯特指数有何不同。为此，我们将介绍一种可用于确认 GHE 分析结果的统计检验方法，即方差比检验（VRT）。接下来，我们将讨论如何应用 GHE 为均值回归交易识别候选外汇交易品种。在此，我们介绍一种生成进场和出场信号的指标。最后，我们将在一个基础的 EA 交易中对其进行测试。

了解广义赫斯特指数

赫斯特指数衡量时间序列的缩放特性。缩放属性是描述系统在大小或时间尺度变化时的行为的基本特征。就时间序列数据而言，缩放特性有助于深入了解不同时间尺度之间的关系以及数据中存在的模式。对于平稳序列，与几何随机游走相比，后续值随时间的变化更为缓慢。为了从数学上量化这种行为，我们分析了序列中的扩散率。方差用作表示其他值偏离序列中第一个值的比率的度量。

在上式中，"K" 代表进行分析的任意滞后期。为了更好地了解序列的性质，我们还必须评估其他滞后的方差。因此，"K" 可以指定任何小于数列长度的正整数值。最大的滞后是可自由支配的，牢记这一点很重要。因此，赫斯特指数与不同滞后下方差的缩放行为有关。根据幂律，其定义如下：

GHE 是对原始公式的概括，其中 2 被替换为一个变量，通常用 "q" 表示。从而将上述公式改为：

和

GHE 通过分析时间序列中连续点之间变化的不同统计特征如何随不同阶矩的变化而变化，从而扩展了最初的赫斯特方法。用数学术语来说，矩是描述分布形状和特征的统计量。qth 阶矩是一种特殊的矩，其中 "q" 是一个决定阶数的参数。对于每个 "q" 值，GHE 都强调时间序列的不同特征。具体来说，当 q=1 时，结果描述了绝对偏差的缩放特性。而 q=2 在研究长程依赖性时最为重要。

在 MQL5 中实现 GHE

本节将介绍 MQL5 中 GHE 的实现。之后，我们将通过分析人工生成的时间序列随机样本对其进行测试。我们的实现方案包含在 GHE.mqh 文件中。文件首先包含 VectorMatrixTools.mqh，其中包含用于初始化常见类型向量和矩阵的各种函数的定义。该文件的内容如下所示。

//+------------------------------------------------------------------+
//|                                            VectorMatrixTools.mqh |
//|                                  Copyright 2023, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2023, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
//+------------------------------------------------------------------+
//|Vector arange initialization                                      |
//+------------------------------------------------------------------+
template<typename T>
void arange(vector<T> &vec,T value=0.0,T step=1.0)
  {
   for(ulong i=0; i<vec.Size(); i++,value+=step)
      vec[i]=value;
  }
//+------------------------------------------------------------------+
//| Vector sliced initialization                                     |
//+------------------------------------------------------------------+
template<typename T>
void slice(vector<T> &vec,const vector<T> &toCopyfrom,ulong start=0,ulong stop=ULONG_MAX, ulong step=1)
  {
   start = (start>=toCopyfrom.Size())?toCopyfrom.Size()-1:start;
   stop  = (stop>=toCopyfrom.Size())?toCopyfrom.Size()-1:stop;
   step  = (step==0)?1:step;

   ulong numerator = (stop>=start)?stop-start:start-stop;
   ulong size = (numerator/step)+1;
   if(!vec.Resize(size))
     {
      Print(__FUNCTION__ " invalid slicing parameters for vector initialization");
      return;
     }

   if(stop>start)
     {
      for(ulong i =start, k = 0; i<toCopyfrom.Size() && k<vec.Size() && i<=stop; i+=step, k++)
         vec[k] = toCopyfrom[i];
     }
   else
     {
      for(long i = long(start), k = 0; i>-1 && k<long(vec.Size()) && i>=long(stop); i-=long(step), k++)
         vec[k] = toCopyfrom[i];
     }
  }
//+------------------------------------------------------------------+
//| Vector sliced initialization  using array                        |
//+------------------------------------------------------------------+
template<typename T>
void assign(vector<T> &vec,const T &toCopyfrom[],ulong start=0,ulong stop=ULONG_MAX, ulong step=1)
  {
   start = (start>=toCopyfrom.Size())?toCopyfrom.Size()-1:start;
   stop  = (stop>=toCopyfrom.Size())?toCopyfrom.Size()-1:stop;
   step  = (step==0)?1:step;

   ulong numerator = (stop>=start)?stop-start:start-stop;
   ulong size = (numerator/step)+1;

   if(size != vec.Size() &&  !vec.Resize(size))
     {
      Print(__FUNCTION__ " invalid slicing parameters for vector initialization");
      return;
     }

   if(stop>start)
     {
      for(ulong i =start, k = 0; i<ulong(toCopyfrom.Size()) && k<vec.Size() && i<=stop; i+=step, k++)
         vec[k] = toCopyfrom[i];
     }
   else
     {
      for(long i = long(start), k = 0; i>-1 && k<long(vec.Size()) && i>=long(stop); i-=long(step), k++)
         vec[k] = toCopyfrom[i];
     }
  }
//+------------------------------------------------------------------+
//| Matrix initialization                                            |
//+------------------------------------------------------------------+
template<typename T>
void rangetrend(matrix<T> &mat,T value=0.0,T step=1.0)
  {
   ulong r = mat.Rows();

   vector col1(r,arange,value,step);

   vector col2 = vector::Ones(r);

   if(!mat.Resize(r,2) || !mat.Col(col1,0) || !mat.Col(col2,1))
     {
      Print(__FUNCTION__ " matrix initialization error: ", GetLastError());
      return;
     }

  }
//+-------------------------------------------------------------------------------------+
//| ols design Matrix initialization with constant and first column from specified array|
//+-------------------------------------------------------------------------------------+
template<typename T>
void olsdmatrix(matrix<T> &mat,const T &toCopyfrom[],ulong start=0,ulong stop=ULONG_MAX, ulong step=1)
  {
   vector col0(1,assign,toCopyfrom,start,stop,step);

   ulong r = col0.Size();

   if(!r)
     {
      Print(__FUNCTION__," failed to initialize first column ");
      return;
     }

   vector col1 = vector::Ones(r);

   if(!mat.Resize(r,2) || !mat.Col(col0,0) || !mat.Col(col1,1))
     {
      Print(__FUNCTION__ " matrix initialization error: ", GetLastError());
      return;
     }

  }
//+------------------------------------------------------------------+
//|vector to array                                                   |
//+------------------------------------------------------------------+
bool vecToArray(const vector &in, double &out[])
  {
//---
   if(in.Size()<1)
     {
      Print(__FUNCTION__," Empty vector");
      return false;
     }
//---
   if(ulong(out.Size())!=in.Size() && ArrayResize(out,int(in.Size()))!=int(in.Size()))
     {
      Print(__FUNCTION__," resize error ", GetLastError());
      return false;
     }
//---
   for(uint i = 0; i<out.Size(); i++)
      out[i]=in[i];
//---
   return true;
//---
  }
//+------------------------------------------------------------------+
//| difference a vector                                               |
//+------------------------------------------------------------------+
vector difference(const vector &in)
  {
//---
   if(in.Size()<1)
     {
      Print(__FUNCTION__," Empty vector");
      return vector::Zeros(1);
     }
//---
   vector yy,zz;
//---
   yy.Init(in.Size()-1,slice,in,1,in.Size()-1,1);
//---
   zz.Init(in.Size()-1,slice,in,0,in.Size()-2,1);
//---
   return yy-zz;
  }
//+------------------------------------------------------------------+

GHE.mqh，包含 "gen_hurst()" 函数及其重载的定义。一个方法是以向量形式提供要分析的数据，另一个方法是以数组形式提供数据。该函数还接受一个整数 "q"，以及可选的整数参数 "lower "和 "upper"（它们有默认值）。这就是上一节描述 GHE 时提到的 "q"。最后两个参数是可选的，"lower" 和 "upper" 共同定义了进行分析的滞后范围，类似于上述公式中 "K" 值的范围。

//+--------------------------------------------------------------------------+
//|overloaded gen_hurst() function that works with series contained in vector|
//+--------------------------------------------------------------------------+
double general_hurst(vector &data, int q, int lower=0,int upper=0)
  {
   double series[];

   if(!vecToArray(data,series))
      return EMPTY_VALUE;
   else
      return general_hurst(series,q,lower,upper);
  }

当出现错误时，该函数将返回与内置常量 EMPTY_VALUE 等价的值，并将有用的字符串信息输出到终端的"专家"选项卡。在 "gen_hurst()"中，首先常规检查传给它的参数。确保它们符合以下条件：

"q" 不能小于 1。
"lower" 不能小于 2，也不能大于或等于 "upper"。
而 "upper" 参数不能超过所分析数据序列长度的一半。如果不满足其中任何一个条件，函数将立即标示错误。

if(data.Size()<100)
     {
      Print("data array is of insufficient length");
      return EMPTY_VALUE;
     }

   if(lower>=upper || lower<2 ||  upper>int(floor(0.5*data.Size())))
     {
      Print("Invalid input for lower and/or upper");
      return EMPTY_VALUE;
     }

   if(q<=0)
     {
      Print("Invalid input for q");
      return EMPTY_VALUE;
     }

   uint len = data.Size();

   int k =0;

   matrix H,mcord,lmcord;
   vector n_vector,dv,vv,Y,ddVd,VVVd,XP,XY,PddVd,PVVVd,Px_vector,Sqx,pt;
   double dv_array[],vv_array[],mx,SSxx,my,SSxy,cc1,cc2,N;

   if(!H.Resize(ulong(upper-lower),1))
     {
      Print(__LINE__," ",__FUNCTION__," ",GetLastError());
      return EMPTY_VALUE;
     }

   for(int i=lower; i<upper; i++)
     {
      vector x_vector(ulong(i),arange,1.0,1.0);

      if(!mcord.Resize(ulong(i),1))
        {
         Print(__LINE__," ",__FUNCTION__," ",GetLastError());
         return EMPTY_VALUE;
        }

      mcord.Fill(0.0);

该函数的内部工作从一个从 "lower" 到 "upper" 的 "for" 循环开始，对于每一个 "i"，它使用 "arange" 函数创建一个包含 "i" 个元素的向量 "x_vector"。然后，它会调整矩阵 "mcord" 的大小，使其具有 "i" 行和一列。

 for(int j=1; j<i+1; j++)
        {
         if(!diff_array(j,data,dv,Y))
            return EMPTY_VALUE;

内循环首先使用辅助函数 "diff_array()" 计算 "data" 数组中的差值，并将其存储到向量 "dv" 和 "Y" 中。

 N = double(Y.Size());

         vector X(ulong(N),arange,1.0,1.0);

         mx = X.Sum()/N;

         XP = MathPow(X,2.0);

         SSxx = XP.Sum() - N*pow(mx,2.0);

         my = Y.Sum()/N;

         XY = X*Y;

         SSxy = XY.Sum() - N*mx*my;

         cc1 = SSxy/SSxx;

         cc2 = my - cc1*mx;

         ddVd = dv - cc1;

         VVVd = Y - cc1*X - cc2;

         PddVd = MathAbs(ddVd);

         PddVd = pow(PddVd,q);

         PVVVd = MathAbs(VVVd);

         PVVVd = pow(PVVVd,q);

         mcord[j-1][0] = PddVd.Mean()/PVVVd.Mean();
        }

这里计算的是指定滞后的方差。其结果存储在矩阵 "mcord" 中。

 Px_vector = MathLog10(x_vector);

      mx = Px_vector.Mean();

      Sqx = MathPow(Px_vector,2.0);

      SSxx = Sqx.Sum() - i*pow(mx,2.0);

      lmcord = log10(mcord);

      my = lmcord.Mean();

      pt = Px_vector*lmcord.Col(0);

      SSxy = pt.Sum() - i*mx*my;

      H[k][0]= SSxy/SSxx;

      k++;

在内循环之外，外循环的最后一段将更新主 "H" 矩阵值。最后，函数会返回 "H" 矩阵除以 "q" 的平均值。

 return H.Mean()/double(q);

为了测试我们的 GHE 功能，我们准备了作为 EA 交易实现的应用程序 GHE.ex5。通过它，人们可以直观地看到具有预定特征的随机序列，并观察 GHE 是如何工作的。通过完全交互式操作，可以调整 GHE 的所有参数，并在一定范围内调整序列的长度。一个有趣的功能是，在应用 GHE 之前，可以对序列进行对数变换，以测试以这种方式预处理数据是否有任何好处。

GHE 互动应用程序

我们都知道，在实际应用中，数据集会受到过多噪音的困扰。由于 GHE 得出的估计值对样本量很敏感，我们需要检验结果的显著性。要做到这一点，可以进行一个名为 "方差比（VR）检验" 的假设检验。

方差比检验

方差比检验（Variance Ratio Test）是用于评估时间序列随机性的统计检验，方法是检验序列的方差是否随时间间隔的延长而成正比增加。该检验所依据的理念是，如果要检验的序列遵循随机漫步，那么在给定时间间隔内序列变化的方差应随时间间隔的长度线性增加。如果方差增加的速度较慢，则可能表明序列变化中存在序列相关性，说明序列是可预测的。方差比测试的是：

VRT

等于 1，其中：
- X() 是感兴趣的时间序列。
- K 是一个任意的滞后值。
- Var() 表示方差。

检验的零假设是时间序列遵循随机漫步，因此方差比应等于 1。如果方差比明显不同于 1，则可能拒绝零假设，这表明时间序列中存在某种形式的可预测性或序列相关性。

方差比检验的实现

VR 检验是通过 VRT.mqh 中定义的 CVarianceRatio 类实现的。有两个方法可用于进行 VR 测试，它们名为 "Vrt()"，一个用于矢量，另一种用于数组。方法参数说明如下：

"lags" 指定方差计算中使用的周期数或滞后期数，在我们希望使用 VR 检验来评估 GHE 估计值的显著性时，我们可以将 "lags" 设置为 "gen_hurst()" 中相应的 "lower" 或 "upper" 参数。该值不能小于 2。
"trend" 是一个枚举，允许指定我们要检验的随机游走类型。只有两个选项有作用，即 TREND_CONST_ONLY 和 TREND_NONE。
"debiased" 表示是否使用测试的去偏差版本，只有当 "overlap" 为真时才适用。设置为 "true" 时，该函数采用偏差校正技术来调整方差比估计值，从而更准确地反映方差之间的真实关系。这在处理小样本序列时非常有用。
"overlap" 表示是否使用所有重叠的区块。如果为假，则序列长度减一后必须是 "lags" 的整数倍。如果不满足这一条件，输入序列末尾的一些值将被丢弃。
"robust" 选择是否考虑异方差（true）或只考虑同方差（false）。在统计分析中，异方差过程的方差不恒定，而同方差序列的方差恒定。

成功执行后，"Vrt()" 方法会返回 true，然后可以调用任何获取（getter）方法来获取测试结果的所有方面。

//+------------------------------------------------------------------+
//| CVarianceRatio  class                                            |
//| Variance ratio hypthesis test for a random walk                  |
//+------------------------------------------------------------------+
class CVarianceRatio
  {
private:
   double            m_pvalue;     //pvalue
   double            m_statistic;  //test statistic
   double            m_variance;   //variance
   double            m_vr;         //variance ratio
   vector            m_critvalues; //critical values

public:
                     CVarianceRatio(void);
                    ~CVarianceRatio(void);

   bool              Vrt(const double &in_data[], ulong lags, ENUM_TREND trend = TREND_CONST_ONLY, bool debiased=true, bool robust=true, bool overlap = true);
   bool              Vrt(const vector &in_vect, ulong lags, ENUM_TREND trend = TREND_CONST_ONLY, bool debiased=true, bool robust=true, bool overlap = true);

   double            Pvalue(void) { return m_pvalue;}
   double            Statistic(void) { return m_statistic;}
   double            Variance(void) { return m_variance;}
   double            VRatio(void) { return m_vr;}
   vector            CritValues(void) { return m_critvalues;}
  };

在 "Vrt()" 中，如果 "overlap" 为 false，我们将检查输入序列的长度是否能被 "lags" 整除。如果没有，我们会修剪系列的末尾，并就数据长度发出警告。然后，我们根据更新后的序列长度重新分配 "nobs"，并计算趋势项 "mu"。在这里，我们计算序列中相邻元素的差值，并将其保存到 "delta_y" 中。使用 "delta_y" 计算方差，并保存在变量 "sigma2_1" 中。如果没有重叠，我们就计算非重叠区块的方差。否则，我们计算重叠区块的方差。如果 "debiased" 和 "overlap" 同时启用，我们就会调整差异。这里，"m_varianced" 的计算取决于 "overlap" 和 "robust"。最后，计算方差比、检验统计量和 p 值。

//+------------------------------------------------------------------+
//| main method for computing Variance ratio test                    |
//+------------------------------------------------------------------+
bool CVarianceRatio::Vrt(const vector &in_vect,ulong lags,ENUM_TREND trend=1,bool debiased=true,bool robust=true,bool overlap=true)
  {
   ulong nobs = in_vect.Size();

   vector y = vector::Zeros(2),delta_y;

   double mu;

   ulong nq = nobs - 1;

   if(in_vect.Size()<1)
     {
      Print(__FUNCTION__, "Invalid input, no data supplied");
      return false;
     }

   if(lags<2 || lags>=in_vect.Size())
     {
      Print(__FUNCTION__," Invalid input for lags");
      return false;
     }

   if(!overlap)
     {
      if(nq % lags != 0)
        {
         ulong extra = nq%lags;
         if(!y.Init(5,slice,in_vect,0,in_vect.Size()-extra-1))
           {
            Print(__FUNCTION__," ",__LINE__);
            return false;
           }
         Print("Warning:Invalid length for input data, size is not exact multiple of lags");
        }
     }
   else
      y.Copy(in_vect);

   nobs = y.Size();

   if(trend == TREND_NONE)
      mu = 0;
   else
      mu = (y[y.Size()-1] - y[0])/double(nobs - 1);

   delta_y = difference(y);

   nq = delta_y.Size();

   vector mudiff = delta_y - mu;

   vector mudiff_sq = MathPow(mudiff,2.0);

   double sigma2_1 = mudiff_sq.Sum()/double(nq);

   double sigma2_q;

   vector delta_y_q;

   if(!overlap)
     {
      vector y1,y2;
      if(!y1.Init(3,slice,y,lags,y.Size()-1,lags) ||
         !y2.Init(3,slice,y,0,y.Size()-lags-1,lags))
        {
         Print(__FUNCTION__," ",__LINE__);
         return false;
        }

      delta_y_q = y1-y2;

      vector delta_d = delta_y_q - double(lags) * mu;

      vector delta_d_sqr = MathPow(delta_d,2.0);

      sigma2_q = delta_d_sqr.Sum()/double(nq);
     }
   else
     {
      vector y1,y2;
      if(!y1.Init(3,slice,y,lags,y.Size()-1) ||
         !y2.Init(3,slice,y,0,y.Size()-lags-1))
        {
         Print(__FUNCTION__," ",__LINE__);
         return false;
        }

      delta_y_q = y1-y2;

      vector delta_d = delta_y_q - double(lags) * mu;

      vector delta_d_sqr = MathPow(delta_d,2.0);

      sigma2_q = delta_d_sqr.Sum()/double(nq*lags);
     }


   if(debiased && overlap)
     {
      sigma2_1 *= double(nq)/double(nq-1);
      double mm = (1.0-(double(lags)/double(nq)));
      double m = double(lags*(nq - lags+1));// * (1.0-double(lags/nq));
      sigma2_q *= double(nq*lags)/(m*mm);
     }

   if(!overlap)
      m_variance = 2.0 * (lags-1);
   else
      if(!robust)
         m_variance = double((2 * (2 * lags - 1) * (lags - 1)) / (3 * lags));
      else
        {
         vector z2, o, p;
         z2=MathPow((delta_y-mu),2.0);
         double scale = pow(z2.Sum(),2.0);
         double theta = 0;
         double delta;
         for(ulong k = 1; k<lags; k++)
           {
            if(!o.Init(3,slice,z2,k,z2.Size()-1) ||
               !p.Init(3,slice,z2,0,z2.Size()-k-1))
              {
               Print(__FUNCTION__," ",__LINE__);
               return false;
              }
            o*=double(nq);
            p/=scale;
            delta = o.Dot(p);
            theta+=4.0*pow((1.0-double(k)/double(lags)),2.0)*delta;
           }
         m_variance = theta;
        }
   m_vr = sigma2_q/sigma2_1;

   m_statistic = sqrt(nq) * (m_vr - 1)/sqrt(m_variance);

   double abs_stat = MathAbs(m_statistic);

   m_pvalue = 2 - 2*CNormalDistr::NormalCDF(abs_stat);

   return true;
  }

为了测试该类，我们修改了用于演示 "gen_hurst()" 函数的应用程序 GHE.ex5。因为 GHE 是由一系列滞后期定义的，而分析集中于这些滞后期。我们可以校准 VRT，在相同的滞后期范围内检验 GHE 结果的显著性。通过在最小和最大滞后期运行 VRT，我们应该可以获得足够的信息。在 GHE.ex5 中，先显示 "lower" 滞后点的方差比，再显示 "upper" 滞后点的方差比。
请记住，差异显著的方差比表明数据具有可预测性。方差比接近 1，表明该序列与随机漫步相差无几。通过使用应用程序，测试不同的参数组合，我们注意到 GHE 和 VRT 结果都受到样本量的影响。

趋势的错误分类

对于长度小于 1000 的序列，这两种方法有时都会得到意想不到的结果。

原始值

此外，在比较使用原始值和对数转换值进行的测试时， GHE 的结果也会出现显著差异。

对数转换

既然我们已经熟悉了 VRT 和 GHE，就可以将它们应用到我们的均值回归策略中。如果已知价格序列是均值回归的，我们就可以根据其当前与均值的偏差来粗略估计价格的走势。我们策略的基础是分析指定时期价格序列的特征。通过这一分析，我们建立了一个模型，可以估算出价格在偏离常态过远后可能会回升的点位。我们需要某种方法来衡量和量化这种分化，以生成入场和退出信号。

Z 分数

Z 分数（z-score）衡量的是价格与其平均值的标准差数量。通过对价格进行归一化处理，z-score 在零附近波动。让我们看看将 z-score 作为指标来实现后的曲线图。完整代码如下所示。

//+------------------------------------------------------------------+
//|                                                       Zscore.mq5 |
//|                                  Copyright 2023, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2023, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#property version   "1.00"
#include<VectorMatrixTools.mqh>
#property indicator_separate_window
#property indicator_buffers 1
#property indicator_plots   1
//--- plot Zscore
#property indicator_label1  "Zscore"
#property indicator_type1   DRAW_LINE
#property indicator_color1  clrBlue
#property indicator_style1  STYLE_SOLID
#property indicator_width1  1
//--- input parameters
input int      z_period = 10;
//--- indicator buffers
double         ZscoreBuffer[];
vector vct;
//+------------------------------------------------------------------+
//| Custom indicator initialization function                         |
//+------------------------------------------------------------------+
int OnInit()
  {
//--- indicator buffers mapping
   SetIndexBuffer(0,ZscoreBuffer,INDICATOR_DATA);
//----
   PlotIndexSetDouble(0,PLOT_EMPTY_VALUE,0);
//---
   PlotIndexSetInteger(0,PLOT_DRAW_BEGIN,z_period-1);
//---
   return(INIT_SUCCEEDED);
  }
//+------------------------------------------------------------------+
//| Custom indicator iteration function                              |
//+------------------------------------------------------------------+
int OnCalculate(const int rates_total,
                const int prev_calculated,
                const datetime &time[],
                const double &open[],
                const double &high[],
                const double &low[],
                const double &close[],
                const long &tick_volume[],
                const long &volume[],
                const int &spread[])
  {
//---
   if(rates_total<z_period)
     {
      Print("Insufficient history");
      return -1;
     }
//---
   int limit;
   if(prev_calculated<=0)
      limit = z_period - 1;
   else
      limit = prev_calculated - 1;
//---
   for(int i = limit; i<rates_total; i++)
     {
      vct.Init(ulong(z_period),assign,close,ulong(i-(z_period-1)),i,1);
      if(vct.Size()==ulong(z_period))
         ZscoreBuffer[i] = (close[i] - vct.Mean())/vct.Std();
      else
         ZscoreBuffer[i]=0.0;
     }

//--- return value of prev_calculated for next call
   return(rates_total);
  }
//+------------------------------------------------------------------+

从图中可以看出，指标值现在看起来更符合正态分布。

Zscore 指标

当 Z 值严重偏离 0，超过某个历史阈值时，就会产生交易信号。如果 Z 值为极负，则表示此时适合做多，反之则表示此时适合做空。这意味着我们需要两个阈值来发出买入和卖出信号。一个负数（用于买入）和一个正数（用于卖出）。有了入场之后，我们就可以考虑出场。一种方法是推导出另一套在特定位置（多头或空头）起作用的阈值。做空时，我们可以在 Z 值回到 0 时平仓。同样，当 Z 值从我们买入时的极端水平上升到 0 时，我们就会关闭多头仓位。

带阈值的指标

现在，我们使用指标 Zscore.ex5 定义了入场点和出场点。让我们把所有这些都整合到一个 EA 中。代码如下所示。

//+------------------------------------------------------------------+
//|                                                MeanReversion.mq5 |
//|                                  Copyright 2024, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2024, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#property version   "1.00"
#resource "\\Indicators\\Zscore.ex5"
#include<ExpertTools.mqh>
//---Input parameters
input int  PeriodLength  = 10;
input double LotsSize = 0.01;
input double  LongOpenLevel = -2.0;
input double  ShortOpenLevel = 2.0;
input double  LongCloseLevel = -0.5;
input double  ShortCloseLevel = 0.5;
input ulong  SlippagePoints = 10;
input ulong  MagicNumber    = 123456;
//---
 int indi_handle;
//---
 double zscore[2];
//+------------------------------------------------------------------+
//| Expert initialization function                                   |
//+------------------------------------------------------------------+
int OnInit()
  {
//---
   if(PeriodLength<2)
    {
     Print("Invalid parameter value for PeriodLength");
     return INIT_FAILED;
    }
//---
   if(!InitializeIndicator())
    return INIT_FAILED;
//---        
   return(INIT_SUCCEEDED);
  }
//+------------------------------------------------------------------+
//| Expert deinitialization function                                 |
//+------------------------------------------------------------------+
void OnDeinit(const int reason)
  {
//---
   
  }
//+------------------------------------------------------------------+
//| Expert tick function                                             |
//+------------------------------------------------------------------+
void OnTick()
  {
//---    
      int signal = GetSignal();
//---      
      if(SumMarketOrders(MagicNumber,_Symbol,-1))
       {
        if(signal==0)
         CloseAll(MagicNumber,_Symbol,-1);
        return; 
       }
      else
        OpenPosition(signal);
//---   
  }
//+------------------------------------------------------------------+
//+------------------------------------------------------------------+
//| Initialize indicator                                             |
//+------------------------------------------------------------------+
bool InitializeIndicator(void)
{
 indi_handle = INVALID_HANDLE;
//---
 int try = 10;
//---
 while(indi_handle == INVALID_HANDLE && try>0)
  {
   indi_handle = (indi_handle==INVALID_HANDLE)?iCustom(NULL,PERIOD_CURRENT,"::Indicators\\Zscore.ex5",PeriodLength):indi_handle;
   try--;
  }
//---
 if(try<0)
  {
   Print("Failed to initialize Zscore indicator ");
   return false;
  }
//---
 return true;
}
//+------------------------------------------------------------------+
//|Get the signal to trade or close                                  |
//+------------------------------------------------------------------+
int GetSignal(const int sig_shift=1)
{
//---
 if( CopyBuffer(indi_handle,int(0),sig_shift,int(2),zscore)<2)
   {
    Print(__FUNCTION__," Error copying from indicator buffers: ", GetLastError());
    return INT_MIN;
   } 
//---   
 if(zscore[1]<LongOpenLevel && zscore[0]>LongOpenLevel)
     return (1);
//---   
 if(zscore[1]>ShortOpenLevel && zscore[0]<ShortOpenLevel)
     return (-1);          
//---   
 if((zscore[1]>LongCloseLevel && zscore[0]<LongCloseLevel) ||
    (zscore[1]<ShortCloseLevel && zscore[0]>ShortCloseLevel))
     return (0);
//---
 return INT_MIN;
//--- 
}
//+------------------------------------------------------------------+
//|  Go long or short                                                |
//+------------------------------------------------------------------+
bool OpenPosition(const int sig)
{

 long pid;
//--- 
 if(LastOrderOpenTime(pid,NULL,MagicNumber)>=iTime(NULL,0,0))
   return false;
//---   
 if(sig==1)
   return SendOrder(_Symbol,0,ORDER_TYPE_BUY,LotsSize,SlippagePoints,0,0,NULL,MagicNumber);
 else
  if(sig==-1)
    return SendOrder(_Symbol,0,ORDER_TYPE_SELL,LotsSize,SlippagePoints,0,0,NULL,MagicNumber); 
//--- 
  return false;                   
}

它非常简单，没有定义止损或止盈水平。我们的目标是首先优化 EA，以获得 Zscore 指标的最佳周期数，以及最佳入场和退出阈值。我们将对几年的数据进行优化，并在样本外测试最优参数，但在此之前，我们要绕个小弯，介绍另一个有趣的工具。在《算法交易》一书中：作者欧内斯特.陈（Ernest Chan）在“制胜策略及其原理”一章中介绍了一种用于制定均值回复策略的有趣工具，即 "均值回归的半衰期"。

均值回归的半衰期

均值回归的半衰期是指偏离均值的程度减半所需的时间。就资产价格而言，均值回归的半衰期表示价格在偏离历史平均值后回归历史平均值的速度。它是衡量均值回归过程发生速度的指标。在数学上，半衰期与均值回归的速度可以用等式来表示：

均值回归的半衰期

其中：
- HL 是半衰期。
- log() 是自然对数。
- lambda 是均值回归的速度。

实际上，半衰期越短，意味着均值回归过程越快，而半衰期越长，意味着均值回归过程越慢。半衰期概念可用于微调均值回归交易策略中的参数，有助于根据历史数据和观察到的均值回归速度优化进场和退出点。均值回归的半衰期是从均值回归过程的数学表示中推导出来的，通常被模拟为 Ornstein-Uhlenbeck 过程。Ornstein-Uhlenbeck 过程是一个随机微分方程，它描述了均值回归行为的连续时间版本。

Chan 认为，可以通过计算均值回归的半衰期来确定均值回归是否是一种合适的策略。首先，如果 lambda 为正值，则根本不应采用均值回归法。即使 lambda 为负值且非常接近于零，也不鼓励应用均值回归，因为这表明半衰期会很长。只有当半衰期较短时，才应采用均值回归法。

均值回归的半衰期在 MeanReversionUtilities.mqh 中以函数形式实现，代码如下。其计算方法是将价格序列与后续值之间的差值序列进行回归。Lambda 等于回归模型的 beta 参数，半衰期的计算方法是用 -log(2) 除以 lambda。

//+------------------------------------------------------------------+
//|Calculate Half life of Mean reversion                             |
//+------------------------------------------------------------------+
double mean_reversion_half_life(vector &data, double &lambda)
  {
//---
   vector yy,zz;
   matrix xx;
//---
   OLS ols_reg;
//---
   yy.Init(data.Size()-1,slice,data,1,data.Size()-1,1);
//---
   zz.Init(data.Size()-1,slice,data,0,data.Size()-2,1);
//---
   if(!xx.Init(zz.Size(),2) || !xx.Col(zz,0) || !xx.Col(vector::Ones(zz.Size()),1) || !ols_reg.Fit(yy-zz,xx))
     {
      Print(__FUNCTION__," Error in calculating half life of mean reversion ", GetLastError());
      return 0;
     }
//---
   vector params = ols_reg.ModelParameters();
   lambda = params[0];
//---
   return (-log(2)/lambda);
//---
  }

我们将把它与 GHE 和 VRT 结合使用，以测试选定年份中几个外汇交易品种的价格样本。我们将利用测试结果来选择一个合适的交易品种，并在该交易品种上应用我们之前构建的 EA。它将在同一时期进行优化，并最终进行抽样测试。下面的脚本接受候选交易品种列表，这些交易品种将使用 GHE、VRT 和半衰期进行测试。

//+------------------------------------------------------------------+
//|                                                 SymbolTester.mq5 |
//|                                  Copyright 2024, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2024, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#property version   "1.00"
#property script_show_inputs
#include<MeanReversionUtilities.mqh>
#include<GHE.mqh>
#include<VRT.mqh>

//--- input parameters
input string   Symbols = "EURUSD,GBPUSD,USDCHF,USDJPY";//Comma separated list of symbols to test
input ENUM_TIMEFRAMES TimeFrame = PERIOD_D1;
input datetime StartDate=D'2020.01.02 00:00:01';
input datetime StopDate=D'2015.01.18 00:00:01';
input int Q_parameter = 2;
input int MinimumLag = 2;
input int MaximumLag = 100;
input bool ApplyLogTransformation = true;
//---
CVarianceRatio vrt;
double ghe,hl,lb,vlower,vupper;
double prices[];
//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
//---Check Size input value
   if(StartDate<=StopDate)
     {
      Print("Invalid input for StartDater or StopDate");
      return;
     }
//---array for symbols
   string symbols[];
//---process list of symbols from user input
   int num_symbols = StringSplit(Symbols,StringGetCharacter(",",0),symbols);
//---incase list contains ending comma
   if(symbols[num_symbols-1]=="")
      num_symbols--;
//---in case there are less than two symbols specified
   if(num_symbols<1)
     {
      Print("Invalid input. Please list at least one symbol");
      return;
     }
//---loop through all paired combinations from list
   for(uint i=0; i<symbols.Size(); i++)
     {

      //--- get prices for the pair of symbols
      if(CopyClose(symbols[i],TimeFrame,StartDate,StopDate,prices)<1)
        {
         Print("Failed to copy close prices ", ::GetLastError());
         return;
        }
      //---
      if(ApplyLogTransformation && !MathLog(prices))
        {
         Print("Mathlog error ", GetLastError());
         return;
        }
      //---
      if(!vrt.Vrt(prices,MinimumLag))
         return;
      //---
      vlower = vrt.VRatio();
      //---
      if(!vrt.Vrt(prices,MaximumLag))
         return;
      //---
      vupper = vrt.VRatio();
      //---
      ghe = general_hurst(prices,Q_parameter,MinimumLag,MaximumLag);
      //---
      hl = mean_reversion_half_life(prices,lb);
      //--- output the results
      Print(symbols[i], " GHE:  ", DoubleToString(ghe)," | Vrt: ",DoubleToString(vlower)," ** ",DoubleToString(vupper)," | HalfLife ",DoubleToString(hl)," | Lambda: ",DoubleToString(lb));
     }
  }
//+------------------------------------------------------------------+

运行脚本会产生以下结果：

19:31:03.143    SymbolTester (USDCHF,D1)        EURUSD GHE:  0.44755644 | Vrt: 0.97454284 ** 0.61945905 | HalfLife 85.60548208 | Lambda: -0.00809700
19:31:03.326    SymbolTester (USDCHF,D1)        GBPUSD GHE:  0.46304381 | Vrt: 1.01218672 ** 0.82086185 | HalfLife 201.38001205 | Lambda: -0.00344199
19:31:03.509    SymbolTester (USDCHF,D1)        USDCHF GHE:  0.42689382 | Vrt: 1.02233286 ** 0.47888803 | HalfLife 28.90550869 | Lambda: -0.02397976
19:31:03.694    SymbolTester (USDCHF,D1)        USDJPY GHE:  0.49198795 | Vrt: 0.99875744 ** 1.06103587 | HalfLife 132.66433924 | Lambda: -0.00522482

在所选日期期间，USDCHF 交易品种的测试结果最为理想。因此，我们将优化交易 USDCHF 的 EA 参数。一个有趣的做法是选择 Zscore 周期进行优化，看看它是否与计算出的半衰期不同。

抽样测试设置

样本内参数设置

在这里，我们可以看到最佳 Zscore 周期，它与计算出的均值回归半衰期非常接近，这是令人鼓舞的。当然，要确定半衰期是否有用，还需要进行更广泛的测试。

优化结果

样本内图形

样本回溯测试

最后，我们用最优参数对 EA 进行样本外测试。

样本外设置

结果看起来并不好。这可能是由于市场在不断变化，因此在 EA 优化期间观察到的特征已不再适用。我们需要更动态的进场和退出阈值，以考虑基本市场动态的变化。

样本外性能

我们可以将这里学到的东西作为进一步开发的基础。我们可以探索的一个途径是应用本文介绍的工具来实现配对交易策略。Zscore 指标不是基于单一价格序列，而是基于两个协整或相关工具的价差。

结论

在本文中，我们演示了广义赫斯特指数在 MQL5 中的实现，并展示了如何使用它来确定价格序列的特征。我们还研究了方差比检验的应用以及均值回归的半衰期。下表是文章所附所有文件的说明。

文件	描述
Mql5\include\ExpertTools.mqh	包含用于 MeanReversion EA 中进行交易操作的函数定义
Mql5\include\GHE.mqh	包含实现广义赫斯特指数函数的定义
Mql5\include\OLS.mqh	包含实现普通最小二乘法回归的 OLS 类的定义
Mql5\include\VRT.mqh	包含 CVarianceRatio 类的定义，该类封装了方差比检验
Mql5\include\VectorMatrixTools.mqh	具有各种函数定义，可快速初始化普通向量和矩阵
Mql5\include\TestUtilities.mqh	在 OLS 类定义中使用的一些声明
Mql5\include\MeanReversionUtilities.mqh	包含各种函数定义，包括一个实现均值回归半衰期的函数定义
Mql5\Indicators\Zscore.mq5	MeanReversion EA 中使用的指标
Mql5\scripts\SymbolTester.mq5	可用于测试交易品种均值回归的脚本
Mql5\Experts\GHE.ex5	可用于探索和实验 GHE 和 VRT 工具的 EA 交易应用程序
Mql5\scripts\MeanReversion.mq5	演示简单的均值回归策略的 EA