Причинно-следственный анализ временных рядов с помощью энтропии переноса

Введение

Энтропия переноса — это статистический инструмент, который дает количественную оценку объему информации, передаваемой из одного временного ряда в другой, обеспечивая понимание характера и поведения целевой переменной. В этой статье мы углубимся в концепцию статистических причинно-следственных связей, рассчитываемой с помощью энтропии переноса. Мы исследуем, как этот метод может выявить направление причинно-следственного взаимовлияния различных процессов. Кроме того, мы приводим подробное описание реализации MQL5 для измерения энтропии переноса, демонстрируя, как этот метод можно применять на практике для анализа потенциально взаимосвязанных временных рядов. Применяя энтропию переноса, постараемся выявить переменные которые могут эффективно дополнить задачи прогнозирования.

Причинно-следственные связи

Эмпирические данные могут быть обманчивы. Тот факт, что две переменные, по-видимому, взаимосвязаны, не означает, что одна из них вызывает изменения в другой, поэтому высказывание «корреляция не есть каузация» звучит справедливо. Корреляция просто измеряет степень связи двух переменных, но не причину такой связи. Например, представьте себе сильную корреляцию между продажами мороженого и стоимостью акций в течение лета. Это не означает, что покупка мороженого приводит к росту акций! Более вероятной причиной является некий скрытый фактор, например само время года, независимо влияющий на обе переменные. Аналогичным образом может существовать связь между ценами на акции компании и ценами на золото, но реальная причина может заключаться в чем-то совсем другом, например в общих настроениях на рынке или в инфляции, влияющей на обе цены. Следующие примеры показывают, что коррелирующие данные могут вводить в заблуждение. Они отмечают связь, но не раскрывают ее причину. Чтобы по-настоящему понять, является ли одно явление причиной другого, нужны более продвинутые инструменты.

Понятие причинности, представление о том, что одно событие приводит к возникновению другого, является фундаментальным в научном исследовании. Однако точное определение причинности представляет собой многогранную задачу, требующую сложных философских, физических и статистических рассуждений. В идеале причина должна неизменно приводить к единственному следствию. Но может быть нелегко выделить один причинный фактор из зачастую сложного переплетения воздействий, влияющих на результат. Например, резкий рост объема торгов может коррелировать с ростом цен акций, но важную роль могут играть и другие факторы, такие как настроения рынка и публикации экономических данных. В таких ситуациях для установления причинно-следственных связей исследователи прибегают к статистическим методам.

Характер как детерминированной (с гарантированным результатом), так и вероятностной (влияющей на вероятность наступления последующего события) причинности зависит от лежащего в ее основе процесса. В детерминированных системах первое событие очевидным образом приводит ко второму, как это наблюдается в случае предсказуемого падения предмета. Напротив, в вероятностных системах акцент смещается на вопрос о том, повышает ли первое событие нашу способность предсказывать наступление второго. Например, хотя недавно прошедший дождь может быть связан с последующим распусканием цветов, этому могут способствовать и другие факторы. В таких случаях возникает вопрос, улучшает ли знание о первом событии нашу способность предсказать второе событие.

Экономист Клайв Грейнджер, опираясь на работу Норберта Винера, разработал концепцию причинно-следственных связей, постулируя, что первый сигнал вызывает второй сигнал, если будущие значения второго сигнала можно лучше объяснить с помощью прошлой информации как из первого, так и из второго сигнала, а не просто с использованием запаздывающих значений только второго сигнала. В своем определении причинности Грейнджер основывался на двух принципах. Во-первых, следствие не может проявиться раньше его причины. Во-вторых, в причине будет содержаться уникальная информация, переданная следствию. Эти принципы предполагают, что для количественной оценки причинности нужно понимать временные характеристики задействованных переменных, а также некоторую меру их информационного содержания. Благодаря этому для анализа причинно-следственных связей хорошо подходят временные ряды.

