股票量化交易软件：利用卡尔曼 (Kalman) 滤波器预测价格方向

导言

货币和股票价格的图表总是包含价格波动, 其频率和幅度有所不同。我们的任务是判断基于这些短期和长期走势的主要趋势。一些交易者在图表上绘制趋势线, 而另一些人则使用指标。在这两种情况下, 我们的目的是将真正的价格走势从受到次要因素影响而导致的噪音中分离出来, 因为噪音只会产生短期效果。在本文中, 赫兹量化交易软件提议利用卡尔曼滤波器将主要走势与市场噪音分开。

在交易中使用数字滤波器的思路并不鲜见。例如, 我曾 描述过 运用低通滤波器。但追求完美是无止境的, 所以我们再考察一个策略, 比较一下结果。

1. 卡尔曼滤波器原理

那么, 什么是卡尔曼滤波器, 为什么我们感兴趣呢？以下过滤器定义来自 维基百科 :

卡尔曼滤波器 是一种使用一系列随时间观测到的测量值的算法, 包含统计噪声和其它不准确性。

这意味着该滤波器最初是为处理噪声数据而设计的。还有, 它能够处理不完整的数据。另一个优点, 它是为动态系统设计并应用的; 我们的价格图表恰好属于这样的系统。

滤波器算法的工作在两个步骤中处理:

1. 外推 (预测)
2. 更新 (校正)

1.1. 外推, 系统数值的预测

滤波器操作算法的第一阶段是利用已分析过程的基础模型。在此模型基础上, 形成单步前瞻预测。

(1.1)

其中:

• xk 是动态系统在第 k 步的外推值,
• Fk 是状态转换模型, 展体现当前系统状态对先前状态的依赖性,
• x^k-1 是系统的前一个状态 (前一步中的滤波值),
• Bk 是控制输入模型, 展现控制对系统的影响,
• uk 是系统上的控制向量。

例如, 控制效果可以是新闻因素。不过, 实际当中效果是未知的, 且被忽略, 而其影响是指噪声。

之后预测系统的协方差误差:

(1.2)

其中:

• Pk 是动态系统状态向量的外推协方差矩阵,
• Fk 是状态转换模型, 展体现当前系统状态对先前状态的依赖性,
• P^k-1 是状态向量的协方差矩阵在前一步的更新,
• Qk 是过程的协方差噪声矩阵。

1.2. 系统值的更新

滤波器算法的第二步从测量实际系统的状态 zk 开始。考虑到真实系统状态和测量误差, 指定系统状态的实际测量值。在赫兹量化交易软件的案例中, 测量误差是噪声对动态系统的影响。

此刻, 我们已有两个不同的数值代表单个动态过程的状态。它们包括第一步计算的动态系统外推值, 和实际的测量值。这些具有一定的几率度的数值, 当中的每一个均表征我们过程的真实状态, 因此, 该值介于这两个值之间。因此, 我们的目标是确定信任度, 即此值或彼值的信任程度。为此目的, 执行卡尔曼滤波器第二阶段的迭代。

利用已有数据, 我们判断实际系统状态自外推值的偏差。

(2.1)

此处:

• yk 是外推之后系统实际状态在第 k 步的偏差,
• zk 是第 k 步中系统的实际状态,
• Hk 是显示实际系统状态对于所计算数据依赖性的测量矩阵 (在实际中经常取值一),
• xk 是动态系统在第 k 步的外推值。

在下一步中, 计算误差向量的协方差矩阵:

(2.2)

此处:

• Sk 是在第 k 步的误差矢量的协方差矩阵,
• Hk 是显示实际系统状态对于计算数据依赖性的测量矩阵,
• Pk 是动态系统状态向量的外推协方差矩阵,
• Rk 是测量噪声的协方差矩阵。

然后检测优化增益。增益反映了计算值和经验值的置信度。

(2.3)

此处:

• Kk 是卡尔曼增益值的矩阵,
• Pk 是动态系统状态向量的外推协方差矩阵,
• Hk 是显示实际系统状态对于计算数据依赖性的测量矩阵,
• Sk 是在第 k 步的误差矢量的协方差矩阵。

现在, 我们使用卡尔曼增益来更新系统状态值和状态向量的协方差矩阵估值。

(2.4)

其中:

• x^k 和 x^k-1 在在第 k 和 k-1 步的更新值,
• Kk 是卡尔曼增益值的矩阵,
• yk 是外推后在第 k 步的系统实际状态的偏差。

(2.5)

其中:

• P^k 是更新后的动态系统状态向量的协方差矩阵,
• I 是标识符矩阵,
• Kk 是卡尔曼增益值的矩阵,
• Hk 是显示实际系统状态对于计算数据依赖性的测量矩阵,
• Pk 是外推的动态系统状态向量的协方差矩阵。

以上所有可以概括为以下规划

2. 卡尔曼滤波器的实现

