Двумерные копулы в MQL5 (Часть 1): Реализация гауссовой копулы и t-копулы Стьюдента для моделирования зависимостей

Введение

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

С учётом этого данная статья открывает серию материалов, посвящённых реализации инструментов для торговых стратегий на основе копул. В первой части мы рассмотрим основы статистических копул в контексте парного трейдинга и моделирования зависимостей. Также будет представлен код на MQL5 для предварительной обработки наборов данных перед подгонкой моделей копул. В дальнейшем серия будет сосредоточена на библиотеке, реализующей часто используемые модели копул. Для начала в этой статье представлены реализации гауссовой копулы и t-копулы Стьюдента.

Ключевые основы теории копул

Для работы с копулами необходимо понимать несколько связанных концепций.

Функция плотности вероятности. Функция плотности вероятности (PDF) характеризует плотность распределения непрерывной случайной величины X в точке x; при этом она не задаёт вероятность для отдельного значения напрямую. В отличие от дискретной функции вероятностей, PDF не задаёт вероятности непосредственно для отдельных значений. Вместо этого она представляет плотность, так что вероятность попадания X в конкретный интервал [a,b] определяется интегралом PDF по этому диапазону.

Функция распределения. На основе PDF строится функция распределения (CDF), которая задаёт вероятность того, что случайная величина X примет значение меньшее или равное x. CDF является неубывающей функцией в диапазоне от 0 до 1 и полностью описывает вероятностное распределение случайной величины.

Маргинальные и совместные распределения. При работе с несколькими случайными величинами необходимо различать маргинальные и совместные распределения. Маргинальное распределение описывает вероятностное распределение одной переменной в составе большего набора, фактически игнорируя влияние остальных переменных. Совместное распределение, напротив, характеризует общее вероятностное распределение всех переменных в наборе данных.  

Преобразование вероятностного интеграла (PIT). PIT — фундаментальная концепция теории вероятностей и важнейший строительный блок для понимания копул. Она утверждает, что если X — непрерывная случайная величина с CDF = F​(X), то преобразованная переменная U=F​(X) имеет равномерное распределение на интервале [0,1]. Иными словами, независимо от того, имеет ли X нормальное, экспоненциальное или любое другое непрерывное распределение, применение её собственной CDF преобразует её в переменную с равномерным распределением на [0,1].

PIT позволяет отобразить любую непрерывную случайную величину в единичный интервал, фактически стандартизируя её при сохранении вероятностной структуры. Преобразовав маргинальные распределения к равномерному виду, можно анализировать связь между переменными без влияния их индивидуальных распределений. Именно в этом равномерном пространстве и работают копулы: они описывают структуру зависимости между переменными после такого преобразования маргинальных распределений.

Эмпирическая функция распределения. Когда истинное базовое распределение переменной неизвестно, можно использовать эмпирическую функцию распределения (ECDF) как оценку истинной CDF. Для выборки из N наблюдений x1​,x2​,…,xN​ ECDF определяется как:

где задаётся функция, равная 1, если x_i <= x, и 0 в противном случае.

Использование ECDF для выполнения PIT часто называют непараметрическим подходом к оцениванию, поскольку он не предполагает конкретной формы распределения маргиналей. Преобразованные точки данных используются как входные данные для подгонки самой копулы.

Файл linear_cdf.mqh содержит определение класса CLinearCDF, используемого для оценки ECDF набора данных. Класс строит ECDF с линейной интерполяцией между ступенями стандартной ECDF. Это эффективно сглаживает результат, сохраняя совпадение со стандартной ECDF в точках выборки.

//+------------------------------------------------------------------+
//|                                                   linear_cdf.mqh |
//|                                  Copyright 2024, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2024, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#include<Arrays/ArrayDouble.mqh>
#include<Arrays/ArrayObj.mqh>
#include "ecdf.mqh"
//+------------------------------------------------------------------+
//|CArray objects                                                    |
//+------------------------------------------------------------------+
class CArrayDoubles:public CArrayObj
  {
   virtual bool      CreateElement(const int index)
     {
      m_data[index] = new CArrayDouble();
      return m_data[index]!=NULL;
     }
  };