Опираясь на характер данных временных рядов, мы можем проанализировать, как информация из одного ряда в определенный момент влияет на предсказуемость другого ряда в более поздний момент. Грейнджер определяет причинность как снижение неопределенности. Если знание значений серии X улучшает прогноз будущего значения серии Y по сравнению с простым использованием прошлых значений самой серии Y, то считается, что серия X прогнозирует серию Y. Основываясь на этой идее, Грейнджер разработал тесты для определения причинно-следственных связей с помощью запаздывающих временных рядов и авторегрессионных моделей. Он предположил, что между X и Y нет причинно-следственной связи, пока включение прошлых значений X не улучшит существенно прогноз будущих значений Y. Такое улучшение обычно измеряют по уменьшению ошибки прогнозирования, например, по лучшему соответствию регрессионной модели. Причинность по Грейнджеру позволяет статистически выявлять взаимосвязи между временными рядами. Однако важно отметить ее ограничения. Причинность по Грейнджеру определяет только направленную связь, а не обязательно окончательный механизм причин и следствий. Помимо тех причин, о которых у нас имеются данные, могут действовать и другие причины. Кроме того, поскольку причинность по Грейнджеру по сути основана на модели авторегрессии, наиболее эффективна она для выявления линейных причинно-следственных связей. Нелинейные причинно-следственные связи требуют иного подхода.

С математической точки зрения, причинность по Грейнджеру можно выразить, рассмотрев два временных ряда, X и Y. Запаздывающие значения от каждого из них обозначаются как X(t−k) и Y(t−k), представляя собой отставание в точке k. Максимальное рассматриваемое отставание обозначается как p. Применяя модель авторегрессии, будущее значение Y регрессируется по своим собственным прошлым значениям.

Это выражение рассматривает будущее значение Y только с точки зрения его прошлых значений. При введении ряда X в модель будущие значения выражаются через прошлые значения из обоих рядов, X и Y.

Если включение прошлых значений X значительно улучшает прогноз ряда Y по сравнению с моделью, использующей только прошлые значения Y, то считается, что X является причиной ряда Y по Грейнджеру. Обычно это оценивают путем проверки нулевой гипотезы о том, что коэффициенты совместно равны нулю. Если опровергнуть эту нулевую гипотезу, это означает, что ряд X предоставляет значимую прогностическую информацию о ряде Y, выходящую за рамки содержимого только прошлых значений Y. Чтобы проверить, является ли ряд X причиной ряда Y по Грейнджеру, сравним две модели в рамках проверки гипотезы:

• Нулевая гипотеза: X не является причиной Y по Грейнджеру.

• Альтернативная гипотеза: X является причиной Y по Грейнджеру.

F-тест используется для сравнения соответствия ограниченной модели (без X) и неограниченной модели (с X) путем изучения остатков соответствующих моделей. Ограниченная сумма квадратов остатков — это остатки от модели без X, а неограниченная сумма квадратов остатков получается из модели с X. F-статистика рассчитывается так:

Где n — количество наблюдений. Рассчитанная F-статистика сравнивается с критическим значением из F-распределения со степенями свободы p и n − 2p − 1. Если F-статистика превышает критическое значение, мы отказываемся от нулевой гипотезы и заключаем, что X является причиной по Грейнджеру для Y. Альтернативным способом расчета F-статистики может быть тест на основе однофакторного дисперсионного анализа (ANOVA). Его формула приведена ниже.

Энтропия переноса

На заре теории информации ученые использовали взаимную информацию для понимания принципов взаимодействия связанных процессов. Эта концепция, основанная на энтропии Клода Шеннона, информирует, накладывается ли информация из одного временного ряда на данные другого. Проще говоря, она показывает, возможно ли закодировать оба ряда вместе, используя меньше информации, чем при кодировании каждого в отдельности. По этой причине взаимную информацию иногда называют избыточностью. Один процесс обменивается информацией с другим, что позволяет эффективно описывать второй процесс путем повторного использования уже полученной из первого процесса информации.

Формально взаимная информация между двумя стохастическими процессами, X(t) и Y(t), устанавливается при превышении суммой их предельных энтропий совместную энтропию объединенной системы. Такая математическая зависимость отражает снижение неопределенности в объединенной системе по сравнению с отдельными процессами. Иными словами, она отражает степень, в которой информация об одном процессе может быть использована для уменьшения внутренней энтропии, связанной с другим процессом. Поскольку энтропия определяется исключительно распределением лежащих в ее основе вероятностей, любое такое распределение можно охарактеризовать соответствующим значением энтропии. Это значение количественно определяет уровень неожиданности, связанной с конкретным результатом, учитывая известное распределение вероятностей.