现在, 我们已知晓了卡尔曼滤波器的工作原理。我们进入到实际实现。以上滤波器公式的矩阵形式允许接收若干个来源的数据。我建议在柱线收盘价基础上构建一个滤波器, 并将矩阵形式简化为离散的。

2.1. 输入数据初始化

在开始编写代码之前, 赫兹量化交易软件先定义输入数据。

如上所述, 卡尔曼滤波器的基础是一个动态过程模型, 用于预测过程的下一个状态。该滤波器最初旨在协同线性系统一起使用的, 其当前状态可轻易地通过前一状态乘以一个系数来定义。我们的情况有很更困难: 赫兹量化交易软件的动态系统是非线性的, 且比率在每一步都有变化。更甚的是, 我们不清楚该系统相邻状态之间的关系。这个任务似乎难以解决。这是一个棘手的解决方案: 我们将利用在这几篇文章中 [1],[2],[3] 描述的自回归模型。

我们开始吧。首先, 我们在这个类中声明 CKalman 类和所需的变量

class CKalman
  {
private:
//---
   uint              ci_HistoryBars;               //用于分析的柱线数量
   uint              ci_Shift;                     //计算自回归的偏移
   string            cs_Symbol;                    //品名
   ENUM_TIMEFRAMES   ce_Timeframe;                 //时间帧
   double            cda_AR[];                     //自回归系数
   int               ci_IP;                        //自回归系数的数量
   datetime          cdt_LastCalculated;           //最后计算时间;
   
   bool              cb_AR_Flag;                   //自回归计算的标志
//--- 卡尔曼滤波器的数值
   double            cd_X;                         // X
   double            cda_F[];                      // F 数组
   double            cd_P;                         // P
   double            cd_Q;                         // Q
   double            cd_y;                         // y
   double            cd_S;                         // S
   double            cd_R;                         // R
   double            cd_K;                         // K
   
public:
                     CKalman(uint bars=6240, uint shift=0, string symbol=NULL, ENUM_TIMEFRAMES period=PERIOD_H1);
                    ~CKalman();
   void              Clear_AR_Flag(void)  {  cb_AR_Flag=false; }
  };

我们在类的初始化函数中为变量分配初值。

CKalman::CKalman(uint bars, uint shift, string symbol, ENUM_TIMEFRAMES period)
  {
   ci_HistoryBars =  bars;
   cs_Symbol      =  (symbol==NULL ? _Symbol : symbol);
   ce_Timeframe   =  period;
   cb_AR_Flag     =  false;
   ci_Shift       =  shift;
   cd_P           =  1;
   cd_K           =  0.9;
  }

我使用了来自文章 [1] 的算法创建一个自回归模型。为此目的, 需要在类中添加两个私有函数。

   bool              Autoregression(void);
   bool              LevinsonRecursion(const double &R[],double &A[],double &K[]);

LevinsonRecursion 函数按原样使用。Autoregression 函数有略微修改, 所以我们来仔细考察这个函数。在函数伊始, 我们检查分析所需的历史数据的可用性。如果没有足够的历史数据, 则返回 false。