//+------------------------------------------------------------------+
//| ecdf quantile transformer                                        |
//+------------------------------------------------------------------+
class CLinearCDF
  {
protected:
   CSpline1DInterpolantShell* m_shells[];
   CArrayDoubles      m_internal;
   CArrayDouble      *m_arrays[];
   double            m_ub,m_lb;
   ulong             m_cols;
   void              clean_up(void)
     {

      for(uint i = 0; i<m_shells.Size(); i++)
         delete m_shells[i];
      for(uint i = 0; i<m_arrays.Size(); i++)
         delete m_arrays[i];
      m_internal.Clear();
     }
public:
                     CLinearCDF(void)
     {
      m_cols = 0;
     }
                    ~CLinearCDF(void)
     {
      clean_up();
     }
   bool              save(const int filehandle)
     {
      if(!m_cols)
         return false;
      return m_internal.Save(filehandle);
     }
   bool              load(const int filehandle)
     {
      if(m_cols)
        {
         clean_up();
         m_cols = 0;
        }

      if(m_internal.Load(filehandle))
        {
         int n_arrays = m_internal.Total();
         m_cols = ulong((n_arrays-1)/2);
         ArrayResize(m_shells,int(m_cols));
         ArrayResize(m_arrays,int(m_cols*2)+1);
         double xa[],ya[];
         m_arrays[0] = m_internal.At(0);
         m_ub = m_arrays[0].At(0);
         m_lb = m_arrays[0].At(1);
         for(int i = 0; i<int(m_cols); i++)
           {
            m_arrays[i*2+1] = m_internal.At(i*2+1);
            m_arrays[i*2+2] = m_internal.At(i*2+2);
            ArrayResize(xa,m_arrays[i*2+1].Total());
            for(int j = 0; j<m_arrays[i*2+1].Total(); j++)
               xa[j] = m_arrays[i*2+1][j];
            ArrayResize(ya,m_arrays[i*2+2].Total());
            for(int j = 0; j<m_arrays[i*2+2].Total(); j++)
               ya[j] = m_arrays[i*2+2][j];
            m_shells[i] = new CSpline1DInterpolantShell();
            CAlglib::Spline1DBuildLinear(xa,ya,m_shells[i]);
           }
         return true;
        }
      else
         return false;
     }
   bool              fit(matrix &data,double upper_bound = 1.0-1.e-5,double lower_bound = 1.e-5)
     {
      if(upper_bound<=lower_bound)
        {
         Print(__FUNCTION__, " invalid inputs ");
         return false;
        }

      if(m_cols)
        {
         clean_up();
         m_cols = 0;
        }

      m_ub = upper_bound;
      m_lb = lower_bound;
      m_cols = data.Cols();
      vector v;
      double xa[],ya[];
      vector slopes;
      vector slope_samples;
      CECDF m_ecdf;
      ArrayResize(m_shells,int(m_cols));
      ArrayResize(m_arrays,int(m_cols*2)+1);
      m_arrays[0] = new CArrayDouble();
      if(!m_arrays[0].Add(m_ub) || !m_arrays[0].Add(m_lb) || !m_internal.Add(m_arrays[0]))
        {
         Print(__FUNCTION__, " error ", GetLastError());
         return false;
        }
      for(uint i = 0; i<m_shells.Size(); i++)
        {
         v = data.Col(i);
         if(!m_ecdf.fit(v))
           {
            Print(__FUNCTION__, " error ");
            return false;
           }

         np::quickSort(v,true,0,long(v.Size()-1));
         slopes = np::unique(v);
         slope_samples = m_ecdf.ecdf(slopes);
         m_shells[i] = new CSpline1DInterpolantShell();
         if(!np::vecAsArray(slopes,xa) || !np::vecAsArray(slope_samples,ya))
           {
            Print(__FUNCTION__, " error ");
            return false;
           }
         m_arrays[i*2+1] = new CArrayDouble();
         m_arrays[i*2+2] = new CArrayDouble();
         if(!m_arrays[i*2+1].AddArray(xa) || !m_arrays[i*2+2].AddArray(ya) ||!m_internal.Add(m_arrays[i*2+1])||!m_internal.Add(m_arrays[i*2+2]))
           {
            Print(__FUNCTION__, " error ", GetLastError());
            return false;
           }
         CAlglib::Spline1DBuildLinear(xa,ya,m_shells[i]);
        }

      return true;
     }

   matrix            to_quantile(matrix &in)
     {
      if(in.Cols()!=m_cols)
        {
         Print(__FUNCTION__, " invalid input ");
         return matrix::Zeros(0,0);
        }

      vector temp;
      matrix out(in.Rows(),m_cols);
      for(ulong i = 0; i<m_cols; ++i)
        {
         temp  = in.Col(i);
         for(ulong j = 0; j<in.Rows(); ++j)
            out[j][i] = MathMax(MathMin(CAlglib::Spline1DCalc(m_shells[i],temp[j]),m_ub),m_lb);
        }

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

Метод fit() принимает матрицу с любым количеством столбцов, которые предполагаются переменными многомерного набора данных. Необязательные входные параметры upper_bound и lower_bound задают ограничения, предотвращающие появление значений, строго равных 0 или 1. Внутри значения стандартной ECDF берутся в уникальных отсортированных точках выборки. Затем на основе этих точек и соответствующих значений ECDF создаётся одномерная модель линейной интерполяции. 

//+------------------------------------------------------------------+
//|                                                    ECDF_Demo.mq5 |
//|                                  Copyright 2025, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2025, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#property version   "1.00"
#property script_show_inputs
#include<ECDF/linear_cdf.mqh>
//--- input parameters
input string   FirstSymbol = "NZDUSD";
input string   SecondSymbol = "AUDUSD";
input ENUM_TIMEFRAMES TimeFrame = PERIOD_D1;
input datetime StartDate=D'2024.01.01 00:00:01';
input ulong HistoryLength = 1000;
input bool Use_Returns = false;
input int DisplayDuration = 30;//display duration in seconds
//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
//---download symbol data
   vector syA,syB;
   if(!syA.CopyRates(FirstSymbol,TimeFrame,COPY_RATES_CLOSE,StartDate,HistoryLength) ||
      !syB.CopyRates(SecondSymbol,TimeFrame,COPY_RATES_CLOSE,StartDate,HistoryLength) ||
      syA.Size()!=syB.Size())
     {
      Print("Error downloading price data ", ::GetLastError());
      return;
     }
//---optionally transform to returns
   if(Use_Returns)
     {
      syA = np::diff(log(syA));
      syB = np::diff(log(syB));
     }
//---create matrix container and fill it
   matrix xy(syA.Size(),2);
   if(!xy.Col(syA,0) || !xy.Col(syB,1))
     {
      Print(" column insertion error ", GetLastError());
      return ;
     }
//---build ECDF model
   CLinearCDF qd;
   if(!qd.fit(xy))
     {
      Print(" CLinearCDF failure ");
      return;
     }
//---transform raw values to quantiles
   matrix yx = qd.to_quantile(xy);
//---display QQ plot
   np::scatter_plot(yx.Col(0),yx.Col(1),"Q-Q plot",FirstSymbol,SecondSymbol,true,0,0,0,0,500,500,true,30);
//---
  }
//+------------------------------------------------------------------

Скрипт ECDF_Demo.ex5 демонстрирует использование класса CLinearCDF, применяя его к двумерному набору ценовых данных. Затем преобразованные значения используются для отображения Q-Q-графика (квантиль-квантильный график) данных. 

Определение копулы

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

Для случайных величин, каждая из которых имеет собственное маргинальное распределение, копула связывает эти маргинальные распределения в совместное распределение, описывая, как переменные соотносятся друг с другом независимо от формы их индивидуальных распределений. Такое разделение маргиналей и зависимости является центральным элементом современного многомерного анализа, особенно когда предположения о нормальности или линейной зависимости неприменимы. Математически это формализуется теоремой Склара. Теорема Склара утверждает, что для любой непрерывной совместной функции распределения F существует единственная копула C такая, что:

где F1​(x1​) и F2​(x2​) — маргинальные CDF величин X1​ и X2​ соответственно.

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

Подобно тому как коэффициент корреляции ограничен диапазоном от −1 до 1, копулы характеризуются аналогичными границами, называемыми границами Фреше—Хёффдинга. Нижняя и верхняя границы Фреше—Хёффдинга — это функции, определяющие максимально отрицательную и положительную зависимость, которую может описать конкретная копула. Сами границы Фреше—Хёффдинга также являются копулами, где d соответствует размерности, то есть количеству переменных.

Хвостовая зависимость

Уникальная природа отдельных функций копул позволяет им отражать разные типы структур зависимости. Свойство копулы, которое обычно представляет интерес в финансовых приложениях, — это способность моделировать хвостовую зависимость. Хвостовая зависимость является мерой совместного движения случайных величин, когда они принимают экстремальные значения. Тип хвостовой зависимости, который может моделировать копула, определяет, как структура ведёт себя при совместных экстремальных событиях. Верхний хвост относится к крайне высоким значениям распределения, а нижний хвост — к крайне низким. Хвостовая зависимость количественно выражается коэффициентом хвостовой зависимости (λ), который вычисляется как предел условной вероятности. Без выражения через копулу хвостовую зависимость трудно количественно оценить.

Нижняя хвостовая зависимость измеряет вероятность того, что одна переменная окажется крайне низкой при условии, что другая переменная также крайне низка.

Верхняя хвостовая зависимость измеряет вероятность того, что одна переменная окажется крайне высокой при условии, что другая переменная также крайне высока.

Типы копул

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

Эллиптические копулы — это класс многомерных копул, получаемых из эллиптических распределений, таких как многомерное нормальное (гауссово) распределение и распределение Стьюдента. Они в основном параметризуются некоторой мерой связи, которая отражает линейную структуру зависимости между переменными. Две наиболее распространённые эллиптические копулы — это гауссова копула и t-копула Стьюдента. В следующем разделе мы перейдём к реализации двумерных эллиптических копул, начиная с бивариантной гауссовой копулы.

Двумерная гауссова копула

Двумерная гауссова копула выводится из стандартного двумерного нормального распределения. Она характеризуется одним параметром: линейным коэффициентом корреляции (ρ) базового двумерного нормального распределения. Этот параметр определяет силу и направление зависимости между двумя случайными величинами. Гауссова копула радиально симметрична и популярна благодаря вычислительной эффективности, хотя обладает слабой хвостовой зависимостью. Это означает, что она недостаточно хорошо описывает сильные одновременные экстремальные события в переменных.

Заголовочный файл gaussian.mqh содержит класс CGaussian, реализующий двумерную гауссову копулу. CGaussian наследуется от абстрактного класса CBivariateCopula, определённого в base.mqh, и служит основой для всех реализаций двумерных копул, которые мы будем описывать.

//+------------------------------------------------------------------+
//| base copula class                                                |
//+------------------------------------------------------------------+
class CBivariateCopula:public CObject
  {
private:
   int               m_len;
protected:
   double            m_eps;
   double            m_theta;
   double            m_rho;
   double            m_nu;
   double            m_threshold;
   ENUM_COPULA_TYPE  m_copula_type;
   vector            m_bounds;
   double            m_tau;
   matrix            m_cov;
   virtual bool      check_theta(void) 
     {
       return (m_theta>m_bounds.Min() < m_theta < m_bounds.Max());  
     }
   virtual double    pdf(double u,double v) { return EMPTY_VALUE; }
   virtual double    cdf(double u,double v) { return EMPTY_VALUE; }
   virtual double    condi_cdf(double u,double v) { return EMPTY_VALUE; }
   virtual vector    pdf(vector& u,vector& v){ return vector::Zeros(0); }
   virtual vector    cdf(vector& u,vector& v){ return vector::Zeros(0); }
   virtual vector    condi_cdf(vector& u,vector& v){ return vector::Zeros(0); }
   virtual double    theta_hat(const double tau){ return EMPTY_VALUE; }
   void              preprocess(double &u, double &v)
     {
      u = MathMin(MathMax(m_eps,u),1.0-m_eps);
      v = MathMin(MathMax(m_eps,v),1.0-m_eps);
     }
   void              preprocess(vector &u, vector &v)
     {
      u.Clip(m_eps,1.0-m_eps);
      v.Clip(m_eps,1.0-m_eps);
     }
   double            kendalTau(const vector &vec1,const vector &vec2)
     {
      double tau = double("nan");
      ulong size=vec1.Size();
      if(size==0 || vec2.Size()!=size)
        {
         Print(__FUNCTION__, " size of input vectors donot match ");
         return(tau);
        }
      //---
      long cnt1=0,cnt2=0,cnt=0;
      //---
      for(long i=0; i<long(size); i++)
        {
         for(long j=i+1; j<long(size); j++)
           {
            double delta1=vec1[i]-vec1[j];
            double delta2=vec2[i]-vec2[j];
            double delta=delta1*delta2;
            if(delta==0)
              {
               if(delta1!=0)
                  cnt1++;
               if(delta2!=0)
                  cnt2++;
              }
            else
              {
               cnt1++;
               cnt2++;
               if(delta>0.0)
                  cnt++;
               else
                  cnt--;
              }
           }
        }
      //--- calculate Kendall tau
      long den=cnt1*cnt2;
      if(den==0)
        {
         Print(__FUNCTION__, " failed zero check at line 76");
         return tau;
        }
      tau=double(cnt)/MathSqrt(den);
      //---
      return(tau);
     }

public:
                     CBivariateCopula(void)
     {
      m_len = 6;
      m_eps = 1.e-5;
      m_theta = m_rho = m_nu = EMPTY_VALUE;
      m_copula_type = WRONG_VALUE;
      m_cov  = matrix::Zeros(0,0);
      m_bounds = vector::Zeros(2);
      m_bounds[0] = -DBL_MIN;
      m_bounds[1] = DBL_MAX;
     }
                    ~CBivariateCopula(void)
     {
     }
   double            Get_theta(void)  { return m_theta; }
   double            Get_tau(void)    { return m_tau;   }
   double            Get_rho(void)    { return m_rho;   }
   matrix            Get_covar(void)  { return m_cov;   }
   double            Get_nu(void)     { return m_nu;    }
   double            Get_threshold(void) { return m_threshold; }
   double            Get_eps(void)    { return m_eps;   }
   virtual void              Set_theta(double theta) {  m_theta = theta; }
   virtual void              Set_tau(double tau) { m_tau = tau; }
   virtual void              Set_threshold(double threshold) { m_threshold = threshold; }
   virtual void              Set_nu(double nu) { m_nu = nu; }
   virtual void              Set_rho(double rho) { m_rho = rho; }
   virtual void              Set_eps(double eps) { m_eps = eps; }
   virtual matrix            Sample(ulong num_samples) { return matrix::Zeros(0,0); }
   virtual void              Set_covariance(matrix &cov)
     {
      m_cov = cov;
      m_rho = m_cov[0][1] / (sqrt(m_cov[0][0]) * sqrt(m_cov[1][1]));
     }
   virtual double    Fit(vector &u, vector&v)
     {
      if(u.Max()>1.0 || v.Max()>1.0 || v.Min()<0.0 ||u.Min()<0.0)
       {
        Print(__FUNCTION__, " Invalid input variable(s) ");
        return EMPTY_VALUE;
       }
      m_tau = kendalTau(u,v);
      m_theta = theta_hat(m_tau);
      if(!check_theta())
        Print(__FUNCTION__, " Invalid theta " );
      return  m_theta;
     }
   double            Log_likelihood(vector &u, vector &v)
     {
      if(u.Size()!=v.Size())
        {
         Print(__FUNCTION__, " vectors are not of equal length ");
         return EMPTY_VALUE;
        }
      
      vector ll = pdf(u,v);

      return (log(ll)).Sum();
     }
   double            Copula_PDF(double u, double v)
     {
      preprocess(u,v);
      return pdf(u,v);
     }
   double            Copula_CDF(double u, double v)
     {
      preprocess(u,v);
      return cdf(u,v);
     }
   double            Conditional_Probability(double u, double v)
     {
      preprocess(u,v);
      return condi_cdf(u,v);
     }
   vector            Copula_PDF(vector& u, vector& v)
     {
      preprocess(u,v);
      return pdf(u,v);
     }
   vector            Copula_CDF(vector& u, vector& v)
     {
      preprocess(u,v);
      return cdf(u,v);
     }
   vector            Conditional_Probability(vector& u, vector& v)
     {
      preprocess(u,v);
      return condi_cdf(u,v);
     }
   double            Theta(double tau = 0.0)
     {
      if(!tau && m_tau)
         return theta_hat(m_tau);
      else
         if(tau)
            return theta_hat(tau);
         else
           {
            Print(__FUNCTION__ " invalid input parameter 'tau' ");
            return EMPTY_VALUE;
           }
     }

   virtual bool      Save(const int file_handle)
     {
      if(file_handle!=INVALID_HANDLE)
        {
         //---
         if(FileWriteLong(file_handle,-1)==sizeof(long))
           {
            //---
            if(FileWriteInteger(file_handle,int(m_copula_type),INT_VALUE)!=INT_VALUE)
               return(false);
            //---
            if(FileWriteInteger(file_handle,m_len,INT_VALUE)!=INT_VALUE)
               return(false);
            //---
            if(FileWriteDouble(file_handle,m_eps)!=sizeof(double) ||
               FileWriteDouble(file_handle,m_theta)!=sizeof(double) ||
               FileWriteDouble(file_handle,m_tau)!=sizeof(double) ||
               FileWriteDouble(file_handle,m_rho)!=sizeof(double) ||
               FileWriteDouble(file_handle,m_nu)!=sizeof(double) ||
               FileWriteDouble(file_handle,m_threshold)!=sizeof(double))
               return false;
            return true;
           }
        }
      return false;
     }
   virtual bool      Load(const int file_handle)
     {
      if(file_handle!=INVALID_HANDLE)
        {
         //---
         if(FileReadLong(file_handle)==-1)
           {
            //---
            m_copula_type = ENUM_COPULA_TYPE(FileReadInteger(file_handle,INT_VALUE));

            if(FileReadInteger(file_handle,INT_VALUE)!=m_len)
               return false;

            ResetLastError();
            m_eps=FileReadDouble(file_handle);
            m_theta=FileReadDouble(file_handle);
            m_tau=FileReadDouble(file_handle);
            m_rho=FileReadDouble(file_handle);
            m_nu=FileReadDouble(file_handle);
            m_threshold=FileReadDouble(file_handle);
            if(GetLastError())
              {
               Print(__FUNCTION__, " possible read error ", GetLastError());
               return false;
              }
            return true;
           }
        }
      return false;
     }
   virtual int       Type(void)
     {
      return int(m_copula_type);
     }
  };

Он позволяет задавать модели копул двумя способами. Модель можно определить явно, используя методы доступа для установки соответствующих параметров копулы. Альтернативно модель можно построить, подогнав её к двумерному набору равномерных переменных (маргинальных CDF исходных значений) с помощью метода Fit(). В случае успеха метод вернёт любое значение, отличное от константы EMPTY_VALUE.

Плотность копулы оценивается вызовом метода Copula_PDF(), а Copula_CDF() вычисляет совместную функцию распределения копулы. Условная CDF копулы рассчитывается вызовом метода Conditional_Probability(). Все перечисленные методы перегружены для приёма скалярных и векторных входных данных. Метод Sample() используется для генерации синтетических наборов данных, соответствующих параметрам модели копулы. Наконец, Load() и Save() обеспечивают сериализацию и десериализацию моделей копул для постоянного хранения.

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

virtual void Set_covariance(matrix& cov) override
     {
      if(cov.Rows()!=2 || cov.Cols()!=2)
       {
        Print(__FUNCTION__," invalid input: expecting 2x2 matrix ");
        return;
       }
      m_cov = cov;
      m_rho = m_cov[0][1] / (sqrt(m_cov[0][0]) * sqrt(m_cov[1][1]));
      m_tau = (2.0/M_PI) * MathArcsin(m_rho);
      m_theta = theta_hat(m_tau);
     }

Скрипт SampleEllipticCopula.ex5 проверяет, соответствуют ли синтетические данные, сгенерированные из эллиптической копулы, заданным параметрам модели. Настраиваемые пользователем входные параметры с префиксом Covar задают элементы ковариационной матрицы. Параметр Degree_Of_Freedom не относится к тестированию гауссовой копулы. Параметр Size определяет количество генерируемых выборок.

//+------------------------------------------------------------------+
//|                                         SampleEllipticCopula.mq5 |
//|                                  Copyright 2025, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2025, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#property version   "1.00"
#property script_show_inputs
#include<Copulas/Bivariate/gaussian.mqh>
#include<Copulas/Bivariate/t.mqh>
//+------------------------------------------------------------------+
//|  enum                                                            |
//+------------------------------------------------------------------+
enum ENUM_ELPT_COPULA
 {
  GAUSSIAN=0,//Gaussian
  STUDENT//Student
 };
//--- input parameters
input ENUM_ELPT_COPULA copula_to_sample_from = GAUSSIAN;
input double Covar_0_0_ = 2.0;
input double Covar_0_1_ = 0.5;
input double Covar_1_0_ = 0.5;
input double Covar_1_1_ = 2.0;
input double Degrees_Of_Freedom = 5.0;
input ulong Size = 1000;

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

//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
//---
   matrix covar(2,2);
//---
   covar[0,0] = Covar_0_0_;
   covar[0,1] = Covar_0_1_;
   covar[1,0] = Covar_1_0_;
   covar[1,1] = Covar_1_1_;
//---
   CBivariateCopula *copula = NULL;
//---
   double set_theta = 0.0;
   switch(copula_to_sample_from)
    {
     case GAUSSIAN:
       copula = (CBivariateCopula*)new CGaussian();
       copula.Set_covariance(covar);
       break;
     case STUDENT:
       copula = (CBivariateCopula*)new CStudent();
       copula.Set_covariance(covar);
       copula.Set_nu(Degrees_Of_Freedom);
       break;
    }
//---
   Print("Specified correlation parameter of ", EnumToString(copula_to_sample_from)," copula :",copula.Get_rho());
   matrix samples = copula.Sample(Size);
   matrix inv_samples = samples;
   
   for(ulong i = 0; i<inv_samples.Rows(); ++i)
     for(ulong j =0; j<inv_samples.Cols(); ++j)
          inv_samples[i,j] = copula_to_sample_from == GAUSSIAN?CAlglib::InvNormalCDF(inv_samples[i,j]):studentTQuantile(inv_samples[i,j],Degrees_Of_Freedom);
          
   Print("Correlation of matrix of variable sampled from  ", EnumToString(copula_to_sample_from)," copula :\n", inv_samples.CorrCoef(false));      
   string cname = EnumToString((ENUM_COPULA_TYPE)copula.Type());
//---
   delete copula;
//---display scatter plot
   np::scatter_plot(inv_samples.Col(0),inv_samples.Col(1),cname+" Scatter plot","X","Y",true,0,0,0,0,500,300,true,10);  
//---
  }

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

