移動エントロピーを用いた時系列の因果分析

はじめに

移動エントロピーは、ある時系列から別の時系列へ伝達される情報量を定量化するための統計的ツールであり、対象となる変数の性質や挙動に関する深い洞察を提供します。本記事では、移動エントロピーを用いた統計的因果関係の概念を掘り下げ、さまざまなプロセス間の因果的影響の方向性がどのように明らかになるかを探ります。また、移動エントロピーを測定するためのMQL5の実装について詳述し、この手法が潜在的に結合している時系列の分析にどのように応用できるかを実証します。移動エントロピーを活用することで、予測タスクの向上に役立つ変数を特定することを目指します。

因果関係

実証的なデータは誤解を招くことがあります。2つの変数が連動して動いているように見えるからといって、一方が他方を引き起こしているわけではありません。だからこそ、「相関関係は因果関係ではない」という言葉は真実なのです。相関は単に2つの変数がどの程度結びついているかを測るものであり、なぜ結びついているかを測るものではありません。例えば、夏のアイスクリームの売上と株価に強い相関関係があるとしましょう。アイスクリームを買えば株が上がるということではありません。より可能性が高いのは、シーズンそのものが、両方の変数に独立して影響を与えているような、隠れた要因です。同様に、ある企業の株価と金価格の間には関連性があるかもしれませんが、本当の原因は、市場全体の感情やインフレが両方の価格に影響を与えるなど、まったく別のものかもしれません。これらの例は、相関データが誤解を招く可能性があることを浮き彫りにしています。つながりを示していますが、その背後にある理由は示していません。あることが別のことを引き起こすかどうかを本当に理解するには、より高度なツールが必要です。

ある事象が別の事象を引き起こすという因果関係の概念は、科学的探求の基本です。しかし、因果関係を正確に定義することは、哲学的、物理学的、統計学的に深い考察を要する多面的な課題です。理想的なのは、原因が必ず特異な結果を生むことです。しかし、結果に影響を及ぼす複雑な影響の網の目から、単一の原因因子を分離することは難しいです。例えば、取引量の急増は株価の上昇と相関するかもしれませんが、市場心理や経済データの発表など、他の要因も重要な役割を果たす可能性があります。このような場合、研究者は因果関係を推測するために統計的手法を用います。

因果関係の性質は、決定論的（結果が保証されている）であれ、確率論的（その後の出来事の可能性に影響を与える）であれ、その根底にあるプロセスにかかっています。決定論的なシステムでは、落下する物体の予測可能な落下で観察されるように、最初の事象は明らかに2番目の事象につながります。逆に、確率的システムでは、最初の事象が2番目の事象の発生を予測する能力を高めるかどうかに焦点が移ります。例えば、最近の降雨はその後の開花に関係しているかもしれませんが、他の環境要因も関係している可能性があります。このような場合、最初の事象に関する知識が2番目の事象を予測する能力を向上させるかどうかが問題となります。

経済学者のクライヴ・グレンジャーは、ノーバート・ウィーナーの研究を基に因果関係の概念を発展させ、第1の信号が第2の信号を引き起こすのは、第2の信号の遅延値のみを使用するのではなく、第1の信号と第2の信号の両方からの過去の情報を使用して第2の信号の将来の値をより適切に説明できる場合であると仮定しました。グレンジャーは因果関係の定義を2つの原則に基づいておこないました。第一に、結果はその原因の前に現れることはありません。第二に、原因には結果に伝達される固有の情報が含まれます。これらの原則は、因果関係を定量化するためには、関係する変数の時間的特性を理解する必要があること、また、それらの情報量の尺度を理解する必要があることを示唆しています。このため、時系列は因果分析に適しています。

時系列データの性質上、ある時点のある系列の情報が、後の時点の別の系列の予測可能性にどのような影響を与えるかを分析することができます。グレンジャーは因果関係を不確実性の低減と定義しています。系列Xの過去の値を知ることで、系列Y自体の過去の値だけを使用する場合と比較して系列Yの将来の値の予測が改善される場合は、XはYを予測できると言われます。この考えに基づいて、グレンジャーは遅延時系列と自己回帰モデルを使用して因果関係のテストを開発しました。彼は、Xの過去の値を含めることでYの将来の値の予測が著しく改善されない限り、XとYの間に因果関係はないと提唱しました。この改善は、通常、回帰モデルにおけるより良い適合のような予測誤差の減少によって測定されます。グレンジャー因果関係は、時系列間の関係を統計的に検出することができます。しかし、その限界に注意することは重要です。グレンジャー因果関係は、方向性のある関係を特定するだけであり、必ずしも決定的な因果メカニズムを特定するものではありません。私たちがデータとして持っているもの以外にも原因があるかもしれません。さらに、グレンジャー因果関係は基本的に自己回帰の枠組みに基づいているため、線形因果関係を明らかにするのに最も効果的です。非線形の因果関係は、別のアプローチを必要とします。

数学的には、グレンジャー因果関係は、2つの時系列XとYを考慮することで表現できます。それぞれの遅延値は、kでの遅延を表すX(t−k)とY(t−k) で表されます。考慮される最大ラグはpとして示されます。自己回帰モデルを適用すると、Yの将来の値は自身の過去の値に回帰されます。

この式は、Yの将来の値を過去の値のみで考えるものです。Xをモデルに導入することで、将来の値はXとYの両方の過去の値で表現されます。

Xの過去値を含めることで、Yの過去値のみを用いたモデルと比較してYの予測値が有意に改善される場合、XはYに対してグレンジャー因果があると言ええます。これは通常、係数が共同でゼロであるという帰無仮説を検定することで評価されます。この帰無仮説が棄却されれば、XはYの過去の値だけに含まれる以上の、Yに関する有意な予測情報を提供していることになります。XにYに対してのグレンジャー因果があるかどうかを検証するために、仮説検定の下で2つのモデルを比較します。

