# Relativistic Ermakov-Milne-Pinney Systems and First Integrals

Fernando Haas  Physics Institute, Federal University of Rio Grande do Sul, Av.
Bento Gonçalves 9500, 91501-970 Porto Alegre, RS, Brazil

###### Abstract

The Ermakov-Milne-Pinney equation is ubiquitous in many areas of physics that
have an explicit time-dependence, including quantum systems with time
dependent Hamiltonian, cosmology, time-dependent harmonic oscillators,
accelerator dynamics, etc. The Eliezer and Gray physical interpretation of the
Ermakov-Lewis invariant is applied as a guiding principle for the derivation
of the special relativistic analog of the Ermakov-Milne-Pinney equation and
associated first integral. The special relativistic extension of the Ray-Reid
system and invariant is obtained. General properties of the relativistic
Ermakov-Milne-Pinney are analyzed. The conservative case of the relativistic
Ermakov-Milne-Pinney equation is described in terms of a pseudo-potential,
reducing the problem to an effective Newtonian form. The non-relativistic
limit is considered as well. A relativistic nonlinear superposition law for
relativistic Ermakov systems is identified. The generalized Ermakov-Milne-
Pinney equation has additional nonlinearities, due to the relativistic
effects.

###### pacs:

0.30.Hq; 03.30.+p; 05.45.-a; 45.20.-d

##  I Introduction

The Ermakov-Milne-Pinney (EMP) equation Ermakov ; Milne ; Pinney

$\ddot{x}+\kappa^{2}(t)x=C/x^{3}\,,$ (1)

where $C$ is a real constant usually taken as positive, is an ubiquitous
nonlinear non-autonomous ordinary differential equation with many
applications, in particular in problems related to the time-dependent harmonic
oscillator or in connection with exact solutions of the one-dimensional time-
independent Schrödinger equation. In more generality, applications of the EMP
equation appear in cosmological models Hawkins ; Rosu ; Riemann , Bose-
Einstein condensates Haas1 ; Herring , photonic lattices Rod , accelerator
dynamics Courant ; Qin , gravitational wave propagation Ger , higher order
spin models Fring , quantum plasmas EPL , limit cycles Llibre , dynamical
symmetries Haas3 , magneto-gasdynamics Rogers , time-dependent non-Hermitian
quantum system Dey ; Fring2 , supersymmetric systems Cen , noncommutative
quantum mechanics Dey2 , etc. Historical notes can be found e.g. in Cari ;
Morris . We note that the name of the EMP equation is not yet a consensus in
the literature. For instance, sometimes it can be referred just as Pinney
equation.

Diverse generalizations of the EMP equation have been proposed, as for
instance allowing $C(t)$ to have a dependence on time Herring ; Morris , the
inclusion of dissipation Haas2 ; Guha , unbalanced systems of EMP equations
with different frequency functions Ath , modified nonlinearities Reid ,
stochastic differential equations with additive noise Cer . Following the
generalization trend, it would be relevant to extend the EMP equation into the
special relativity domain, with potential applications for the dynamics of
charged particles in high energy density fields Hart . To date only few
results weakly related to special relativistic EMP equations are available, in
connection with the Dirac equation Thi  and the relativistic isotropic two-
dimensional time-dependent harmonic oscillator Petrov . Naturally, the
relativistic extension is not straightforward due to the intrinsic extra
nonlinearity imposed by the presence of the Lorentz factor. Moreover, the
intimate connection between the EMP equation and the time-dependent harmonic
oscillator raises the question of what defines the relativistic harmonic
oscillator. In this regard, the relativistic motion immediately induces
anharmonicities even in the case of a quadratic potential with a frequency-
amplitude dependence Guerrero ; Belendez . Here we follow the approach of many
authors Edery ; Poszwa ; Znojil ; MK ; Gold ; Harvey ; They ; Char2 ; Llibre2
; Li , adopting the spinless Salpeter equation Sal1  with a quadratic external
potential as our definition of relativistic harmonic oscillator. The choice is
due to its simplicity, where just the Newtonian kinetic energy is replaced by
its special relativistic counterpart. This version of the relativistic
harmonic oscillator model has been recently experimentally probed Fuji .

The present work proposes a systematic approach towards the relativistic EMP
equation and beyond, using as a simple guiding principle the Eliezer and Gray
Eliezer  physical interpretation of the Ermakov-Lewis associated first
integral for the isotropic two-dimensional time-dependent harmonic oscillator.
The reason is that the Eliezer and Gray method based on the conservation of
the angular momentum of an auxiliary planar motion provides the standard
physical reasoning behind the conservation of the Ermakov invariant RayReid ;
Mancas . Moreover, the extension of the Ray-Reid Ermakov systems RayReid  to
the relativistic domain will be also obtained, in a certain sense to be
explained. In brief, a relativistic extension of the celebrated EMP equation
and extensions is proposed (which obviously has no direct relation to special
solutions of Klein-Gordon, Dirac or similar relativistic partial differential
equations).

This work is organized as follows. In Section II, we review the Eliezer and
Gray interpretation. In Section III, we consider the relativistic isotropic
two-dimensional time-dependent harmonic oscillator and follow the Eliezer and
Gray method, in order to identify a relativistic EMP equation together with
the associated conservation law. From the structure of the dynamical
equations, the appropriate relativistic Ray-Reid system will be also derived,
as well as with the corresponding first integral. Section IV provides an
alternative derivation based on a dynamical rescaling of time parameter.
Sections V and VI deal with basic properties of the relativistic EMP equation,
mainly in the autonomous case. Section VII is dedicated to a nonlinear
superposition law relating the solutions of the relativistic EMP system.
Finally, Section VIII is reserved to the conclusions and extra remarks.

##  II The Eliezer and Gray physical interpretation

We briefly reproduce the Eliezer and Gray physical interpretation Eliezer  of
the Ermakov-Lewis invariant, which will serve us as a guidance for a
relativistic generalization of the EMP system. Consider the Lagrangian

$L=\frac{1}{2}(\dot{x}^{2}+\dot{y}^{2})-V(x,y,t)\,,\quad
V(x,y,t)=\frac{\kappa^{2}(t)(x^{2}+y^{2})}{2}$ (2)

the corresponding auxiliary plane motion

$\ddot{\bf r}+\kappa^{2}(t){\bf r}=0\,,$ (3)

where ${\bf r}$ has the Cartesian components $(x,y)$, and the EMP equation

$\ddot{\rho}+\kappa^{2}(t)\rho=\frac{J^{2}}{\rho^{3}}\,,$ (4)

where $J$ is a real constant. In terms of polar coordinates $(\rho,\theta)$,
where $\rho=|{\bf r}|,x=\rho\cos\theta,y=\rho\sin\theta$, the equations of
motion become

$\displaystyle\ddot{\rho}-\rho\dot{\theta}^{2}+\kappa^{2}\rho$
$\displaystyle=$ $\displaystyle 0\,,$ (5)
$\displaystyle\frac{1}{\rho}\frac{d}{dt}(\rho^{2}\dot{\theta})$
$\displaystyle=$ $\displaystyle 0\,.$ (6)

Equation (6) implies the constancy of the angular momentum

$J=\rho^{2}\dot{\theta}={\rm const.}$ (7)

assuming for simplicity an unit mass. Combining Eqs. (5) and (7) we derive Eq.
(4). Considering the equation for the $x-$component of the auxiliary motion
together with Eq. (4), the constancy of the Ermakov-Lewis invariant $I$ given
by

$I=\frac{1}{2}\left[(\rho\dot{x}-x\dot{\rho})^{2}+\frac{J^{2}x^{2}}{\rho^{2}}\right]$
(8)

is directly verified, $dI/dt=0$ along trajectories. Expressing in terms of
polar coordinates, one has

$I=\frac{1}{2}(J^{2}\sin^{2}\theta+J^{2}\cos^{2}\theta)=\frac{J^{2}}{2}\,.$
(9)

Therefore, the invariance of $I$ is equivalent to the invariance of the
angular momentum of the auxiliary plane motion.

For reference, it is worth to consider the Ray-Reid (RR) generalization
RayReid  of the EMP system and invariant, namely

$\displaystyle\ddot{x}+\kappa^{2}x$ $\displaystyle=$
$\displaystyle\frac{f(y/x)}{yx^{2}}\,,$ (10)
$\displaystyle\ddot{y}+\kappa^{2}y$ $\displaystyle=$
$\displaystyle\frac{g(x/y)}{xy^{2}}\,,$ (11)

where $f,g$ are arbitrary functions of the indicated arguments. The RR first
integral for Eqs. (10) and (11) is

$I_{RR}=\frac{1}{2}(x\dot{y}-y\dot{x})^{2}+\int^{y/x}f(s)ds+\int^{x/y}g(s)ds\,.$
(12)

One can directly verify that $dI_{RR}/dt=0$ along trajectories.

##  III A relativistic Ermakov-Milne-Pinney system

Proceeding in strict analogy with the NR case, consider the equations of
motion for the 2D relativistic unit rest mass time-dependent harmonic
oscillator, which can be derived Goldstein ; Hand  from the Lagrangian

$L=-\frac{c^{2}}{\gamma}-V(x,y,t)\,,\quad
V(x,y,t)=\frac{\kappa^{2}(t)(x^{2}+y^{2})}{2}\,,$ (13)

where $c$ is the speed of light and
$\gamma=[1-(\dot{x}^{2}+\dot{y}^{2})/c^{2}]^{-1/2}$. The Euler-Lagrange
equations are

$\displaystyle\left(1-\frac{\dot{y}^{2}}{c^{2}}\right)\ddot{x}+\frac{\dot{x}\dot{y}}{c^{2}}\ddot{y}$
$\displaystyle=$ $\displaystyle-\frac{\kappa^{2}x}{\gamma^{3}}\,,$ (14)
$\displaystyle\left(1-\frac{\dot{x}^{2}}{c^{2}}\right)\ddot{y}+\frac{\dot{x}\dot{y}}{c^{2}}\ddot{x}$
$\displaystyle=$ $\displaystyle-\frac{\kappa^{2}y}{\gamma^{3}}\,,$ (15)

which can be disentangled as

$\displaystyle\ddot{x}+\frac{\kappa^{2}}{\gamma}\frac{x}{\gamma_{x}^{2}}$
$\displaystyle=$
$\displaystyle\frac{\kappa^{2}}{\gamma}\frac{\dot{x}\dot{y}}{c^{2}}y\,,$ (16)
$\displaystyle\ddot{y}+\frac{\kappa^{2}}{\gamma}\frac{y}{\gamma_{y}^{2}}$
$\displaystyle=$
$\displaystyle\frac{\kappa^{2}}{\gamma}\frac{\dot{x}\dot{y}}{c^{2}}x\,,$ (17)

where
$\gamma_{x}=(1-\dot{x}^{2}/c^{2})^{-1/2},\,\gamma_{y}=(1-\dot{y}^{2}/c^{2})^{-1/2}$.
The equations of motion are in agreement with Petrov . Similarly to the NR
case, employing cylindrical coordinates $x=\rho\cos\theta,y=\rho\sin\theta$,
we have the conserved angular momentum expressed as

$J=\frac{\partial L}{\partial\dot{\theta}}=\gamma\rho^{2}\dot{\theta}\,.$ (18)

Moreover, we get the Lorentz factor expressible in terms of $\rho,\dot{\rho}$
as

$\gamma=\left(\frac{1+J^{2}/c^{2}\rho^{2}}{1-\dot{\rho}^{2}/c^{2}}\right)^{1/2}\,.$
(19)

The NR formal limit $c\rightarrow\infty$ yields $\gamma=1$ also from Eq. (19).
Using the angular momentum to eliminate the angular velocity, we obtain

$\displaystyle\ddot{x}+\frac{\kappa^{2}}{\gamma}\left(x-\frac{\rho\dot{\rho}\dot{x}}{c^{2}}\right)$
$\displaystyle=$ $\displaystyle 0\,,$ (20)
$\displaystyle\ddot{\rho}+\frac{\kappa^{2}}{\gamma}\left(1-\frac{\dot{\rho}^{2}}{c^{2}}\right)\rho$
$\displaystyle=$ $\displaystyle\frac{J^{2}}{\gamma^{2}\rho^{3}}\,.$ (21)

From the identity
$\gamma^{2}(\rho\dot{x}-x\dot{\rho})^{2}=J^{2}\sin^{2}\theta$, it becomes
self-evident that the quantity

$I_{R}=\frac{1}{2}\left[\gamma^{2}(\rho\dot{x}-x\dot{\rho})^{2}+\frac{J^{2}x^{2}}{\rho^{2}}\right]$
(22)