CQ      0       14:06:16.963    SampleEllipticCopula (XAUEUR,D1)        Specified correlation parameter of GAUSSIAN copula :0.24999999999999994
DJ      0       14:06:17.021    SampleEllipticCopula (XAUEUR,D1)        Correlation of matrix of variable sampled from  GAUSSIAN copula :
DG      0       14:06:17.022    SampleEllipticCopula (XAUEUR,D1)        [[1,0.2497651923653427]
RK      0       14:06:17.022    SampleEllipticCopula (XAUEUR,D1)         [0.2497651923653427,1]]

В следующем разделе мы обратимся к другой эллиптической копуле — t-копуле Стьюдента.

Двумерная t-копула Стьюдента

Двумерная t-копула Стьюдента является одной из наиболее заметных альтернатив гауссовой копуле, поскольку она учитывает хвостовую зависимость, которую гауссова копула описать не может. t-копула Стьюдента определяется двумя параметрами: коэффициентом корреляции и числом степеней свободы, обычно обозначаемыми ρ и ν соответственно. В отличие от гауссовой копулы, t-копула имеет ненулевую хвостовую зависимость, то есть экстремальное событие в одной переменной повышает вероятность экстремального события в другой.

где:

• t ρ , ν tρ,ν​ = CDF двумерного t-распределения Стьюдента с корреляцией ρ ρ и ν ν степенями свободы,
• t ν − 1 tν−1​ = обратная CDF (квантильная функция) одномерного t-распределения с ν ν степенями свободы,
• u , v ∈ [ 0 , 1 ] u,v∈[0,1].

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

