In the previous article we presented the MQL5 implementations of bivariate elliptical copulae, specifically the Gaussian copula and Student’s t-copula. When it comes to the application of these copulae, there are at least two drawbacks: their algebraic expression is complicated, and their level of symmetry can be problematic. These problems have provided the impetus for the search for copulae with convenient algebraic expressions that can also model asymmetric dependency in datasets. One of the most popular families of copulae which satisfy both requirements is the class of Archimedean copulae. Therefore, in this article, we discuss common bivariate Archimedean copulae and their implementation in MQL5. We will see the expansion of the bivariate copulas library by adding implementations of Frank, Joe, Gumbel, Clayton, N13 and N14 copulae.
Introducing Archimedean copulae
A bivariate Archimedean copula is a specific type of copula C(u,v) used in statistics to model the dependence structure between two random variables with uniform marginal distributions. Its defining property is its generative form, which can be expressed using a single, continuous, strictly decreasing, and convex function called the generator function ϕ.
The generator ϕ must satisfy ϕ(1)=0. This structure gives Archimedean copulae a high degree of symmetry (C(u,v)=C(v,u)) and allows for a wide range of dependence structures to be modeled simply by choosing different generator functions. Archimedean copulae are fundamentally more flexible than common copulae, but that flexibility is usually reduced for practical use. Theoretically, there are an infinite number of choices for the generator function. Every slight change in the function creates a new, unique copula. In practice, one usually chooses a parametric family governed by a unique generator function, defined by a single scalar input variable. This makes estimation and modeling much easier. Therefore, most common Archimedean families are characterized by a single parameter embedded within the generator function, which governs the strength of the dependence.
The generator function defines an Archimedean copula and captures the entire dependence structure between the random variables. In simple terms, its role is to translate the marginal probabilities into a dependence structure that can be easily combined using simple addition. The function transforms the marginal variable. Because the generator is strictly decreasing, low probabilities get mapped to large numbers, and high probabilities get mapped to 0. This is an inversion of the probability scale. The specific shape and parameters of the generator function entirely determine the resulting copula family and, therefore, the exact way the two variables depend on each other. We begin our exploration of Archimedean copulae with the bivariate Frank copula in the next section.
Bivariate Frank copula
The bivariate Frank copula is a symmetric Archimedean copula that can model both positive and negative dependencies without exhibiting tail dependence. It is defined by the generator function:
where θ (theta) determines the strength and direction of association. θ values of zero are undefined for the Frank copula’s generator. Positive values of θ indicate positive dependence, while negative values correspond to negative dependence; as θ tends to 0, the copula approaches the independence copula. The Frank copula is radially symmetric, meaning it treats both tails equally, and does not exhibit upper or lower tail dependence. This makes it well-suited for modeling moderate dependence structures that are symmetric around the center of the distribution, such as in applications involving continuous variables where extreme co-movements are unlikely.
The class CFrank defined in frank.mqh, inherits from CBivariateCopula to represent a Frank copula.
#property copyright “Copyright 2025, MetaQuotes Ltd.” #property link “https://www.mql5.com” #include “base.mqh” class CFrank:public CBivariateCopula { private: class CInt_Function_1_Func : public CIntegrator1_Func { private: double m_th; public: CInt_Function_1_Func(void) {} ~CInt_Function_1_Func(void) {} void set_theta(double theta) { m_th = theta; } virtual void Int_Func(double x,double xminusa,double bminusx,double &y,CObject &obj) { y = x/m_th/(exp(x)-1.0); } }; class CBrent1:public CBrentQ { private: double debyel(double theta) { CInt_Function_1_Func fint; fint.set_theta(theta); CObject obj; CAutoGKStateShell s; double integral; CAutoGKReportShell rep; CAlglib::AutoGKSmooth(0.0,theta,s); CAlglib::AutoGKIntegrate(s,fint,obj); CAlglib::AutoGKResults(s,integral,rep); CAutoGKReport report = rep.GetInnerObj(); if(report.m_terminationtype<0.0) Print(__FUNCTION__, " integration error ",report.m_terminationtype); return integral; } public: CBrent1(void) { } ~CBrent1(void) { } virtual double objective(double u, double v,double y) { return (1.0 – 4.0/u+4.0*debyel(u)/u) – y; } }; virtual double theta_hat(const double tau) override { CBrent1 fun; return fun.minimize(0,tau,-100.0,100.0,2.e-12,4.0*2.e-16,100); } vector generate_pair(double v1, double v2) { vector out(2); out[0] = v1; out[1] = -1.0/m_theta*log(1.0+(v2*(1.0-exp(-m_theta)))/(v2*(exp(-m_theta*v1)-1.0)-exp(-m_theta*v1))); return out; } protected: virtual double pdf(double u,double v) override { double et,eut,evt,pd; et = exp(m_theta); eut = exp(u*m_theta); evt = exp(v*m_theta); pd = et*eut*evt*(et – 1.0)*m_theta/pow(et+eut*evt-et*eut-et*evt,2.0); return pd; } virtual double cdf(double u, double v) override { return (-1.0 / m_theta * log(1.0 + (exp(-1.0 * m_theta * u) – 1.0) * (exp(-1.0 * m_theta * v) – 1.0)/ (exp(-1 * m_theta) – 1.0))); } virtual double condi_cdf(double u, double v) override { double enut = exp(-u*m_theta); double envt = exp(-v*m_theta); double ent = exp(-m_theta); double denominator = ((ent – 1.0) + (enut – 1.0) * (envt – 1.0)); if(denominator) return (envt * (enut – 1.0)/ (denominator)); else return EMPTY_VALUE; } virtual vector pdf(vector &u,vector &v) override { vector eut,evt,pd; double et = exp(m_theta); eut = exp(u*m_theta); evt = exp(v*m_theta); pd = et*eut*evt*(et – 1.0)*m_theta/pow(et+eut*evt-et*eut-et*evt,2.0); return pd; } virtual vector cdf(vector &u, vector &v) override { return (-1.0 / m_theta * log(1.0 + (exp(-1.0 * m_theta * u) – 1.0) * (exp(-1.0 * m_theta * v) – 1.0)/ (exp(-1 * m_theta) – 1.0))); } virtual vector condi_cdf(vector &u, vector &v) override { vector enut = exp(-1.0*u*m_theta); vector envt = exp(-1.0*v*m_theta); double ent = exp(-m_theta); return (envt * (enut – 1.0)/ ((ent – 1.0) + (enut – 1.0) * (envt – 1.0))); } public: CFrank(void) { m_threshold = 1.e-10; m_copula_type = FRANK_COPULA; m_invalid_theta = 0.0; } ~CFrank(void) { } virtual matrix Sample(ulong num_samples) override { double unf_v[],unf_c[]; vector v,c,u; if(!MathRandomUniform(0.0,1.0,int(num_samples),unf_v) || !MathRandomUniform(0.0,1.0,int(num_samples),unf_c) || !v.Assign(unf_v) || !c.Assign(unf_c)) { Print(__FUNCTION__, " failed to get uniform random numbers ", GetLastError()); return matrix::Zeros(0,0); } matrix out(num_samples,2); for(ulong irow = 0; irow