is a first integral, since

$I_{R}=\frac{1}{2}(J^{2}\sin^{2}\theta+J^{2}\cos^{2}\theta)=\frac{J^{2}}{2}\,,$
(23)

in complete analogy with the Ermakov-Lewis invariant for the NR problem. It
can be also directly verified that $dI_{R}/dt=0$ along the trajectories of the
system (20)-(21). The invariant has the same physical interpretation of the NR
Ermakov invariant in terms of the conservation of the angular motion of the
associated auxiliary 2D motion. In this context, it is justified to interpret
Eqs. (20) and (21) as the (special) relativistic EMP system. Equation (21) is
a relativistic EMP equation (REMP), and the first integral in Eq. (22) is the
relativistic Ermakov-Lewis invariant of the problem. Obviously, in the formal
limit $c\rightarrow\infty$ one recovers the NR case. Understanding the Lorentz
factor in the sense of Eq. (19), the relativistic EMP system is a pair of
nonlinear second-order ordinary differential equations for $x,\rho$. Unlike
the NR case, the equation for $x$ is not uncoupled.

To recapitulate, we have just used the standard Lagrangian for the special
relativistic two-dimensional motion Landau ; Jackson ; Corben ; Arya , in the
particular case where the particle is under the influence of a force linear in
position (a quadratic potential) as seen in a fixed local inertial frame with
respect to which the motion is observed Goldstein ; Hand .

In passing, we write the Hamiltonian ${\cal H}$ associated with (13),

${\cal
H}=c^{2}\left(1+\frac{p_{x}^{2}+p_{y}^{2}}{c^{2}}\right)^{1/2}+V(x,y,t)\,,\quad
p_{x}=\gamma\dot{x}\,,\,\,p_{y}=\gamma\dot{y}\,,$ (24)

which is not a constant of motion in the non-conservative case where
$\kappa(t)$ is a time-dependent function. Even in the autonomous case where
$d{\cal H}/dt=0$, the nature of ${\cal H}$ and the invariant $I_{R}$ in Eq.
(22) is different: the former is energy-like, while the latter is angular
momentum-like.

From the structure of Eqs. (20) and (21) and after some trial and error, it is
possible to identify a relativistic Ray-Reid (RRR) system, namely

$\displaystyle\ddot{x}+\frac{\kappa^{2}}{\gamma}\left(x-\frac{\rho\dot{\rho}\dot{x}}{c^{2}}\right)$
$\displaystyle=$
$\displaystyle\left(1-\frac{\dot{x}^{2}}{c^{2}}\right)\frac{f(y/x)}{\gamma^{2}yx^{2}}-\frac{\dot{x}\dot{y}}{\gamma^{2}c^{2}}\frac{g(x/y)}{xy^{2}}\,,$
(25)
$\displaystyle\ddot{y}+\frac{\kappa^{2}}{\gamma}\left(y-\frac{\rho\dot{\rho}\dot{y}}{c^{2}}\right)$
$\displaystyle=$
$\displaystyle-\,\frac{\dot{x}\dot{y}}{\gamma^{2}c^{2}}\frac{f(y/x)}{yx^{2}}+\left(1-\frac{\dot{y}^{2}}{c^{2}}\right)\frac{g(x/y)}{\gamma^{2}xy^{2}}\,,$
(26)

where $\rho\dot{\rho}=x\dot{x}+y\dot{y}$ and the Lorentz factor is understood
as a function of $\dot{x},\dot{y}$. The invariant for Eqs. (25) and (26) is

$I_{RRR}=\frac{\gamma^{2}}{2}(x\dot{y}-y\dot{x})^{2}+\int^{y/x}f(s)ds+\int^{x/y}g(s)ds\,.$
(27)

It can be verified that $dI_{RRR}/dt=0$ along trajectories. It is apparent
that Eqs. (25), (26) and (27) define a RRR system and its invariant, showing a
complete symmetry between the $x$ and $y$ variables and recovering the RR
system and invariant in the formal NR limit $c\rightarrow\infty$, as shown by
comparison with Eqs. (10)-(12). A derivation provided by a dynamical rescaling
of time will be described in the next Section.

Although our treatment has as motivation the relativistic time-dependent
harmonic oscillator, it is clear that the invariants shown in Eqs. (22) and
(27) do not depend on $\kappa$. Therefore, one is authorized to allow for more
general functional dependencies of $\kappa$, e.g.
$\kappa=\kappa(t,x,y,\dot{x},\dot{y},\dots)$ in Eqs. (25) and (26),
maintaining the constancy of the RRR invariant. Similar remarks apply to the
NR case ReidRay ; HG ; GH . Notice that Eqs. (25) and (26) do not have a
Lagrangian structure in general.

##  IV Derivation from a dynamical rescaling of time

Start from the dynamical system

$x^{\prime\prime}+\omega^{2}x=\frac{f(y/x)}{yx^{2}}\,,\quad
y^{\prime\prime}+\omega^{2}y=\frac{g(x/y)}{xy^{2}}\,,$ (28)

where a prime denotes derivative with respect to the independent variable
$\tau$, or $x^{\prime}=dx/d\tau,y^{\prime}=dy/d\tau$, and where $\omega$ can
be anything in the same spirit of the last Section. In the same way as for the
RR system, it is obvious that (28) possess the invariant

$I=\frac{1}{2}(xy^{\prime}-yx^{\prime})^{2}+\int^{y/x}f(s)ds+\int^{x/y}g(s)ds\,,\quad\frac{dI}{d\tau}=0\,.$
(29)

If we now move to a new independent variable $t$ defined by

$\frac{dt}{d\tau}=\gamma\,,\quad\gamma=\left(1-\frac{\dot{x}^{2}+\dot{y}^{2}}{c^{2}}\right)^{-1/2}\,,\quad\dot{x}=\frac{dx}{dt}\,,\quad\dot{y}=\frac{dy}{dt}\,,$
(30)

then it can be checked after some algebra that the RRR system (25) and (26) is
recovered, with $\kappa^{2}=\omega^{2}/\gamma$, and that the invariant (29)
transforms into (27). The procedure gives a more transparent derivation of the
RRR system from a dynamical rescaling of time starting from the RR system.

It is worthwhile to have a discussion about Lorentz invariance. In the same
sense as the Newtonian harmonic oscillator is obviously not Galilean
invariant, the relativistic harmonic oscillator is not Lorentz invariant. In
this respect we remember that the symmetry group of the NR Ermakov system is
the $SL(2,\Re)$ group Leach , which has no relation with the Galilean group.
Moreover notice that the force equation $dp_{i}/dt=F_{i}$ where $p_{i}$ are
the components of the relativistic momentum is relativistic only “in a certain
sense”, quoting the words from Golsdtein’s book Goldstein , since time has
been keep as entirely distinct from the spatial coordinates. However, a fully
covariant interpretation is possible, regarding (13) as the Lagrangian for a
charged particle under a scalar electrostatic potential which is a quadratic
function of the spatial coordinates, in the instant inertial frame with
respect the motion is observed. In this context, a fully covariant formalism
could be easily performed following the standard approach Goldstein ; Hand ;
Landau ; Jackson ; Corben ; Arya  starting with the Lagrangian
$L=-mc^{2}/\gamma+q{\bf A}\cdot{\bf r}-q\phi$ (the well-known notation is
employed) in the two-spatial dimensions case, with a scalar potential
$\phi\sim(x^{2}+y^{2})$ together with ${\bf A}=0$. This is the Lagrangian (13)
assuming unit rest mass without loss of generality. In addition, the
relativistic harmonic oscillator and REMP equations reproduce the Newtonian
systems in the limit $u/c\rightarrow 0$, where $u$ is a measure of the maximal
velocity of the problem. For instance if $A_{0}$ is the amplitude of the
motion under a linear force $F=-\kappa^{2}x$, the relativistic effects become
negligible provided $\kappa A_{0}/c\rightarrow 0$. The case of a very strong
external field acting on a charged particle [29] is a suitable system to probe
the relativistic effects. In this context since the non-relativistic Ermakov
system has found applications in accelerator physics Qin , the relativistic
version has potential applications for charged particle motions under high-
intensity external fields.

Notice that the derivation of the equations (20) and (21) containing the REMP
equation, together with the relativistic Ermakov-Lewis invariant (22), follows
a very different route compared to the RRR system (25)-(26) and associated
invariant (27). To obtain the REMP equation, the starting point was the
relativistic Lagrangian for the 2D time-dependent harmonic oscillator,
together with polar coordinates, exactly the same as the procedure for the NR
EMP system but then with the special relativistic Lagrangian. On the other
hand, the RRR system can be considered “relativistic” only in the ad hoc sense
that it reduces to the NR Ray-Reid system in the limit $c\rightarrow\infty$
where $c$ is some reference velocity, since its derivation does not starts
from any physical relativistic Lagrangian. In the same context, the use of a
dynamical rescaling of time clearly has not a relativistic spirit. First of
all, the $\gamma$ factor in Eq. (30) is defined in terms of the particle
velocity, not the relative velocity between two inertial frames. Second, it is
obviously not a Lorentz transformation since the space variables are kept the
same. Third, we had not the objective of setting a sort of symmetry
transformation for the equations (28). Rather, the goal was the use of a non-
invariance transformation towards the derivation of the RRR system. Although
at the moment there is the lacking of a a physical interpretation of the RRR
system, which as discussed is relativistic only in a formal sense, we think it
is worthwhile to present it since it has a striking analogy with the
traditional RR system, being a non-trivial pair of coupled nonlinear second-
order ordinary differential equations possessing an exact invariant.

##  V The $J=0$ case

For a vanishing angular momentum, one has $\dot{\theta}=0$ so that it can be
chosen $\theta=0$, without loss of generality. In this case, $y=0$ and Eq.
(20) becomes

$\frac{d}{dt}(\gamma\dot{x})=-\kappa^{2}x\,,\quad\gamma=\left(1-\frac{\dot{x}^{2}}{c^{2}}\right)^{-1/2}\,,$
(31)

which is Bouquet  the equation for an one-dimensional (1D) relativistic time-
dependent harmonic oscillator.

In the case of a constant frequency $\kappa$, obviously the energy is
conserved. In terms of rescaled variables $\bar{x}=\kappa x/c,\,\bar{t}=\kappa
t,\,\bar{v}=d\bar{x}/d\bar{t}$, the corresponding first integral $H_{1D}\geq
1$ is

$H_{1D}=(1-\bar{v}^{2})^{-1/2}+\frac{\bar{x}^{2}}{2}\,,$ (32)

with phase-space contour plots shown in Fig. 1. Clearly, only bounded
trajectories are admissible, as expected. The return points are located at
$\bar{x}=\pm\sqrt{2}(H_{1D}-1)^{1/2}$. In rescaled coordinates, the NR limit
$H_{1D}\simeq 1$ corresponds to circular trajectories in phase space,
$\bar{v}^{2}+\bar{x}^{2}\simeq 2(H_{1D}-1)$. For larger values of $H_{1D}$,
the return points are also larger, implying more relativistic effects and
increasing anharmonicity, as seen in Fig. 1.

Figure 1: Phase-space contour plots of the energy first integral (32) for the
1D conservative relativistic harmonic oscillator described by Eq. (31) with
constant $\kappa$, for $\bar{x}=\kappa x/c,\,\bar{v}=\dot{x}/c$ and different
values of $H_{1D}$, as indicated.

It is instructive to rewrite the conservation law in a Newtonian form,

$\frac{\bar{v}^{2}}{2}+V_{1D}(\bar{x})=0\,,$ (33)

where the pseudo-potential $V_{1D}(\bar{x})$ is defined by

$V_{1D}(\bar{x})=-\,\frac{1}{2}+\frac{1}{2(H_{1D}-\bar{x}^{2}/2)^{2}}\,,$ (34)

shown in Fig. 2. In the NR limit, the variable $\bar{x}$ is limited to small
values, so that $V_{1D}={\rm cte.}+\bar{x}^{2}/(2H_{1D}^{3})+\dots$ with
$H_{1D}\simeq 1$, while larger values of $H_{1D}$ correspond to enhanced
relativistic effects and anharmonicity. The quantity $V_{1D}$ provides a kind
of transfer of the nonlinearity from the kinetic energy term to an anharmonic
(pseudo) potential term, but seemingly has not a physical interpretation.

Figure 2: Pseudo-potential $V_{1D}(\bar{x})$ from Eq. (34) and different
values of $H_{1D}$, as indicated.

The conserved energy can be used for the quadrature of the motion in terms of
elliptic functions. For this purpose, we set