• 帰無仮説：XにYに対してのグレンジャー因果がない

• 代替仮説：XにYに対してのグレンジャー因果がある

F検定は、それぞれのモデルの残差を調べることによって、制限付きモデル（Xなし）と制限なしモデル（Xあり）の適合度を比較するために使用されます。制限付き残差2乗和は、Xを含まないモデルからの残差であり、制限なし残差2乗和は、Xを含むモデルからの残差です。F検定は次のように計算されます。

ここで、nは観測の数です。計算されたF統計量は、自由度がpと「n - 2p - 1」のF分布の臨界値と比較されます。F統計量が臨界値より大きければ、帰無仮説を棄却し、XにYに対してのグレンジャー因果があると結論づけます。あるいは、F統計量は、一元配置分散分析（ANOVA）検定を用いて計算することもできます。計算式は以下の通りです。

移動エントロピー

情報理論の黎明期、科学者たちは相互情報を使って、結合プロセスがどのように相互作用するかを理解しました。クロード・シャノンのエントロピーに基づくこの概念は、ある時系列内の情報が別の時系列と重複しているかどうかを教えてくれます。もっと簡単に言えば、両方の系列を別々に符号化するよりも少ない情報で一緒に符号化できるかどうかを明らかにします。このため、相互情報は冗長性と呼ばれることもあります。あるプロセスは別のプロセスと情報を共有し、最初のプロセスからすでに取り込まれた情報を再利用することで、2番目のプロセスを効率的に記述することができます。

形式的には、2つの確率過程X(t)とY(t)の間の相互情報は、それらの限界エントロピーの和が結合系の結合エントロピーを超えるときに成立します。この数学的関係は、個々のプロセスに比べ、複合システムに関する不確実性の減少を反映しています。言い換えれば、一方のプロセスに関する情報が、もう一方のプロセスに関連する固有のエントロピーを低減するためにどの程度利用できるかを捉えるものです。エントロピーは基礎となる確率分布によってのみ決まるので、そのような分布はすべて、関連するエントロピー値によって特徴付けることができます。この値は、既知の確率分布が与えられた場合に、特定の結果に関連する意外性のレベルを定量化したものです。

この概念は、グレンジャー因果関係の文脈で特に重要になります。時系列間の潜在的な因果関係を調査する場合、その目的は、潜在的なソースプロセスからの情報を組み込むことによって、ターゲットプロセスに関連する不確実性を低減することです。副次的な時系列が含まれることで、対象プロセス分布のエントロピーが減少することが証明された場合、ソース系列から対象系列への統計的因果関係の影響が存在することを示唆します。この減少は、移動エントロピー（英語）と呼ばれます。

移動エントロピー(TE: Transfer Entropy)は、2つの時系列間の情報伝達の方向を測定するために、カルバック・ライブラー情報量の概念に基づいています。具体的には、TEは条件付き相互情報の考え方に基づいており、カルバック・ライブラー情報量(KL)（または相対エントロピーとしても知られる）を用いて表現することができます。KL情報量は、2つの確率分布の差を測定します。TEでは、KL距離は、Yの現在の状態とXとYの過去の状態の結合確率分布と、これらの状態の周辺分布の積との差を測定します。数学的には、時系列Xから時系列Yへの情報の伝達は次のように表すことができます。

ここで、y(t+1) はYの将来の状態、 y(t)はYの過去の状態、x(t)はXの過去の状態です。この定式化は、移動エントロピーが、y(t)に加えてx(t)からの情報を考慮したときに y(t+1)​の確率分布がどれだけ変化するかを測定することを強調しています。

2009年、Lionel Barnett、Adam Barrett、Anil Sethは共著で論文「Granger Causality and Transfer Entropy Are Equivalent for Gaussian Variables」を発表し、時系列がガウス分布に従う場合、移動エントロピーは グレンジャー因果性のF統計量の半分に相当することを実証しました。

この結果は、後にコードで実装する線形移動エントロピーの定義となります。非線形の因果関係を説明するために、Thomas Schreiberの研究に従って不確実性の削減の概念を拡張し、時系列を異なる遷移確率分布を持つマルコフ過程として扱います。

不確実性の低減をモデル化するSchreiberのアプローチは、時系列X(t)とY(t)を、既知の遷移確率p(x)とq(x)を持つマルコフ過程として扱うことで、情報理論を活用しています。線形モデルに依存するグレンジャーの自己回帰モデルとは異なり、このアプローチでは条件付き相互情報を用いて情報伝達を記述します。相互情報はエントロピーの差から得られるので、条件付き相互情報は、各エントロピー項を付加情報に条件付けることによって得られます。そして、条件付き相互情報量の式にラグ付き変数を代入することで移動エントロピーが計算され、条件付き相互エントロピーを用いて特定のラグkにおけるX(t)からY(t)への情報伝達を分析することができます。

計算上、共同エントロピーは確率分布しか必要としないので、この方法は魅力的です。一つのラグkに対する移動エントロピーは、4つの別々の結合エントロピー項として表すことができ、これはデータから正確な確率分布で簡単に計算できます。この式の利点は、より多くの遅行次元を扱えることです。しかし、ラグが1つ増えるごとに状態空間の次元が2つ増え、確率密度の推定に関連する有限データの問題が指数関数的に増大するため、移動エントロピーを正確に定量化する能力に大きな影響を与えます。

このアプローチの重要な強みは、そのノンパラメトリックな性質にあります。他の手法と異なり、定常性以上の基礎的なデータ分布に関する仮定をしないため、データ生成プロセスに関する予備知識がなくても適用できます。しかし、この利点には注意点があります。結果は、基礎となる分布の正確な推定に大きく依存します。移動エントロピーは、4つのエントロピーの項を計算するために、限られたデータを用いて、関係する確率過程の真の確率分布を近似する必要があります。この推定の精度は、移動エントロピーの所見の信頼性に大きく影響します。このことを念頭に置くと、計算されたエントロピー値がスプリアスである可能性を考慮しなければなりません。したがって、結果の頑健性を判断する何らかの方法が必要です。