t-копула Стьюдента реализована в t.mqh как класс CStudent, который также является подклассом CBivariateCopula.

//+------------------------------------------------------------------+
//|                                                            t.mqh |
//|                                  Copyright 2025, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2025, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#include "base.mqh"
#ifndef _NP_
#include<np.mqh>
#endif
#include<Math/Stat/T.mqh>
//+------------------------------------------------------------------+
//|utilities for Student's t distribution: Cheng Algorithm           |
//+------------------------------------------------------------------+
double            chengFuAlgo(double probability,double df)
  {
   double k = ceil(df*0.5);
   double a = 1.0 - probability;
   double qi = EMPTY_VALUE;
   if(a!=0.5)
     {
      qi = sqrt(2.0*pow(1.0-2.0*a,2.0)/(1.0 - pow(1.0-2.0*a,2.0)));
      double i = 0.0;
      double gy,j;
      while(i<20.0)
        {
         gy = 0.0;
         j = 0.0;
         while(j<=k-1.0)
           {
            gy+=MathFactorial(int(2*j))/pow(2.0,2.0*j)/pow(MathFactorial(int(j)),2.0)*pow(1.0+pow(qi,2.0)/(2.0*k),-j);
            j+=1.0;
           }
         qi = 1.0/sqrt(1.0/(2.0*k)*(pow(gy/1.0-2.0*a,2.0) - 1.0));
         i+=1.0;
        }
      if(a>0.5)
         return -qi;
      else
         return qi;
     }
   else
      return 0.0;
  }