$H_{1D}=1+\frac{A^{2}}{2}\,,\quad-A<\bar{x}<A\,,$ (35)

where $A$ is the amplitude of the motion, and define

$\bar{x}=A\sin\phi\,,\quad\phi=\phi(\bar{t})\,.$ (36)

After some algebra, Eq. (33) is rewritten as

$\left(\frac{d\phi}{d\bar{t}}\right)^{2}=\frac{1+\frac{A^{2}}{4}\cos^{2}\phi}{(1+\frac{A^{2}}{2}\cos^{2}\phi)^{2}}\,,$
(37)

which can be immediately integrated yielding

$\sqrt{4+A^{2}}\,E\left(\phi,\frac{A}{\sqrt{4+A^{2}}}\right)-\frac{2}{\sqrt{4+A^{2}}}\,F\left(\phi,\frac{A}{\sqrt{4+A^{2}}}\right)=\bar{t}-\bar{t_{0}}\,,$
(38)

where $\bar{t}_{0}$ is an integration constant and $F,E$ are incomplete
elliptic integrals of the first and second kind, respectively defined MilneT
according to

$F(\phi,k)=\int_{0}^{\phi}\frac{d\phi^{\prime}}{(1-k^{2}\sin^{2}\phi^{\prime})^{1/2}}\,,\quad
E(\phi,k)=\int_{0}^{\phi}(1-k^{2}\sin^{2}\phi^{\prime})^{1/2}d\phi^{\prime}\,.$
(39)

Using MATHEMATICA it is easy to produce a series solution,

$\phi(\bar{t})=\bar{t}-\bar{t}_{0}-\frac{3A^{2}}{32}(2(\bar{t}-\bar{t}_{0})+\sin[2(\bar{t}-\bar{t}_{0})])+{\cal
O}(A^{4})\,,$ (40)

recovering the Newtonian result in the limit of small amplitude (which is also
the non-relativistic limit). Similarly the period $T$ for which $\phi=2\pi$ is

$T=2\pi\left(1+\frac{3A^{2}}{16}\right)+{\cal O}(A^{4})\,.$ (41)

Approximate periodic solutions for the 1D conservative relativistic harmonic
oscillator can also be found in Mickens .

##  VI Relativistic conservative Ermakov-Milne-Pinney equation

The REMP defined in Eq. (21) with $J\neq 0$ does not have collapsing
($\rho\rightarrow 0$) solutions due to the inverse cubic term. For instance,
suppose $\kappa={\rm cte.}$ and the rescaling
$\bar{\rho}=\kappa\rho/c,\bar{t}=\kappa
t,\bar{v}=d\bar{\rho}/d\bar{t},\bar{J}=\kappa J/c^{2}$, so that

$\frac{d\bar{v}}{d\bar{t}}+\frac{1}{\gamma}(1-\bar{v}^{2})\bar{\rho}=\frac{\bar{J}^{2}}{\gamma^{2}\bar{\rho}^{3}}\,.$
(42)

The energy first integral is

$H=\gamma+\frac{\bar{\rho}^{2}}{2}\geq
1\,,\quad\gamma=\left(\frac{1+\bar{J}^{2}/\bar{\rho}^{2}}{1-\bar{v}^{2}}\right)^{1/2}\,.$
(43)

It is immediate to conclude that $0<\bar{\rho}<\sqrt{2H}$.

Rewriting in a Newtonian form yields

$\bar{v}^{2}/2+V(\bar{\rho})=0\,,$ (44)

with a pseudo-potential defined by

$V(\bar{\rho})=-\frac{1}{2}+\frac{1+\bar{J}^{2}/\bar{\rho}^{2}}{2(H-\bar{\rho}^{2}/2)^{2}}\,,$
(45)

shown in Fig. 3. It is possible to show that $V(\bar{\rho}_{*})<0$, where
$0<\bar{\rho}_{*}<2H$ is the equilibrium point of $V$, viz.

$\bar{\rho}_{*}=\frac{1}{2}\left(-3\bar{J}^{2}+\sqrt{9\bar{J}^{4}+16\bar{J}^{2}H}\right)^{1/2}\,.$
(46)

Figure 3: Pseudo-potential $V$ from Eq. (45).

The return points correspond to $V=0$, similarly to Fig. 2, with stable
oscillations around $\bar{\rho}_{*}$. Elementary algebra shows that the
condition for periodic motions ($V(\bar{\rho}_{*})<0$) is equivalent to

$\bar{J}^{2}<F(H)=\frac{4}{27}\left(-9H+H^{3}+(3+H^{2})^{3/2}\right)\,,$ (47)

or, sufficiently small angular momentum. In the NR limit one has
${\bar{J}}^{2}<(H-1)^{2}$ disregarding ${\cal O}(H-1)^{3}$ terms, so that the
periodicity condition can be shown to be automatically satisfied. The
characteristic function $F$ is shown in Fig. 4. A numerical investigation
shows that the periodicity condition is always meet, for meaningful initial
conditions $0<\bar{\rho}(0)<\sqrt{2H},\,-1<\bar{v}(0)<1$, as expected.

Figure 4: Characteristic function $F(H)$ from Eq. (47) for $H\geq 1$.

The quadrature of Eq. (44) can be made in terms of elliptic functions, after
determining the return points corresponding to the potential well from Eq.
(45). However, the result is exceedingly complicated.

##  VII Nonlinear superposition law

Suppose $\rho=\rho(t)$ a particular solution of the REMP equation (21).
Introducing the new variables

$Q=\frac{x}{\rho}\,,\quad T=\int\frac{dt}{\gamma\rho^{2}}\,,$ (48)

where $\gamma=\gamma(t)$ is given by Eq. (19) converts the relativistic
Ermakov-Lewis invariant (22) into

$I_{R}=\frac{1}{2}\left(\frac{dQ}{dT}\right)^{2}+\frac{J^{2}Q^{2}}{2}\,,$ (49)

formally the same as the energy first integral for a 1D conservative NR
harmonic oscillator. A quadrature yields

$Q=\frac{\sqrt{2I_{R}}}{|J|}\sin(JT+\delta)\,,$ (50)

or

$x=\rho\sin\left(J\int\frac{dt}{\gamma\rho^{2}}+\delta\right)\,,$ (51)

since $I_{R}=J^{2}/2$, adopting $J>0$ and where $\delta$ is a constant phase.
The nonlinear superposition law (51) generalizes the Newtonian result ReidRay
to the relativistic context. In concrete applications, typically the
particular solution $\rho$ should be numerically found.

##  VIII Conclusion

In this work, considerable progress was achieved, in the generalization of
Ermakov systems towards the special relativity domain. The Eliezer and Gray
physical interpretation of the Ermakov-Lewis invariant, was used as a guide
for the derivation of the relativistic analog of the EMP equation, together
with the corresponding first integral for the relativistic planar time-
dependent harmonic oscillator. General aspects of the relativistic EMP
equation have been addressed, and a nonlinear superposition law was derived.
In spite of the successful results, it is still possible to derive other
classes of relativistic Ermakov systems, not arising from the correspondence
with the Eliezer and Gray physical interpretation. For instance, symmetry
principles can be a guiding principle, although the $SL(2,\Re)$ group
structure of non-relativistic Ermakov systems obviously tends to be broken in
the relativistic domain. In the same footing, to carry to the relativistic
domain the linearization of standard Ermakov systems would be a probably
unfeasible task. In future works it is worthwhile to see how one can actually
solve the new relativistic Ermakov systems with an explicit time-dependence.
For this purpose the search for quasi-exact solutions fring5  can be a
fruitful approach.

###### Acknowledgements.

This research was funded by Conselho Nacional de Desenvolvimento Científico e
Tecnológico (CNPq).