Эта концепция становится особенно актуальной в контексте причинности по Грейнджеру. При исследовании потенциальных причинно-следственных связей между временными рядами целью является уменьшение неопределенности, связанной с целевым процессом, путем включения информации из потенциального исходного процесса. Если включение вторичного временного ряда доказуемо уменьшает энтропию распределения целевого процесса, это предполагает наличие статистически значимого причинного воздействия со стороны исходного ряда на целевой. Такое уменьшение называют энтропией переноса.

Энтропия переноса (ЭП) основана на концепции расстояния Кульбака-Лейблера для измерения направления переноса информации между двумя временными рядами. В частности, ЭП основана на идее условной взаимной информации и может быть выражена с помощью расстояния Кульбака-Лейблера (РКЛ), известного также как расхождение Кульбака-Лейблера или относительная энтропия. РКЛ измеряет разницу между двумя распределениями вероятностей. В ЭП это расстояние измеряет разницу между совместным распределением вероятностей текущего состояния Y и прошлых состояний X и Y и произведением предельных распределений этих состояний. Математически перенос информации из временного ряда X во временной ряд Y можно выразить следующим образом:

где y(t+1)​ — будущее состояние Y, y(t) ​— прошлое состояние Y, а x(t) — прошлое состояние X. Такая формулировка подчеркивает, что энтропия переноса измеряет, насколько изменится распределение вероятностей y(t+1) при рассмотрении информации из x(t) в дополнение к y(t)​.

В 2009 году Лайонел Барнетт (Lionel Barnett), Адам Барретт (Adam Barrett) и Анил Сет (Anil Seth) опубликовали статью «Причинность по Грейнджеру и энтропия переноса эквивалентны для гауссовых переменных» ("Granger Causality and Transfer Entropy Are Equivalent for Gaussian Variables"), демонстрируя, что при следовании временных рядов гауссовскому распределению энтропия переноса эквивалентна половине F-статистики для причинности по Грейнджеру.

Этот результат дает определение линейной энтропии переноса, которую мы реализуем в коде позднее. Для учета нелинейной причинности расширяем концепцию снижения неопределенности, следуя работе Томаса Шрайбера (Thomas Schreiber), в которой временной ряд рассматривают как марковский процесс с различными распределениями вероятностей перехода.

Подход Шрайбера к моделированию снижения неопределенности пользуется теорией информации, рассматривая временные ряды X(t) и Y(t) как марковские процессы с известными вероятностями перехода p(x) и q(x). В отличие от авторегрессионной модели, опирающейся на линейные модели, этот подход применяет для описания передачи данных условную взаимную информацию. Поскольку взаимная информация выводится из разницы энтропий, условную энтропию получают путем добавления к каждому члену энтропии дополнительной информации в качестве условия. Затем вычисляют энтропию переноса подстановкой запаздывающих переменных в уравнение условной взаимной информации, что позволяет анализировать перенос данных из X(t) в Y(t) при определенном запаздывании k с помощью условной взаимной энтропии.

С вычислительной точки зрения этот метод привлекателен, поскольку для совместной энтропии требуется только распределение вероятностей. Энтропию переноса для отдельного запаздывания k можно выразить как четыре отдельных члена совместной энтропии, которые легко вычислить с помощью точного распределения вероятностей на основе данных. Преимуществом этой формулы является ее способность обрабатывать больше запаздывающих измерений. Однако каждое дополнительное отставание увеличивает размерность пространства состояния на два, что оказывает значительное влияние на возможность точного количественного определения энтропии переноса из-за экспоненциального роста проблем с конечными данными, связанных с оценкой плотности вероятностей.

Главное преимущество этого подхода состоит в его непараметрическом характере. В отличие от других методов, он не делает каких-либо предположений о базовом распределении данных вне стационарности, что позволяет применять его, не имея предварительных знаний о процессах генерации данных. Однако у этого преимущества есть другая сторона: результаты в значительной степени зависят от точной оценки базового распределения. Расчет энтропии переноса требует аппроксимации инстинного распределения вероятностей соответствующих стохастических процессов с помощью ограниченных данных для вычисления четырех членов энтропии. Точность такой оценки в значительной мере влияет на надежность результатов переноса энтропии. Имея это в виду, необходимо учитывать возможность получения при расчетах ложных значений энтропии. Поэтому нам нужен какой-то способ определения надежности результатов.