移動エントロピーの推定にノンパラメトリックアプローチを採用するという私たちの主張には、その結果がゴミではなく真実を伝えるものであることを保証するという大きな課題が伴います。したがって、移動エントロピーの解釈には、推定値の統計的有意性を評価するという、より有益なアプローチを検討するのが得策です。一般的な有意性検定では、時系列データをあらかじめ定義された回数シャッフルします。その後、シャッフルされた各バージョンについて移動エントロピーが計算されます。p値はその後、移動エントロピーが元の値より低いシャッフルされたデータの割合として計算されます。

別のアプローチでは、シャッフルされたデータの平均値から結果が何標準偏差離れているかを計算する必要があります。シャッフルは時間構造を破壊するので、シャッフルされた移動エントロピー値の平均はゼロに近くなると予想されます。この平均値を中心としたデータの広がりは、元の結果の重要性を反映しています。算出された値はz得点と呼ばれます。z得点は通常、p値に比べてシャッフルの回数が少なく、計算効率が高いです。 

p値の場合、目標は可能な限りゼロに近い確率を得ることです。統計的有意性を示すz得点は3.0以上でなければなりません。

MQL5での実装

移動エントロピーの定量化とその重要性を決定するツールを実装したコードは、transfer_entropy.mqhに収められています。このファイルには、CTransEntropyというクラスの定義と、その他のヘルパークラスや関数が含まれています。この授業では、時系列データの統計分析の枠組みを提供し、特に変数間の因果関係を評価することを目的とします。線形グレンジャー因果関係（線形移動エントロピー）と非線形移動エントロピーを定量化するための2つの異なるメソッドが公開されています。これは双方向で計算され、変数間の情報の流れをより完全に把握することができます。

データの潜在的な非定常性に対処するために、このクラスはウィンドウ処理プロシージャを組み込んでいます。ユーザーは、ウィンドウサイズとストライドを定義することができ、より小さく重なり合ったセグメントでデータを分析することができます。このアプローチは、各ウィンドウに固有の結果をもたらし、因果関係の強さの時間的変動を特定することを容易にします。さらに、非定常データの分析に伴う課題も軽減されます。このクラスはまた、統合された有意性検定メカニズムも提供します。ユーザーは、周辺分布を保持しながら実行するデータシャッフルの回数を指定することができます。これらのシャッフルされたデータセットに基づいて、このクラスは各方向の移動エントロピーのp値とz得点を計算します。これらの統計値は、観察された因果関係や情報伝達が偶然によるものである可能性についての本質的な洞察を提供し、分析の頑健性を高めます。

クラスのインスタンスは、パラメータなしのデフォルトコンストラクタを使ってインスタンス化されます。

public:
                     CTransEntropy(void)
     {
      if(!m_transfer_entropies.Resize(2))
         Print(__FUNCTION__, " error ", GetLastError());

     }

このメソッドでは、与えられたデータセットでオブジェクトを初期化し、分析のためのさまざまなパラメータを設定します。

bool              Initialize(matrix &in, ulong endog_index, ulong exog_index, ulong lag, bool maxLagOnly=true, ulong winsize=0,ulong winstride=0)
     {
      if(!lag || lag>in.Rows()/2)
        {
         Print(__FUNCTION__, " Invalid parameter(s) : lag must be > 0  and < rows/2");
         return false;
        }

      if(endog_index==exog_index)
        {
         Print(__FUNCTION__, " Invalid parameter(s) : endog cannot be = exog ");
         return false;
        }

      if(!m_dataset.Resize(in.Rows(),2))
        {
         Print(__FUNCTION__, " error ", GetLastError());
         return false;
        }

      if(!m_dataset.Col(in.Col(endog_index),0) || !m_dataset.Col(in.Col(exog_index),1))
        {
         Print(__FUNCTION__, " error ", GetLastError());
         return false;
        }

      if(!m_wins.Initialize(m_dataset,lag,maxLagOnly,winsize,winstride))
         return false;

      m_tlag = lag;
      m_endog = endog_index;
      m_exog = exog_index;
      m_maxlagonly = maxLagOnly;

      return true;
     }

最初に必要なパラメータは、少なくとも2つの列を持つ行列であり、分析対象の時系列は入力行列の列に配置されなければなりません。非定常データを扱う場合は、事前にデータを差分することを推奨します。2番目と3番目のパラメータは入力データ行列の列添字で、それぞれ内生時系列（従属時系列）と外生時系列（独立時系列）を示します。

第4のパラメータであるlagは、分析で考慮されるラグパラメータを定義します。次のブーリアンパラメータ、maxLagOnlyは、ラグが単一の項を定義するか（trueの場合）、ラグまでのすべてのラグ付き値を定義するか（falseの場合）を決定します。最後から2番目のパラメータwinsizeはウィンドウの長さを表します。0に設定すると、データにウィンドウ処理は適用されません。最後に、winstrideはオプションでウィンドウ処理のウィンドウのストライドを設定し、時系列データ上を通過する連続ウィンドウ間のステップを定義します。

この方法は、内生的指数と外生的指数が同じでないことを確認することから始まります。もしそうなら、エラーメッセージを表示してfalseを返します。内部行列m_datasetは、解析対象の二変量データセットを格納するためにサイズ変更されます。そして、入力行列からendog_indexとexog_indexで指定された列を、それぞれm_datasetの1列目と2列目にコピーします。ウィンドウ処理が要求された場合、ヘルパークラスCDataWindowsがm_dataset行列のウィンドウ処理に使用されます。これが完了すると、メソッドは、後で使用するために、提供されたパラメータで内部変数を設定します。