## References

  * (1) Ermakov, V.P. Differentzial’nyya uravneniya vtorago poryadka. Usloviya integriruemosti v konechnom vide. Univ. Isz. Kiev Series III 1880, 9, 1–25 (English translation: Harin, A.O. Second-order differential equations. Conditions of complete integrability. Appl. Anal. Discrete Math. 2008, 2, 123–145).
  * (2) Milne, W.E. The numerical determination of characteristic numbers. Phys. Rev. 1930, 35, 863–867.
  * (3) Pinney, E. The nonlinear equation $y^{\prime\prime}(x)+p(x)y(x)+cy^{-3}(x)=0$. Proc. Amer. Math. Soc. 1950, 1, 681.
  * (4) Hawkins, R.M.; Lidsey, J.E. Ermakov-Pinney equation in scalar field cosmologies. Phys. Rev. D 2002, 66, 023523–023530.
  * (5) Rosu, H.C.; Mancas, S.C.; Chen, P. Barotropic FRW Cosmologies with Chiellini damping. Phys. Lett. A 2015, 379, 882–887.
  * (6) Tsamparlis, M.; Paliathanasis, A.; Basilakos, S.; Capozziello, S. Scalar-tensor gravity cosmology: Noether symmetries and analytical solutions. Phys. Rev. D 2014, 89, 063532–063540.
  * (7) Haas, F. Anisotropic Bose-Einstein condensates and completely integrable dynamical systems. Phys. Rev. A 2002, 65, 033603–033608.
  * (8) Herring, G.; Kevrekidis, P.G.; Williams, F.; Christodoulakis, T.; Frantzeskakis, D.F. From Feshbach-resonance managed Bose-Einstein condensates to anisotropic universes: applications of the Ermakov-Pinney equation with time-dependent nonlinearity. Phys. Lett. A 2007, 367, 140–148.
  * (9) Rodríguez-Lara, B.M.; Aleahmad, P.; Moya-Cessa, H.M.; Christodoulides, D.N. Ermakov–Lewis symmetry in photonic lattices. Opt. Lett. 2014, 39, 2083–2085.
  * (10) Courant, E.D.; Snyder, H.D. Theory of the alternating-gradient synchrotron. Ann. Phys. 1958, 3, 1–48.
  * (11) Qin, H.; Davidson, R.C. Symmetries and invariants of the oscillator and envelope equations with time-dependent frequency. Phys. Rev. ST Accel. Beams 2006, 9, 054001–054005.
  * (12) Geralico, A.; Landolfi, G.; Ruggeri, C.; Soliani, G. Novel approach to the study of quantum effects in the early Universe. Phys. Rev. D 2004, 69, 043504–043513.
  * (13) Fring, A.; Frith, T. Metric versus observable operator representation, higher spin models. Eur. Phys. J. Plus 2018, 133, 57–65.
  * (14) Haas, F.; Bret, A. Nonlinear low-frequency collisional quantum Buneman instability. Europhys. Lett. 2012, 97, 26001, 7 pages.
  * (15) Llibre, J.; Pérez-Chavela, E. Limit cycles for a class of second order differential equations. Phys. Lett. A 2011, 375, 1080–1083.
  * (16) Haas, F.; Goedert, J. On the linearization of the generalized Ermakov systems. J. Phys. A: Math. Gen. 1999, 32, 2835–2844.
  * (17) Rogers, C.; Schief, W. K. The pulsrodon in 2+1-dimensional magneto-gasdynamics: Hamiltonian structure and integrability. J. Math. Phys. 2011, 52, 083701–083721.
  * (18) Dey, S.; Fring, A.; Gouba, L. Milne quantization for non-Hermitian systems. J. Phys. A: Math. Theor. 2015, 48, 40FT01, 10 pages.
  * (19) Fring, A.; Frith, T. Mending the broken PT-regime via an explicit time-dependent Dyson map. Phys. Lett. A 2017, 381, 2318–2323.
  * (20) Cen, J.; Fring, A.; Frith, T. Time-dependent Darboux (supersymmetric) transformations for non-Hermitian quantum systems. J. Phys. A: Math. Theor. 2019, 52, 115302, 23 pages.
  * (21) Dey, S.; Fring, A. Noncommutative quantum mechanics in a time-dependent background. Phys. Rev. D 2014, 90, 084005, 14 pages.
  * (22) Cariñena, J.F.; Lucas, J.; Rañada, M.F. A nonlinear superposition rule for solutions of the Milne-Pinney equation. SIGMA 2008, 4, 031, 18 pages.
  * (23) Morris, R.M.; Leach, P.G.L. The effects of symmetry-breaking functions on the Ermakov-Pinney equation. Appl. Anal. Discrete Math. 2017, 11, 62–73.
  * (24) Haas, F. The damped Pinney equation and its applications to dissipative quantum mechanics. Physica Scripta 2010, 81, 025004, 8 pages.
  * (25) Guha, P.; Ghose-Choudhury, A. Nonlocal transformations of the generalized Liénard type equations and dissipative Ermakov-Milne-Pinney systems. Int. J. Geom. Meth. Mod. Phys. 2019, 16, 1950107, 18 pages.
  * (26) Athorne, C. On solving a class of unbalanced Ermakov- Pinney systems. J. Phys. A: Math. Gen. 2001, 34, L563–L566.
  * (27) Reid, J.L. An exact solution of the nonlinear differential equation $\ddot{y}+p(t)y=qm(t)/y^{2m-1}$. Proc. Am. Math. Soc. 1971, 27, 601–602.
  * (28) Cervantes-López, E.; Espinoza, P.B.; Gallegos, A.; Rosu, H.C. Ermakov–Ray–Reid systems with additive noise. Physica A 2015, 439, 44–47.
  * (29) Hartemann, F.V. High-field Electrodynamics. CRC Press: Boca Raton, USA, 2002.
  * (30) Thylwe, K.-E.; McCabe, P. Bohr-Sommerfeld quantization condition for Dirac states derived from an Ermakov-type invariant. J. Math. Phys. 2013, 54, 052301–052308.
  * (31) Petrov, S.V. Classical dynamics of the relativistic oscillator. Eur. J. Phys. 2016, 37, 065605, 11 pages.
  * (32) Guerrero, J.; Aldaya, V. A perturbative approach to the relativistic harmonic oscillator. Mod. Phys. Lett. A 1999, 14, 1689–1699.
  * (33) Beléndez, A.; Pascual, C.; Méndez, D.I.; Neipp, C. Solution of the relativistic (an)harmonic oscillator using the harmonic balance method. J. Sound Vib. 2008, 311, 1447–1456.
  * (34) Edery, A.; Laporte, P. First and second-order relativistic corrections to the two and higher-dimensional isotropic harmonic oscillator obeying the spinless Salpeter equation. J. Phys. Commun. 2018, 2, 025024, 16 pages.
  * (35) Poszwa, A. Relativistic generalizations of the quantum harmonic oscillator Acta Physica Polonica A 2014, 126, 1226–1234.
  * (36) Znojil, M. Harmonic oscillations in quasi-relativistic regime. J. Phys. A: Math. Gen 1996, 29, 2905–2918.
  * (37) Balasubramanya, M.K. The action variable and frequency of a relativistic harmonic oscillator. Theor. Math. Phys. 2010, 162, 341–346.
  * (38) Gold, L. Note on the relativistic harmonic oscillator. J. Franklin Inst. 1957, 264, 25–27.
  * (39) Harvey, A.L. Relativistic harmonic oscillator. Phys. Rev. D 1972, 6, 1474–1475.
  * (40) Teychenné, D.; Bonnaud, G.; Bobin, J. Oscillatory relativistic motion of a particle in a power-law or sinusoidal-shaped potential well. Phys. Rev. E 1994, 49, 3253–3263.
  * (41) Chargui, Y. On an approximation of the two-body spinless Salpeter equation. Eur. Phys. J. Plus 2018, 133, 543, 6 pages.
  * (42) Llibre, J.; Makhlouf, A. On the periodic solutions of the relativistic driven harmonic oscillator. J. Math. Phys. 2020, 61, 012501–012507.
  * (43) Li, Z.; Liu, J.; Lucha, W.; Ma, W.; Schöberl, F.F. Relativistic harmonic oscillator. J. Math. Phys. 2005, 46, 103514–103525.
  * (44) Salpeter, E.E. Mass corrections to the fine structure of hydrogen-like atoms. Phys. Rev. 1952, 87, 328–343.
  * (45) Fujiwara, K.M.; Geiger, Z.A.; Singh, K.; Senaratne, R.; Rajagopal, S.V.; Lipatov, M.; Shimasaki, T.; Weld, D.M. Experimental realization of a relativistic harmonic oscillator. New J. Phys. 2018, 20, 063027, 7 pages.
  * (46) Eliezer, C.J.; Gray, A. A note on the time-dependent harmonic oscillator. SIAM J. Appl. Math. 1976, 30, 463–468.
  * (47) Ray, J.R.; Reid, J.L. More exact invariants for the time-dependent harmonic oscillator. Phys. Lett. A 1979, 71, 317–318.
  * (48) Mancas, S.C.; Rosu, H.C. Ermakov-Lewis invariants and Reid systems. Phys. Lett. A 2014, 378, 2113–2117.
  * (49) Goldstein, H.; Poole, C.; Safko, J. Classical Mechanics, 3rd ed. Addison Wesley: San Francisco, USA, 2002; pp. 316.
  * (50) Hand, L.N.; Finch, J.D. Analytical Mechanics. Cambridge University Press: Cambridge, UK, 1998; p. 551.
  * (51) Landau, L.D..; Lifshitz, E.M. The Classical Theory of Fields, 4rd revised English ed. Butterworth-Heinemann: Oxford, UK, 1975; pp. 48.
  * (52) Jackson, J.D. Classical Electrodynamics. John Wiley: New York, USA, 1962; pp. 408.
  * (53) Corben, H.C.; Stehle, P. Classical Mechanics. John Wiley: New York, USA, 1960; pp. 297.
  * (54) Arya, A. Introduction to Classical Mechanics, 2nd ed. Benjamin Cummings: San Francisco, USA, 1997; p. 694.
  * (55) Reid, J.L.; Ray, J.R. Ermakov systems, nonlinear superposition, and solutions of nonlinear equations of motion. J. Math. Phys. 1980, 21, 1583–1587.
  * (56) Haas, F.; Goedert, J. On the Hamiltonian structure of Ermakov systems. J. Phys. A: Math. Gen. 1996, 29, 4083–4092.
  * (57) Goedert, J.; Haas, F. On the Lie symmetries of a class of generalized Ermakov systems. Phys. Lett. A 1998, 239, 348–352.
  * (58) Leach, P.G.L. Generalized Ermakov systems. Phys. Lett. A 1991, 158, 102–106.
  * (59) Martin, J.; Bouquet, S. Search for first integrals in relativistic time‐dependent Hamiltonian systems. J. Math. Phys. 1994, 35, 181–189.
  * (60) Milne-Thomson, L.M. Elliptic integrals, in: Abramowitz, M.; Stegun, I.A. (Eds.), Handbook of Mathematical Functions, Dover Publications: New York, USA, 1972; pp. 589.
  * (61) Mickens, R.E. Periodic solutions of the relativistic harmonic oscillator. J. Sound Vib. 1998, 212, 905–908.
  * (62) Fring, A.; Tenney, R. Time-independent approximations for time-dependent optical potentials. Eur. Phys. J. Plus 2020, 135, 163, 19 pages.


# Nonlinear Granger Causality using Kernel Ridge Regression

Wojciech “Victor” Fulmyk<EMAIL_ADDRESS>
###### Abstract

I introduce a novel algorithm and accompanying Python library, named
mlcausality, designed for the identification of nonlinear Granger causal
relationships. This novel algorithm uses a flexible plug-in architecture that
enables researchers to employ any nonlinear regressor as the base prediction
model. Subsequently, I conduct a comprehensive performance analysis of
mlcausality when the prediction regressor is the kernel ridge regressor with
the radial basis function kernel. The results demonstrate that mlcausality
employing kernel ridge regression achieves competitive AUC scores across a
diverse set of simulated data. Furthermore, mlcausality with kernel ridge
regression yields more finely calibrated $p$-values in comparison to rival
algorithms. This enhancement enables mlcausality to attain superior accuracy
scores when using intuitive $p$-value-based thresholding criteria. Finally,
mlcausality with kernel ridge regression exhibits significantly reduced
computation times compared to existing nonlinear Granger causality algorithms.
In fact, in numerous instances, this innovative approach achieves superior
solutions within computational timeframes that are an order of magnitude
shorter than those required by competing algorithms.

Keywords: machine learning, Granger causality, time-series, nonlinearity,
kernel ridge regression, mlcausality, synthetic network, causal discovery,
nonparametric methods.

##  1 Introduction

Identifying the direction of relationships amongst simultaneously observed
time-series continues to generate significant interest among researchers from
many different fields. Since Nobel prize-winning economist Clive Granger first
introduced the concept of Granger causality (Granger, 1969), significant
theoretical contributions towards the development of Granger-causal methods
have been made by computer scientists, biostatisticians, physicists,
mathematicians, and many others. Much of the recent research has focused on
identifying nonlinear Granger causality, possibly in the presence of other
confounding time-series (Wismüller et al., 2021; Rosoł et al., 2022; Lacasa et
al., 2015; Gao et al., 2017; DSouza et al., 2017). The novel approach proposed
herein, which I call mlcausality, tackles the problem of identifying nonlinear
relationships using a regressor-agnostic non-parametric test. Consequently,
mlcausality has a plug-in architecture that allows the researcher to use any
nonlinear regressor, such as a kernel ridge regressor, a support vector
regressor, a random forest regressor, or a gradient boosting regressor, as the
base prediction model. Due to mlcausality’s plug-in architecture and the large
number of regressors with which it can be used, the rest of this paper focuses
exclusively on analyzing mlcausality’s performance when the kernel ridge
regressor with the radial basis function kernel is used, with an analysis of
the performance of other regressors deferred to subsequent papers. For many
nonlinear networks mlcausality with kernel ridge regression outperforms rival
algorithms in terms of both network recovery performance as well as execution
speed.

##  2 Methodology

###  2.1 An introduction to Granger causality

Fundamentally, Granger causality describes predictability: given two time-
series $X$ and $Y$, if the histories of both $X$ and $Y$ predict future values
of $Y$ better than the history of $Y$ alone, then $X$ is said to Granger cause
$Y$. Before continuing, it is important to note that Granger causality itself
is a misnomer because no truly causal or mechanical link between the time-
series in question is implied. In other words, if $X$ Granger causes $Y$, then
it could be the case that $X$ actually causes $Y$, or it could be the case
that $X$ simply leads $Y$ and therefore provides useful information about the
future value(s) of $Y$ without actually causing $Y$.

Classical Granger causality is formulated as follows. Given two time-series
$X$ and $Y$ and all available information $U$, $X$ Granger causes $Y$ if the
variance of the prediction error when the history of $X$ is included is lower
than the variance of the prediction error when $X$ is excluded:

$\sigma^{2}(Y|U)<\sigma^{2}(Y|(U-X))$ (1)

where $(U-X)$ denotes all available information except for time series $X$.

To evaluate equation 1 two linear models on lagged variables are estimated, an
“unrestricted” model that includes the lags of $X$ and a “restricted” model
that excludes the lags of $X$. Although it is theoretically possible to
include different amounts of lags for each time-series, in practice an equal
amount of lags for all time-series is typically used. Given an equal amount of
lags for all variables, the classical Granger causality model evaluates:

$\displaystyle Y(t)=\alpha_{r,0}+\sum_{l=1}^{L}{\alpha_{r,l}Y(t-l)}+E_{r}$ (2)
$\displaystyle
Y(t)=\alpha_{u,0}+\sum_{l=1}^{L}{\alpha_{u,l}Y(t-l)}+\sum_{l=1}^{L}{\beta_{u,l}X(t-l)}+E_{u}$
(3)

where $L$ is the total number of lags; $\alpha_{m,0}$ are the intercepts for
the restricted $r$ and unrestricted $u$ models $m$; $\alpha_{m,l}$ are the
coefficients for lags $l$ of $Y$; $\beta_{m,l}$ are the coefficients for lags
$l$ of $X$, and $E_{m}$ are the error terms.

Note that, under classical Granger causality, the restricted model is a nested
model of the unrestricted model. For linear regression this guarantees that
the fit of the unrestricted model will be no worse than that of the restricted
model. In order to test whether the variance of the unrestricted model is
significantly lower, in the statistical sense, than the variance of the
restricted model, an F-test is performed:

$F=\frac{RSS_{r}-RSS_{u}/L}{RSS_{u}/(N-2L-1)}$ (4)

where $N$ is the total number of training samples for which all lag terms
exist, $RSS_{m}$ is the residual sum of squares for model $m$, and $L$ is the
total number of lags. Under the null hypothesis, if the fit from the
unrestricted model is not significantly better than that of the restricted
model, the $F$-statistic from equation 4 follows an $F$-distribution with
$(L,N-2L-1)$ degrees of freedom.

The primary limitation of the classical Granger causality model is that, by
construction, the classical model is only able to identify Granger causality
if $X$ and $Y$ are linearly related. Identifying nonlinear relationships using
classical Granger causality would require performing nonlinear data
transformations, such as by taking natural logs of the data. However, making
such explicit data transformations is not always a trivial task, especially in
cases where the data generating process is not well known to the researcher.