Наше стремление применить непараметрический подход к оценке энтропии переноса сопряжено с серьезной проблемой: необходимо гарантировать, что результаты отразят истину, а не мусор. Поэтому нам стоит рассмотреть более информативный подход к интерпретации энтропии переноса, который подразумевает оценку статистической значимости рассчитанного значения. Обычные тесты на значимость включают в себя перемешивание данных временного ряда заданное количество раз. Затем для каждой перемешанной версии рассчитывается энтропия переноса. Вслед за этим p-значение рассчитывается как доля перемешанных данных с энтропией переноса ниже исходного значения.

Другой подход требует расчета количества стандартных отклонений результата от среднего значения перемешанных данных. Поскольку перемешивание нарушает временную структуру, ожидается, что среднее значение при перемешивании энтропии переноса будет близко к нулю. Разброс данных вокруг этого среднего значения отражает значимость первоначального результата. Рассчитанное значение называют z-оценкой. Как правило, z-оценки требуют меньшего количества перемешиваний по сравнению с p-значениями, что делает их эффективнее с точки зрения вычислений. 

В случае p-значения цель — получить вероятность, максимально близкую к нулю. В то время как указывающая на статистическую значимость z-оценка должна превышать 3.0.

Реализация MQL5

Код, реализующий инструменты для количественной оценки энтропии переноса и определения ее значимости, включен в transfer_entropy.mqh. Файл содержит определение класса CTransEntropy, а также других вспомогательных классов и функций. Этот класс предлагает структуру для статистического анализа данных временных рядов, специально ориентированную на оценку причинно-следственных связей между переменными. Он представляет два различных метода количественной оценки линейной причинности по Грейнджеру (линейной энтропии переноса) и нелинейной энтропии переноса. Рассчитывается в обоих направлениях, предоставляя более полную картину потока информации между переменными.

Для устранения потенциальной нестационарности данных в этот класс включена оконная процедура. Пользователи могут определять размер и шаг окон, что позволяет анализировать данные в меньших, перекрывающихся сегментах. Такой подход дает специфичные для каждого окна результаты, что упрощает выявление временных изменений стойкости причинно-следственных связей. Кроме того, он сглаживает проблемы, связанные с анализом нестационарных данных. Класс также предоставляет интегрированный механизм проверки значимости. Пользователи могут указать количество перемешиваний данных, которые необходимо выполнить, сохраняя при этом предельные распределения. На основе таких перемешанных наборов данных класс вычисляет p-значения и z-оценки для энтропии переноса в каждом направлении. Эти статистические значения позволяют понять, какова вероятность обусловленности наблюдаемых причинно-следственных связей или переноса информации случайностью, что повышает надежность анализа.

Экземпляр класса создается с помощью конструктора по умолчанию без параметров.

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

     }

Затем пользователи должны вызвать метод Initialize(), который инициализирует объект с определенным набором данных и настраивает различные параметры для анализа.

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;
     }

Первый требуемый параметр — матрица, в которой не менее двух столбцов; анализируемый временной ряд нужно поместить в столбцы входной матрицы. При работе с нестационарными данными рекомендуется предварительно провести их дифференциацию. Второй и третий параметры — индексы столбцов матрицы входных данных, указывающие на эндогенный (зависимый) и экзогенный (независимый) временные ряды соответственно.

Четвертый параметр, запаздывание, определяет параметр lag, который рассматривается в анализе. Следующий логический параметр, maxLagOnly, определяет, задает lag один элемент (если true) или все запаздывающие значения вплоть до lag (если false). Предпоследний параметр, winsize, обозначает длину окна. Если он установлен на 0, к данным не будет применяться разбиение на окна. Наконец, параметр winstride опционально устанавливает размерность окна для операций с окнами, определяя шаг между последовательными окнами при их прохождении по данным временного ряда.

Метод начинается с необходимости убедиться, что эндо- и экзогенные индексы не совпадают. Если же они совпадают, он выводит сообщение об ошибке и возвращает значение false. Размер внутренней матрицы m_dataset изменяется для хранения анализируемого двумерного набора данных. Затем он копирует столбцы, указанные endog_index и exog_index, из входной матрицы в первый и второй столбцы m_dataset соответственно. Если это необходимо, для оконного отображения матрицы m_dataset используется вспомогательный класс CDataWindows. После этого метод устанавливает внутренние переменные с предоставленными параметрами для дальнейшего использования.