//+------------------------------------------------------------------+
//|class that generates windows of the dataset to be analyzed        |
//+------------------------------------------------------------------+
class CDataWindows
  {
private:
   matrix m_dwins[],
          m_data;
   ulong  m_lag,
          m_win_size,
          m_stride_size;

   bool m_max_lag_only,
        m_has_windows;

   matrix            applylags(void)
     {
      matrix out=np::sliceMatrixRows(m_data,m_lag);

      if(m_max_lag_only)
        {
         if(!out.Resize(out.Rows(),m_data.Cols()+2))
           {
            Print(__FUNCTION__, " error ", GetLastError());
            return matrix::Zeros(1,1);
           }

         for(ulong i = 2; i<4; i++)
           {
            vector col = m_data.Col(i-2);
            col = np::sliceVector(col,0,col.Size()-m_lag);

            if(!out.Col(col,i))
              {
               Print(__FUNCTION__, " error ", GetLastError());
               return matrix::Zeros(1,1);
              }
           }
        }
      else
        {
         if(!out.Resize(out.Rows(),m_data.Cols()+(m_lag*2)))
           {
            Print(__FUNCTION__, " error ", GetLastError());
            return matrix::Zeros(1,1);
           }

         for(ulong i = 0,k = 2; i<2; i++)
           {
            for(ulong t = 1; t<(m_lag+1); t++,k++)
              {
               vector col = m_data.Col(i);
               col = np::sliceVector(col,m_lag-t,col.Size()-t);

               if(!out.Col(col,k))
                 {
                  Print(__FUNCTION__, " error ", GetLastError());
                  return matrix::Zeros(1,1);
                 }
              }

           }
        }

      return out;

     }

   bool              applywindows(void)
     {
      if(m_dwins.Size())
         ArrayFree(m_dwins);

      for(ulong i = (m_stride_size+m_win_size); i<m_data.Rows(); i+=ulong(MathMax(m_stride_size,1)))
        {
         if(ArrayResize(m_dwins,int(m_dwins.Size()+1),100)<0)
           {
            Print(__FUNCTION__," error ", GetLastError());
            return false;
           }
         m_dwins[m_dwins.Size()-1] = np::sliceMatrixRows(m_data,i-m_win_size,(i-m_win_size)+m_win_size);
        }

      return true;
     }


public:
                     CDataWindows(void)
     {

     }

                    ~CDataWindows(void)
     {

     }

   bool              Initialize(matrix &data, ulong lag, bool max_lag_only=true, ulong window_size=0, ulong window_stride =0)
     {
      if(data.Cols()<2)
        {
         Print(__FUNCTION__, " matrix should contain at least 2 columns ");
         return false;
        }

      m_data = data;

      m_max_lag_only = max_lag_only;

      if(lag)
        {
         m_lag = lag;
         m_data = applylags();
        }

      if(window_size)
        {
         m_win_size = window_size;
         m_stride_size = window_stride;
         m_has_windows = true;
         if(!applywindows())
            return false;
        }
      else
        {
         m_has_windows = false;

         if(m_dwins.Size())
            ArrayFree(m_dwins);

         if(ArrayResize(m_dwins,1)<0)
           {
            Print(__FUNCTION__," error ", GetLastError());
            return false;
           }

         m_dwins[0]=m_data;
        }

      return true;
     }

   matrix            getWindowAt(ulong ind)
     {
      if(ind < ulong(m_dwins.Size()))
         return m_dwins[ind];
      else
        {
         Print(__FUNCTION__, " Index out of bounds ");
         return matrix::Zeros(1,1);
        }
     }

   ulong             numWindows(void)
     {
      return ulong(m_dwins.Size());
     }

   bool              hasWindows(void)
     {
      return m_has_windows;
     }
  };

Initialize()メソッドが正常に完了したら、Calculate_Linear_TE()または Calculate_NonLinear_TE()を呼び出して、それぞれ線形および非線形の移動エントロピーをテストすることができます。どちらのメソッドも、完了するとブール値を返します。メソッドCalculate_Linear_TE()は、オプションのパラメータ n_shufflesを1つ取ることができます。n_shuffleが0（デフォルト）の場合、有意差検定はおこなわれません。

bool              Calculate_Linear_TE(ulong n_shuffles=0)
     {
      ulong c = m_wins.numWindows();

      matrix TE(c,2);
      matrix sTE(c,2);
      matrix pvals(c,2);
      matrix zscores(c,2);

      for(ulong i=0; i<m_wins.numWindows(); i++)
        {
         matrix df = m_wins.getWindowAt(i);

         m_transfer_entropies[0] = linear_transfer(df,0,1);

         m_transfer_entropies[1] = linear_transfer(df,1,0);


         if(!TE.Row(m_transfer_entropies,i))
           {
            Print(__FUNCTION__, " error ", GetLastError());
            return false;
           }

         SigResult rlts;

         if(n_shuffles)
           {
            significance(df,m_transfer_entropies,m_endog,m_exog,m_tlag,m_maxlagonly,n_shuffles,rlts);

            if(!sTE.Row(rlts.mean,i) || !pvals.Row(rlts.pvalue,i) || !zscores.Row(rlts.zscore,i))
              {
               Print(__FUNCTION__, " error ", GetLastError());
               return false;
              }

           }

        }

      m_results.TE_XY = TE.Col(0);
      m_results.TE_YX = TE.Col(1);
      m_results.p_value_XY = pvals.Col(0);
      m_results.p_value_YX = pvals.Col(1);
      m_results.z_score_XY = zscores.Col(0);
      m_results.z_score_YX = zscores.Col(1);
      m_results.Ave_TE_XY = sTE.Col(0);
      m_results.Ave_TE_YX = sTE.Col(1);

      return true;
     }

このメソッドでは、グランジャーズ法を用いて線形移動エントロピーを計算します。これはprivateメソッドlinear_transfer()で実装されています。このルーチンの最後の2つのパラメータは、入力行列の従属変数と独立変数（列）を特定します。列のインデックスを入れ替えてメソッドを2回呼び出すことで、両方向の移動エントロピーを得ることができます。