On a positive note, the addition of superfluous confounding time-series into
the classical model is possible as long as the denominator degrees of freedom
remain strictly positive. In other words, if there was an additional
confounding time-series $Z$ that one would want to include in the analysis,
then one could simply add $\sum_{l=1}^{L}{\gamma_{m,l}Z(t-l)}$ to both the
restricted and unrestricted models in equations 2 and 3 and alter the
denominator degrees of freedom in equation 4 to $N-3L-1$. The null hypothesis
with confounding time-series $Z$ is that
$\sigma^{2}(Y|X_{t-1},...,X_{t-L},Y_{t-1},...,Y_{t-L},Z_{t-1},...,Z_{t-L})\nless\sigma^{2}(Y|Y_{t-1},...,Y_{t-L},Z_{t-1},...,Z_{t-L})$.
Additional confounding time-series can be added analogously.

###  2.2 Kernel Ridge Regression

The following provides a brief review of kernel ridge regression. For a more
thorough treatment, the reader is encouraged to review Murphy (2012); Vovk
(2013); and Exterkate et al. (2016).

Let $\bm{x}\in\mathbb{R}^{D}$ represent a feature vector for $D$ features, and
$\bm{X}\in\mathbb{R}^{N\times D}$ represent a design matrix for $N$
observations. Given a target vector $\bm{y}\in\mathbb{R}^{N}$, ridge
regression minimizes the following objective function:

$\displaystyle\min_{w}\
||\bm{y}-\bm{X}\bm{w}||_{2}^{2}+\lambda||\bm{w}||_{2}^{2}$ (5)
$\displaystyle=$ $\displaystyle\min_{w}\
(\bm{y}-\bm{X}\bm{w})^{T}(\bm{y}-\bm{X}\bm{w})+\lambda\bm{w}^{T}\bm{w}$

where $\bm{w}\in\mathbb{R}^{D}$ are the regression coefficients and $\lambda$
is a penalty term that penalizes the magnitudes of those coefficients. Note
that $\lambda$ acts as a regularization term: if $\lambda=0$, then objective
function 5 collapses to a linear regression, while values of $\lambda$ greater
than zero encourage smaller coefficient magnitudes at the optimum.

The solution is found by taking the derivative of the objective function and
setting that derivative equal to zero:

$\displaystyle\bm{w}=(\bm{X}^{T}\bm{X}+\lambda\bm{I}_{D})^{-1}\bm{X}^{T}\bm{y}$
(6)
$\displaystyle\bm{w}=\bm{X}^{T}(\bm{X}\bm{X}^{T}+\lambda\bm{I}_{N})^{-1}\bm{y}$
(7)

where equation 7 is due to the Sherman-Morrison-Woodbury formula (Woodbury,
1950).

The ridge regression model, as described above, produces nothing more than a
linear solution that penalizes the magnitudes of the coefficients in relation
to $\lambda$ and is therefore incapable of efficiently modelling nonlinear
relationships by itself. In order to introduce nonlinearity, the ridge
regression model is kernalized. Intuitively, kernel methods map feature
vectors into a higher, possibly infinitely dimensional space. The goal of
kernel ridge regression is to then apply (inherently linear) ridge regression
on the mappings of the feature vectors in that higher dimensional space.
However, from a computational perspective, it turns out that one does not need
to actually map feature vectors into a different space at all; rather, it is
sufficient to just calculate relationships between all feature vectors using a
kernel function. This process is commonly referred to as the “kernel trick.”

Formally, define a symmetric kernel function
$\kappa(\bm{x},\bm{x}^{\prime})=\kappa(\bm{x}^{\prime},\bm{x})\in\mathbb{R}$
as a real-valued function for feature vectors $\bm{x},\bm{x}^{\prime}\in\chi$,
with $|\chi|=N$. If the Gram matrix $\bm{K}$ with elements
$\kappa(\bm{x}_{i},\bm{x}_{j})_{ij}\ \forall\ i,j\in\\{1,...,N\\}$ is positive
definite, then by Mercer’s theorem (Mercer, 1909) there exists a function
$\phi$ mapping $\bm{x}$ to $\mathbb{R}^{M}$ such that

$\kappa(\bm{x},\bm{x}^{\prime})=\phi(\bm{x})^{T}\phi(\bm{x}^{\prime})$ (8)

Moreover, Mercer’s theorem guarantees that $M$, the dimension of the space to
which $\phi$ maps feature vectors into, can be arbitrarily large, potentially
infinite. A kernel that satisfies the above Gram matrix condition is called a
Mercer kernel.

Ridge regression can be kernalized as follows. Suppose that one wishes to
evaluate ridge regression not on the feature vectors $\bm{x}$ directly but
rather on some mapping $\phi(\bm{x})$ into a higher, possibly infinitely
dimensional, space. Moreover, let $\bm{\phi}(\cdot)$ represent the multi-
observational analogue of $\phi(\cdot)$ that admits $\bm{X}$, the design
matrix that stores feature vectors $\bm{x}$ for multiple observations, as an
input. Then ridge regression could be partially kernalized by replacing
$\bm{X}\bm{X}^{T}$ in equation 7 with the Gram matrix $\bm{K}$ for some Mercer
kernel $\kappa(\bm{x},\bm{x}^{\prime})=\phi(\bm{x})^{T}\phi(\bm{x}^{\prime})$,
and by replacing $\bm{X}^{T}$ with $\bm{\phi}(\bm{X})^{T}$. The partially
kernalized solution then becomes:

$\displaystyle\bm{w}=\bm{\phi}(\bm{X})^{T}(\bm{K}+\lambda\bm{I}_{N})^{-1}\bm{y}$
(9)

Now, let $\bm{\alpha}=(\bm{K}+\lambda\bm{I}_{N})^{-1}\bm{y}$. Equation 9 then
becomes:

$\displaystyle\bm{w}$ $\displaystyle=\bm{\phi}(\bm{X})^{T}\bm{\alpha}$ (10)
$\displaystyle=\sum_{i=1}^{N}\alpha_{i}\phi(\bm{x}_{i})$ (11)

Note that equation 10 is not yet usable because the kernalization is only
partial: we still have to deal with $\bm{\phi}(\bm{X})^{T}$, which is
potentially difficult to evaluate because all feature vectors in $\bm{X}$
would have to be transformed using $\phi(\cdot)$ directly. However, if just
the in-sample predictions from the kernel ridge regression are needed, then
for each feature vector $\bm{x}$, one just has to evaluate

$\displaystyle\hat{f}(\bm{x})$ $\displaystyle=\bm{w}^{T}\phi(\bm{x})$ (12)
$\displaystyle=\sum_{i=1}^{N}\alpha_{i}\phi(\bm{x}_{i})^{T}\phi(\bm{x})$ (13)
$\displaystyle=\sum_{i=1}^{N}\alpha_{i}\kappa(\bm{x}_{i},\bm{x})$ (14)

where the last line completes the kernalization. It is thus possible to
generate nonlinear in-sample predictions by fitting an (inherently linear)
ridge regression on a mapping of the feature vectors into a higher dimensional
space, without actually performing the mapping into that higher dimensional
space. In order to perform out-of-sample predictions, one would first fit
equation 14 using the training data and then use the obtained in-sample
predictions to precompute $\bm{\alpha}$; out-of-sample predictions could then
be made directly using equation 14 with that precomputed $\bm{\alpha}$ and a
previously unseen feature vector $\bm{x}$ that comes from the test set.

Using kernel ridge regression in practice involves choosing a suitable kernel
function and, by design, the mlcausality Python library does not provide any
restrictions with regards to this: the user is free to choose any kernel
function supported by the scikit-learn Python library. However, all of the
subsequent analysis presented herein was performed using the radial basis
function (RBF) kernel, also known as the Gaussian kernel:

$\displaystyle\kappa(\bm{x},\bm{x}^{\prime})=e^{-\gamma||\bm{x}-\bm{x}^{\prime}||_{2}^{2}}$
(15)

where $\gamma\in(0,\infty)$ is a hyperparameter. The RBF kernel is infinitely
dimensional, as can be seen in the multinomial expansion in Shashua (2009).

###  2.3 The mlcausality Algorithm

Suppose there is a set $\vmathbb{X}$ of equally-spaced time-series that share
a common time index and do not contain any gaps. Suppose that each of these
time series is represented by a column vector $\vmathbb{x}\in\vmathbb{X}$ and
that each of these vectors has length $N+L$. Suppose that $\vmathbb{X}$
contains at least 2 such time series, that is, that $|\vmathbb{X}|=G\geq 2$.
Let $L$ represent the number of lags. Moreover, let
$\bm{X}_{u}\in\mathbb{R}^{N\times LG}$ be a matrix of $L$ lags $\forall$ $G$
time-series in $\vmathbb{X}$, as depicted in the matrix below:

$\mathbf{\bm{X}_{u}}=\begin{bmatrix}\vmathbb{x}_{1,1\phantom{-1}}&\cdots&\vmathbb{x}_{1,L\phantom{-1}\phantom{-1}}&\cdots\cdots&\vmathbb{x}_{G,1\phantom{-1}}&\cdots&\vmathbb{x}_{G,L\phantom{-1}\phantom{-1}}\\\
\vmathbb{x}_{1,2\phantom{-1}}&\cdots&\vmathbb{x}_{1,L+1\phantom{-1}}&\cdots\cdots&\vmathbb{x}_{G,2\phantom{-1}}&\cdots&\vmathbb{x}_{G,L+1\phantom{-1}}\\\
\vmathbb{x}_{1,3\phantom{-1}}&\cdots&\vmathbb{x}_{1,L+2\phantom{-1}}&\cdots\cdots&\vmathbb{x}_{G,3\phantom{-1}}&\cdots&\vmathbb{x}_{G,L+2\phantom{-1}}\\\
\vdots&\cdots&\vdots&\ddots&\vdots&\cdots&\vdots\\\
\vmathbb{x}_{1,N-2}&\cdots&\vmathbb{x}_{1,N+L-3}&\cdots\cdots&\vmathbb{x}_{G,N-2}&\cdots&\vmathbb{x}_{G,N+L-3}\\\
\vmathbb{x}_{1,N-1}&\cdots&\vmathbb{x}_{1,N+L-2}&\cdots\cdots&\vmathbb{x}_{G,N-1}&\cdots&\vmathbb{x}_{G,N+L-2}\\\
\vmathbb{x}_{1,N\phantom{-1}}&\cdots&\vmathbb{x}_{1,N+L-1}&\cdots\cdots&\vmathbb{x}_{G,N\phantom{-1}}&\cdots&\vmathbb{x}_{G,N+L-1}\\\
\end{bmatrix}$

where the first coordinate in the subscript of an element indicates the time-
series and the second coordinate represents the time index. Note that a single
row in the above matrix represents a concatenation of ${1,...,L}$ lags
$\forall$ $G$ time-series in $\vmathbb{X}$, and that only rows that have no
missing lags are kept. Similarly, let $\bm{X}_{r}\in\mathbb{R}^{N\times
L(G-1)}$ be a matrix of $L$ lags $\forall$ $G-1$ time-series in $\vmathbb{X}$
other than time-series $\bm{q}\in\vmathbb{X}$. Then, in order to test whether
time-series $\bm{q}\in\vmathbb{X}$ Granger-causes time-series
$\bm{y}\in\vmathbb{X}$ with $\bm{q}\neq\bm{y}$, mlcausality evaluates

$\displaystyle\hat{\bm{y}}_{r}=f_{r}(\bm{X}_{r},\bm{y}_{true},*_{r})$ (16)
$\displaystyle\hat{\bm{y}}_{u}=f_{u}(\bm{X}_{u},\bm{y}_{true},*_{u})$ (17)

where the subscripts $r$ and $u$ represent the restricted and unrestricted
models, respectively; $f$ represents any nonlinear regressor, such as the
kernel ridge regressor described above; and $*$ is a placeholder for any
additional hyperparameters the chosen regressor accepts.

For kernel ridge regression with the RBF kernel $*_{r}=(\lambda,\gamma_{r})$
and $*_{u}=(\lambda,\gamma_{u})$. Note that, in theory, $\lambda$ can be
varied between the restricted and unrestricted models, however, different
$\lambda$ values for the restricted and unrestricted models would make it
difficult to assess whether the restricted and unrestricted models perform
differently solely on the bases of the exclusion or inclusion of time-series
$\bm{q}$; as such, $\lambda$ will be kept the same for both models. It is not
clear which value of $\lambda$ should be chosen for the nonlinear Granger
causality identification task, although simulation results suggest that the
performance of mlcausality with the kernel ridge regressor and the RBF kernel
is not very sensitive to the choice of $\lambda$. Consequently, all results
presented herein use $\lambda=1$, which is the scikit-learn Python library
default value for the penalty term in kernel ridge regression.