//+------------------------------------------------------------------+
//| Hills algorithm                                                  |
//+------------------------------------------------------------------+
double            hillsAlgo(double probability,double df)
  {
   double z,a,b,c,x,y,d;
   z = a = b = c = x = y = d = EMPTY_VALUE;
   bool neg = false;
   if(probability>0.5)
     {
      z = 2.0*(1.0-probability);
     }
   else
     {
      neg = true;
      z = 2.0*probability;
     }
   a = 1.0/(df - 0.5);
   b = 48.0/(a*a);
   c = ((20700.0*a/b - 98.0)*a - 16.0)*a + 96.36;
   d = ((94.5/(b + c) - 3.0)/b + 1.0)*sqrt(a*M_PI/2.0)*df;
   x = z*d;
   y = pow(x,(2.0/df));

   if(y>0.05 + a)
     {
      x = CAlglib::InvNormalCDF(z*0.5);
      y = x*x;
      if(df<5.0)
         c = c+0.3*(df-4.5)*(x+0.6);
      c = c + (((0.05*d*x - 5.0)*x - 7.0)*x - 2.0)*x + b;
      y = (((((0.4*y + 6.3)*y + 36.0)*y + 94.5)/c - y - 3.0)/b + 1.0)*x;
      y = a*y*y;
      if(y > 0.002)
         y = exp(y) - 1.0;
      else
         y = y + 0.5*y*y;
     }
   else
      y = ((1.0/(((df + 6.0)/(df*y) - 0.089*d - 0.822)*(df + 2.0)*3.0) + 0.5/(df + 4.0))*y - 1.0)*(df + 1.0)/(df + 2.0) + 1.0/y;
   double q = sqrt(df*y);
   return (neg)?-1.0*q:q;
  }
//+------------------------------------------------------------------+
//| Student's t Inverse CDF                                          |
//+------------------------------------------------------------------+
double            studentTQuantile(double probability, double df)
  {
   if(df == 1.0)
      return tan(M_PI*(probability-1.0/2.0));
   else
      if(df == 2.0)
         return (2.0*probability-1.0)*sqrt(2.0/(4.0*probability*(1.0-probability)));
      else
         if(df == 4.0)
           {
            double a = 4.0*probability*(1.0 - probability);
            double q = cos(1.0/3.0*acos(sqrt(a)))/sqrt(a);
            return np::sign(probability - 0.5)*2.0*sqrt(q-1.0);
           }
         else
            if(!MathMod(df,2.0))
               return chengFuAlgo(probability,df);
            else
               return hillsAlgo(probability,df);
  }
//+------------------------------------------------------------------+
//| Student T Probability Density Function                           |
//+------------------------------------------------------------------+
double studentTDensity(double x, double df)
  {
   return CAlglib::GammaFunction((df + 1.0)/2.0)/(sqrt(df*M_PI)*CAlglib::GammaFunction(df/2.0))*pow((1.0 + pow(x,2.0)/df),(-((df + 1.0)/2.0)));
  }