double            linear_transfer(matrix &testdata,long dep_index, long indep_index)
     {
      vector joint_residuals,independent_residuals;
      double entropy=0.0;

      OLS ols;

      double gc;
      vector y;
      matrix x,xx;

      matrix joint;
      if(m_maxlagonly)
         joint = np::sliceMatrixCols(testdata,2);
      else
        {
         if(!joint.Resize(testdata.Rows(), testdata.Cols()-1))
           {
            Print(__FUNCTION__, " error ", GetLastError());
            return entropy;
           }
         matrix sliced = np::sliceMatrixCols(testdata,2);
         if(!np::matrixCopyCols(joint,sliced,1) || !joint.Col(testdata.Col(indep_index),0))
           {
            Print(__FUNCTION__, " error ", GetLastError());
            return entropy;
           }
        }
      matrix indep = (m_maxlagonly)?np::sliceMatrixCols(testdata,dep_index+2,dep_index+3):np::sliceMatrixCols(testdata,(dep_index==0)?2:dep_index+m_tlag+1,(dep_index==0)?2+m_tlag:END);

      y = testdata.Col(dep_index);

      if(dep_index>indep_index)
        {
         if(m_maxlagonly)
           {
            if(!joint.SwapCols(0,1))
              {
               Print(__FUNCTION__, " error ", GetLastError());
               return entropy;
              }
           }
         else
           {
            for(ulong i = 0; i<m_tlag; i++)
              {
               if(!joint.SwapCols(i,i+m_tlag))
                 {
                  Print(__FUNCTION__, " error ", GetLastError());
                  return entropy;
                 }
              }
           }
        }

      if(!addtrend(joint,xx))
         return entropy;

      if(!ols.Fit(y,xx))
         return entropy;

      joint_residuals = ols.Residuals();

      if(!addtrend(indep,x))
         return entropy;

      if(!ols.Fit(y,x))
         return entropy;

      independent_residuals = ols.Residuals();

      gc = log(independent_residuals.Var()/joint_residuals.Var());

      entropy = gc/2.0;

      return entropy;

     }

Calculate_NonLinear_TE()メソッドは、n_shuffleと一緒にnumBinsという追加パラメータを取ります。このパラメータは、変数の確率密度の推定に使用するビンの数を定義します。 

bool              Calculate_NonLinear_TE(ulong numBins, ulong n_shuffles=0)
     {
      ulong c = m_wins.numWindows();

      matrix TE(c,2);
      matrix sTE(c,2);
      matrix pvals(c,2);
      matrix zscores(c,2);

      for(ulong i=0; i<m_wins.numWindows(); i++)
        {
         matrix df = m_wins.getWindowAt(i);

         m_transfer_entropies[0] = nonlinear_transfer(df,0,1,numBins);

         m_transfer_entropies[1] = nonlinear_transfer(df,1,0,numBins);


         if(!TE.Row(m_transfer_entropies,i))
           {
            Print(__FUNCTION__, " error ", GetLastError());
            return false;
           }

         SigResult rlts;

         if(n_shuffles)
           {
            significance(df,m_transfer_entropies,m_endog,m_exog,m_tlag,m_maxlagonly,n_shuffles,rlts,numBins,NONLINEAR_TE);

            if(!sTE.Row(rlts.mean,i) || !pvals.Row(rlts.pvalue,i) || !zscores.Row(rlts.zscore,i))
              {
               Print(__FUNCTION__, " error ", GetLastError());
               return false;
              }

           }

        }

      m_results.TE_XY = TE.Col(0);
      m_results.TE_YX = TE.Col(1);
      m_results.p_value_XY = pvals.Col(0);
      m_results.p_value_YX = pvals.Col(1);
      m_results.z_score_XY = zscores.Col(0);
      m_results.z_score_YX = zscores.Col(1);
      m_results.Ave_TE_XY = sTE.Col(0);
      m_results.Ave_TE_YX = sTE.Col(1);

      return true;


     }

確率密度の推定にはヒストグラム法を用います。この方法が選ばれたのは、実装が最も簡単だからです。一般化された移動エントロピーの計算の責任は、privateメソッドであるnonlinear_entropy()とget_entropy()に委ねられています。 

double            get_entropy(matrix &testdata, ulong num_bins)
     {

      vector hist;
      vector bounds[];
      hist=vector::Ones(10);

      if(!np::histogramdd(testdata,num_bins,hist,bounds))
        {
         Print(__FUNCTION__, " error ");
         return EMPTY_VALUE;
        }

      vector pdf = hist/hist.Sum();
      vector lpdf = pdf;

      for(ulong i = 0; i<pdf.Size(); i++)
        {
         if(lpdf[i]==0.0)
            lpdf[i] = 1.0;
        }

      vector ent = pdf*log(lpdf);

      return -1.0*ent.Sum();

     }

結合条件付きエントロピーと独立条件付きエントロピーを計算するために使用された4つの成分値は、 nonlinear_transfer()で結合され、最終的な推定値を得ます。