On the other hand, $\gamma$ is subscripted with the model type because it is
natural to use $\gamma=1/D$, where $D$ is the number of features, i.e. the
number of columns in either $\bm{X}_{r}$ or $\bm{X}_{u}$, for the problem at
hand. This is because parameter $\gamma$ in the RBF kernel acts as a
multiplier for the square of the Euclidean norm, where the number of elements
in the sum of that norm is equal to the number of elements in a feature
vector. Using an identical $\gamma$ for both models would not account or
correct for the greater amount of features in the unrestricted model stemming
from its inclusion of time-series $\bm{q}$. As such, all results presented
herein use $\gamma=1/D$ which acts as a weight that scales the squared norm in
the RBF kernel by the number of features in that norm.

####  2.3.1 The Statistical Test

A notable distinction between classical Granger causality and mlcausality lies
in the substitution of linear regression with a nonlinear regressor, such as
the kernel ridge regressor or the support vector regressor. The task now is to
identify an appropriate approach for assessing whether the forecasts generated
by the unrestricted model exhibit superior performance compared to those
produced by the restricted model in the context of the usage of that nonlinear
regressor.

Recall that classical Granger causality relies on an $F$-test to determine if
the variance of the unrestricted model significantly differs from that of the
restricted model. Nevertheless, the incorporation of nonlinear regressors
within the framework of mlcausality can potentially lead to suboptimal
outcomes due to the $F$-test’s sensitivity to normality, or even to an
outright violation of some of the assumptions underpinning the $F$-test.
Furthermore, depending on the specific regressor employed, it may prove
challenging or even infeasible to precisely calculate the required degrees of
freedom for executing the $F$-test. Consequently, instead of resorting to the
$F$-test, mlcausality opts for the sign test (Dixon and Mood, 1946), a non-
parametric test that imposes minimal assumptions.

The sign test compares the counts of positive and negative values and checks
whether they follow a binomial distribution with the probability of success
set to 0.5. The null hypothesis of the sign test is that the median of the
distribution is equal to zero, while the one-sided alternative hypothesis is
that the the median is greater than zero.

The data subjected to the sign test is formulated in the following manner.
First, the absolute values of the errors from the restricted and unrestricted
models are calculated:

$\displaystyle|\bm{E}_{r}|=|\hat{\bm{y}}_{r}-\bm{y}_{true}|$ (18)
$\displaystyle|\bm{E}_{u}|=|\hat{\bm{y}}_{u}-\bm{y}_{true}|$ (19)

where $|\bm{v}|$ indicates, through an abuse of notation quite common in
computer science, the absolute value of every element in vector $\bm{v}$.
Then, a difference between the absolute values of the errors is found:
$\bm{\delta}=|\bm{E}_{r}|-|\bm{E}_{u}|$. The sign test is subsequently
employed on $\bm{\delta}$.

There are several issues related to the usage of the sign test in this context
that must be addressed. First off, note that, although the sign test does not
come with very restrictive assumptions, there is one key assumption that must
be met: that of independence of all elements in $\bm{\delta}$. For the context
presented herein independence can be assured by splitting the original set of
time-series data into distinct train and test sets on the time domain. Once
the split is formed, both the restricted and unrestricted models are trained
on the train data only, with the difference vector $\bm{\delta}$ being
constructed solely from the restricted and unrestricted models’ predictions of
the test data. Note that the train-test split has to account for the time-
series characteristics of the underlying data; in other words, the training
data has to precede the test data, and there has to be a gap of at least $L$
observations between the train and test sets in order to avoid data leakage.
In order to satisfy the above, by default, the mlcausality Python library uses
the first 70% of the observations as the training data and the remaining 30%
of the observations as the testing data, with a gap of $L$ observations
between these two datasets. All results presented in this paper are with
respect to this default data split, although superior results could be
obtained with different splits in some cases. Note that the need to split the
data into distinct train and test sets represents a departure from classical
Granger causality where the train and test sets are identical and equal to all
the data that is available. The main implication of this departure is that, in
order for the mlcausality test to work as intended, the train and test data
have to come from the same distribution.

Secondly, it is important to acknowledge that the sign test allocates
identical weights to all time periods and does not take into account the
extent to which one model outperforms the other at a specific time point. An
alternative to the sign test that assigns more weight to time periods in
$|\bm{\delta}|$ that have greater values is the Wilcoxon signed rank test
(Wilcoxon, 1945). The null hypothesis of the Wilcoxon signed rank test is that
the difference between two related paired samples is symmetric around a value
less than or equal to zero, while the alternative hypothesis is that the
difference is symmetric around a value greater than zero. The Wilcoxon signed
rank test evaluates the following test statistic:

$\displaystyle T=\sum_{i=1}^{N}sgn(\delta_{i})R_{i}$ (20)

where $sgn(\delta_{i})$ is the sign of $\delta_{i}$, and $R_{i}$ is the rank
of the magnitude of $\delta_{i}$ in vector $\bm{\delta}$, with the smallest
magnitude being assigned rank one, the second smallest magnitude being
assigned rank two, and so on.

The Wilcoxon signed rank test would provide greater power than the sign test
if the assumptions of the Wilcoxon signed rank test were fully met; in
particular, if all elements in $\bm{\delta}$ were independent, and if the
distribution of $\bm{\delta}$ was symmetric. The independence assumption can
be satisfied with a suitable train-test split, but the symmetry assumption
required by the Wilcoxon signed rank test is more restrictive than the
minimalist assumptions of the sign test, which only requires the satisfaction
of independence. Consequently, the mlcausality Python library uses the sign
test by default, although there are overrides that allow the usage of the
Wilcoxon signed rank test for those that wish to use it. In practice, at least
for the simulated networks presented herein, the sign test tends to perform
very well, and there is typically no need to depart from the sign test in
favour of the Wilcoxon signed rank test.

####  2.3.2 Data Preprocessing

Prior to feeding the raw data into the mlcausality algorithm, as defined
above, some preprocessing steps may be needed or desired.

First off, as is typical for these types of analysis, one should ensure that
all time-series fed into the mlcausality algorithm are stationary. Moreover,
depending on the regressor that is used, additional transformations may be
desired. When the kernel ridge regressor is used all time-series should be, at
a minimum, scaled to similar magnitudes. This is necessary because, if there
are large discrepancies in the magnitudes of the features, then the effect of
penalty parameter $\lambda$ will differ substantially from one feature to the
next, which could cause some features to either under-penalized or over-
penalized compared to others based on their scale alone.

In practice, for the mlcausality algorithm exclusively, the time-series data
for analysis presented herein underwent a transformation using a quantile
transformer. To be more specific, each individual time-series was partitioned
into 1000 quantiles and then subjected to a transformation that mapped them
onto a uniform distribution scaled within the range of 0 to 1. This particular
transformation yields outstanding results when employed in conjunction with
the sign test utilized by mlcausality. With this quantile transformation, the
criterion for better prediction is one that forecasts a quantile closer to the
quantile of the actual value, regardless of whether the prediction is closer
in terms of the units of the outcome variable. In essence, the quantile
transformation effectively compels the sign test to assess models based on
their capacity to predict the quantile of the actual value, rather than their
capability to predict the actual value itself.

##  3 Results

In what follows, the performance of mlcausality (MLC in the following plots)
is compared to the following three algorithms: the large-scale nonlinear
granger causality (lsNGC) algorithm (Wismüller et al., 2021); the mutual
nonlinear cross-mapping methods using local models (LM) algorithm (Sugihara et
al. (2012), using the software implementation provided by Javier (2021)); and
the PC-momentary conditional independence (PCMCI) algorithm (Runge et al.,
2019). For lsNGC, $c_{f}=25$ and $c_{g}=5$, which are suggested in Wismüller
et al. (2021), unless the lsNGC algorithm threw an error, at which point
$c_{f}$ was progressively decreased until all simulated data ran successfully.
For all other algorithms the software defaults were kept.

All models were run on a computer with an AMD Ryzen 5 3600 6-Core processor
with 12 threads in total. In order to ensure that the comparison between
mlcausality and other competing algorithms is as fair as possible, all
algorithms were parallelized to run on 12 processes, one for each thread.

The results presented herein are for lag orders selected using Cao’s minimum
embedding dimension selection method (Cao, 1997). In particular, for every
network and time-series length combination, Cao’s algorithm is run on all
time-series in that combination, with the largest identified minimum embedding
dimension used as the basis for construction that combination’s lag order.

###  3.1 The Simulated Networks

The following briefly describes the data generating processes for the
simulated networks for which performance results are presented. All networks
were initialized using normally distributed white noise $w(t)$ with mean = 0
and variance = 1. Every network was analyzed for time-series lengths of 500,
1000, 1500, and 2000 time-points. For every network and every time-series
length, 50 independent sets were generated after discarding a 500 time-point
burn-in. The network plots for all tested networks are available in figure 1.

(a) 5 Linear and 5 Nonlinear

(b) 7 Nonlinear

(c) 9 Nonlinear

(d) 11 Nonlinear

(e) 34 Zachary1

(f) 34 Zachary2

Figure 1: Network plots

5-node linear network: This network was first proposed as example 3 in Baccalá
and Sameshima (2001). The network is generated using the following
multivariate autoregressive model:

$\displaystyle x_{1}(t)=0.95\sqrt{2}x_{1}(t-1)-0.9025x_{1}(t-2)+w_{1}(t)$ (21)
$\displaystyle x_{2}(t)=0.5x_{1}(t-2)+w_{2}(t)$ $\displaystyle
x_{3}(t)=-0.4x_{1}(t-3)+w_{3}(t)$ $\displaystyle
x_{4}(t)=-0.5x_{1}(t-2)+0.25\sqrt{2}x_{4}(t-1)+0.25\sqrt{2}x_{5}(t-1)+w_{4}(t)$
$\displaystyle
x_{5}(t)=-0.25\sqrt{2}x_{4}(t-1)+0.25\sqrt{2}x_{5}(t-1)+w_{5}(t)$

This network has the following causal connections: $x_{1}\rightarrow x_{2}$,
$x_{1}\rightarrow x_{3}$, $x_{1}\rightarrow x_{4}$, $x_{4}\rightarrow x_{5}$,
and $x_{5}\rightarrow x_{4}$. This network structure has the potential to
present several difficulties for a Granger causality recovery algorithm. For
instance, although there is no direct causal link between $x_{2}$, $x_{3}$,
and $x_{4}$, they may all be correlated because of the causal effect of
$x_{1}$ on all of them.

5-node nonlinear network: This network was introduced in Wismüller et al.
(2021). The causal connections in this 5-node nonlinear network are identical
to those of the 5-node linear network, except some of those connections are
converted into nonlinear ones:

$\displaystyle x_{1}(t)=0.95\sqrt{2}x_{1}(t-1)-0.9025x_{1}(t-2)+w_{1}(t)$ (22)
$\displaystyle x_{2}(t)=0.5x_{1}^{2}(t-2)+w_{2}(t)$ $\displaystyle
x_{3}(t)=-0.4x_{1}(t-3)+w_{3}(t)$ $\displaystyle
x_{4}(t)=-0.5x_{1}^{2}(t-2)+0.5\sqrt{2}x_{4}(t-1)+0.25\sqrt{2}x_{5}(t-1)+w_{4}(t)$
$\displaystyle x_{5}(t)=-0.5\sqrt{2}x_{4}(t-1)+0.5\sqrt{2}x_{5}(t-1)+w_{5}(t)$

7-node nonlinear network: This network introduces many very complicated
nonlinear interactions:

$\displaystyle x_{1}(t)=0.95\sqrt{2}x_{1}(t-1)-0.9025x_{1}(t-2)+w_{1}(t)$ (23)
$\displaystyle x_{2}(t)=-0.04x_{1}^{3}(t-3)+0.04x_{1}^{3}(t-1)+w_{2}(t)$
$\displaystyle
x_{3}(t)=-0.04\sqrt{2}x_{2}^{3}(t-1)+0.04\sqrt{2}x_{2}^{3}(t-2)+w_{3}(t)$
$\displaystyle
x_{4}(t)=ln(1+|x_{3}(t-1)|)*sgn(x_{3}(t-1))+0.001x_{7}^{3}(t-2)-0.001x_{7}^{3}(t-3)+w_{4}(t)$
$\displaystyle x_{5}(t)=0.04*clip(w_{5a}(t),-1,1)*x_{6}^{5}(t-2)+w_{5b}(t)$
$\displaystyle x_{6}(t)=0.04*x_{1}^{3}(t-2)+0.04*x_{3}^{3}(t-1)+w_{6}(t)$
$\displaystyle
x_{7}(t)=clip(w_{7a}(t),-0.5,0.5)*(0.04x_{1}^{3}(t-2)+0.1x_{6}^{2}(t-1)-0.1x_{6}^{2}(t-2))+w_{7b}(t)$