//+------------------------------------------------------------------+
//| Bivariate Student Copula                                         |
//+------------------------------------------------------------------+
class CStudent : public CBivariateCopula
  {
private:
   virtual double    theta_hat(const double tau) override
     {
      return sin(tau*M_PI/2.0);
     }
   matrix            generate_corr_student(ulong num, matrix& cov, double nu_)
     {
      vector means = vector::Zeros(2);
      matrix normal = np::multivariate_normal(means,cov,num);

      double chi[];
      if(!MathRandomChiSquare(nu_,int(num),chi))
        {
         Print(__FUNCTION__, " Math random chisquare error ", GetLastError());
         return matrix::Zeros(0,0);
        }

      matrix out = matrix::Zeros(num,2);

      for(ulong i = 0; i<out.Rows(); i++)
        {
         double chisqrt = sqrt(chi[i]/nu_);
         out[i,0] = normal[i,0]/chisqrt;
         out[i,1] = normal[i,1]/chisqrt;
        }
      return out;
     }
   double            bvtdist(double z1,double z2, vector& mu, matrix& cov, double df)
     {
      double x1 = z1 - mu[0];
      double x2 = z2 - mu[1];
      double det_cov = cov[0,0]*cov[1,1]-cov[0,1]*cov[1,0];

      double xt = (-2.0*cov[0,1]*x1*x2+cov[0,0]*(pow(x1,2.0)+pow(x2,2.0)))/det_cov;
      double numerator = CAlglib::GammaFunction((2.0+df)/2.0);
      double denominator = (CAlglib::GammaFunction(df/2.0)*df*M_PI*sqrt(det_cov)*pow(1.0+xt/df,(2.0+df)/2.0));
      return numerator/denominator;
     }
   //+------------------------------------------------------------------+
   //|  inner integral function object                                  |
   //+------------------------------------------------------------------+
   class CInnerInt : public CIntegrator1_Func
     {
   private:
      double            m_x;
      double            m_rho;
      matrix            m_corr;
      vector            m_mu;
      double            m_nu;
   public:
      //---
                     CInnerInt(void)
        {
         m_mu = vector::Zeros(2);
        }
                    ~CInnerInt(void) {}
      void              Set(double x, double rho, double nu)
        {
         m_x = x;
         m_nu = nu;
         m_rho = rho;
         matrix crr = {{1.0, m_rho}, {m_rho, 1.0}};
         m_corr = crr;
        }
      virtual void      Int_Func(double x,double xminusa,double bminusx,double &y,CObject &obj)
        {
         double x1 = m_x - 0.0;
         double x2 = x - 0.0;
         double det_cov = m_corr[0,0]*m_corr[1,1]-m_corr[0,1]*m_corr[1,0];

         double xt = (-2.0*m_corr[0,1]*x1*x2+m_corr[0,0]*(pow(x1,2.0)+pow(x2,2.0)))/det_cov;
         double numerator = CAlglib::GammaFunction((2.0+m_nu)/2.0);
         double denominator = (CAlglib::GammaFunction(m_nu/2.0)*m_nu*M_PI*sqrt(det_cov)*pow(1.0+xt/m_nu,(2.0+m_nu)/2.0));
         y =  numerator/denominator;
        }
     };
   //+------------------------------------------------------------------+
   //|  outer integral function object                                  |
   //+------------------------------------------------------------------+
   class COuterInt : public CIntegrator1_Func
     {
   public:
      double            ay_limit;
      double            by_limit;
      double            nu;
      double            rho;
      //---
                     COuterInt(void) {}
                    ~COuterInt(void) {}
      virtual void      Int_Func(double x,double xminusa,double bminusx,double &y,CObject &obj)
        {
         CInnerInt fint;
         //---
         fint.Set(x,rho,nu);
         //---
         CAutoGKStateShell s;
         //---
         double integral;
         //---
         CObject ob;
         //---
         CAutoGKReportShell rep;
         //---
         CAlglib::AutoGKSmooth(ay_limit,by_limit,s);
         //---
         CAlglib::AutoGKIntegrate(s,fint,ob);
         //---
         CAlglib::AutoGKResults(s,integral,rep);
         //---
         CAutoGKReport report = rep.GetInnerObj();
         //---
         if(report.m_terminationtype<0.0)
            Print(__FUNCTION__, " integration error ",report.m_terminationtype);
         //---
         y = integral;
        }
     };
protected:
   virtual double    pdf(double u,double v) override
     {
      double y1,y2;
      y1 = studentTQuantile(u,m_nu);
      y2 = studentTQuantile(v,m_nu);
      vector means = vector::Zeros(2);
      matrix corr = {{1.0,m_rho},{m_rho,1.0}};
      double numerator = bvtdist(y1,y2,means,corr,m_nu);
      double denominator,pdf1,pdf2;
      pdf1 = studentTDensity(y1,m_nu);
      pdf2 = studentTDensity(y2,m_nu);
      denominator = pdf1*pdf2;
      return numerator/denominator;
     }
   virtual double    cdf(double u,double v) override
     {
      double uu,vv;
      uu = MathMax(u,1.e-6);
      vv = MathMax(v,1.e-6);
      //---
      int errcode = 0;
      //---
      COuterInt fout;
      fout.rho = m_rho;
      fout.nu = m_nu;
      fout.ay_limit = double(LONG_MIN);
      fout.by_limit = studentTQuantile(vv,m_nu);
      double ax_limit = double(LONG_MIN);
      double bx_limit = studentTQuantile(uu,m_nu);
      //---
      CObject obj;
      CAutoGKStateShell ss;
      //---
      double outer_integral;
      CAutoGKReportShell repp;
      //---
      CAlglib::AutoGKSmooth(ax_limit,bx_limit,ss);
      //---
      CAlglib::AutoGKIntegrate(ss,fout,obj);
      //---
      CAlglib::AutoGKResults(ss,outer_integral,repp);
      //---
      CAutoGKReport report = repp.GetInnerObj();
      //---
      if(report.m_terminationtype<0.0)
         Print(__FUNCTION__, " integration error ",report.m_terminationtype);
      //---
      return MathMax(MathMin(outer_integral,1.0),0.0);
     }
   virtual double    condi_cdf(double u,double v) override
     {
      double inv_u,inv_v;
      int errcode = 0;
      inv_u = studentTQuantile(u,m_nu);
      inv_v = studentTQuantile(v,m_nu);
      double numerator,denominator;
      numerator = (inv_u-m_rho*inv_v)*sqrt(m_nu+1.0);
      denominator = sqrt((1.0-pow(m_rho,2.))*(pow(inv_v,2.0)+m_nu));
      double tcdf = MathCumulativeDistributionT(numerator/denominator,m_nu+1.0,errcode);
      if(errcode)
        {
         Print(__FUNCTION__, " mathcdf error ", errcode);
         return EMPTY_VALUE;
        }
      return tcdf;
     }
   virtual vector    pdf(vector &u,vector &v) override
     {
      vector out(u.Size());
      for(ulong i = 0; i<u.Size(); ++i)
         out[i] = pdf(u[i],v[i]);
      return out;
     }
   virtual vector    cdf(vector &u,vector &v) override
     {
      vector out(u.Size());
      for(ulong i = 0; i<u.Size(); ++i)
         out[i] = cdf(u[i],v[i]);
      return out;
     }
   virtual vector    condi_cdf(vector &u,vector &v) override
     {
      vector out(u.Size());
      for(ulong i = 0; i<u.Size(); ++i)
         out[i] = condi_cdf(u[i],v[i]);
      return out;
     }
public:
                     CStudent(void)
     {
      m_copula_type = STUDENT_COPULA;
     }
                    ~CStudent(void)
     {
     }

   virtual matrix    Sample(ulong num_samples) override
     {
      matrix spairs = generate_corr_student(num_samples,m_cov,m_nu);
      matrix out = spairs;
      int err1,err2;
      err1 = err2 = 0;
      for(ulong i = 0; i<spairs.Rows(); ++i)
        {
         out[i,0] = MathCumulativeDistributionT(spairs[i,0],m_nu,err1);
         out[i,1] = MathCumulativeDistributionT(spairs[i,1],m_nu,err2);
         if(err1 || err2)
           {
            Print(__FUNCTION__, " mathcdf error ", err1?err1:err2);
            return matrix::Zeros(0,0);
           }
        }
      return out;
     }
   virtual double    Fit(vector &u, vector&v) override
     {
      if(u.Max()>1.0 || v.Max()>1.0 || v.Min()<0.0 ||u.Min()<0.0)
        {
         Print(__FUNCTION__, " Invalid input variable(s) ");
         return EMPTY_VALUE;
        }
      m_tau = kendalTau(u,v);
      m_theta = theta_hat(m_tau);
      if(m_nu == EMPTY_VALUE)
         m_nu = fit_nu(u,v);
      matrix vals = matrix::Zeros(u.Size(),2);
      for(ulong i = 0; i<vals.Rows(); ++i)
        {
         vals[i,0] = studentTQuantile(u[i],m_nu);
         vals[i,1] = studentTQuantile(v[i],m_nu);
        }
      m_cov = vals.Cov(false);
      m_rho = m_cov[0,1]/(sqrt(m_cov[0,0])*sqrt(m_cov[1,1]));
      return m_rho;
     }
  };