double            nonlinear_transfer(matrix &testdata,long dep_index, long indep_index, ulong numbins)
     {
      double entropy=0.0;

      matrix one;
      matrix two;
      matrix three;
      matrix four;

      if(m_maxlagonly)
        {
         if(!one.Resize(testdata.Rows(),3) || !two.Resize(testdata.Rows(),2) || !three.Resize(testdata.Rows(),2) || !four.Resize(testdata.Rows(),1) ||
            !one.Col(testdata.Col(dep_index),0) || !one.Col(testdata.Col(dep_index+2),1) || !one.Col(testdata.Col(indep_index+2),2) ||
            !two.Col(testdata.Col(indep_index+2),0) || !two.Col(testdata.Col(dep_index+2),1) ||
            !three.Col(testdata.Col(dep_index),0) || !three.Col(testdata.Col(dep_index+2),1) ||
            !four.Col(testdata.Col(dep_index),0))
           {
            Print(__FUNCTION__, " error ", GetLastError());
            return entropy;
           }
        }
      else
        {

         if(!one.Resize(testdata.Rows(), testdata.Cols()-1) || !two.Resize(testdata.Rows(), testdata.Cols()-2) ||
            !three.Resize(testdata.Rows(), m_tlag+1))
           {
            Print(__FUNCTION__, " error ", GetLastError());
            return entropy;
           }

         matrix deplag = np::sliceMatrixCols(testdata,dep_index?dep_index+m_tlag+1:2,dep_index?END:2+m_tlag);
         matrix indlag = np::sliceMatrixCols(testdata,indep_index?indep_index+m_tlag+1:2,indep_index?END:2+m_tlag);
         //one
         if(!np::matrixCopyCols(one,deplag,1,1+m_tlag) || !np::matrixCopyCols(one,indlag,1+m_tlag) || !one.Col(testdata.Col(dep_index),0))
           {
            Print(__FUNCTION__, " error ", GetLastError());
            return entropy;
           }
         //two
         if(!np::matrixCopyCols(two,indlag,indlag.Cols()) || !np::matrixCopyCols(two,deplag,indlag.Cols()))
           {
            Print(__FUNCTION__, " error ", GetLastError());
            return entropy;
           }
         //three
         if(!np::matrixCopyCols(three,deplag,1) || !three.Col(testdata.Col(dep_index),0))
           {
            Print(__FUNCTION__, " error ", GetLastError());
            return entropy;
           }
         //four
         four = deplag;
        }

      double h1=get_entropy(one,numbins);
      double h2=get_entropy(two,numbins);
      double h3=get_entropy(three,numbins);
      double h4=get_entropy(four,numbins);

      // entropy = independent conditional entropy (h3-h4)  - joint conditional entropy (h1-h2)
      entropy = (h3-h4) - (h1-h2);

      return entropy;

     }

テストの総合的な結果には、get_results()メソッドを使ってアクセスできます。この構造体の各メンバーは、結果の異なる側面を参照し、各ベクトルの長さは、Initialize()メソッドで設定されたパラメータと、実行された移動エントロピー解析のタイプに依存します。

//+------------------------------------------------------------------+
//| Transfer entropy results struct                                  |
//+------------------------------------------------------------------+
struct TEResult
  {
   vector            TE_XY;
   vector            TE_YX;
   vector            p_value_XY;
   vector            p_value_YX;
   vector            z_score_XY;
   vector            z_score_YX;
   vector            Ave_TE_XY;
   vector            Ave_TE_YX;
  };

結果構造体の特性を以下に示します。

構造体プロパティ	詳細
TE_XY	外生変数から内生変数へのエントロピーの 移動
TE_YX	内生変数から外生変数へのエントロピーの 移動
z_score_XY	外生変数から内生変数への移動エントロピーの 重要度
z_score_YX	内生変数から外生変数への移動エントロピーの 重要度
p_value_XY	外生変数から内生変数への移動エントロピーの p値有意性
p_value_YX	外生変数から内生変数への移動エントロピーの p値有意性
Ave_TE_XY	外生変数から内生変数への平均移動 エントロピー
Ave_TE_YX	内生変数から外生変数への平均移動 エントロピー

get_transfer_entropies()を呼び出すと、データセットの最後のウィンドウについて、両方向に測定された推定移動エントロピーのベクトルが返されます。結果の順序は、クラスに渡された元のデータの列の順序に従います。つまり、ベクトルの最初のエントロピー値は、最初の列の系列に対応します。

例

クラスの機能性は、あらかじめ決められた特徴を持つランダムに生成されたシリーズに対してテストを実行することでテストされます。この系列は、generate_time_series.mqhで定義されている以下の関数を使用して生成されます。

//+------------------------------------------------------------------+
//|Generate a random walk time series under Geometric Brownian Motion|
//+------------------------------------------------------------------+
vector random_series(double initial_val, ulong steps, ulong len, double muu, double sgma)
  {
   vector out(len);

   out[0] = initial_val;

   int err=0;

   for(ulong i=1; i<len; i++)
     {
      out[i] = out[i-1]*(1.0+(muu*(double(steps)/double(len)))+(MathRandomNormal(muu,sgma,err)*sqrt(double(steps)/double(len))));
      if(err)
        {
         Print(__FUNCTION__, " MathRandonNormal() ", GetLastError());
         return vector::Zeros(1);
        }
     }

   return out;
  }

関数random_series()は、幾何学的ブラウン運動に特徴的なランダムウォーク時系列を生成します。パラメータは次の通りです。

• initial_val ：時系列の初期値
• steps：ランダムウォークの総ステップ数
• len：生成する時系列の長さ
• muu：GBMのドリフト項（平均）
• sgma ：GBMのボラティリティ（標準偏差）

//+-----------------------------------------------------------------------------------------------+
//|Generate two time series under Geometric Brownian Motion with S2 dependent in part on S1-lagged|
//+-----------------------------------------------------------------------------------------------+
matrix coupled_random_series(double init_1,double init_2,ulong steps, ulong len, double muu_1, double muu_2, double sgma_1, double sgma_2,
                            double alpha, double epsilon, ulong lag)
  {

   vector gbm1 = random_series(init_1,steps,len,muu_1,sgma_1);
   vector gbm2 = random_series(init_2,steps,len,muu_2,sgma_2);

   if(gbm1.Size()!=gbm2.Size())
     {
      return matrix::Zeros(1,1);
     }

   matrix out(gbm2.Size()-lag,2);

   for(ulong i = lag; i<gbm2.Size(); i++)
     {
      gbm2[i]=(1.0-alpha)*(epsilon*gbm2[i-lag] + (1.0-epsilon) * gbm2[i]) + (alpha) * gbm1[i-lag];
      out[i-lag][0] = gbm2[i];
      out[i-lag][1] = gbm1[i];
     }

   return out;
  }