where $clip(v,a,b)$ is a function that limits $v$ to the range $[a,b]$ and
$sgn(v)$ is the sign of $v$. The connections for this network are as follows:
$x_{1}\rightarrow x_{2}$, $x_{1}\rightarrow x_{6}$, $x_{1}\rightarrow x_{7}$,
$x_{2}\rightarrow x_{3}$, $x_{3}\rightarrow x_{4}$, $x_{3}\rightarrow x_{6}$,
$x_{6}\rightarrow x_{5}$, $x_{6}\rightarrow x_{7}$, and $x_{7}\rightarrow
x_{4}$.

9-node nonlinear network: This network combines many nonlinear interactions
with autoregressive terms for all variables:

$\displaystyle x_{1}(t)=0.95\sqrt{2}x_{1}(t-1)-0.9025x_{1}(t-2)+w_{1}(t)$ (24)
$\displaystyle
x_{2}(t)=0.5x_{1}^{2}(t-2)+0.5x_{2}^{2}(t-1)-0.4x_{2}^{2}(t-2)+w_{2}(t)$
$\displaystyle
x_{3}(t)=-0.4x_{1}(t-3)+0.5x_{3}^{2}(t-1)-0.4x_{3}^{2}(t-2)+w_{3}(t)$
$\displaystyle
x_{4}(t)=-0.5x_{1}^{2}(t-2)+0.5x_{4}^{2}(t-1)-0.4x_{4}^{2}(t-2)+0.5\sqrt{2}x_{4}(t-1)+0.25\sqrt{2}x_{5}(t-1)+w_{4}(t)$
$\displaystyle x_{5}(t)=-0.5\sqrt{2}x_{4}(t-1)+0.5\sqrt{2}x_{5}(t-1)+w_{5}(t)$
$\displaystyle
x_{6}(t)=sgn(x_{4}(t-1))*ln(|x_{4}(t-1)|+1)+0.5x_{6}^{2}(t-1)-0.4x_{6}^{2}(t-2)+w_{6}(t)$
$\displaystyle
x_{7}(t)=0.04*clip(w_{7a}(t),-1,1)*x_{6}^{5}(t-2)+0.5x_{7}^{2}(t-1)-0.4x_{7}^{2}(t-2)+w_{7b}(t)$
$\displaystyle
x_{8}(t)=0.4x_{1}(t-2)+0.25x_{3}^{3}(t-1)+0.5x_{8}^{2}(t-1)-0.4x_{8}^{2}(t-2)+w_{8}(t)$
$\displaystyle
x_{9}(t)=clip(w_{9a}(t),-0.5,0.5)*(0.2x_{1}(t-2)+0.1x_{8}^{2}(t-1)-0.1x_{8}^{2}(t-2))+0.5x_{9}^{2}(t-1)-0.4x_{9}^{2}(t-2)+w_{9b}(t)$

where $clip(v,a,b)$ is a function that limits $v$ to the range $[a,b]$ and
$sgn(v)$ is the sign of $v$. The connections of this network are as follows:
$x_{1}\rightarrow x_{2}$, $x_{1}\rightarrow x_{3}$, $x_{1}\rightarrow x_{4}$,
$x_{1}\rightarrow x_{8}$, $x_{1}\rightarrow x_{9}$, $x_{3}\rightarrow x_{8}$,
$x_{4}\rightarrow x_{5}$, $x_{4}\rightarrow x_{6}$, $x_{5}\rightarrow x_{4}$,
$x_{6}\rightarrow x_{7}$, and $x_{8}\rightarrow x_{9}$.

11-node nonlinear network: Yet another network with many nonlinear
interactions:

$\displaystyle x_{1}(t)=0.25x_{1}^{2}(t-1)-0.25x_{1}^{2}(t-2)+w_{1}(t)$ (25)
$\displaystyle x_{2}(t)=ln(1+|x_{1}(t-2)|)*sgn(x_{1}(t-2))+w_{2}(t)$
$\displaystyle x_{3}(t)=-0.1x_{2}^{3}(t-3)+w_{3}(t)$ $\displaystyle
x_{4}(t)=-0.5x_{2}^{2}(t-2)+0.5\sqrt{2}x_{4}(t-1)+0.25\sqrt{2}x_{5}(t-1)+w_{4}(t)$
$\displaystyle x_{5}(t)=-0.5\sqrt{2}x_{4}(t-1)+0.5\sqrt{2}x_{5}(t-1)+w_{5}(t)$
$\displaystyle x_{6}(t)=ln(1+|x_{4}(t-1)|)*sgn(x_{4}(t-1))+w_{6}(t)$
$\displaystyle x_{7}(t)=0.04*clip(w_{7a}(t),-1,1)*x_{6}^{5}(t-2)+w_{7b}(t)$
$\displaystyle x_{8}(t)=0.4x_{1}(t-2)+0.25x_{3}^{3}(t-1)+w_{8}(t)$
$\displaystyle
x_{9}(t)=clip(w_{9a}(t),-0.5,0.5)*(0.2x_{1}(t-2)+0.1x_{8}^{2}(t-1)-0.1x_{8}^{2}(t-2))+w_{9b}(t)$
$\displaystyle
x_{10}(t)=0.25x_{1}^{2}(t-3)-0.01x_{2}^{2}(t-3)+0.15x_{3}^{3}(t-3)+w_{10}(t)$
$\displaystyle
x_{11}(t)=0.1x_{2}^{4}(t-1)-0.1x_{2}^{4}(t-2)+0.1x_{6}^{3}(t-3)+w_{11}(t)$

where $clip(v,a,b)$ is a function that limits $v$ to the range $[a,b]$ and
$sgn(v)$ is the sign of $v$. The connections of this network are as follows:
$x_{1}\rightarrow x_{2}$, $x_{1}\rightarrow x_{8}$, $x_{1}\rightarrow x_{9}$,
$x_{1}\rightarrow x_{10}$, $x_{2}\rightarrow x_{3}$, $x_{2}\rightarrow x_{4}$,
$x_{2}\rightarrow x_{10}$, $x_{2}\rightarrow x_{11}$, $x_{3}\rightarrow
x_{8}$, $x_{3}\rightarrow x_{10}$, $x_{4}\rightarrow x_{5}$, $x_{4}\rightarrow
x_{6}$, $x_{5}\rightarrow x_{4}$, $x_{6}\rightarrow x_{7}$, and
$x_{8}\rightarrow x_{9}$.

34-node Zachary1 and Zachary2 networks: These networks are identical to the
Zachary1 and Zachary2 networks in Wismüller et al. (2021) and are constructed
using the undirected connections in the Zachary karate club dataset (Zachary,
1977). The Zachary dataset is a social network composed of 34 members of a
karate club that lists links between pairs of members that interacted outside
the club. The nodal interactions, adapted from Marinazzo et al. (2008), are as
follows:

$\displaystyle
x_{i}(t)=\bigg{(}1-\sum_{j=1}^{n}c_{ij}\bigg{)}(1-ax_{i}^{2}(t-1))+\sum_{j=1}^{n}c_{ij}(1-ax_{j}^{2}(t-1))+sw_{i}(t)$
(26)

where $c_{ij}$ indicates the coupling $j\rightarrow i$. For the Zachary1
network, all 78 edges linking the 34 nodes in the network are assumed to be
bidirectional and $a=1.8$, $s=0.01$, and $c=0.025$. For the Zachary2 network,
5 of the 78 edges linking the 34 nodes in the network are randomly selected to
be bidirectional, while for the rest of the links, the direction is assigned
randomly. Moreover, for the Zachary2 network $c=0.05$. Note that each of the
50 independent sets of Zachary2 networks could have a different network
structure because the directions of the links are randomly assigned for each
of the 50 independent sets.

###  3.2 Evaluating mlcausality’s Network Recovery Performance

####  3.2.1 Non-thresholded Metrics

Figure 2 shows the area under the receiver operating characteristic curve
(AUC) for all model, network, and time-series length combinations. For
nonlinear networks up to around 11 nodes in size mlcausality with kernel ridge
regression exhibits leading, joint-leading, or near-leading AUC performance,
with AUC performance significantly declining for the 34-node Zachary1 and
Zachary2 networks. For networks with 9 nodes or less substantial performance
improvements are not observed as the time-series length increases. For the
11-node nonlinear, Zachary1, and Zachary2 networks mlcausality’s performance
steadily improves as the time-series length increases: this indicates that for
networks greater than around 10 nodes in length at least 2000 time-points are
needed to achieve peak performance with the default 70-30 train-test split. A
similar pattern is observed for rival algorithms when recovering the Zachary1
and Zachary2 networks but not the 11-node nonlinear network. This implies that
mlcausality is more “data hungry” than competing algorithms, a fact that
should not be surprising given that mlcausality splits the data into separate
train and test sets. It can therefore be concluded that for networks of up to
10 nodes in size that have at least 500 time-points mlcausality achieves
leading or near-leading AUC performance. Moreover, given the data-hungry
nature of mlcausality and the trajectory of improving AUC scores as the number
of time-points increases for the 34-node Zachary1 and Zachary2 networks, it is
plausible to assume that mlcausality might be capable of matching or exceeding
the AUC performance of other algorithms given sufficently long time-series
data.

Figure 2: AUC boxplots for time-series of different lengths. In all cases the
median of the distribution is represented by a black square; the upper and
lower edges of a box represent the 25th and 75th percentiles; and the whiskers
extend to the distribution maximums and minimums. MLC indicates the
mlcaulsality model; lsNGC indicates the large-scale nonlinear granger
causality algorithm; LM indicates the mutual nonlinear cross-mapping methods
using local models approach; and PCMCI indicates the PC-momentary conditional
independence algorithm.

Figure 3 compares the Brier scores for all network, time-series length, and
model combinations. Here, we see that mlcausality exhibits significantly lower
Brier scores than competing algorithms. The implication of lower Brier scores
is that mlcausality’s $p$-values are better calibrated than rival algorithms
and are a truer reflection of actual probabilities. Furthermore, for most of
the analyzed networks mlcausality’s Brier scores tend to decrease as the time-
series length increases: this implies that the $p$-values tend to converge
towards true probabilities as time-series lengths increase. The same is not
true for any of the other tested networks, with lsNGC in particular generally
exhibiting higher Brier scores with increased time-series lengths. The above
suggests that $p$-value-based thresholding rules are likely to perform better
for mlcausality than competing algorithms, especially with long time-series
data.

Figure 3: Brier score boxplots for time-series of different lengths. In all
cases the median of the distribution is represented by a black square; the
upper and lower edges of a box represent the 25th and 75th percentiles; and
the whiskers extend to the distribution maximums and minimums. MLC indicates
the mlcaulsality model; lsNGC indicates the large-scale nonlinear granger
causality algorithm; LM indicates the mutual nonlinear cross-mapping methods
using local models approach; and PCMCI indicates the PC-momentary conditional
independence algorithm.

####  3.2.2 Metrics at G-mean Optimal Thresholds

Figure 4 shows the accuracy at optimal geometric mean of sensitivity and
specificity (G-mean) $p$-value thresholds. Specifically, for every model,
network, and time-series combination, 100 equally-spaced $p$-value thresholds
were tested, with the $p$-value threshold that generated the highest median
G-mean for the 50 independent sets in that combination being chosen. Figure 4
indicates that mlcausality exhibits leading or joint-leading accuracy at
optimal $p$-value thresholds for networks up to and including 11 nodes in
size, and near-leading accuracy for the 34-node Zachary networks when the
number of observations is high (2000+ time-points).

Moreover, figure 4 suggests that relying on AUC scores alone can be a somewhat
misleading indicator of true model performance when a thresholding criteria
must be implemented. For instance, despite near-perfect AUC scores for the
lsNGC model on the 5-node linear network, at optimal G-mean thresholds, the
accuracies for those same model and network combinations are well below 1.
This seemingly disturbing discrepancy is not the result of an error. The lsNGC
AUC scores indicate that, for nearly every independent instance of the 5-node
linear model, there exists some split that leads to perfect classification. On
the other hand, the accuracy scores in figure 4 are evaluated at thresholds
that maximize the median G-mean for _all_ 50 independent sets in a
network–time-series combination. Although almost every single instance of the
5-node linear network can be perfectly classified by some $p$-value split when
using the lsNGC model, that split is not stable from one independent set of
data to the next. Furthermore, as the number of time-points increases, the
accuracy of the lsNGC model tends to decline, which is entirely consistent
with the increasing Brier scores in time-series lengths observed in figure 3.