//+------------------------------------------------------------------+

Ниже приведён вывод двух запусков SampleEllipticCopula.ex5 с выбранным соответствующим вариантом t-копулы.

IM      0       14:07:00.097    SampleEllipticCopula (XAUEUR,D1)        Specified correlation parameter of STUDENT copula :0.24999999999999994
NQ      0       14:07:00.295    SampleEllipticCopula (XAUEUR,D1)        Correlation of matrix of variable sampled from  STUDENT copula :
GF      0       14:07:00.295    SampleEllipticCopula (XAUEUR,D1)        [[1,0.2441795870726779]
MJ      0       14:07:00.295    SampleEllipticCopula (XAUEUR,D1)         [0.2441795870726779,1]]

Тестирование и проверка

Для тестирования остальных методов обеих реализаций копул предусмотрены скрипты TestBivariateGaussianCopula.ex5 и TestBivariateStudentTCopula.ex5. Оба скрипта создают конкретные модели копул и проверяют код PDF, CDF и условной CDF. Эти тесты были воспроизведены из модульных тестов ArbitrageLab Python-пакета, на котором основан представленный код MQL5. Это позволяет проверить, работает ли код MQL5 так же, как исходная реализация на Python.

//+------------------------------------------------------------------+
//|                                  TestBivariateGaussianCopula.mq5 |
//|                                  Copyright 2025, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2025, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#property version   "1.00"
//--- input parameters
#include<Copulas/Bivariate/gaussian.mqh>
//---
 double Covar_0_0_ = 2.0;
 double Covar_0_1_ = 0.5;
 double Covar_1_0_ = 0.5;
 double Covar_1_1_ = 2.0;
//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
//---
   matrix covar(2,2);
//---
   covar[0,0] = Covar_0_0_;
   covar[0,1] = Covar_0_1_;
   covar[1,0] = Covar_1_0_;
   covar[1,1] = Covar_1_1_;
//---   
   CGaussian norm_copula;
   norm_copula.Set_covariance(covar);
//---   
   Print("Rho ", norm_copula.Get_rho());
   Print("Joint cdf(0.7,0.0001) ",norm_copula.Copula_CDF(0.7,0.0001));
   Print("Joint cdf(0.0001,0.7) ",norm_copula.Copula_CDF(0.0001,0.7));
   Print("Joint cdf(0.7,1.0) ",norm_copula.Copula_CDF(0.7,1.0));
   Print("Joint cdf(0.5,0.7) ",norm_copula.Copula_CDF(0.5,0.7));
//---   
   Print("Joint pdf(0.5,0.7) ",norm_copula.Copula_PDF(0.5,0.7));
   Print("Joint pdf(0.7,0.5) ",norm_copula.Copula_PDF(0.7,0.5));
   Print("Joint pdf(0.6,0.7) ",norm_copula.Copula_PDF(0.6,0.7));
//---   
   Print("Conditional CDF (U<=0.5|V=0.7) ",norm_copula.Conditional_Probability(0.5,0.7));
  }
//+------------------------------------------------------------------+

Ниже приведены результаты теста методов CGaussian.

KR      0       14:40:26.403    TestBivariateGaussianCopula (XAUEUR,D1) Rho 0.24999999999999994
JR      0       14:40:26.403    TestBivariateGaussianCopula (XAUEUR,D1) Joint cdf(0.7,0.0001) 0.00009407442939459683
JG      0       14:40:26.403    TestBivariateGaussianCopula (XAUEUR,D1) Joint cdf(0.0001,0.7) 0.00009407442939459683
KK      0       14:40:26.403    TestBivariateGaussianCopula (XAUEUR,D1) Joint cdf(0.7,1.0) 0.6999973036455064
FL      0       14:40:26.403    TestBivariateGaussianCopula (XAUEUR,D1) Joint cdf(0.5,0.7) 0.38494413131861455
IR      0       14:40:26.403    TestBivariateGaussianCopula (XAUEUR,D1) Joint pdf(0.5,0.7) 1.0233716657780745
KG      0       14:40:26.403    TestBivariateGaussianCopula (XAUEUR,D1) Joint pdf(0.7,0.5) 1.0233716657780745
LK      0       14:40:26.403    TestBivariateGaussianCopula (XAUEUR,D1) Joint pdf(0.6,0.7) 1.0580116369285744
EP      0       14:40:26.403    TestBivariateGaussianCopula (XAUEUR,D1) Conditional CDF (U<=0.5|V=0.7) 0.4461479582632463

Для сравнения ниже приведён Python-код модульного теста, относящегося к гауссовой копуле.