関数coupled_random_series()は、2つの連成ランダムウォーク時系列を生成し、2つ目の系列(gbm2)は1つ目の系列(gbm1)のラグ値に部分的に依存します。この関数は，2列の行列を返し，従属系列は最初の列に配置されます。関数のパラメータは以下の通りです。

• init_1 ：最初の時系列の初期値
• init_2 ：2番目の時系列の初期値
• steps：ランダムウォークの総ステップ数
• len：生成する時系列の長さ
• muu_1：第1系列のドリフト項
• muu_2：第2系列のドリフト項
• sgma_1：第1シリーズのボラティリティ
• sgma_2：第2シリーズのボラティリティ
• alpha：独立系列が従属系列に与える影響を表す混合パラメータ
• epsilon：遅行従属系列値の影響を調整するパラメータ
• lag：従属系列の独立系列に対する依存性を表すラグ

CTransEntropyクラスの機能を実証するために、2つのMetaTrader 5スクリプトが準備されました。両方のスクリプトは、データセットを分析し、従属変数（時系列）で観察される従属性を最もよく特徴付ける独立変数（時系列）のラグを検出するために、このクラスがどのように使用できるかを示しています。第一の方法は、異なるラグでの移動エントロピーの分析によって得られた一連の結果から、最も重要な方向性エントロピーの値を決定するために、目視検査に頼るものです。この方法はスクリプトLagDetection.ex5で実装されています。

//+------------------------------------------------------------------+
//|                                                 LagDetection.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<transfer_entropy.mqh>
#include<generate_time_series.mqh>
//--- input parameters

input double   Init1=100.0;
input double   Init2=90.0;
input ulong    Steps=1;
input ulong    Len=500;
input double   Avg1=0;
input double   Avg2=0;
input double   Sigma1=1;
input double   Sigma2=1;
input double   Alph=0.5;
input double   Epsilon=0.3;
input ulong    Lag=3;
input bool     UseSeed = true;
input ulong    Bins=3;
input ENUM_TE_TYPE testtype=NONLINEAR_TE;
input ulong    NumLagsToTest = 10;
input int      PlotViewTimeInSecs = 20;
//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
   if(UseSeed)
    {
     MathSrand(256);
    }
//---
   if(!NumLagsToTest)
     {
      Print(" Invalid input parameter value for \'NumLagsToTest\'. It must be > 0 ");
      return;
     }

   matrix series = coupled_random_series(Init1,Init2,Steps,Len,Avg1,Avg2,Sigma1,Sigma2,Alph,Epsilon,Lag);

   series = log(series);

   series = np::diff(series,1,false);
   
   matrix entropies(NumLagsToTest,2);

   for(ulong k = 0; k<NumLagsToTest; k++)
     {
      CTransEntropy ote;

      if(!ote.Initialize(series,0,1,k+1))
         return;

      if((testtype==NONLINEAR_TE && !ote.Calculate_NonLinear_TE(Bins)) ||
         (testtype==LINEAR_TE && !ote.Calculate_Linear_TE()))
         return;

      vector res = ote.get_transfer_entropies();

      entropies.Row(res,k);
     }

   Print(" entropies ", entropies);

   CGraphic* g = np::plotMatrix(entropies,"Transfer Entropies","Col 0,Col 1","Lag","TE");

   if(g==NULL)
      return;
   else
     {
      Sleep(int(MathAbs(PlotViewTimeInSecs))*1000);
      g.Destroy();
      delete g;
     }

   return;
  }
//+------------------------------------------------------------------+

スクリプトの最初の11個のユーザーがアクセス可能な入力パラメータは、生成された系列のプロパティを制御します。最後の4つの入力パラメータは、分析の様々な側面を設定します。

• Bins：データの確率密度の推定に使用するヒストグラム法のビンの数を設定します。
• testtype：線形または非線形の移動エントロピー分析を選択できるようにします。
• NumLagsToTest：テストを実施するラグの最大数を1から順に設定します。
• PlotViewTimeInSecs：プログラムが終了するまでのグラフの表示時間を指定します。
• UseSeed：trueの場合、テスト結果の再現性を保証するための乱数発生器のシードを有効にします。
  

このスクリプトは、所定の依存関係を持つ2つの時系列を生成し、異なるラグで移動エントロピーを推定します。データは分析前に差分されています。おそらくこの場合は不要でしょうが、良い練習になります。その結果（移動エントロピー）はグラフに可視化され、横軸の対応するラグに対して縦軸に移動エントロピーがプロットされます。テストに成功すると、ランダム系列を生成するために選択されたラグに明確なピークを持つプロットが生成されるはずです。

プログラムを実行すると、線形検定が系列の生成に使用されたラグ依存性を首尾よく特定したことがわかります。従属系列は、ランダムに生成されたデータセットの最初の列にあることを思い出してください。

非線形テストオプションを使用して再度テストをおこなうと、同様の結果が得られます。この例では、エントロピー値の大きさは著しく小さくなります。これは、データの確率分布を推定するためにヒストグラム法を使用することの限界によるものと考えられます。また、選択されたビンの数は、推定された移動エントロピーに影響することに注意すべきです。 

次のデモンストレーションでは、特定のラグで得られたエントロピー値の有意性を検証します。これはスクリプトLagDetectionUsingSignificance.ex5で実装されています。