Figure 4: Accuracy boxplots for a threshold at the G-mean max for time-series
of different lengths. In all cases the median of the distribution is
represented by a black square; the upper and lower edges of a box represent
the 25th and 75th percentiles; and the whiskers extend to the distribution
maximums and minimums. MLC indicates the mlcaulsality model; lsNGC indicates
the large-scale nonlinear granger causality algorithm; LM indicates the mutual
nonlinear cross-mapping methods using local models approach; and PCMCI
indicates the PC-momentary conditional independence algorithm.

Figure 5 shows the balanced accuracy scores at the maximal G-mean threshold.
Balanced accuracy differs from accuracy in that balanced accuracy measures the
average of accuracy scores for both the minority and majority classes (in
effect, this is equivalent to taking the average of the sensitivity and the
specificity). The balanced accuracy scores confirm that mlcaulsality exhibits
leading performance at optimal $p$-value thresholds for networks up to and
including 9 nodes in size. For the 11-node and Zachary2 networks mlcaulsality
achieves near-leading balanced accuracy only when the number of observations
is very high (2000+ time-points). For the Zachary1 network near-leading
balanced accuracy is never achieved, but mlcaulsality does exhibit rapidly
rising balanced accuracy scores in time-series length which suggests that
leading or near-leading performance may be possible with a sufficiently long
time-series data.

Figure 5: Balanced accuracy boxplots for a threshold at the G-mean max for
time-series of different lengths. In all cases the median of the distribution
is represented by a black square; the upper and lower edges of a box represent
the 25th and 75th percentiles; and the whiskers extend to the distribution
maximums and minimums. MLC indicates the mlcaulsality model; lsNGC indicates
the large-scale nonlinear granger causality algorithm; LM indicates the mutual
nonlinear cross-mapping methods using local models approach; and PCMCI
indicates the PC-momentary conditional independence algorithm.

####  3.2.3 Metrics at the $p$-value = 0.05 Threshold

When Granger causality algorithms are used in practice, the true nature of the
Granger causal relationships in the studied networks are unknown to the
researcher, and therefore the researcher will not be able to use an optimal
threshold derived from unobserved sensitivity and specificity such as G-mean
or Youden’s index. However, the researcher will have to pick and settle on
some threshold anyway in order to recover the relationships and make the
analysis meaningful. Figures 6 and 7 below show the accuracy and balanced
accuracy scores when the $p$-value threshold is set to 0.05, a commonly used
significance level. Comparing accuracy and balanced accuracy scores in figures
6 and 7 to their counterparts at G-mean optimal thresholds in figures 4 and 5
reveals that the 0.05 $p$-value threshold yields excellent, near-optimal
results for the mlcausality algorithm for the analyzed networks. Moreover, at
the 0.05 threshold, mlcausality achieves superior accuracy to all competing
networks, and leading or joint-leading balanced accuracy scores for networks
of 9 nodes or less. For networks with 11 nodes or more mlcausality’s balanced
accuracy at the 0.05 threshold increases with time-series length, which
suggests that, for sufficiently long time-series, mlcausality may be able to
match or exceed the balanced accuracy performance of competing algorithms.

Figure 6: Accuracy boxplots for a threshold at $p$-value=0.05 for time-series
of different lengths. In all cases the median of the distribution is
represented by a black square; the upper and lower edges of a box represent
the 25th and 75th percentiles; and the whiskers extend to the distribution
maximums and minimums. MLC indicates the mlcaulsality model; lsNGC indicates
the large-scale nonlinear granger causality algorithm; LM indicates the mutual
nonlinear cross-mapping methods using local models approach; and PCMCI
indicates the PC-momentary conditional independence algorithm. Figure 7:
Balanced accuracy boxplots for a threshold at $p$-value=0.05 for time-series
of different lengths. In all cases the median of the distribution is
represented by a black square; the upper and lower edges of a box represent
the 25th and 75th percentiles; and the whiskers extend to the distribution
maximums and minimums. MLC indicates the mlcaulsality model; lsNGC indicates
the large-scale nonlinear granger causality algorithm; LM indicates the mutual
nonlinear cross-mapping methods using local models approach; and PCMCI
indicates the PC-momentary conditional independence algorithm.

###  3.3 Evaluating mlcausality’s Runtime Performance

Table 1 presents the runtimes for all model, network, and time-series length
combinations. In all cases, mlcausality with kernel ridge regression and the
RBF kernel runs significantly faster than rival algorithms, in some cases more
than 10 times faster than the second fastest algorithm. Moreover, mlcausality
appears to scale much better, in terms of runtime, to increasing network sizes
and time-series lengths than competing algorithms. When coupled with improving
AUC, accuracy, and balanced accuracy performance in time-series length,
mlcausality’s excellent runtime performance provides strong arguments in favor
of its usage to handle exceptionally large and complex time-series networks.

Network Details | Runtime for 50 iterations using multiprocessing (seconds)
---|---
Network | Num. obs. | mlcausality | lsNGC | LM | PCMCI
5-linear | 500 |  <1 |  23 |  4 |  5
1000 |  <1 |  47 |  11 |  6
1500 |  2 |  72 |  22 |  7
2000 |  5 |  95 |  36 |  7
5-nonlinear | 500 |  <1 |  24 |  4 |  4
1000 |  <1 |  48 |  11 |  5
1500 |  2 |  72 |  22 |  6
2000 |  5 |  95 |  35 |  6
7-nonlinear | 500 |  <1 |  36 |  9 |  9
1000 |  1 |  73 |  24 |  11
1500 |  3 |  111 |  45 |  14
2000 |  7 |  147 |  72 |  17
9-nonlinear | 500 |  <1 |  60 |  16 |  11
1000 |  1 |  120 |  42 |  13
1500 |  4 |  180 |  77 |  15
2000 |  9 |  237 |  120 |  17
11-nonlinear | 500 |  <1 |  74 |  25 |  17
1000 |  2 |  147 |  67 |  20
1500 |  6 |  220 |  124 |  23
2000 |  11 |  294 |  198 |  27
34-Zachary1 | 500 |  5 |  245 |  261 |  344
1000 |  15 |  479 |  674 |  319
1500 |  33 |  719 |  1262 |  403
2000 |  62 |  942 |  2012 |  473
34-Zachary2 | 500 |  5 |  245 |  259 |  329
1000 |  16 |  492 |  680 |  438
1500 |  33 |  730 |  1267 |  532
2000 |  61 |  976 |  2006 |  632
Table 1: Runtimes (in seconds) for different models and networks for a
computer with an AMD Ryzen 5 3600 6-Core dual-thread processor and 32 Gb of
ram. 50 independent sets of each network were evaluated. All algorithms were
parallelized to run 12 processes, one for each thread on the Ryzen processor.

##  4 Conclusion

In this paper I presented a new method and associated Python library,
mlcausality, for identifying nonlinear Granger causal relationships. Despite
mlcausality containing a plug-in architecture that allows for the usage of any
nonlinear regressor, the analysis presented herein focused specifically on the
kernel ridge regressor with the radial basis function kernel. On simulated
networks mlcausality with the kernel ridge regressor and the RBF kernel
achieves excellent performance at runtimes that are, in many cases, up to 10
times lower than the second fastest competing algorithm. For moderately sized
networks of up to 11 nodes mlcausality achieves leading or highly competitive
recovery performance as measured by AUC, and tends to achieve superior
accuracy and balanced accuracy when thresholded using $p$-value-based
criteria. Furthermore, mlcausality’s improving AUC, accuracy, balanced
accuracy and Brier scores in time-series length for large networks, when
coupled with mlcausality’s excellent runtime performance, position this
algorithm well to handle large networks with many time-steps. In short,
mlcausality achieves leading nonlinear Granger causality performance at a
fraction of the computational time of rival algorithms.

##  5 Code Availability

The mlcausality library itself is publicly available at
https://github.com/WojtekFulmyk/mlcausality. Replication codes for all results
in this paper are available at https://github.com/WojtekFulmyk/mlcausality-
krr-paper-replication.

## References

  * Baccalá and Sameshima (2001) Luiz A Baccalá and Koichi Sameshima.  Partial directed coherence: a new concept in neural structure determination.  _Biological cybernetics_ , 84(6):463–474, 2001\.
  * Cao (1997) Liangyue Cao.  Practical method for determining the minimum embedding dimension of a scalar time series.  _Physica D: Nonlinear Phenomena_ , 110(1-2):43–50, 1997.
  * Dixon and Mood (1946) Wilfrid J Dixon and Alexander M Mood.  The statistical sign test.  _Journal of the American Statistical Association_ , 41(236):557–566, 1946.
  * DSouza et al. (2017) Adora M DSouza, Anas Z Abidin, Lutz Leistritz, and Axel Wismüller.  Exploring connectivity with large-scale granger causality on resting-state functional mri.  _Journal of neuroscience methods_ , 287:68–79, 2017.
  * Exterkate et al. (2016) Peter Exterkate, Patrick JF Groenen, Christiaan Heij, and Dick van Dijk.  Nonlinear forecasting with many predictors using kernel ridge regression.  _International Journal of Forecasting_ , 32(3):736–753, 2016.
  * Gao et al. (2017) Zhong-Ke Gao, Michael Small, and Juergen Kurths.  Complex network analysis of time series.  _Europhysics Letters_ , 116(5):50001, 2017.
  * Granger (1969) Clive WJ Granger.  Investigating causal relations by econometric models and cross-spectral methods.  _Econometrica: journal of the Econometric Society_ , pages 424–438, 1969.
  * Javier (2021) Prince Joseph Erneszer Javier.  causal-ccm a Python implementation of Convergent Cross Mapping, 6 2021\.  URL https://github.com/PrinceJavier/causal_ccm.
  * Lacasa et al. (2015) Lucas Lacasa, Vincenzo Nicosia, and Vito Latora.  Network structure of multivariate time series.  _Scientific reports_ , 5(1):1–9, 2015.
  * Marinazzo et al. (2008) Daniele Marinazzo, Mario Pellicoro, and Sebastiano Stramaglia.  Kernel-granger causality and the analysis of dynamical networks.  _Physical review E_ , 77(5):056215, 2008.
  * Mercer (1909) James Mercer.  Xvi. functions of positive and negative type, and their connection the theory of integral equations.  _Philosophical transactions of the royal society of London. Series A, containing papers of a mathematical or physical character_ , 209(441-458):415–446, 1909.
  * Murphy (2012) Kevin P Murphy.  _Machine learning: a probabilistic perspective_.  MIT press, 2012.
  * Rosoł et al. (2022) Maciej Rosoł, Marcel Młyńczak, and Gerard Cybulski.  Granger causality test with nonlinear neural-network-based methods: Python package and simulation study.  _Computer Methods and Programs in Biomedicine_ , 216:106669, 2022.
  * Runge et al. (2019) Jakob Runge, Peer Nowack, Marlene Kretschmer, Seth Flaxman, and Dino Sejdinovic.  Detecting and quantifying causal associations in large nonlinear time series datasets.  _Science advances_ , 5(11), 2019.
  * Shashua (2009) Amnon Shashua.  Introduction to machine learning: Class notes 67577.  _arXiv preprint arXiv:0904.3664_ , 2009.
  * Sugihara et al. (2012) George Sugihara, Robert May, Hao Ye, Chih-hao Hsieh, Ethan Deyle, Michael Fogarty, and Stephan Munch.  Detecting causality in complex ecosystems.  _science_ , 338(6106):496–500, 2012.
  * Vovk (2013) Vladimir Vovk.  Kernel ridge regression.  In _Empirical Inference: Festschrift in Honor of Vladimir N. Vapnik_ , pages 105–116. Springer, 2013.
  * Wilcoxon (1945) F Wilcoxon.  Individual comparisons by ranking methods. biom. bull., 1, 80–83, 1945\.
  * Wismüller et al. (2021) Axel Wismüller, Adora M Dsouza, M Ali Vosoughi, and Anas Abidin.  Large-scale nonlinear granger causality for inferring directed dependence from short multivariate time-series data.  _Scientific reports_ , 11(1):7817, 2021.
  * Woodbury (1950) Max A Woodbury.  _Inverting modified matrices_.  Department of Statistics, Princeton University, 1950.
  * Zachary (1977) Wayne W Zachary.  An information flow model for conflict and fission in small groups.  _Journal of anthropological research_ , 33(4):452–473, 1977.