bool CKalman::Autoregression(void)
  {
   //--- 检查数据不足
   if(Bars(cs_Symbol,ce_Timeframe)<(int)ci_HistoryBars)
      return false;

现在, 我们加载所需的历史数据并填充实际状态转移模型系数的数组。

//---
   double   cda_QuotesCenter[];                          //计算的数据

//--- 令所有价格可用
   double close[];
   int NumTS=CopyClose(cs_Symbol,ce_Timeframe,ci_Shift+1,ci_HistoryBars+1,close)-1;
   if(NumTS<=0)
      return false;
   ArraySetAsSeries(close,true);
   if(ArraySize(cda_QuotesCenter)!=NumTS)
     {
      if(ArrayResize(cda_QuotesCenter,NumTS)<NumTS)
         return false;
     }
   for(int i=0;i<NumTS;i++)
      cda_QuotesCenter[i]=close[i]/close[i+1];           // 计算系数

在准备操作之后, 我们检测自回归模型的系数个数, 并计算它们的值。

   ci_IP=(int)MathRound(50*MathLog10(NumTS));
   if(ci_IP>NumTS*0.7)
      ci_IP=(int)MathRound(NumTS*0.7);                         // 自回归模型的顺序
  
   double cor[],tdat[];
   if(ci_IP<=0 || ArrayResize(cor,ci_IP)<ci_IP || ArrayResize(cda_AR,ci_IP)<ci_IP || ArrayResize(tdat,ci_IP)<ci_IP)
      return false;
   double a=0;
   for(int i=0;i<NumTS;i++)
      a+=cda_QuotesCenter[i]*cda_QuotesCenter[i];    
   for(int i=1;i<=ci_IP;i++)
     {  
      double c=0;
      for(int k=i;k<NumTS;k++)
         c+=cda_QuotesCenter[k]*cda_QuotesCenter[k-i];
      cor[i-1]=c/a;                                            // 自回归
     } 
  
   if(!LevinsonRecursion(cor,cda_AR,tdat))                     // 莱文森-德宾 (Levinson-Durbin) 递归
      return false;

现在我们将自回归系数的总和降低到 "1", 并将计算执行的标志设置为 'true'。

   double sum=0;
   for(int i=0;i<ci_IP;i++)
     {
      sum+=cda_AR[i];
     }
   if(sum==0)
      return false;
  
   double k=1/sum;
   for(int i=0;i<ci_IP;i++)
      cda_AR[i]*=k;

   cb_AR_Flag=true;

接下来, 我们初始化滤波器所需的变量。为了计算噪声协方差, 我们使用所在分析周期的 Close 值的偏差的均方根值。

   cd_R=MathStandardDeviation(close);

为了确定过程噪声协方差的值, 赫兹量化交易软件首先计算自回归模型值的数组, 并找出模型值的均方根偏差。

   double auto_reg[];
   ArrayResize(auto_reg,NumTS-ci_IP);
   for(int i=(NumTS-ci_IP)-2;i>=0;i--)
     {
      auto_reg[i]=0;
      for(int c=0;c<ci_IP;c++)
        {
         auto_reg[i]+=cda_AR[c]*cda_QuotesCenter[i+c];
        }
     }
   cd_Q=MathStandardDeviation(auto_reg);

然后, 我们将实际的状态转换系数复制到 cda_F 数组, 从其可以进一步使用它们来计算新的系数。

   ArrayFree(cda_F);
   if(ArrayResize(cda_F,(ci_IP+1))<=0)
      return false;
   ArrayCopy(cda_F,cda_QuotesCenter,0,NumTS-ci_IP,ci_IP+1);

对于我们系统的初始值, 我们使用的是最后 10 个值的算术平均值。

   cd_X=MathMean(close,0,10);

2.2. 价格走势预测

在收到滤波器操作所需的所有初始数据之后, 我们可以继续进行实际实现。卡尔曼滤波器操作的第一步是 单步前瞻系统状态预测。我们创建一个 Forecast 公有函数, 在其内我们会实现函数 1.1 和 1.2。

double            Forecast(void);

在函数伊始, 我们检查回归模型是否已经被计算。应在必要时调用其计算函数。在模型重新计算出错的情况下返回 EMPTY_VALUE,

double CKalman::Forecast()
  {
   if(!cb_AR_Flag)
     {
      ArrayFree(cda_AR);
      if(Autoregression())
        {
         return EMPTY_VALUE;
        }
     }

之后, 赫兹量化交易软件计算状态转换系数并将其保存到 cda_F 数组的 "0" 单元中, 随后数值依次顺移一个单元。

   Shift(cda_F);
   cda_F[0]=0;
   for(int i=0;i<ci_IP;i++)
      cda_F[0]+=cda_F[i+1]*cda_AR[i];

然后我们重新计算系统状态和出错概率。

   cd_X=cd_X*cda_F[0];
   cd_P=MathPow(cda_F[0],2)*cd_P+cd_Q;

该函数在最后返回预测的系统状态。在我们的情况中, 这是一根新柱线的预测收盘价。

   return cd_X;
  }

2.3. 校正系统状态

在下一个阶段, 在收到实际的柱线收盘价之后, 我们校正系统状态。为此目的, 我们来创建 Correction 公有函数。在函数参数中, 我们将传递实际的系统状态值, 即实际的柱线收盘价格。

double            Correction(double z);

在该函数中实现了文章给出的 理论部分1.2。其完整代码可在附件中找到。在操作结束时, 该函数返回系统状态的更新 (校正) 值。

3. 卡尔曼滤波器的实际演示

我们来测试这个基于卡尔曼滤波器的类是如何工作的。我们创建一个基于这个类的指标。在新的蜡烛条开盘时, 指标调用系统更新函数, 然后调用函数来预测当前柱线的收盘价。类函数的调用顺序是相反的, 因为我们要为前一根收盘柱线调用更新 (校正) 函数, 再来预测当前的新柱线, 而其收盘价尚未得知。

指标有两个缓冲区。系统状态的预测值将被添加到第一个缓冲区之中, 更新的数值将被添加到第二个缓冲区中。我有意使用两个缓冲区, 令指标不会重绘, 我们可以看到系统在第二个滤波器的操作阶段是如何更新 (校正) 地。指标代码很简单, 可以在下面的附件中找到。这是指标操作的结果。

图表上显示三条虚线:

• 黑线显示实际的柱线收盘价
• 红线显示预测值
• 蓝线是由卡尔曼滤波器更新的系统状态

如您所见, 两条线都接近实际的收盘价, 且其显示的翻转点几率优良。请注意, 指标不会重绘, 并且在收盘价格尚不明了时, 红线在柱线开盘时即已绘制。

此图表显示了此滤波器的一致性, 以及利用此滤波器创建交易系统的可行性。