def test_gaussian(self):
        """
        Test Gaussian copula class.
        """

        cov = [[2, 0.5], [0.5, 2]]
        cop = GaussianCopula(cov=cov)
        # Check describe
        descr = cop.describe()
        self.assertEqual(descr['Descriptive Name'], 'Bivariate Gaussian Copula')
        self.assertEqual(descr['Class Name'], 'Gaussian')
        self.assertEqual(descr['cov'], cov)
        self.assertAlmostEqual(descr['rho'], 0.25, delta=1e-5)

        # Check copula joint cumulative density C(U=u,V=v)
        self.assertAlmostEqual(cop.C(0.7, 1e-4), cop.C(1e-4, 0.7), delta=1e-4)
        self.assertAlmostEqual(cop.C(0.7, 1), cop.C(1, 0.7), delta=1e-4)
        self.assertAlmostEqual(cop.C(0.7, 1e-8), 0, delta=1e-4)
        self.assertAlmostEqual(cop.C(0.7, 1), 0.7, delta=1e-4)
        self.assertAlmostEqual(cop.C(0.5, 0.7), 0.384944, delta=1e-4)

        # Check copula joint probability density c(U=u,V=v)
        self.assertAlmostEqual(cop.c(0.5, 0.7), cop.c(0.7, 0.5), delta=1e-8)
        self.assertAlmostEqual(cop.c(0.5, 0.7), 1.023371665778, delta=1e-4)
        self.assertAlmostEqual(cop.c(0.6, 0.7), 1.058011636928, delta=1e-4)

        # Check copula conditional cdf Prob(U<=u|V=v)
        self.assertAlmostEqual(cop.condi_cdf(0.5, 0.7), 0.446148, delta=1e-4

Далее приведены результаты для класса CStudent.

ER      0       14:40:34.138    TestBivariateStudentTCopula (XAUEUR,D1) Rho 0.24999999999999994
JG      0       14:40:35.369    TestBivariateStudentTCopula (XAUEUR,D1) Joint cdf(0.0,0.2) 0.000006839679521344817
JD      0       14:40:36.613    TestBivariateStudentTCopula (XAUEUR,D1) Joint cdf(0.2,0.0) 0.000006839679521344818
HK      0       14:40:38.041    TestBivariateStudentTCopula (XAUEUR,D1) Joint cdf(0.7,1.0) 0.6999969941893684
PO      0       14:40:39.584    TestBivariateStudentTCopula (XAUEUR,D1) Joint cdf(1.0,0.7) 0.6999969941893683
OR      0       14:40:41.280    TestBivariateStudentTCopula (XAUEUR,D1) Joint cdf(1.0,1.0) 0.9999810861796908
RF      0       14:40:42.452    TestBivariateStudentTCopula (XAUEUR,D1) Joint cdf(0.0,0.0) 0.0000010861794979038194
HJ      0       14:40:43.757    TestBivariateStudentTCopula (XAUEUR,D1) Joint cdf(0.3,0.7) 0.23534922600994124
DL      0       14:40:43.757    TestBivariateStudentTCopula (XAUEUR,D1) Joint pdf(0.5,0.7) 1.0915055449624917
FP      0       14:40:43.757    TestBivariateStudentTCopula (XAUEUR,D1) Joint pdf(0.7,0.5) 1.0915055449624917
MD      0       14:40:43.757    TestBivariateStudentTCopula (XAUEUR,D1) Joint pdf(0.6,0.7) 1.1416004955531887
DD      0       14:40:43.757    TestBivariateStudentTCopula (XAUEUR,D1) Conditional CDF (U<=0.5|V=0.7) 0.441518430406228
MG      0       14:40:43.757    TestBivariateStudentTCopula (XAUEUR,D1) Loglikelihood 2.135711666815534

Модульный тест t-копулы Стьюдента.

def test_student(self):
        """
        Test Student copula class (Student-t).
        """

        cov = [[2, 0.5], [0.5, 2]]
        nu = 5
        cop = StudentCopula(cov=cov, nu=nu)
        # Check describe
        descr = cop.describe()
        self.assertEqual(descr['Descriptive Name'], 'Bivariate Student-t Copula')
        self.assertEqual(descr['Class Name'], 'Student')
        self.assertEqual(descr['cov'], cov)
        self.assertEqual(descr['nu (degrees of freedom)'], nu)
        self.assertAlmostEqual(descr['rho'], 0.25, delta=1e-5)
        # More to be added here for test on C(U<=u, V<=v)

        # Check copula joint probability density c(U=u,V=v)
        self.assertAlmostEqual(cop.c(0.5, 0.7), cop.c(0.7, 0.5), delta=1e-8)
        self.assertAlmostEqual(cop.c(0.5, 0.7), 1.09150554, delta=1e-4)
        self.assertAlmostEqual(cop.c(0.6, 0.7), 1.1416005, delta=1e-4)

        # Check copula joint cumulative density C(U=u,V=v)
        self.assertAlmostEqual(cop.C(0, 0.2), cop.C(0.2, 0), delta=1e-4)
        self.assertAlmostEqual(cop.C(0.7, 1), cop.C(1, 0.7), delta=1e-4)
        self.assertAlmostEqual(cop.C(1, 1), 1, delta=1e-4)
        self.assertAlmostEqual(cop.C(0, 0), 0, delta=1e-4)
        self.assertAlmostEqual(cop.C(0.3, 0.7), 0.23534923332657925, delta=1e-4)

        # Check copula conditional cdf Prob(U<=u|V=v)
        self.assertAlmostEqual(cop.condi_cdf(0.5, 0.7), 0.4415184293094455, delta=1e-4)

        # Check theta(tau)
        self.assertAlmostEqual(cop.theta_hat(2 * np.arcsin(0.2) / np.pi), 0.2, delta=1e-4)

        # log_ml function in ccalc
        u = np.array([0.1, 0.21, 0.5, 0.8])
        v = np.array([0.01, 0.25, 0.4, 0.7])
        new_cop = StudentCopula(nu=5)
        new_cop.fit(u, v)
        ll = new_cop.get_log_likelihood_sum(u=u,
                                            v=v)
        self.assertAlmostEqual(ll, 2.1357117471178584, delta=1e-5

Как видно из вывода, библиотека, по-видимому, работает корректно.

Заключительные замечания

Если сравнивать реализованные к настоящему моменту эллиптические копулы, t-копула может лучше подходить для моделирования совместных экстремальных событий. Однако она не способна описывать асимметричную хвостовую зависимость. Такая ситуация возникает, когда зависимость в одном хвосте сильнее, чем в другом. Для этого более полезными будут архимедовы копулы, такие как Clayton или Gumbel.

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

Файл или папка	Описание
MQL5/include/np.mqh	Заголовочный файл с различными служебными функциями для векторов и матриц.
MQL5/include/Copulas/Bivariate	Эта папка содержит все заголовочные файлы реализаций копул.
MQL5/include/ECDF	Эта папка содержит заголовочные файлы реализации эмпирической CDF, описанной в статье.
MQL5/script/ECDF_Demo.mq5	Этот скрипт демонстрирует применение реализации эмпирической CDF, описанной в статье.
MQL5/script/TestBivariateGaussianCopula.mq5	Этот скрипт демонстрирует использование класса CGaussian.
MQL5/script/TestBivariateStudentTCopula.mq5	Этот скрипт демонстрирует использование класса CStudent.
MQL5/script/SampleEllipticCopula.mq5	Этот скрипт демонстрирует генерацию выборки из эллиптической копулы.