//+------------------------------------------------------------------+
//|                                LagDetectionUsingSignificance.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<transfer_entropy.mqh>
#include<generate_time_series.mqh>
//--- input parameters
input double   Init1=100.0;
input double   Init2=90.0;
input ulong    Steps=1;
input ulong    Len=500;
input double   Avg1=0;
input double   Avg2=0;
input double   Sigma1=1;
input double   Sigma2=1;
input double   Alph=0.5;
input double   Epsilon=0.3;
input ulong    Lag=3;
input bool     UseSeed = true;
input ulong    Bins=3;
input ENUM_TE_TYPE testtype=LINEAR_TE;
input ulong    LagToTest = 3;
input ulong    NumIterations = 100;
//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
   if(UseSeed)
    {
     MathSrand(256);
    }
//---
   if(!LagToTest)
     {
      Print(" Invalid input parameter value for \'LagToTest\'. It must be > 0 ");
      return;
     }

   matrix series = coupled_random_series(Init1,Init2,Steps,Len,Avg1,Avg2,Sigma1,Sigma2,Alph,Epsilon,Lag);

   series = log(series);

   series = np::diff(series,1,false);

   matrix entropies(1,2);


   CTransEntropy ote;

   if(!ote.Initialize(series,0,1,LagToTest))
      return;

   if((testtype==NONLINEAR_TE && !ote.Calculate_NonLinear_TE(Bins,NumIterations)) ||
      (testtype==LINEAR_TE && !ote.Calculate_Linear_TE(NumIterations)))
      return;

   vector res = ote.get_transfer_entropies();

   entropies.Row(res,0);

   TEResult alres = ote.get_results();

   Print(" significance: ", " pvalue 1->0 ",alres.p_value_XY, " pvalue 0->1 ",alres.p_value_YX);
   Print(" zscore 1->0 ",alres.z_score_XY, " zscore 0->1 ",alres.z_score_YX);

   return;
  }
//+------------------------------------------------------------------+

このスクリプトには、最後の2つを除いて、ユーザーが調整可能な同様の入力パラメータがあります。

• LagToTest：テストが実施される特定のラグを設定します。
• NumIterations：有意性検定のためにデータをシャッフルする回数を定義します。

このスクリプトは、従属系列のペアを生成し、選択されたラグで検定を実行します。移動エントロピーは、対応するp値とz得点とともに、端末のExpertsタブに書き込まれます。

最初の実行では、LagToTestとLagパラメータを同じ値に設定してスクリプトを実行します。結果は以下の通りです。彼らは、1列目の系列が行列の2列目の系列に依存していることを示しています。

JS      0       21:33:43.464    LagDetectionUsingSignificance (Crash 1000 Index,M10)     significance:  pvalue 1->0 [0] pvalue 0->1 [0.66]
LE      0       21:33:43.464    LagDetectionUsingSignificance (Crash 1000 Index,M10)     zscore 1->0 [638.8518379295961] zscore 0->1 [-0.5746565128024472]

2回目の実行では、LagToTestパラメータの値のみを変更し、その結果を前回の実行結果と比較します。

p値とz得点の違いに注目してください。この場合、p値もz得点も重要ではありません。

RQ      0       21:33:55.147    LagDetectionUsingSignificance (Crash 1000 Index,M10)     significance:  pvalue 1->0 [0.37] pvalue 0->1 [0.85]
GS      0       21:33:55.147    LagDetectionUsingSignificance (Crash 1000 Index,M10)     zscore 1->0 [-0.2224969673139822] zscore 0->1 [-0.6582062358345131]

テストの結果は、CTransEntropyクラスが良好に動作することを示していますが、より大きなラグで分析をおこなう場合、特に複数のラグ項のオプションが有効な場合（maxLagOnlyがfalse）には、大きな制限があります。これは特に非線形テストで問題となります。これは、データの分布を推定するためにヒストグラム法を使用することに起因します。確率密度の推定にヒストグラム法を使うことには、重大な欠点があります。ビン幅（またはビン番号）の選択は、ヒストグラムの外観と精度に大きく影響します。ビン幅が小さすぎると、ノイズが多く断片的なヒストグラムになり、ビン幅が大きすぎると、重要な詳細が不明瞭になったり、特徴が平滑化されたりします。最大の問題は、ヒストグラムが主に一次元データに有効であるという事実に関連しています。高次元のデータでは、ビンの数は指数関数的に増加します。考慮すべきラグが多数ある場合、利用可能な計算リソースへの要求は著しく伸びる可能性があります。したがって、一般化移動エントロピーを用いて複数のラグを考慮した分析をおこなう場合は、ラグの最大数を小さくしておくことが推奨されます。 

結論

総括すると、CTransEntropyクラスは、線形および非線形の両方のコンテキストにおいて移動エントロピーの分析を可能にします。実用的なデモンストレーションにより、ある時系列が別の時系列に及ぼす影響を検出・定量化する能力が示され、その結果は目視検査および有意差検定によって検証されました。このクラスは様々なシナリオに柔軟に対応し、時系列アプリケーションにおける因果関係の解明に貴重な洞察を提供します。ただし、特に非線形手法を使用する際には、複数のラグを分析することに伴う計算負荷に注意が必要です。効率的なパフォーマンスと正確な結果を得るためには、考慮するラグ数を適切に制限することが推奨されます。全体として、CTransEntropyクラスは複雑な依存関係を明らかにし、動的システムの理解を深めるための強力なツールです。

ファイル	詳細
Mql5\include\generate_time_series.mqh	ランダムな時系列を生成する関数を含む
Mql5\include\ np.mqh	ベクトルと行列の効用関数のコレクション
Mql5\include\ OLS.mqh	通常の最小二乗回帰を実装するOLSクラスの定義を含む
Mql5\include\ TestUtilities.mqh	OLS評価のためのデータセットを準備するために使用されるツール群を提供する
Mql5\include\ transfer_entropy	移動エントロピー解析を実装するCTransEntropyクラスの定義を含む
Mql5\scripts\LagDetection.mq5	CTransEntropyクラスの機能を示すスクリプト
Mql5\scripts\LagDetectionUsingSignificance.mq5	CTransEntropyを使用した結果を解釈するための異なるアプローチを示す2つ目のスクリプト