//+------------------------------------------------------------------+
//|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. Если n_shuffles равен нулю (по умолчанию), тесты на значимость не проводятся.

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;
     }

Метод вычисляет линейную энтропию переноса с помощью метода Грейнджера. Это реализовано в частном методе — linear_transfer(). Два последних параметра этой процедуры определяют зависимую и независимую переменную (столбец) во входной матрице. Дважды вызвав метод при переключенных индексах столбцов, можно получить энтропию переноса в обоих направлениях.

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() принимает параметр numBins в дополнение к параметру n_shuffles. Этот параметр определяет количество интервалов, применяемых при оценке плотности вероятностей переменных. 

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;


     }

Для оценки плотности вероятностей используют метод гистограмм. Он был выбран, потому что его проще всего реализовать. Ответственность за вычисление обобщенной версии энтропии переноса передана частным методам 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();

     }

Четыре значения компонентов, используемые для расчета совместной и независимой условной энтропий, объединены в 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() генерирует два связанных временных ряда случайного блуждания, где второй ряд (gbm2) частично зависит от запаздывающих значений первого ряда (gbm1). Функция возвращает матрицу с двумя столбцами, в первом столбце которой размещается зависимый ряд. Функция имеет следующие параметры:

• init_1 : Начальное значение первого временного ряда.
• init_2 : Начальное значение второго временного ряда.
• steps : Общее количество шагов в случайном блуждании.
• len : Длина генерируемого временного ряда.
• muu_1 : Время ухода первого ряда.
• muu_2 : Время ухода второго ряда.
• sgma_1 : Волатильность первого ряда.
• sgma_2 : Волатильность второго ряда.
• alpha : Параметр смешивания для влияния независимого ряда на зависимый ряд.
• epsilon : Параметр, регулирующий влияние запаздывающих значений зависимого ряда.
• lag : Задержка для подчиненного отношения зависимого ряда к независимому ряду.

Для демонстрации возможностей класса CTransEntropy были приготовлены два скрипта для 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, включает начальное значение для генератора случайных чисел, чтобы обеспечить воспроизводимость результатов теста.
  

Скрипт генерирует два временных ряда с предопределенной зависимостью и оценивает энтропию переноса при запаздываниях разной продолжительности. Обратите внимание: перед анализом данные были дифференцированы. Возможно, в данном случае в этом нет необходимости, но это обычная практика. Затем результаты (полученные энтропии переноса) визуализируются на графике, где энтропия переноса откладывается по вертикальной оси, а соответствующее запаздывание — по горизонтальной. Если тест прошел успешно, должен получиться график с четким пиком на запаздывании, выбранном для генерации случайного ряда.

Запуск программы показывает, что линейный тест успешно выявил зависимость запаздывания, используемую для генеррации ряда. Напомним, что зависимый ряд находится в первом столбце случайно сгенерированного набора данных.

Повторная проверка с использованием возможности нелинейного тестирования дает сходные результаты. В этом случае величина значения энтропии заметно меньше. Это может быть связано с ограничениями использования метода гистограмм для оценки распределения вероятностей данных. Следует также отметить, что количество выбранных интервалов повлияет на расчетную энтропию переноса.  

В следующем примере протестируем значимость показателей энтропии, полученных при определенных запаздываниях. Это реализовано в скрипте 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;
  }
//+------------------------------------------------------------------+

Скрипт имеет аналогиченые настраиваемые пользователем входные параметры, кроме двух последних:

• LagToTest : Устанавливает определенное запаздывание, с которым будет проводиться тест.
• NumIterations : Определяет количество повторов перемешивания данных для проверки значимости.

Скрипт генерирует пару зависимых рядов и выполняет тест с выбранным запаздыванием. Энтропия переноса вместе с соответствующими p-значением и z-оценкой записываются во вкладку «Эксперты» терминала.

При первом запуске скрипт выполняется с параметрами LagToTest и Lag, для которых задано одно и то же значение. Результаты показаны ниже. Они демонстрируют, что ряд в первом столбце матрицы зависит от ряда во втором столбце.

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]

При втором запуске модифицируем только значение параметра 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