∎

11institutetext: †S N Bose National Centre For Basic Sciences 22institutetext:
Block-JD, Sector-III, Salt Lake, Kolkata 700106, India
22email<EMAIL_ADDRESS>(Corresponding author∗) 33institutetext: ‡Indian
Association for the Cultivation of Science 44institutetext: Jadavpur, Kolkata
700032, India

#  Reduction of kinetic equations to Liénard-Levinson-Smith form: Counting
limit cycles

Sandip Saha† Gautam Gangopadhyay†∗ Deb Shankar Ray‡

###### Abstract

We have presented an unified scheme to express a class of system of equations
in two variables into a Liénard - Levinson - Smith(LLS) oscillator form. We
have derived the condition for limit cycle with special reference to Rayleigh
and Liénard systems for arbitrary polynomial functions of damping and
restoring force. Krylov-Boguliubov(K-B) method is implemented to determine the
maximum number of limit cycles admissible for a LLS oscillator atleast in the
weak damping limit. Scheme is illustrated by a number of model systems with
single cycle as well as the multiple cycle cases.

##  1 Introduction

Various open kinetic systemsmurraybio ; murraynld ; epstein ; strogatz ;
goldbook ; remickens  in physics, chemistry and biology, are generically
described by a minimal model of autonomous coupled differential
equationsslross ; arnold ; birkhoff ; jjstoker ; smale1967 ; len0  of two
variables. They exhibit self sustained oscillation in the form of stable limit
cycle in a phase plane in many examples, such as, chemical reactionsepstein ;
brusselator2009 ; len4 , biological rhythmsmurraynld ; murraybio ; kaiser83 ;
kaiser91 ; k-dsr ; strogatz , vibrations in mechanicalrand2012 , optical
system and musical instrumentstrogatz ; dsrrayleigh , to name a few. A
Rayleighdsrrayleigh  equation in violin string and van der Pol oscillation in
electric circuit are the classic examples, in this context. More generally,
Liénardperko ; len4 ; strogatz ; lls1 ; lls2 ; lls3  equation underlines the
concrete criteria for the existence of at least one limit cycle for a general
class of such systems of which van der Pol is a special case of the form
$\ddot{x}+f(x)\dot{x}+x=0$ where $f(x)=\epsilon(x^{2}-1)$ and Liénard
transformation is $\dot{x}=y-F(x)$ and $\dot{y}=-x$ with
$F(x)=\int_{0}^{x}f(\tau)d\tau$. A further generalisation of Liénard equation
is the LLS equationlls1 ; lls2 ; lls3 ; remickens ,
$\ddot{x}+F(x,\dot{x})\dot{x}+G(x)=0$, sometimes called the generalised
Liénard equation. Casting a general system of kinetic equations in two
variables which describe a variety of scenarios in physical, chemical, bio-
chemical and ecological sciences into LLS formremickens  is often not
straight-forwardremickens ; len4 ; limiso ; strogatz . To this end we have
provided a scheme for a wide class of open nonlinear equations, cast in the
LLS form so that the later becomes amenable to several techniques in nonlinear
dynamics.

Our next objective is to find the nature and the number of limit cycles for a
given LLS equation thereby addressing the second part of the Hilbert’s
$16^{th}$ problem. The problem of counting limit cycle has a long legacy since
Hilbert, Smale and many others and still continues it without complete
understandingbirkhoff ; jjstoker ; smale1998 ; perko ; gaiko2008 ;
countinglcjkb ; infdampinglcjkb ; lcbounestjkb . Our scheme is based on the KB
method of averagingkbbook ; strogatz ; slross ; len0 , a variant of multi-
scale perturbation techniquestrogatz ; chen1 ; chen2 ; sarkar2012  to derive
amplitude equation with considering the polynomial forms of the nonlinear
damping and restoring force functions. We have illustrated our results on a
variety of known model systemsperko ; gaiko2008 ; giacomini  with single and
multiple limit cyclesperko ; gaiko2008 ; giacomini .

##  2 Reduction of Kinetic Equations to Liénard - Levinson - Smith(LLS) form:
Conditions for limit cycle

We consider here a set of two dimensional autonomous kinetic equations for an
open system. Our aim is to cast the equations into a form of a variant of LLS
oscillatorstrogatz ; len4 ; limiso  or LLS oscillatorstrogatz ; len4 ; limiso
; len0 ; remickens  which can further be reduced to Rayleigh and Liénard form.
Let us begin with the system of autonomous kinetic equations

$\displaystyle\frac{dx}{dt}$ $\displaystyle=a_{0}+a_{1}x+a_{2}y+f(x,y),$
$\displaystyle\frac{dy}{dt}$ $\displaystyle=b_{0}+b_{1}x+b_{2}y+g(x,y),$ (1)

where $x(t)$ and $y(t)$ are, for example, field variables or populations of
species of chemical, biological or ecological process strogatz ; goldbook ;
murraynld ; murraybio  with $a_{i},b_{i}$ for $i=0,1,2$ are all real
parameters expressed in terms of the appropriate kinetic constants. Let,
($x_{s},y_{s}$) be the fixed point of the system and $f(x,y)$ and $g(x,y)$ are
the non-linear functions of $x$ and $y$. A first step is shifting the steady
state ($x_{s},y_{s}$) to the origin ($0,0$) with the help of a linear
transformation as LLS system is a second order homogeneous ordinary
differential equation.

The linear transformation can be chosen by introducing a new pair of variables
$(\xi,u)$, both of which are functions of $x$ and $y$ where
$\xi=\beta_{0}+\beta_{1}x+\beta_{2}y$ with
$\beta_{0}=-(\beta_{1}x_{s}+\beta_{2}y_{s})$ i.e.
$\xi=\beta_{1}(x-x_{s})+\beta_{2}(y-y_{s})$ such that $\dot{\xi}=u$.
$\beta_{1},\beta_{2}$ are weighted constants such that it makes the new steady
state at the origin, $\xi_{s}=0$, $u_{s}=0$. $u$ is expressed as
$u=\alpha_{0}+\alpha_{1}x+\alpha_{2}y$, with $\beta_{i}$, $\alpha_{i}$ for
$i=0,1,2$ are all real constants which can be expressed in terms of system
parameters. From the inverse transformation we can easily obtain the
expressions for $x$ and $y$ as given by

$\displaystyle x$
$\displaystyle=\frac{\alpha_{2}\left(\beta_{0}-\xi\right)+\beta_{2}\left(u-\alpha_{0}\right)}{\alpha_{1}\beta_{2}-\alpha_{2}\beta_{1}}=L(\xi,u),$
$\displaystyle y$
$\displaystyle=\frac{\alpha_{1}\left(\xi-\beta_{0}\right)+\beta_{1}\left(\alpha_{0}-u\right)}{\alpha_{1}\beta_{2}-\alpha_{2}\beta_{1}}=K(\xi,u),$
(2)

provided that $\alpha_{1}\beta_{2}-\alpha_{2}\beta_{1}\neq 0$. Differentiating
again, $\dot{\xi}=u$ with respect to the independent variable $t$ we get,

$\displaystyle\ddot{\xi}$
$\displaystyle=\dot{u}=\alpha_{1}\dot{x}+\alpha_{2}\dot{y}$
$\displaystyle=\alpha_{1}\\{a_{0}+a_{1}L(\xi,\dot{\xi})+a_{2}K(\xi,\dot{\xi})+\varphi(\xi,\dot{\xi})\\}+\alpha_{2}\\{b_{0}+b_{1}L(\xi,\dot{\xi})+b_{2}K(\xi,\dot{\xi})+\phi(\xi,\dot{\xi})\\},$
(3)

where, $L(\xi,\dot{\xi})=c_{1}\xi+c_{2}\dot{\xi}+c_{L}$ and
$K(\xi,\dot{\xi})=c_{3}\xi+c_{4}\dot{\xi}+c_{K}$ with
$\begin{bmatrix}c_{1}&c_{2}&c_{L}\\\ c_{3}&c_{4}&c_{K}\end{bmatrix}$
=$\frac{1}{\alpha_{1}\beta_{2}-\alpha_{2}\beta_{1}}$
$\begin{bmatrix}-\alpha_{2}&\beta_{2}&\alpha_{2}\beta_{0}-\alpha_{0}\beta_{2}\\\
\alpha_{1}&-\beta_{1}&\alpha_{0}\beta_{1}-\alpha_{1}\beta_{0}\end{bmatrix}$.
The functions $\varphi$ and $\phi$ can be expressed as a power series
expansion as,

$\varphi(\xi,\dot{\xi})=\sum_{n,m=0}^{\infty}\varphi_{nm}\xi^{n}\dot{\xi}^{m}\hskip
14.22636ptand\hskip
28.45274pt\phi(\xi,\dot{\xi})=\sum_{n,m=0}^{\infty}\phi_{nm}\xi^{n}\dot{\xi}^{m},$
(4)

with, $\phi(\xi,\dot{\xi})=\mu\varphi(\xi,\dot{\xi})$, as the functions $f$
and $g$ are related through $\mu$ by $g=\mu f$, $\mu\in\mathbb{R}$. So, after
putting the above form in equation (3) one can find,

$\displaystyle\ddot{\xi}$
$\displaystyle=\alpha_{1}a_{0}+\alpha_{1}a_{1}(c_{1}\xi+c_{2}\dot{\xi}+c_{L})+\alpha_{1}a_{2}(c_{3}\xi+c_{4}\dot{\xi}+c_{K})+(\alpha_{1}+\mu\alpha_{2})\sum_{n,m=0}^{\infty}\varphi_{nm}\xi^{n}\dot{\xi}^{m}$
$\displaystyle+\alpha_{2}b_{0}+\alpha_{2}b_{1}(c_{1}\xi+c_{2}\dot{\xi}+c_{L})+\alpha_{2}b_{2}(c_{3}\xi+c_{4}\dot{\xi}+c_{K}),i.e.,$
$\displaystyle\ddot{\xi}$
$\displaystyle=A_{00}+\left(A_{10}+\sum_{n>1}A_{n0}\xi^{n-1}\right)\xi+\left(A_{01}+\sum_{n>0}A_{n1}\xi^{n}+\sum_{n\geq
0}\sum_{m>1}A_{nm}\xi^{n}\dot{\xi}^{m-1}\right)\dot{\xi},$ (5)

where,
$\alpha_{1}a_{0}+\alpha_{2}b_{0}+(\alpha_{1}+\mu\alpha_{2})\varphi_{00}+(\alpha_{1}a_{1}+\alpha_{2}b_{1})c_{L}+(\alpha_{1}a_{2}+\alpha_{2}b_{2})c_{K}=A_{00}=0$
(by definition of a zero fixed point of $\xi$),
$A_{10}=\alpha_{1}(a_{1}c_{1}+a_{2}c_{3})+\alpha_{2}(b_{1}c_{1}+b_{2}c_{3})+(\alpha_{1}+\mu\alpha_{2})\varphi_{10}$,
$A_{01}=\alpha_{1}(a_{1}c_{2}+a_{2}c_{4})+\alpha_{2}(b_{1}c_{2}+b_{2}c_{4})+(\alpha_{1}+\mu\alpha_{2})\varphi_{01}$,
$A_{n0}=(\alpha_{1}+\mu\alpha_{2})\varphi_{n0}$,
$A_{n1}=(\alpha_{1}+\mu\alpha_{2})\varphi_{n1}$ and
$A_{nm}=(\alpha_{1}+\mu\alpha_{2})\varphi_{nm}$, where indices follow the
values as given in the summation over $m,n\in\mathbb{Z}^{+}$. Finally, the
above equation looks like,

$\ddot{\xi}+F(\xi,\dot{\xi})\dot{\xi}+G(\xi)=0,$ (6)

where, the functions $F(\xi,\dot{\xi})$ and $G(\xi)$ are given by

$\displaystyle F(\xi,\dot{\xi})$
$\displaystyle=-[A_{01}+\sum_{n>0}A_{n1}\xi^{n}+\sum_{n\geq
0}\sum_{m>1}A_{nm}\xi^{n}\dot{\xi}^{m-1}],$ $\displaystyle G(\xi)$
$\displaystyle=-[A_{10}+\sum_{n>1}A_{n0}\xi^{n-1}]\xi.$ (7)

Equation (6) is a well known equation of generalised Liénard form called LLS
equation. The condition for existence of having at least a locally stable
limit cycle of the dynamical system is $F(0,0)<0\implies A_{01}>0$. It can be
shown from the linear stability analysis that there is a relation between
$F(0,0)$ and eigenvalues ($\lambda_{\pm}$) with, $F(0,0)=-2\
Re(\lambda_{\pm})$. For a LLS system, there are six conditions to have a limit
cycle are given in lls2 ; lls3 ; remickens ; len4 . Out of these six
conditions, the condition $F(0,0)$ plays an important role to have a locally
stable or unstable limit cycle for such kind of systemlen4 ; limiso ; len0 ;
remickens  depending upon the sign of $F(0,0)$ is $<0$ or $>0$, respectively.
In particular, two situations may arise:

I: For $A_{nm}=0$, with $n\geq 2,\forall m$ i.e. there be an unique steady
state ($\xi_{s}=0$) with restoring force linear in $\xi$, then the above form
of (6) looks like

$\ddot{\xi}+F_{R}(\dot{\xi})\dot{\xi}+G_{R}(\dot{\xi})\xi=0,$ (8)

where,

$\displaystyle F_{R}(\dot{\xi})$
$\displaystyle=-[A_{01}+\sum_{m>1}A_{0m}\dot{\xi}^{m-1}],\hskip
10.0ptG_{R}(\dot{\xi})=-[A_{10}+\sum_{m>0}A_{1m}\dot{\xi}^{m}],$ (9)

which is in the form of generalised Rayleigh oscillatordsrrayleigh , the limit
cycle condition modifies to, $F_{R}(0)<0$.

II: For $A_{nm}=0$, with $m\geq 2,\forall n$, which corresponds to Liénard
equation with an unique steady state ($\xi_{s}=0$). This is of the form

$\ddot{\xi}+F_{L}(\xi)\dot{\xi}+G_{L}(\xi)=0,$ (10)

where,

$\displaystyle F_{L}(\xi)$
$\displaystyle=-[A_{01}+\sum_{n>0}A_{n1}\xi^{n}],\hskip
10.0ptG_{L}(\xi)=-[A_{10}+\sum_{n>1}A_{n0}\xi^{n-1}]\xi,$ (11)

where the limit cycle condition is $F_{L}(0)<0$. We know that, for a Liénard
system, the damping force function, $F_{L}(\xi)$ and the restoring force
function, $G_{L}(\xi)$ are even and odd functions of $\xi$, respectively.

However, for generalised Liénard or LLS system the odd-even properties of
$G(\xi)$ and $F(\xi,\dot{\xi})$ have complex ramificationsremickens  for
practical systems. Here, we have examined the properties with the help of
Krylov-Boguliubov averaging method.

##  3 Maximum Number of Limit Cycles

We now restrict ourselves to the case of LLS systems where the
$F(\xi,\dot{\xi})$ and $G(\xi)$ are the polynomial functions of $\xi$ and
$\dot{\xi}$. It is well known that linear functional forms of $F$ and $G$
preclude the existence of limit cycle. This can be readily seen by considering
the typical examples, e.g., a Harmonic oscillator or a weakly nonlinear
oscillator with a potential $\frac{1}{2}\omega_{0}^{2}x^{2}+\frac{1}{3}\lambda
x^{4}$, $0<\lambda<1$ or a Lotka-Volterra model, where one encounters a
center. We therefore consider the polynomial form of nonlinear damping
function $F(\xi,\dot{\xi})$ and restoring force function $G(\xi)$ for our
analysis of limit cycle. In what follows we employ K-B method of averaging to
show that the characteristic even/odd powers of polynomials play crucial role
in determining the behaviour of the associated amplitude and phase equations.

To begin with we consider some fixed values of $m,n$ of equation (7) to
truncate the series at $M,N$, for the highest power of $\dot{\xi}$ and $\xi$,
respectively. For explicit structure of a prototypical example of an amplitude
equation we choose upto $M=N=3$ for illustration. This includes all possible
cases for the even and odd nature of $F(\xi,\dot{\xi})$ and $G(\xi)$,
respectively. Then the above form of $F(\xi,\dot{\xi})$ and $G(\xi)$ will be
in the following reduced forms,

$\displaystyle F(\xi,\dot{\xi})$
$\displaystyle=-[A_{01}+A_{11}\xi+A_{21}\xi^{2}+A_{31}\xi^{3}+A_{02}\dot{\xi}+A_{12}\xi\dot{\xi}+A_{22}\xi^{2}\dot{\xi}+A_{32}\xi^{3}\dot{\xi}$
$\displaystyle+A_{03}\dot{\xi}^{2}+A_{13}\xi\dot{\xi}^{2}+A_{23}\xi^{2}\dot{\xi}^{2}+A_{33}\xi^{3}\dot{\xi}^{2}],$
$\displaystyle G(\xi)$
$\displaystyle=-[A_{10}\xi+A_{20}\xi^{2}+A_{30}\xi^{3}].$ (12)

Let us take $|F(0,0)|=\sigma\in\mathbb{R}^{+}$, an arbitrary constant with
$F(\xi,\dot{\xi})=\sigma F_{\sigma}(\xi,\dot{\xi})$. Then the LLS equation can
be rewritten as

$\displaystyle\ddot{\xi}+\sigma F_{\sigma}(\xi,\dot{\xi})\dot{\xi}+G(\xi)=0.$
(13)

Therefore the final equation takes the form of a non-linear oscillator after
rescaling $t$ by $\tau$ taking, $\omega t\rightarrow\tau$ as

$\displaystyle\ddot{Z}(\tau)+\epsilon h(Z(\tau),\dot{Z}(\tau))+Z(\tau)=0,$
(14)

where, $0<\epsilon=\frac{\sigma}{\omega^{2}}\ll 1$, $\omega^{2}=-A_{10}>0$ and
$Z(\tau)\equiv\xi(t)$ and $\omega\dot{Z}(\tau)\equiv\dot{\xi}(t)$. Equation
(14) is now ready for the treatment using K-B method with

$h(Z,\dot{Z})=-\left[\\{B_{01}+B_{11}Z+B_{21}Z^{2}+B_{31}Z^{3}+B_{02}\omega\dot{Z}+B_{12}Z\omega\dot{Z}+B_{22}Z^{2}\omega\dot{Z}+B_{32}Z^{3}\omega\dot{Z}+B_{03}\omega^{2}\dot{Z}^{2}+B_{13}Z\omega^{2}\dot{Z}^{2}+B_{23}Z^{2}\omega^{2}\dot{Z}^{2}+B_{33}Z^{3}\omega^{2}\dot{Z}^{2}\\}\omega\dot{Z}+B_{20}Z^{2}+B_{30}Z^{3}\right],$
(15)

where $B_{ij}=\frac{A_{ij}}{\sigma}$, $i,j=0,1,2,3$ with $B_{00}=0$ and
$B_{01}$ will take the fixed value, -1, 0, or 1 depending upon the nature of
the fixed point is stable focus, center/center-type or limit cycle,
respectively. Now choosing, $Z(\tau)\approx r(\tau)\hskip
2.84544pt\cos(\tau+\phi(\tau))$ as a solution of eq. (14) we have
$\dot{Z}(\tau)\approx-r(\tau)\hskip 2.84544pt\sin(\tau+\phi(\tau))$ with
slowly varying radius $r(\tau)=\sqrt{Z^{2}+\dot{Z}^{2}}$ and phase
$\phi(\tau)=-\tau+tan^{-1}(-\frac{\dot{Z}}{Z})$. The function $h(Z,\dot{Z})$
contains all the non-linear terms and $\epsilon$ is the non-linearity
controlling parameter i.e. one has to satisfy $0<\sigma\ll\omega^{2}$. Then
one can obtain $\dot{r}(\tau)=\epsilon h\sin(\tau+\phi(\tau))$ and
$\dot{\phi}(\tau)=\frac{\epsilon h}{r(\tau)}\cos(\tau+\phi(\tau))$ i.e. the
time derivative of amplitude and phase are of $O(\epsilon)$. So, after taking
a running averagestrogatz ; remickens ; slross  of a time dependent function
$U$ defined as, $\overline{U}(\tau)=\frac{1}{2\pi}\int_{0}^{2\pi}U(s)ds$, one
finds, $\dot{\overline{r}}=\langle\epsilon
h\sin(\tau+\phi(\tau))\rangle_{\tau}$ and
$\dot{\overline{\phi}}=\langle\frac{\epsilon
h}{r(\tau)}\cos(\tau+\phi(\tau))\rangle_{\tau}$, which gives,

$\displaystyle\dot{\overline{r}}$
$\displaystyle=\frac{\epsilon\omega\overline{r}}{16}\\{\overline{r}^{2}\left(B_{23}\overline{r}^{2}\omega^{2}+6B_{03}\omega^{2}+2B_{21}\right)+8B_{01}\\}+O(\epsilon^{2}),$
$\displaystyle\dot{\overline{\phi}}$
$\displaystyle=-\frac{\epsilon\overline{r}^{2}}{16}\left(B_{32}\overline{r}^{2}\omega^{2}+2B_{12}\omega^{2}+6B_{30}\right)+O(\epsilon^{2}).$
(16)

Now from a close look at the equation for $\dot{\overline{r}}$, it is apparent
that only even elements of $F(\xi,\dot{\xi})$ appears but none of any elements
of $G(\xi)$ is present due to the zero averages of $\sin^{\mu}\cos^{\nu}$
terms with $\mu=1$ and $\nu\in\mathbb{Z}$. The non-zero averages arise only
when $\mu,\nu$ both are even i.e.
$\mu=2\eta_{1},\nu=2\eta_{2};\eta_{1},\eta_{2}\in\mathbb{Z}$. Thus, the effect
in $\dot{\overline{r}}$ appears only through the even coefficients of
$F(\xi,\dot{\xi})$ i.e. by examining the respective variables in the
$\dot{\overline{r}}$ equation, we find that only some even coefficients appear
for the first order correction. On the other hand $\dot{\overline{\phi}}$
contains only even coefficients of $F(\xi,\dot{\xi})$ which are not in
amplitude equation along with odd coefficients of $G(\xi)$ which shows that
only odd $G(\xi)$ plays a role here.

So, from the equation of $\dot{\overline{r}}$, one finds that there exist at
most $4$ non-zero values of $\overline{r}$. If out of the four roots every
pair appears as conjugate then there are three possibilities. The cases are,
(i) two different sets of complex conjugate roots giving an asymptotically
stable solution, (ii) one pair of complex conjugate roots and two real roots
of equal magnitude with opposite sign implying a limit cycle solution having
only one cycle and (iii) either four real roots of equal magnitude with
opposite sign having double multiplicity gives a limit cycle solution with
only one cycle or two different sets of real roots of equal magnitude with
opposite sign, may give limit cycle solution with two different cycles of
different radius. The unique zero values of the roots of $\overline{r}$ gives
a center or center-typelen3.5  situation. So, in short, the existence of a
nonzero real root will provide the radius of the cycle which will be stable or
unstable depending on the $-ve$ or $+ve$ sign of
$\frac{d\dot{\overline{r}}}{d\overline{r}}$, at
$\overline{r}=\overline{r}_{ss}$ and at $\overline{r}_{ss}=0$
$\frac{d\dot{\overline{r}}}{d\overline{r}}>0$ or $<0$ gives the nature of the
fixed point.

As an example, for Kaiser modelkaiser83 ; kaiser91 ; k-y2007a ; k-y2007b ;
k-y2007c ; k-y2010 ; k-dsr ; k-y2012 , there exist three limit cycles for a
certain range of $\alpha,\beta$. So, if we choose the parameters, $\alpha$ and
$\beta$ from the three limit cycle zone then there exist six real roots with
three different pairs i.e., three different radii exist according to three
cycles. But, slightly away from the three limit cycle zone, there will exist
only a pair of real roots with the same magnitude and other four will appear
as a complex conjugate pairs and together produces only a stable limit cycle.

Note that, to have a stable limit cycle solution, one condition must be
satisfied i.e. $F(0,0)<0$. But, it fails to give how many limit cycles the
system can admit. According to the root finding algorithm one can guess the
maximum number of cycles of a LLS system. The condition $F(0,0)<0$ plays an
important role as a check for the existence of atleast one stable limit cycle.
But for 2-cycle situations one can have at first the locally unstable limit
cycle before locating the outer stable limit cycle and in this situation
$F(0,0)>0$.

Based on these considerations we have prepared a table(Table-I) illustrating
the possible cases for the maximum number of non-zero real roots or the limit
cycles.

Table-I: Maxumum number of limit cycles for LLS system  $N$ $M$ $N+M$ Max. No.
of Non-zero Max. No. of Real Roots (Even) Cycle(s) Even Even Even
$N+M-2=(N-1)+(M-1)$ $\frac{N+M}{2}-1$ Even Odd Odd $N+M-1=(N)+(M-1)$
$\frac{N+M-1}{2}$ Odd Even Odd $N+M-3=(N-2)+(M-1)$ $\frac{N+M-3}{2}$ Odd Odd
Even $N+M-2=(N-1)+(M-1)$ $\frac{N+M}{2}-1$

Now if we denote the non-zero real values of $\overline{r}$ as an existence of
limit cycles as $\overline{r}$ gives the radius of the cycle where at the same
time a pair of conjugate (one $+ve$ and one $-ve$) roots of equal magnitude
exists for such kind of LLS systems then out of these two roots, radius will
be measured by the magnitude and each distinct magnitude counts the number of
cycles. For example, if there exists six roots, say, $(p,-p)$ occurring twice
and $(q,-q)$ occurring once then the number of cycles will be $2$ of radius
$p$ and $q$, respectively. So, if there are all real roots occurring once,
then the number of cycles will be atmost $\frac{N+M-2}{2}$ or
$\frac{N+M-1}{2}$ or $\frac{N+M-3}{2}$. For LLS equation with $N,M$ are the
maximum power of $\xi$ and $\dot{\xi}$ respectively, we have performed the K-B
analysis numerically for $N=10,M=10$. The result is given in table-II. For
Rayleigh system with $N=1$, for all $M\geq 1$, maximum number of limit cycle
will be $\frac{M-1}{2}$ or $\frac{M-2}{2}$ for odd or even $M$, respectively.
For Liénard system with $M=1$ for all $N\geq 1$, the maximum number of limit
cycle becomes $\frac{N-1}{2}$ or $\frac{N}{2}$ for odd or even $N$,
respectively. The above table is valid for an arbitrary finite polynomials of
$F$ and $G$. For the case of arbitrary infinite polynomialcountinglcjkb ;
infdampinglcjkb ; lcbounestjkb  cases maximum number of limit cycles can be
stated for finite truncation.

##  4 Applications to Some Model Systems

Here we have examined three classes of physical models where the analysis of
the maximum number of limit cycles holds. This connection with the general
model system is discussed with polynomial damping and restoring force
function.

###  4.1 One-Cycle Cases: van der Pol Oscillator, Simple Glycolytic
Oscillator, Modified Brusselator Model

Considering the van der Pol oscillatorstrogatz ; limiso ; slross ; kaiser83 ;
epstein ; murraynld ; powerlaw  with equation,
$\ddot{x}+\epsilon(x^{2}-1)\dot{x}+x=0$ having the weak nonlinearity for
$0<\epsilon\ll 1$ produces a locally stable limit cycle with $F(0,0)<0$. So,
if we compare with the general table we have $N=2$ and $M=1$. This gives a
condition for a unique stable limit cycle.

Next considering the Liénard formlimiso ; len3.5  of simple Glycolytic
oscillatorstrogatz ; limiso ; len3.5 ; epstein ; murraynld ; goldbook  as,

$\ddot{\xi}+\left[(1+a+3b^{2})-2b\xi-2bk-3b\dot{\xi}+\xi\dot{\xi}+k\dot{\xi}+\dot{\xi}^{2}\right]\dot{\xi}+(a+b^{2})\xi=0;a,b>0,k=b+\frac{b}{a+b^{2}}$

, has a unique stable limit cycle with $F(0,0)<0$limiso  having $N=1$ and
$M=3$. This gives one limit cycle.

Furthermore, considering the Modified Brusselator model having the Liénard
formlen4 ; limiso ; goldbook ,

$\ddot{\xi}+\left[-\frac{2a_{1}\xi}{\alpha}-b+\frac{a_{1}^{2}}{\alpha^{2}}+\alpha-\frac{2a_{1}\dot{\xi}}{\alpha^{2}}+\frac{b\dot{\xi}}{a_{1}}+\frac{\dot{\xi}^{2}}{\alpha^{2}}+\frac{\xi\dot{\xi}}{\alpha}\right]\dot{\xi}+\frac{a_{1}^{2}\xi}{\alpha}=0;a_{1},b,\alpha>0$

gives a unique stable limit cycle with $F(0,0)<0$len4 ; limiso , where $N=1$
and $M=3$ again giving rise to the same situation.

###  4.2 Two-Cycle Cases

We rewrite the Liénard form according to Ref.perko ; giacomini  as,
$\dot{x(t)}=y(t)-F(x(t)),\dot{y(t)}=-x(t)$, where $F(x(t))$ is an odd
polynomial. After taking derivative it takes the following form
$\ddot{x}+F^{\prime}(x)\dot{x}+x=0$, where $F^{\prime}(x)=\frac{\partial
F(x)}{\partial x}$ now becomes the form of an even polynomial. For
$F(x)=a_{1}x+a_{2}x^{2}+a_{3}x^{3}$, it has been shownperko ; giacomini ; lins
that the system allows a unique limit cycle if $a_{1}a_{3}<0$, which will be
stable if $a_{1}<0$ and unstable if $a_{1}>0$. This corresponds to the table
$N=2,M=1$. Further extension by Rychkovrychkov  shows that for
$F(x)=(a_{1}x+a_{3}x^{3}+a_{5}x^{5})$ the number of limit cycle is atmost two.
Numerical simulation corroborates this observation when $F(x)$ is chosen as in
Ref.perko ; giacomini ; lins , $F(x)=0.32x^{5}-\frac{4}{3}x^{3}+0.8x$. For
this case the inner one is unstable limit cycle as $F(0,0)=0.8>0$ i.e. a
stable fixed point but the outer one is a stable limit cycle. Here, as per
Table-I we have $N=4$ and $M=1$. The above table thus gives the strategies to
find out the number of limit cycles(both stable and unstable) a system can
have. On the other hand our analysis by K-B method provides a hint towards a
choice of the parameter space for search of real roots of the radial equation.

###  4.3 Three-Cycles Case: Kaiser Bi-rythmicity Model

Extending van der Pol oscillator model with a nonlinear function of higher
order polynomial, Kaiserkaiser83 ; kaiser91 ; k-y2007a ; k-y2007b ; k-y2007c ;
k-y2010 ; k-dsr ; k-y2012  has described bi-rythmicity with the nonlinear
equation,

$\displaystyle\ddot{x}-\mu(1-x^{2}+\alpha x^{4}-\beta
x^{6})\dot{x}+x=E\cos\Omega t.$ (17)

Here, $\alpha,\beta,\mu>0$ tune the non linearity. This is a prototype self-
sustained oscillatory system in absence of $E$ and $\Omega$ which are the
amplitude and the frequency of the external excitation, respectively. The
model exhibits an extremely rich bifurcation behaviour and the system actually
produces bi-rhythmicity. It has been emphasized that in the undriven case, the
model is a multi-limit cycle oscillator and has three limit cycles, two of
them are stable and between the two stable limit cycles there is an unstable
one which divides the basins of attraction of the two stable cycles. In
presence of $E$ and $\Omega$, the above system exhibits some interesting
featureskaiser83 ; kaiser91 ; k-y2007a ; k-y2007b ; k-y2007c ; k-y2010 ;
k-y2012 . From table-I one has $N=6$ and $M=1$ with even-odd sub cases, while
$E=0$. Thus, there may have 6 roots for the radial equation if $\mu>0$ and
$\alpha,\beta$(controlling parameters of the radii) are chosen from three
limit cycle zone($\alpha=0.144,\beta=0.005$) and finally, the number of
distinct values will be 3 which implies that the system can have atmost three
limit cycles (but here it is exactly 3). Further, if we choose $\beta=0$ for
the above undriven Kaiser model with $\alpha=0.1$, one can have two limit
cycles with radii $\approx 2.35$ and $\approx 3.80$, respectively, of which
the smaller one will be stable and the larger one will be unstable.

Note that, if there are odd number of limit cycles, say $l$, then out of the
$l-$cycles, $\frac{l+1}{2}$ will be stable limit cycles and the remaining
$\frac{l-1}{2}$ will be unstable limit cycles. For example, for the van der
Pol oscillator, Glycolytic oscillator, Modified Brusselator model etc. only
one limit cycle exists which is stable. For Kaiser model, $l=3$ and one can
observe the situations accordingly. So, for odd number of cycles innermost one
will be locally stable.

###  4.4 k-Cycle Cases:

####  4.4.1 A Model With N=1 and M=2k+1

For counting the number of limit cycles Gaikogaiko2008  has shown, for a
Liénard-type system i.e., LLS equation having the form,

$\displaystyle\ddot{x}-\left(\mu_{1}+\mu_{2}\dot{x}+\mu_{3}\dot{x}^{2}+\dots+\mu_{2k}\dot{x}^{2k-1}+\mu_{2k+1}\dot{x}^{2k}\right)\dot{x}+x$
$\displaystyle=0,$ (18)

can have atmost $k$ limit cycles if and only if, $\mu_{1}>0$. The
resultgaiko2008  correlates with our result. For any value of $k\in\mathbb{Z}$
it fits the odd-odd case of the general table and accordingly, $M$ and $N$ are
$2k+1$ and $1$, respectively, and finally the number of cycles will be atmost
$\frac{N+M}{2}-1=k$.

####  4.4.2 A Model With N=2k and M=1

Blows and Lloydblows ; perko  have stated that “For the Liénard or LLS system
$\dot{x}=y-F(x),\dot{y}=-g(x)$ with $g(x)=x$ and
$F(x)=a_{1}x+a_{2}x^{2}+\dots+a_{2k+1}x^{2k+1}$ has at most $k$ local limit
cycles and there are coefficients with $a_{1},a_{3},\dots,a_{2k+1}$ altering
in sign”. This can be found from the table-I with $N=2k$ and $M=1$ to give the
condition of atmost $k$ limit cycles. For example, taking $k=3$ with
$F(x)=-\epsilon(72x-\frac{392}{3}x^{3}+\frac{224}{5}x^{5}-\frac{128}{35}x^{7})$
has exactly three limit cycles for sufficiently small $\epsilon\neq 0$ which
are circles with radii $1,2$ and $3$. The above statement nicely corresponds
to the Theorem-6, pp-260perko .

Table-II: Table for highest degree polynomial $N+M(\oplus)$ for LLS system
together with maximum number of distinct conjugate roots$(R)$, with $1\leq
N,M\leq 10$

$\oplus,R$ |  |  |  |  | M |  |
---|---|---|---|---|---|---|---
|  | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10
| 1 | $2,0$ | $3,0$ | $4,2$ | $5,2$ | $6,4$ | $7,4$ | $8,6$ | $9,6$ | $10,8$ | $11,8$
| 2 | $3,2$ | $4,2$ | $5,4$ | $6,4$ | $7,6$ | $8,6$ | $9,8$ | $10,8$ | $11,10$ | $12,10$
| 3 | $4,2$ | $5,2$ | $6,4$ | $7,4$ | $8,6$ | $9,6$ | $10,8$ | $11,8$ | $12,10$ | $13,10$
| 4 | $5,4$ | $6,4$ | $7,6$ | $8,6$ | $9,8$ | $10,8$ | $11,10$ | $12,10$ | $13,12$ | $14,12$
N | 5 | $6,4$ | $7,4$ | $8,6$ | $9,6$ | $10,8$ | $11,8$ | $12,10$ | $13,10$ | $14,12$ | $15,12$
6 | $7,6$ | $8,6$ | $9,8$ | $10,8$ | $11,10$ | $12,10$ | $13,12$ | $14,12$ | $15,14$ | $16,14$
| 7 | $8,6$ | $9,6$ | $10,8$ | $11,8$ | $12,10$ | $13,10$ | $14,12$ | $15,12$ | $16,14$ | $17,14$
| 8 | $9,8$ | $10,8$ | $11,10$ | $12,10$ | $13,12$ | $14,12$ | $15,14$ | $16,14$ | $17,16$ | $18,16$
| 9 | $10,8$ | $11,8$ | $12,10$ | $13,10$ | $14,12$ | $15,12$ | $16,14$ | $17,14$ | $18,16$ | $19,16$
| 10 | $11,10$ | $12,10$ | $13,12$ | $14,12$ | $15,14$ | $16,14$ | $17,16$ | $18,16$ | $19,18$ | $20,18$

Counting the number of limit cycles through RG method in first order will give
similar result which was done by Das et. al.countinglcjkb ; infdampinglcjkb ;
lcbounestjkb  for some models. We have verified similar results for (3,3)
polynomial cases for ($F$,$G$) functions using RG method which become
increasingly very difficult and almost impossible upto (10,10) case than K-B
averaging method as tabulated in this work. It is very useful to count the
number of limit cycles from the table by just looking at the LLS form. For
example, the number of limit cycles of all models in Ref.countinglcjkb ;
infdampinglcjkb ; lcbounestjkb  along with the models in our work can be
estimated from our table. The table-I can also be utilized to prepare a model
of a desired number of limit cycles in a systematic way.

##  5 Conclusions

We have presented a scheme to cast a set of a class of coupled nonlinear
equations in two variables into a LLS form. By expressing the nonlinear
damping and forcing functions as polynomial we have implemented K-B method of
averaging to explore the number of admissible limit cycles of the dynamical
systems. Our results can be summarised as follows:

  1. 1.

For a LLS system, the number of limit cycles will be atmost $\frac{N+M-2}{2}$
when $N$ and $M$ degree of the polynomials for damping and restoring force
both are even or odd. Again, $\frac{N+M-1}{2}$ cycles can be found when $N$ is
even and $M$ is odd and finally, $\frac{N+M-3}{2}$ cycles when $N$ is odd and
$M$ is even.

  2. 2.

For a Liénard system, in particular, the formula of counting the number of
limit cycles follows the same with $M=1$ and $N\in\mathbb{Z^{+}}$. Also for
the generalised Rayleigh situation there occurs a linear restoring force so
that $N$ is $1$ and $M\in\mathbb{Z^{+}}$ .

  3. 3.

We have validated our general result with the help of a variety of physical
systems with one, two, three upto arbitrary k-cycles.

  4. 4.

This method stated in our work can also be utilized to prepare a model of a
desired number of limit cycles in a systematic way.

Acknowledgement

Sandip Saha acknowledges RGNF, UGC, India for the partial financial support.

Compliance with ethical standards
Conflicts of interests

The authors declare that there is no conflict of interests regarding the
publication of this paper.

##  6 Appendix: Lotka-Volterra System

To obtain the LLS form of Lotka-Volterra Systemepstein ; strogatz ; goldbook ,
let us set $z=\delta x+\beta y$ then $\dot{z}=\alpha\delta x-\beta\gamma y=u$
$\implies$ $x=\frac{\dot{z}+\gamma z}{(\alpha+\gamma)\delta}$ and
$y=\frac{-\dot{z}+\alpha z}{(\alpha+\gamma)\beta}$. After taking $t$
derivative upon $\dot{z}$ one can have,

$\displaystyle\ddot{z}$ $\displaystyle=(\alpha-\gamma)\dot{z}+\alpha\gamma
z+\frac{\dot{z}^{2}}{\alpha+\gamma}+\frac{\gamma-\alpha}{\alpha+\gamma}z\dot{z}-\frac{\alpha\gamma}{\alpha+\gamma}z^{2}.$

The fixed point $(0,0)$ gives a saddle solution which is not of any interest
in the present context. Choosing the remaining non-zero fixed point for
further investigations, and after taking perturbation $z=\xi+z_{s}$ around the
fixed point $z_{s}=\alpha+\gamma=\delta x_{s}+\beta y_{s}\neq 0$, one can get
the LLS form with $F(\xi,\dot{\xi})=a_{1}\xi+a_{2}\dot{\xi}$ with
$a_{1}=\frac{\alpha-\gamma}{\alpha+\gamma}$ and
$a_{2}=-\frac{1}{\alpha+\gamma}$. It is to be noted that $G(\xi)$ contains
nonlinearity with $G(\xi)=\omega^{2}\xi+a_{3}\xi^{2}$ where
$\omega=\sqrt{\alpha\gamma}=Im(\lambda)$($+ve$ sense) and
$a_{3}=\frac{\alpha\gamma}{\alpha+\gamma}$. After introducing a small
parameter $\epsilon_{1}$ (say) in the constants, $a_{i},b_{i}$ such that
$a_{i}=\epsilon_{1}b_{i},i=1,2,3$ the above equation reduces to
$\ddot{\xi}+\epsilon_{1}(b_{1}\xi+b_{2}\dot{\xi})\dot{\xi}+\omega^{2}\xi+\epsilon_{1}b_{3}\xi^{2}=0$.

## References

  * (1) J. D. Murray, Mathematical Biology.  Berlin: Springer, 1989.
  * (2) J. D. Murray, “1977lectures on non-linear differential equation models in biology,” 1977.
  * (3) I. R. Epstein and J. A. Pojman, An introduction to nonlinear chemical dynamics: oscillations, waves, patterns, and chaos.  Oxford University Press, New York, 1998.
  * (4) S. H. Strogatz, Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering.  Westview Press, USA, 1994.
  * (5) A. Goldbeter and M. J. Berridge, Biochemical Oscillations and Cellular Rhythms: The Molecular Bases of Periodic and Chaotic Behaviour.  Cambridge University Press, 1996.
  * (6) R. E. Mickens, Oscillations in planar dynamic systems, vol. 37.  World Scientific, 1996.
  * (7) S. L. Ross, Differential Equations.  Wiley, 1984.
  * (8) V. I. Arnold and Y. Ilyashenko, Ordinary Differential Equations, Encyclopedia Mathematical Science, 1st edn. Berlin: Springer, 1988.
  * (9) G. D. Birkhoff, Dynamical Systems.  Providence: A. M. S. Publications, 1927.
  * (10) J. Stoker, Nonlinear vibrations in mechanical and electrical systems.  Wiley Classics Library, Wiley, 1992.
  * (11) S. Smale, “Differentiable dynamical systems,” Bull. Am. Math. Soc., vol. 73, pp. 747–817, 1967.
  * (12) D. W. Jordan and P. Smith, Nonlinear Ordinary Differential Equations: An introduction for Scientists and Engineers, 4th edn. Oxford: Oxford University Press, 2007.
  * (13) A. I. Lavrova, E. B. Postnikov, and Y. M. Romanovsky, “Brusselator — an abstract chemical reaction?,” Physics-Uspekhi, vol. 52, no. 12, p. 1239, 2009.
  * (14) S. Ghosh and D. S. Ray, “Liénard-type chemical oscillator,” The European Physical Journal B, vol. 87, p. 65, Mar 2014.
  * (15) F. Kaiser, Theory of Resonant Effects of RF and MW Energy, pp. 251–282.  Boston, MA: Springer US, 1983.
  * (16) F. Kaiser and C. Eichwald, “Bifurcation structure of a driven, multi-limit-cycle van der pol oscillator (i): The superharmonic resonance structure,” International Journal of Bifurcation and Chaos, vol. 01, no. 02, pp. 485–491, 1991.
  * (17) P. Ghosh, S. Sen, S. S. Riaz, and D. S. Ray, “Controlling birhythmicity in a self-sustained oscillator by time-delayed feedback,” Phys. Rev. E, vol. 83, p. 036205, Mar 2011.
  * (18) R. H. Rand, “Lecture notes on nonlinear vibrations,” 2012.
  * (19) S. Ghosh and D. S. Ray, “Rayleigh-type parametric chemical oscillation,” The Journal of Chemical Physics, vol. 143, no. 12, p. 124901, 2015.
  * (20) L. Perko, Differential equations and dynamical systems (Second Edition), vol. 7.  Springer - Verlag, 1996.
  * (21) A. Liénard Rev. Gen. Electricite, vol. 23, p. 901, 1928.
  * (22) N. Levinson and O. K. Smith Duke Math. J., vol. 9, p. 382, 1942.
  * (23) N. Levinson Ann. Math., vol. 45, p. 723, 1944.
  * (24) S. Saha and G. Gangopadhyay, “Isochronicity and limit cycle oscillation in chemical systems,” Journal of Mathematical Chemistry, vol. 55, pp. 887–910, Mar 2017.
  * (25) S. Smale, “Mathematical problems for the next century,” The mathematical intelligencer, vol. 20, no. 2, pp. 7–15, 1998.
  * (26) V. A. Gaiko, “Limit cycles of liénard-type dynamical systems,” Cubo, vol. 10, pp. 115–132, 01 2008.
  * (27) D. Das, D. Banerjee, J. K. Bhattacharjee, and A. K. Mallik, “Counting limit cycles with the help of the renormalization group,” The European Physical Journal D, vol. 61, pp. 443–448, Jan 2011.
  * (28) D. Das, D. Banerjee, and J. K. Bhattacharjee, “Finding limit cycles in self-excited oscillators with infinite-series damping functions,” The European Physical Journal D, vol. 69, p. 85, Mar 2015.
  * (29) A. Dutta, D. Das, D. Banerjee, and J. K. Bhattacharjee, “Estimating the boundaries of a limit cycle in a 2d dynamical system using renormalization group,” Communications in Nonlinear Science and Numerical Simulation, vol. 57, pp. 47 – 57, 2018.
  * (30) N. M. Krylov and N. N. Bogolyubov, Introduction to non-linear mechanics.  Princeton Univ. Press, 1947.
  * (31) L. Y. Chen, N. Goldenfeld, and Y. Oono, “Renormalization group theory for global asymptotic analysis,” Phys. Rev. Lett., vol. 73, pp. 1311–1315, Sep 1994.
  * (32) L. Y. Chen, N. Goldenfeld, and Y. Oono, “Renormalization group and singular perturbations: Multiple scales, boundary layers, and reductive perturbation theory,” Phys. Rev. E, vol. 54, pp. 376–394, Jul 1996.
  * (33) A. Sarkar, P. Guha, A. Ghose-Choudhury, J. K. Bhattacharjee, A. K. Mallik, and P. G. L. Leach, “On the properties of a variant of the riccati system of equations,” Journal of Physics A: Mathematical and Theoretical, vol. 45, no. 41, p. 415101, 2012.
  * (34) H. Giacomini and S. Neukirch, “Number of limit cycles of the liénard equation,” Phys. Rev. E, vol. 56, pp. 3809–3813, Oct 1997.
  * (35) A. Sarkar, J. K. Bhattacharjee, S. Chakraborty, and D. B. Banerjee, “Center or limit cycle: renormalization group as a probe,” The European Physical Journal D, vol. 64, pp. 479–489, Oct 2011.
  * (36) H. G. E. Kadji, R. Yamapi, and J. B. Chabi Orou, “Synchronization of two coupled self-excited systems with multi-limit cycles,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 17, no. 3, p. 033113, 2007\.
  * (37) H. G. E. Kadji, J. B. Chabi Orou, R. Yamapi, and P. Woafo, “Nonlinear dynamics and strange attractors in the biological system,” Chaos, Solitons & Fractals, vol. 32, no. 2, pp. 862 – 882, 2007.
  * (38) R. Yamapi, B. R. Nana Nbendjo, and H. G. E. Kadji, “Dynamics and active control of motion of a driven multi-limit-cycle van der pol oscillator,” International Journal of Bifurcation and Chaos, vol. 17, no. 04, pp. 1343–1354, 2007.
  * (39) R. Yamapi, G. Filatrella, and M. A. Aziz-Alaoui, “Global stability analysis of birhythmicity in a self-sustained oscillator,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 20, no. 1, p. 013114, 2010\.
  * (40) A. Chéagé Chamgoué, R. Yamapi, and P. Woafo, “Dynamics of a biological system with time-delayed noise,” The European Physical Journal Plus, vol. 127, p. 59, May 2012.
  * (41) S. Saha and G. Gangopadhyay, “When an oscillating center in an open system undergoes power law decay,” Journal of Mathematical Chemistry, Nov 2018\.
  * (42) A. Lins, W. de Melo, and C. C. Pugh, “On liénard’s equation, lecture notes in mathematics,” 1977.
  * (43) G. S. Rychkov, “The maximal number of limit cycles of the system $\dot{y}=-x,\dot{x}=y-\sum\limits_{i=0}^{2}a_{i+1}x^{2i+1}$ is equal to two,” Differentsial’nye Uravneniya, vol. 11, no. 2, pp. 390–391, 1975.
  * (44) T. R. Blows and N. G. Lloyd, “The number of small-amplitude limit cycles of liénard equations,” Mathematical Proceedings of the Cambridge Philosophical Society, vol. 95, no. 2, p. 359–366, 1984.


11footnotetext:  MODESTYA Group, Department of Statistics, Mathematical
Analysis and Optimization, Universidade de Santiago de Compostela, Campus
Vida, Santiago de Compostela, Spain.22footnotetext:  Technological Institute
for Industrial Mathematics (ITMATI), Campus Vida, Santiago de Compostela,
Spain.

# Variable selection in Functional Additive Regression Models

Manuel Febrero–Bande1,2, Wenceslao González–Manteiga1,2
and Manuel Oviedo de la Fuente1,2

###### Abstract

This paper considers the problem of variable selection in regression models in
the case of functional variables that may be mixed with other type of
variables (scalar, multivariate, directional, etc.). Our proposal begins with
a simple null model and sequentially selects a new variable to be incorporated
into the model based on the use of distance correlation proposed by Székely et
al. (2007). For the sake of simplicity, this paper only uses additive models.
However, the proposed algorithm may assess the type of contribution (linear,
non linear, …) of each variable. The algorithm has shown quite promising
results when applied to simulations and real data sets.

##  1 Introduction

The variable selection problem in a general regression model tries to find the
subset of covariates that best predicts or explains a response. In the
classical approach, the covariates and the response are scalar (or
multivariate) and the model established among them is linear.

The stepwise regression, the most widely-used model selection technique
throughout the 80’s and the 90’s, is rooted in the classical papers by Akaike
(1973), Mallows (1973), Schwarz (1978) and Stone (1979). The main idea is to
use some diagnostic tools, directly derived from the linear model, to evaluate
the contribution of a new covariate and decide whether it should be included
in the model. The final subset is usually constructed using two main
strategies: the forward selection that begins with a simple null model and
tests, at each step, the inclusion of a new covariate in the model; and the
backward selection that starts with the full model including all candidate
variables and removes the most insignificant one at each step. It is also
possible to mix both strategies, testing at each step which variables can be
included or excluded in the optimal regression subset of covariates. In any
case, the stepwise regression is anchored in the diagnostics of the linear
model; it is therefore blind to detecting contributions other than the linear
one among the covariates and the response.

The work by Tibshirani (1996) proposing the LASSO estimator opens a new
direction on variable selection procedures. The main innovation of the LASSO
estimator is that it includes a $l_{1}$-type constraint for the coefficient
vector $\beta$ to force some parameters (components of $\beta$) to equal zero
and, thereby, obtains the optimal subset of covariates such as those with non-
zero coefficients. The effect of the constraint also helps in the optimization
step and satisfactorily deals with the sparsity phenomenon. See Zou (2006) for
a revision of the oracle properties of LASSO. Interesting examples may be
found in the literature following the same line but using several penalties,
constraints or using different structures in the regression models, such as:
LARS (Efron et al. (2002)), SCAD (Fan and Li (2001)), COSSO (Lin and Zhang
(2006)), additive models (Xue (2009)) and extensions to partial linear,
additive or semiparametric models like PLM (Du et al. (2012)), APLM (Liu et
al. (2011)) or GAPLM (Wang et al. (2011)). All these works have two common
characteristics: each paper is based on a specific model and all the
covariates must be included in the model at the same time. The latter could
lead to highly demanding computing algorithms that are difficult to implement
at times. In particular, for high-dimensional or functional data problems, the
previous steps that commonly include variable standardization and/or variable
representation may notably increase the complexity and the cost of the
algorithms. See, for example, Hastie et al. (2015) for a review of some of the
aforementioned methods.

Another type of strategy is a pure feature selection where the covariate is
selected without a model. This is the approach employed in mRMR (minimum
Redundancy Maximum Relevance, Peng et al. (2005)). To enter into the model, a
new candidate covariate must have a great relevancy with the response while
maintaining a lower redundancy with the covariates already selected in the
model. The main advantage of this approach is that it is an incremental rule;
once a variate has been selected, it cannot be deselected in a later step. On
the other hand, the measures for redundancy and relevancy must be chosen in
function of the regression model applied to ensure good predictive results in
the final model. The FLASH method proposed in Radchenko and James (2011) is a
modification of the LASSO technique that sequentially includes a new variate
changing the penalty at each step. This greedy increasing strategy is also
employed by Boosting (see, for example, Bühlmann and Yu (2003) or Ferraty and
Vieu (2009) in a functional data context). Boosting is not a purely feature
selection method but rather a predictive procedure that selects at each step
the best covariate/model with respect to the unexplained part of the response.
The final prediction is constructed as a combination of the different steps.
All the previous solutions are not completely satisfactory in a functional
data framework, specially when the number of possible covariates can be
arbitrarily large. Specifically, we are interested in an automatic regression
procedure capable of dealing with a large number of covariates of different
nature, possibly very closely related to one another.

Figure 1: Price and Energy Market Demand curves negotiated in the Iberian
Energy Market. Period: Jan, 2014- Dec, 2014 (Source: omie.es).

Our motivating example comes from the energy market. Figure 1 shows the daily
profile of Price and Energy (Electricity) Market Demand (both measured hourly)
at the Iberian Energy Market from Jan, 2014 to Dec, 2014 with a color code,
red and green, that considers, respectively, two groups: M-F (Monday to
Friday) and WEnd (Saturday and Sunday & Sun). We are interested in predicting
the price or the demand at certain hours of the following day. We dispose of
many different sources of information to do this. These sources may include
other variables related with energy market or generation, meteorological
information, calendar effects or any transformation/filter of the preceding.
Figures 2 and 3 show a small sample of selected covariates. Figure 2 includes
some of the energy generation variates (by type) included in the energy
generation pool in the same period. The total demand can be decomposed as a
sum of the different energy generations covering the demand. Figure 3 is a
small sample of calendar and meteorological information. In both figures, the
color code is the same as in Figure 1. Note that the number of variables that
may be included in the model is rather high (over $150$).

Figure 2: Daily profile curves of amount of generated energy by type. Period:
Jan, 2014- Dec, 2014 (Source: demanda.ree.es). Red and green colors
represents, respectively, M-F and W-End days.

The curves in Figure 2 are the components of the energy pool by generation
type, which may be useful for predicting price or demand. Due to the Iberian
Energy Market regulations, the odds of becoming part of the final energy pool
consumed due to its price or availability are unalike among all the types of
energy. For instance, hydroelectric energy is only offered to the market when
the price is high in the presence of founded expectations on refilling the
reservoir (using the weather forecasts).

Figure 3: Examples of available scalar covariates for the energy prediction
problem. Period: Jan, 2014- Dec, 2014 (Source: aemet.es). Red and green colors
represents, respectively, M-F and W-End days.

The scalar covariates for our prediction problem can be linked to calendar
effects (like month or day-of-week), meteorological information (like
temperature or wind speed) and transformations from the functional variables
(such as the market demand at midday
$X(t)\mathbbm{1}_{\left\\{t=12h\right\\}}$). This means that many new
variables can be created from the original ones (for instance, by using
derivatives or considering certain subintervals), many of which may certainly
be closely related.

Our aim is to select significant covariates for a general additive regression
model with scalar response:

$Y_{i}=\alpha+\sum_{j=1}^{J}f_{j}\left(X_{i}^{(j)}\right)+\varepsilon_{i},\quad
i=1,\ldots,N$ (1)

where the covariates are chosen from the set
$S=\left\\{X^{1},X^{2},\cdots,X^{k},\cdots\right\\}$ of different potential
covariates (functional, vectorial, …). The notation $X^{(j)}$ refers to the
$j$-th covariate selected for the model. The number of variates can be
extraordinarily large, so we intent to construct the regression model
sequentially, i.e. from the trivial model up to the one that includes all the
useful information provided by the covariates in the set.

The rest of this paper is structured as follows. Section 2 presents our
proposal, Section 3 presents some artificial examples and simulation studies
and Section 4 includes the applications to our motivating example.

##  2 The procedure

Our proposal borrows some ideas from other variable selection methods. First,
we separate the selection process from the estimation step by picking the
possible candidates using a measure that is not based in a concrete model. The
method for selecting the covariate must be general so that the model may
include scalar, multivariate, directional or functional covariates in the
model on an impartial way. A model is chosen (among those in a catalog) that
can be applicable to the selected covariates and the response, ensuring that
the model is correct, i.e. it fulfills the considered structural hypothesis
and its performance is superior to the previous one. Typically, we must decide
the contribution of the new variable jointly with the selected ones in
precedent steps. The new model is compared with the previous one using an
appropriate test to decide whether to definitely incorporate the new variate
into the model. If the new model is not better than the previous one, the
candidate is discarded and the current model remains invariant. If the new
model is better then the new model becomes the current model. Finally, using
the residuals of the current model, new candidates are trying following the
same steps. This procedure is repeated until no more variables enter the
model. The idea of using residuals makes that, in subsequent steps, the
algorithm attempts to capture information not chosen by previous steps. This
also sharply contrasts with the rules of a pure feature method. In a pure
feature method it is necessary to compute a measure of redundancy among the
covariates in the model and the candidates to be incorporated that can be
quite difficult due to different nature of the variables.

As aforementioned, the key idea for the selection of a feature is to find a
tool that can be homogeneously applied to variates of different nature. This
paper is based on the exhaustive use of the distance correlation,
$\mathcal{R}$, proposed by Székely et al. (2007), and recently reviewed in
Székely and Rizzo (2017). The distance correlation fulfills the following two
conditions:

  * i)

$\mathcal{R}(X,Y)$ is defined for $X$ and $Y$ random vector variables in
arbitrary finite dimension spaces.

  * ii)

$\mathcal{R}(X,Y)=0$ characterizes independence of $X$ and $Y$.

Both conditions that are established in Székely et al. (2007), signal the
presence of a way to measure the relationship among $X$ and $Y$ homogeneously
for arbitrary dimensions of these vectors. Indeed, the work by Lyons (2013)
extends these properties for metric spaces of a certain type (_strong negative
type_). In particular, Hilbert spaces or any space that can be embedded in a
Hilbert space are of strong negative type. Along this paper, we will only use
Hilbert spaces for all covariates (for scalar covariates this means that we
are using the classical Euclidean spaces). As $\mathcal{R}$ characterizes
independence, the distance correlation can detect relationships among
variables other than the linear one. A test for independence using the
distance correlation can be derived following Székely et al. (2007) and
Székely and Rizzo (2013) that allows us to contrast when a variable could be a
reasonable candidate to enter the model.

In Yenigün and Rizzo (2015), a definition of partial distance correlation
(PDC) among $X$ and $Y$ given $Z$ was introduced based on computing the
distance correlation among the residuals of two models: $Y$ respect to $Z$ and
$X$ respect to $Z$. The PDC is the key of that paper for a variable selection
method although has two important drawbacks: First, PDC is constructed under
linear relationship assumptions among variables and, in the original paper, it
is only applied to classical linear models with scalar variates (the authors
called it “linear PDC”). Second, its implementation is done using the function
pdcor of R–package energy (Rizzo and Székely (2014)) which only uses the
distance matrices among elements of $X$, $Y$ and $Z$. Specifically, $Z$ (the
variables already in the model) could be a mix of functional, scalar or
multivariate variables where an appropriate distance using all of them must be
hard to compute. Even restricting ourselves to the scalar case, those
variables should have a similar scale. On the contrary, our proposal can be
applied to any mix of variables given that in every step, the distance
correlation is computed among the residuals of the current model with each
candidate. Taking into account that the residuals have the same nature as the
response variable, the distance correlation can always be computed at each
step, supposing that we were able to compute $\mathcal{R}(X,Y)$ at first
iteration.

The computation of $\mathcal{R}(X,Y)$ for a given sample (empirical distance
correlation) is quite straightforward because it only depends on the distances
among data. See again, Székely et al. (2007) for details. Specifically,
supposing we are interested in computing the distance correlation among
$X^{j}$ and $Y$, let $a_{kl}=d(X_{k}^{j},X_{l}^{j})$, the distance in the
Hilbert space of $X^{j}$ among cases $k$ and $l$, and
$A_{kl}=a_{kl}-\bar{a}_{k\cdot}-\bar{a}_{\cdot l}+\bar{a}_{\cdot\cdot}$ and
respectively, $b_{kl}=d(Y_{k},Y_{l})$, and
$B_{kl}=b_{kl}-\bar{b}_{k\cdot}-\bar{b}_{\cdot l}+\bar{b}_{\cdot\cdot}$.

The distance correlation is simply computed as:

$\mathcal{R}(X^{j},Y)=\frac{\sum_{k}\sum_{l}A_{kl}B_{kl}}{\sqrt{\sum_{k}\sum_{l}A_{kl}^{2}}\sqrt{\sum_{k}\sum_{l}B_{kl}^{2}}}$

The simplicity of the $\mathcal{R}(X^{j},Y)$ has an important drawback in
terms of computational time. As it is based on distances, the memory and the
number of necessary operations is of order $N(N-1)/2$ where $N$ is the sample
size. So, for large sample sizes, the computation must be divide into some
small tasks to avoid memory overruns. See the end of section 4 for details on
how to solve this issue.

###  2.1 The algorithm

Our proposal can be formalized as follows:

  1. 1.

Let $Y$ the response and $S=\left\\{X^{1},\ldots,X^{p}\right\\}$ the set of
all variables that can be included in the model.

  2. 2.

Set $\hat{Y}=\bar{Y}$, and let $M^{(0)}=\emptyset$ the initial set of the
variates included in the model. Set $i=0$.

  3. 3.

Compute the residuals of the current model: $\hat{\varepsilon}=Y-\hat{Y}$.

  4. 4.

Choose $X^{j}\in S\neq\emptyset$ such that: 1)
$\mathcal{R}\left(\hat{\varepsilon},X^{j}\right)\geq\mathcal{R}\left(\hat{\varepsilon},X^{k}\right),\forall
k\neq j\in S$ and 2) the null hypothesis for the test of independence among
$\left\\{X^{j}\right\\}$ and $\hat{\varepsilon}$ is rejected. IF NOT, END.

  5. 5.

Update the sets $M$ and $S$: $M^{(i+1)}=M^{(i)}\cup\left\\{X^{j}\right\\}$,
and $S=S\backslash\left\\{X^{j}\right\\}$.

  6. 6.

Compute the new model for $Y$ using $M^{(i+1)}$ choosing the best contribution
of the new covariate. Typically, there will be a catalog of all possible ways
of constructing correct models with the variates in $M^{(i+1)}$ fixing the
contributions of the variates in $M^{(i)}$ and adding the new one.

  7. 7.

Analyze the contribution of $X^{j}$ in the new model respect to the current:
IF this contribution is not relevant (typically comparing with the current
model)
THEN $M^{(i+1)}=M^{(i+1)}\backslash\left\\{X^{j}\right\\}$ and the current
model remains unalterable
ELSE the new model becomes the current model and provides new predictions
($\hat{Y}$). Along this paper we have employed an additive model:
$\hat{Y}=\bar{Y}+\sum_{m\in M}\hat{f}_{m}\left(X^{(m)}\right)$ where at each
step $\hat{f}_{m}$ could be linear or nonlinear.

  8. 8.

Update the number of iterations: $i=i+1$ and go to 3

  9. 9.

END. The current model is the final model with the variates included in
$M^{(i)}$. $S$ is either the empty set or contains those variables that accept
the null hypothesis of the test of independence respect to the residuals of
the current model.

Steps 1–3 establish the null model as the initial model for beginning the
procedure. Step 4 selects the variable from the set of available ones $S$. A
new candidate is selected when maximizes the distance correlation among
available ones and the test of independence between $\hat{\varepsilon}$ and
the candidate rejects the null hypothesis. For simplicity, every time a
distance correlation is showed, the value is filtered by the test of
independence, i.e. all tables show
$\mathcal{R}\left(X^{k},\hat{\varepsilon}\right)\mathbbm{1}_{\left\\{H_{1}\right\\}}$
that only takes values distinct from zero for those covariates that rejects
the null hypothesis of independence. From the practical point of view, the
algorithm spends most of its time computing the latter quantity for each
variable in this step. This task is not so fast because the distance matrices
for a sample size $N$ require $N(N-1)/2-N$ stored positions; so, both the
number of operations and the amount of consumed memory are quadratic functions
of the sample size. We may alleviate the computational burden using two simple
tricks: first, we may easily parallelize this task to into small tasks that
compute the distance correlation separately for different covariates in $S$ or
even, when the distance matrix is too large, the computation of $\mathcal{R}$
can be done through all the submatrices below the diagonal. Secondly, we can
compute the distance matrices needed for the algorithm in advance for all
covariates avoiding a repetitive calculation of these matrices for each
iteration.

Step 5 updates the sets of variates. The selected variable is included as a
possible candidate in set $M$ and it is removed from set $S$. So, each
candidate is only checked once in a clear forward strategy.

Step 6 tries the possible models that can be conformed with the variables
included in the updated $M$. Here, the previous contributions of the variables
already in the model remain fixed and the catalog is used to check what could
be the contribution of the new candidate. For instance, it may happen that my
catalog of models only contains linear models or models with a limited number
of nonlinear terms. In the latter case, if that limit is reached in previous
steps, the new candidate can only be included with a linear contribution. This
means that the catalog of possible models may not suffice to explain the
relationship shown by the distance correlation. As a diagnostic tool, in
latter iterations the test of independence can be computed for those variates
in $M$ to analyze if the contribution of covariate already in the model is
collecting all the possible information of that covariate. If the test is
rejected for any of the variables of $M$, it is an evidence that the model is
not completely exploiting its information. Please note that this diagnostic
tool is never done in step 5 to avoid an infinite loop. Of course, a final
arrangement is leaved to the user when the procedure ends. For instance, a
certain variate entering the model at an earlier stage, may become irrelevant
some iterations later, when other variates add their contributions to the
model. Of course, a final diagnostic of the model could lead to changes in the
contribution of the covariates. Yet the comparison among the different
variables must be carefully analyzed given its nature. Of course, the
procedure could be modified to include backward steps as in the classical
stepwise method but, also due the distinct nature of the covariates, it could
be difficult to find an homogenize backward step that works in an universal
way. To avoid this type of problems, we have only applied additive models in
this paper because this type of model is available for (nearly) any kind of
variates and is quite flexible taking into account that each contribution may
be linear or nonlinear (smooth). For instance, if we want to include
iterations among variables, we can define a new covariate computing that
interaction to be included as a candidate in the additive model.

Step 7 analyzes the contribution of the new candidate in the model with
respect to the previous one. Typically, depending on the model, this can be
done with a Generalized Likelihood Ratio Test. If the contribution of the new
candidate is not relevant, the variate is removed from set $M$ and the
previous model remains as the current model trying with another candidates. If
the new model is better than the previous one, the new residuals are computed
and the algorithm backs to step 3 to recompute the distance correlation among
these residuals and the variables remaining in set $S$. The procedure ends
when no more variables can be added to the model because $S$ is empty or all
variables in $S$ accepts the independence null hypothesis.

##  3 Numerical results

In the general framework, there is no previous papers that can be considered
direct competitors of our proposal. But, as commented before, Yenigün and
Rizzo (2015) proposes an algorithm valid for scalar variates under linear
relationships with the response. We have compared our proposal in the same
scenarios of that paper. As in the original paper, samples of size $N=100$ are
generated from different regression models and with several distributions for
the $p=8$ candidate covariates (typically combinations of $N(0,\sigma^{2})$
and $U[-a,a]$). This process was repeated $B=500$ times and the number of
times every single covariate was selected is recorded as well as the Root Mean
Squared Prediction Error:

$RMSPE=\frac{1}{B}\sum_{i=1}^{B}\sqrt{\frac{1}{n}\sum_{j=1}^{n}(Y_{j}^{(i)}-\hat{Y}_{j}^{(i)})^{2}}$

for the final model using a test sample also of size $n=100$. The scenarios
are defined as follows:

  1. YR1:

$Y=Z_{1}+Z_{2}+Z_{3}+\varepsilon,\quad Z_{i}\sim N(0,1),\;\varepsilon\sim
N(0,\sigma^{2}=2^{2})$

  2. YR2:

$Y=\log[4+\sin(3Z_{1})+\sin(Z_{2})+Z_{3}^{2}+Z_{4}+0.1\varepsilon],\quad
Z_{1}\sim N(0,1),\;Z_{2}\sim N(0,2^{2}),\\\ Z_{3}\sim
U[-1.5,1.5],\;Z_{4},\ldots,Z_{8}\sim U[-1,1]$ and $\varepsilon\sim N(0,1)$

  3. YR3:

$Y=\left|Z_{1}\right|+Z_{2}^{2}+Z_{3}^{2},\quad Z_{1}\sim
N(0,1.4^{2}),\;Z_{2}\sim U[-1.7,1.7],\;Z_{3}\sim N(0,.8^{2}),\\\
Z_{4},\ldots,Z_{8}\sim N(0,1)$.

  4. YR4:

$Y=Z_{1}+Z_{2}+Z_{3}+\varepsilon,\;Z\sim N_{p}(0,\Sigma)$ with
$\Sigma=\left(\begin{array}[]{cccc}1&\theta&\cdots&\theta\\\
\theta&1&\cdots&\theta\\\ \vdots&\vdots&\ddots&\vdots\\\
\theta&\theta&\cdots&1\end{array}\right)$ where $\theta=0.6$ and
$\varepsilon\sim N(0,2^{2})$.

  5. YR5:

$Y=Z_{1}+Z_{2}+Z_{3}+\varepsilon,\;\\\ Z\sim N_{p}(0,\Sigma)$ with
$\Sigma=\left(\begin{array}[]{ccccc}1&\theta&\theta^{2}&\cdots&\theta^{p-1}\\\
\theta&1&\theta&\cdots&\theta^{p-2}\\\
\theta^{2}&\theta&1&\cdots&\theta^{p-3}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
\theta^{p-1}&\theta^{p-2}&\theta^{p-3}&\cdots&1\end{array}\right)$
where $\theta=0.6$ and $\varepsilon\sim N(0,2^{2})$.

For these scenarios, we have compared the proposal of Yenigün and Rizzo (2015)
(PDC) with ours in two versions: restricting the model to be have only linear
contributions (LM) or letting every contribution to be nonlinear (AM). The
results are provided in Table 1. In the independent linear case YR1, the three
methods show a similar performance both in variable identification and in
terms of prediction error. In the second scenario, YR2, the constraint to deal
only with linear models produces high values of RMSPE for PDC and LM on the
contrary of the AM case that mimics the theoretical $\sigma$. Indeed, $Z_{1}$
is not identified the same proportion of times than $Z_{2},Z_{3}$ or $Z_{4}$
except in the case of AM. The third scenario, YR3, shows similar results than
YR2 although here the difference among linear versions and the nonlinear one
is now more important in terms of RMSPE. Recall that these two scenarios are
nonlinear. Scenarios YR4 and YR5 are linear scenarios with dependence among
covariates. In both cases, it seems that our proposals (LM or AM) identify the
relevant covariates less accuracy than PDC but really, PDC tends to include
more variates than necessary. In YR4, the minimum of proportion of times that
$Z_{4},\ldots,Z_{8}$ enters the model is $0.436$ where for LM and AM the
maximum that a irrelevant covariate enters the model is $0.044$. Something
similar happens in YR5 although due to the dependence structure, only $Z_{4}$
enters the model a high proportion of times with the PDC. The difference is
that YR4 considers a constant correlation among all covariates but in YR5,
only $Z_{4}$ shows an important correlation with $Z_{3}$ (the correlation
decreases respect to the column number). Taking into account the RSMPE, it
seems that LM and AM do not need to select all the relevant variables to
achieve similar prediction levels. On the contrary, PDC needs to include some
of the irrelevant ones as part of the model. Intuitively, a high dependence
among covariates leads to that not all the relevant variates may enter the
model because, and due to the dependence structure, its contribution may have
been incorporated into the model by other highly related. So, as a conclusion
of this small simulation study, our general proposal AM works well in linear
scenarios, clearly better in nonlinear ones and in the case of dependence
among covariates reduces the amount of irrelevant variates that are introduced
in the model without losing predictive performance.

Table 1:  Proportion of times that each variable enters the model for scenarios YR1–YR5. The last column contains the Root Mean Square Prediction Error for $500$ repetitions. | $Z_{1}$ | $Z_{2}$ | $Z_{3}$ | $Z_{4}$ | $Z_{5}$ | $Z_{6}$ | $Z_{7}$ | $Z_{8}$ | RMSPE
---|---|---|---|---|---|---|---|---|---
YR1–PDC | 0.968 | 0.968 | 0.976 | 0.078 | 0.096 | 0.082 | 0.078 | 0.080 | 2.070
YR1–LM | 0.990 | 0.996 | 0.998 | 0.050 | 0.096 | 0.076 | 0.054 | 0.086 | 2.060
YR1–AM | 0.986 | 0.994 | 0.994 | 0.046 | 0.084 | 0.066 | 0.056 | 0.078 | 2.110
YR2–PDC | 0.658 | 0.980 | 0.886 | 0.996 | 0.066 | 0.060 | 0.082 | 0.080 | 0.280
YR2–LM | 0.776 | 0.996 | 0.952 | 0.994 | 0.050 | 0.064 | 0.068 | 0.062 | 0.280
YR2–AM | 1.000 | 1.000 | 1.000 | 1.000 | 0.052 | 0.054 | 0.078 | 0.066 | 0.090
YR3–PDC | 1.000 | 0.996 | 0.948 | 0.068 | 0.080 | 0.064 | 0.084 | 0.072 | 1.560
YR3–LM | 0.994 | 0.994 | 0.934 | 0.074 | 0.064 | 0.060 | 0.076 | 0.068 | 1.560
YR3–AM | 1.000 | 1.000 | 1.000 | 0.076 | 0.054 | 0.056 | 0.048 | 0.078 | 0.060
YR4–PDC | 0.992 | 0.988 | 0.998 | 0.436 | 0.476 | 0.452 | 0.464 | 0.468 | 2.060
YR4–LM | 0.766 | 0.792 | 0.772 | 0.044 | 0.040 | 0.040 | 0.022 | 0.034 | 2.130
YR4–AM | 0.754 | 0.776 | 0.764 | 0.038 | 0.038 | 0.038 | 0.014 | 0.032 | 2.180
YR5–PDC | 0.998 | 1.000 | 0.998 | 0.350 | 0.042 | 0.014 | 0.020 | 0.016 | 2.040
YR5–LM | 0.876 | 0.856 | 0.820 | 0.082 | 0.058 | 0.034 | 0.046 | 0.050 | 2.100
YR5–AM | 0.876 | 0.862 | 0.808 | 0.086 | 0.060 | 0.044 | 0.038 | 0.062 | 2.150

For assessing the situation when the variables are of different nature, we
have developed a simulation study to check the performance of the algorithm in
a mixed scenario with functional and scalar variables. Five functional and
five scalar variables were simulated and the response was constructed as a
function of the two first functional and the two first scalar variables. The
functional variables:
$\left\\{\mathcal{X}_{1},\ldots,\mathcal{X}_{5}\right\\}$ are generated
following independent Ornstein-Uhlenbeck processes in $[0,1]$, and the scalar
variables $\left\\{Z_{1},\ldots,Z_{5}\right\\}$ following, respectively,
$U[0,1]$, $N(0,1)$, $N(0,1)$, $U[0,1]$ and $N(0,1)$ independently of each
other. We constructed the response as follows:

$Y=10+a_{1}\left<\mathcal{X}_{1},\beta_{1}\right>+a_{2}\left|\left|\mathcal{X}_{2}^{2}\right|\right|+3a_{3}Z_{1}+a_{4}Z_{2}^{2}+\varepsilon$

with $\beta_{1}=2t+\sin{4\pi t+0.1},t\in[0,1]$ and $\varepsilon\sim
N(0,.25^{2})$.

The coefficients $a=\left\\{a_{1},a_{2},a_{3},a_{4}\right\\}$ were introduced
to emphasize/mask each part of the model. We estimated the model through a
Functional Additive Model (Müller and Yao (2008)) using the first four
principal components for functional covariates and a standard additive model
for the scalar ones. Samples with $N=200$ were generated and the process was
repeated $B=500$ times to count the proportion that a particular covariate
enters the model. Table 2 shows the excellent results for different
combinations of $\left\\{a_{1},a_{2},a_{3},a_{4}\right\\}$. In all cases,
$\mathcal{X}_{1}$, $Z_{1}$ and $Z_{2}$ were selected all the times in the
$500$ repetitions. $\mathcal{X}_{2}$ was also selected all the times except
for $a_{2}=\frac{1}{8}$ that corresponds to a small contribution of that
variate to the model. This can be explained because the
$\textrm{Var}\left(0.125\left|\left|\mathcal{X}_{2}^{2}\right|\right|\right)\approx
0.01$ that it is about $\frac{1}{6}$ of the residual variance. The non
relevant variables enter the model about 5% of the times.

Table 2: % of times that the variate was included in the model for $B=500$ replications. Independent case. $a$ | $\mathcal{X}_{1}$ | $\mathcal{X}_{2}$ | $\mathcal{X}_{3}$ | $\mathcal{X}_{4}$ | $\mathcal{X}_{5}$ | $Z_{1}$ | $Z_{2}$ | $Z_{3}$ | $Z_{4}$ | $Z_{5}$
---|---|---|---|---|---|---|---|---|---|---
$\left\\{1,1,1,1\right\\}$ | 1.000 | 1.000 | 0.046 | 0.038 | 0.030 | 1 | 1 | 0.030 | 0.040 | 0.054
$\left\\{\frac{1}{2},\frac{1}{2},\frac{1}{2},\frac{1}{2}\right\\}$ | 1.000 | 1.000 | 0.048 | 0.044 | 0.046 | 1 | 1 | 0.060 | 0.060 | 0.046
$\left\\{\frac{1}{4},\frac{1}{4},\frac{1}{4},\frac{1}{4}\right\\}$ | 1.000 | 0.984 | 0.040 | 0.044 | 0.058 | 1 | 1 | 0.058 | 0.040 | 0.060
$\left\\{\frac{1}{4},\frac{1}{4},1,1\right\\}$ | 1.000 | 0.988 | 0.056 | 0.030 | 0.052 | 1 | 1 | 0.058 | 0.044 | 0.052
$\left\\{1,1,\frac{1}{4},\frac{1}{4}\right\\}$ | 1.000 | 1.000 | 0.058 | 0.052 | 0.042 | 1 | 1 | 0.048 | 0.064 | 0.056
$\left\\{\frac{1}{8},\frac{1}{8},\frac{1}{8},\frac{1}{8}\right\\}$ | 0.908 | 0.474 | 0.048 | 0.050 | 0.044 | 1 | 1 | 0.058 | 0.060 | 0.060

We have repeated the simulation study but forcing that some irrelevant
covariates ($\mathcal{X}_{3}$, $\mathcal{X}_{4}$, $\mathcal{X}_{5}$, $Z_{3}$,
$Z_{4}$, $Z_{5}$) may have a strong relationship with any of the important
ones ($\mathcal{X}_{1}$, $\mathcal{X}_{2}$, $Z_{1}$, $Z_{2}$) to check how the
procedure approaches covariates with strong collinearities.

To this end, we computed new covariates $\mathcal{X}_{3}^{*}$,
$\mathcal{X}_{4}^{*}$, $Z_{3}^{*}$ and $Z_{4}^{*}$ in the following way:
$\mathcal{X}_{3}^{*}=0.95\mathcal{X}_{1}+0.05\mathcal{X}_{3}$,
$\mathcal{X}_{4}^{*}=0.95\mathcal{X}_{2}+0.05\mathcal{X}_{4}$,
$Z_{3}^{*}=0.95Z_{1}+0.05Z_{3}$ and $Z_{4}^{*}=0.95Z_{2}+0.05Z_{4}$. This
ensures that $\mathcal{X}_{3}^{*}$ (respectively, $\mathcal{X}_{4}^{*}$,
$Z_{3}^{*}$, $Z_{4}^{*}$) has nearly the same information as $\mathcal{X}_{1}$
(respectively, $\mathcal{X}_{2}$, $Z_{1}$, $Z_{2}$). Table 3 shows the
proportion of times that every variate is included in the model for this new
scenario. Now the proportion of times shown in Table 2 for any of the relevant
variables is shared among that variable and its (almost) copy. The important
message is that when a variate is included in the model, its copy is banned
for consequent steps and the model can be safely estimated without redundant
information. Typically now, the sum of proportion of times of each covariate
and its (almost) copy in Table 3 is the proportion of times that can be seen
in Table 2 for the covariate.

Table 3: % of times that the variate was included in the model for $B=500$ replications under strong collinearities among covariates. $a$ | $\mathcal{X}_{1}$ | $\mathcal{X}_{2}$ | $\mathcal{X}_{3}^{*}$ | $\mathcal{X}_{4}^{*}$ | $\mathcal{X}_{5}$ | $Z_{1}$ | $Z_{2}$ | $Z_{3}^{*}$ | $Z_{4}^{*}$ | $Z_{5}$
---|---|---|---|---|---|---|---|---|---|---
$\left\\{1,1,1,1\right\\}$ | 0.584 | 0.644 | 0.416 | 0.356 | 0.046 | 0.774 | 0.574 | 0.228 | 0.426 | 0.038
$\left\\{\frac{1}{2},\frac{1}{2},\frac{1}{2},\frac{1}{2}\right\\}$ | 0.550 | 0.640 | 0.450 | 0.360 | 0.048 | 0.738 | 0.580 | 0.262 | 0.420 | 0.052
$\left\\{\frac{1}{4},\frac{1}{4},\frac{1}{4},\frac{1}{4}\right\\}$ | 0.580 | 0.594 | 0.420 | 0.396 | 0.032 | 0.724 | 0.538 | 0.276 | 0.462 | 0.060
$\left\\{\frac{1}{4},\frac{1}{4},1,1\right\\}$ | 0.558 | 0.564 | 0.442 | 0.416 | 0.048 | 0.846 | 0.584 | 0.154 | 0.416 | 0.060
$\left\\{1,1,\frac{1}{4},\frac{1}{4}\right\\}$ | 0.592 | 0.684 | 0.408 | 0.316 | 0.040 | 0.662 | 0.540 | 0.338 | 0.460 | 0.072
$\left\\{\frac{1}{8},\frac{1}{8},\frac{1}{8},\frac{1}{8}\right\\}$ | 0.484 | 0.250 | 0.406 | 0.222 | 0.038 | 0.664 | 0.502 | 0.336 | 0.496 | 0.048

Table 4 shows the distance correlation with the residuals for a single run at
each iteration of the procedure under independence among covariates. The first
column of this table is the iteration step. This allows us to analyze the
order in which the covariates enter the model: $Z_{1}$ ($0.177$), $Z_{2}$
($0.141$), $\mathcal{X}_{1}$ ($0.22$), $\mathcal{X}_{2}$ ($0.103$). In all
cases, $\mathcal{R}(X,\hat{\varepsilon})$ increases with the iterations until
that covariate is chosen. This is an expected behaviour in a case in which all
the covariates are independent of each other and the relationship with the
residual becomes stronger as previous effects of other covariates are removed.
Also, if we compute
$\mathcal{R}(X,\hat{\varepsilon})\mathbbm{1}_{\left\\{\mathcal{H}_{1}\right\\}}$
for $X\in M$, we obtain zeros suggesting that the model was able to
incorporate all the important information of that covariates. In this case,
the algorithm ends in iteration five when no more covariates are added to the
model because none of the remaining ($\mathcal{X}_{3}$, $\mathcal{X}_{4}$,
$\mathcal{X}_{5}$, $Z_{3}$, $Z_{4}$, $Z_{5}$) rejects the null hypothesis of
independence with respect to the residuals.

Table 4:  $\mathcal{R}(X,\hat{\varepsilon})\mathbbm{1}_{\left\\{\mathcal{H}_{1}\right\\}}$ for one run with $a=\left\\{1,1,1,1\right\\}$. At each iteration, a variable enters the model when maximizes the row. Empty cells correspond to variates selected in previous steps. | $\mathcal{X}_{1}$ | $\mathcal{X}_{2}$ | $\mathcal{X}_{3}$ | $\mathcal{X}_{4}$ | $\mathcal{X}_{5}$ | $Z_{1}$ | $Z_{2}$ | $Z_{3}$ | $Z_{4}$ | $Z_{5}$
---|---|---|---|---|---|---|---|---|---|---
1 | 0.040 | 0.029 | 0.000 | 0.000 | 0.000 | 0.177 | 0.108 | 0.000 | 0.000 | 0.000
2 | 0.122 | 0.049 | 0.012 | 0.000 | 0.000 |  | 0.141 | 0.000 | 0.000 | 0.000
3 | 0.220 | 0.062 | 0.000 | 0.000 | 0.000 |  |  | 0.000 | 0.000 | 0.000
4 |  | 0.103 | 0.014 | 0.000 | 0.000 |  |  | 0.000 | 0.000 | 0.000
5 |  |  | 0.000 | 0.000 | 0.000 |  |  | 0.000 | 0.000 | 0.000

The following example is related with a classification problem and so, the
response is now binomial. We have included also an important number of
needless variates to check how the algorithm deals with such situations.
Consider $50$ iid variates distributed as $U[-1,1]$ and named, respectively,
$\\{X_{1},X_{2},Z_{1},\ldots,Z_{48}\\}$. The change in the names is to
emphasize that only the first two variates are used in the true classification
rule. Considering only the data points $X_{1}$ and $X_{2}$ belonging to the
circular crown with $r=0.6$ and $R=1$ and its norm in $\mathbb{R}^{2}$,
$\left|\left|(X_{1},X_{2})\right|\right|_{2}=\sqrt{X_{1}^{2}+X_{2}^{2}}$, we
conform the first group with all data points that
$0.6\leq\left|\left|(X_{1},X_{2})\right|\right|_{2}\leq 0.8$ and the second
fulfilling a similar condition
$\left|\left|(X_{1},X_{2})\right|\right|_{2}>0.8$ as shown in Figure 4 that
shows a pairs plot of the first four variates $\\{X_{1},X_{2},Z_{1},Z_{2}\\}$
for a simple run with sample size $N=1000$. Here the red and blue colors
indicate the group for every datum. The choice of the best variates to
construct the classification rule is an impossible task for a classical
correlation coefficient due to the nonlinear nature of the problem.

Figure 4: Example of four covariates for one run in the classification
scenario.

We have applied the same algorithm described before to this case but using a
generalized additive model:
$\mathbb{E}\left(Y\right)=g^{-1}\left(\alpha+\sum_{j=1}^{J}f_{j}\left(X^{(j)}\right)\right)$
where $g$ is the link function. Note that the “true” model is not additive but
a joint function of $X_{1}$ and $X_{2}$. The only change respect to the
regression example is the way of computing the distance among elements of the
categorical response. In this case, the distance matrix among elements of $Y$
only have two values: $0$ when both points belongs to the same group and
$\sqrt{2}$ (the distance between $(1,0)$ and $(0,1)$) on the contrary case. We
repeated this scenario $100$ times and applied the obtained model at each
repetition to $200$ new points to check its prediction ability. The relevant
variates $\\{X_{1},X_{2}\\}$ for the model were selected in all the runs as
the two first choices. The order is important because it is an indicator of
the strength of the information provided by the variates. On the other hand,
the proportion of times that the variates $\\{Z_{1},\ldots,Z_{48}\\}$ enter
the model were in the interval $\left[0,0.07\right]$. The mean of
misclassification error for new data was also excellent: $0.7$%.

##  4 Real data application

We have applied our proposal to the Market Demand Energy in the Iberian Energy
Market and its Price using the available information from different sources.
Indeed, it is possible to create new variables as transformations/functions of
the original ones. For instance, all the information included as a functional
variable was also included as scalar covariates (each discretization point was
a new one). We employed the FSAM/FAM model in all cases given its flexibility
and availability not only for multivariate variables, but also for functional
ones (see, for instance, Müller and Yao (2008) or Febrero-Bande and González-
Manteiga (2013)).

###  4.1 Energy Market Demand

We have collected information about $N=2,459$ days from August, 2008 to
December, 2014 of the Market Demand Energy at day $t$ and hour $H$
($\textrm{En}_{t}(H)$) in the IntraDay Electricity Iberian Market. Our aim is
to forecast this variable with the available information up to day $t-1$
obtained from the sources that measure energy or price. The meteorological
information is considered concurrent in time because the forecasting of these
variates can be done accurately. Here, the meteorological variables are taken
from Madrid-Barajas airport as a coarse information. Detailed contributions of
the meteorological information are expected that can occur when the energy
demand is analyzed at low levels (substations, counties, etc) but no when we
are dealing with national summaries. In any case, the station in Madrid seems
to be a good choice due to two reasons that may have influence in high energy
consumptions: Madrid has the highest population in Spain and it is surrounded
by an important industrial belt. Finally, we consider also the calendar
effects for every day although these variables may have redundant information
with meteorological ones. For instance, temperatures and solar radiations are
closely related with month as a calendar effect.

The following list summarizes all the variates involved:

  * •

Energy Market Variables (source:www.omie.es): Daily profiles of Energy (En)
and Price (Pr) at $t-1$ and $t-7$.

  * •

Total Load and Generated energy type (source:demanda.ree.es): Daily profiles
of Load (every ten minutes) (Lo), Nuclear ($\textbf{Nu})$, Fuel/Gas (Fu), Coal
(Ca), Combined Cycle (Cc), Solar (So), Wind Energy (WE), Hydroelectric (Hy),
Cogeneration (Co), Rest (Re).

  * •

Meteorological information at Madrid-Barajas airport (source:aemet.es):
Temperatures ($T_{Max}$, $T_{Min}$, $T_{Med}$, $T_{A}=T_{Max}-T_{Min}$), Wind
Speed ($WS$), Solar Radiation ($SR$), Precipitation, Pressure, $\ldots$.

  * •

Every discretization value of the functional variates ($t-1$ and $t-7$).

  * •

Categorical: Year($YY$), Month($MM$), Day-of-Week:
$DoW=\left\\{\mathbbm{1}_{\left\\{Mon\right\\}},\ldots,\mathbbm{1}_{\left\\{Sun\right\\}}\right\\}$.

Specifically, we are interested on estimating the model 1 where
$Y_{i}=\textrm{En}_{t}(H)$, the number $J$ and the covariates $X_{i}^{(j)}$
must be selected from the above list but evaluated before day $t$ (except for
meteorological and calendar variables), $f_{j}$ is the contribution of each
variable and $\varepsilon_{i}$ is the error. The sample size is the number of
days: $N=2,459$.

Table 5 summarizes the results of this application for four specific hours
along the day where the variates are listed in order of inclusion into the
model and the functional variables are marked in bold. In all those models,
the deviance is mostly explained by the effect of the three first variables
selected. As an example, the evolution in deviance of the model for
$\textrm{En}_{t}(18)$ was 59.9%, 71.9%, 80.0%, 80.7% and 81.7%.

Table 5: Summary of models for energy with its selected variables (in order of entering) Resp. | Covariates (in order) | Dev. expl. | $\sigma_{\varepsilon}$
---|---|---|---
$\textrm{En}_{t}(6)$ |  $\textrm{En}_{t-1}(6)$, $\textbf{En}_{t-7}$, $WS$, $\textbf{Pr}_{t-1}$, $\textbf{Pr}_{t-7}$, $\textrm{En}_{t-1}(24)$, $\textbf{Re}_{t-1}$ | 84.4% | 1207.6
$\textrm{En}_{t}(12)$ |  $\textbf{En}_{t-7}$, $\textrm{En}_{t-1}(16)$, $DoW$, $\textbf{Cc}_{t-1}$, $\textrm{En}_{t-1}(23)$ | 84.2% | 1540.3
$\textrm{En}_{t}(18)$ |  $\textbf{En}_{t-7}$, $\textrm{En}_{t-1}(17)$, $DoW$, $WS$, $\textbf{Pr}_{t-1}$, | 81.7% | 1662.6
$\textrm{En}_{t}(24)$ |  $\textrm{En}_{t-1}(24)$, $\textbf{Lo}_{t-7}$, $DoW$, $\textbf{Cc}_{t-1}$ | 86.0% | 1221.0

In these models, the information of the first contributor tends to be the same
as that of a week before rather than that of the previous day except for the
consumption in the early morning hours (06:00, 24:00). This suggests a strong
weekly pattern in the market demand profile complemented with the appearance
of the indicators for Day-of-Week as a contributor in the third place. The
different behaviour of the early morning hours models comes from the fact that
the energy demand at that hours is quite stable and has no main dependence on
the activity of a given day of the week. Therefore, it is enough to consider
the value of the day before. The inclusion into the models of variates like
$\textbf{Pr}_{t-1}$ or $\textbf{Cc}_{t-1}$ although shocking at first glance,
may be have a simple explanation. These two variates $\textbf{Pr}_{t-1}$ and
$\textbf{Cc}_{t-1}$ are closely related to each other because, following the
rules of Iberian Energy Market, the final price is fixed as the maximum of all
types of energy required to cover the demand being the combined cycle, one of
the expensive ones. The price is therefore higher when the proportion of
energy produced by combined cycle power stations is also high. These variates
are probably included in the models to improve the prediction on the days
characterized as having an energy demand that cannot be solved by using
renewable sources (or other cheaper ones) although, in any case, its
contribution is marginal.

###  4.2 Energy Price

The second application example corresponds to the negotiated price at day $t$
and hour $H$ ($\textrm{Pr}_{t}(H)$) in the Iberian Energy Market. The set of
possible covariates is the same as in the previous example and the results are
shown in Table 6.

Table 6: Summary of models for price with the selected variables (in order of entering) Resp. | Covariates (in order) | Dev. expl. | $\sigma_{\varepsilon}$
---|---|---|---
$\textrm{Pr}_{t}(6)$ |  $\textbf{Pr}_{t-1}$, $T_{A}$, $T_{Max}$, $T_{Med}$, $SR$, $\textrm{Pr}_{t-1}(24)$, $\textbf{Cc}_{t-1}$ | 79.7% | 6.93
$\textrm{Pr}_{t}(12)$ |  $\textrm{Pr}_{t-1}(23)$, $T_{Max}$, $\textbf{Lo}_{t-7}$, $\textrm{Pr}_{t-1}(9)$, $WS$, $\textbf{Cc}_{t-1}$, $\textbf{Pr}_{t-7}$, $\textbf{Cc}_{t-7}$ | 83.6% | 6.46
$\textrm{Pr}_{t}(18)$ |  $\textrm{Pr}_{t-1}(18)$, $DoW$, $WS$, $\textbf{Pr}_{t-7}$, $T_{Med}$, $\textrm{Pr}_{t-1}(16)$, $\textbf{Cc}_{t-1}$, $\textbf{Cc}_{t-7}$, $\textbf{Pr}_{t-1}$, $\textbf{WE}_{t-1}$ | 82.7% | 6.85
$\textrm{Pr}_{t}(24)$ |  $\textrm{Pr}_{t-1}(24)$, $\textbf{Pr}_{t-7}$, $\textbf{Cc}_{t-1}$, $\textbf{Cc}_{t-7}$, $YY$, $T_{A}$, $T_{Med}$ | 74.1% | 7.05

The main contributor in all price models is the price of the preceding day,
which indicates a strong persistence of this variable. Also, all models
include the combined cycle generation (in positions more or less advanced).
Surprisingly now, the variables related with demand do not become part of the
models. This is contrary to what is expected by the classical economic theory.
The meteorological variates supply this gap. For instance, the model for
$\textrm{En}_{t}(18)$ includes the mean of wind speed and the mean of
temperature as covariates. These two variates can jointly explain the high
price of energy on summer days with high temperatures but no wind, when the
energy system must provide high amounts of energy (for feeding the air
conditioning equipments) with low probability of using wind sources. The
inclusion of other meteorological variates into the models follows similar
rules. As commented before, it is not expected a high influence from
meteorological variables because the information comes from a particular site
(Barajas airport) and it would be hard to explain the price or the consumption
of energy for a whole country using only the information of one site even
though this is an important location. Surely, meteorological information is
useful for predicting demand in small regions but its contribution for price
can only be explained in terms of general rules. The calendar effects only
appear in two models. The Day-of-Week effect appears in a prominent position
in the model for $\textrm{En}_{t}(18)$ telling that difference in price
between Mondays and Sundays (the maximum difference) is about twelve euros
(the biggest differences are among weekend and non-weekend days). The Year
effect only appears in the model for $\textrm{En}_{t}(24)$ showing a clearly
increasing pattern along the years, which was particularly intense in the
period 2009-2011.

An important remark concerns the practical issues of the implementation of the
distance correlation when the sample size grows. The main difficulty is
unrelated to the ease of the implementation rather with the memory consumption
of the method. The best strategy for implementing this procedure making an
extensive use of the distance correlation with respect to the same covariates,
is to compute and store the distance matrix for every covariate in advance.
However, strategy is impossible when the sample size or the number of
covariates grows. The overall consumption of memory when storing the distance
matrices for $p$ covariates with $N$ elements each is
$(p+1)\left(N(N-1)/2-N\right)$. In this example, we have $p=300$ covariates
(many of which correspond to the discretization points of the functional ones)
with a sample size $N=2459$. This amount typically exceeds the available
resources of a desktop computer (even for the powerful ones). To overcome this
difficulty, the obvious response is to use a HPC facility where the
computation of distance correlation for each pair can be distributed along the
available nodes (storing the distance matrices in the shared storage device).
When it is impossible to access to a HPC facility, one may still compute the
distance correlation for an arbitrary sample size by simply dividing all the
computations in blocks using submatrices of dimension $L\times L$. As an
example, Table 7 shows the maximum consumed memory and the execution time (in
seconds) when the test of independence based on distance correlation for
different sample sizes and two different dimensions of the submatrices is
executed in a Intel Core i7-3770 with 32GB of RAM. The last row of the table
shows that with $N=25000$ the computer was unable to reserve enough memory.
The submatrix trick can manage that sample size maintaining a limit for the
memory consumption, but it does so with the cost of larger execution times.
The execution time using $L=500$ seems slightly better than with $L=1000$.
However we cannot conclude a general rule. For instance, for $N=10000$ the
execution times with $L=100,250$ and $2000$ were, respectively, $52.10,39.84$
and $64.93$.

Table 7: Maximum memory consumption and computation times (s.) for computing the test of independence using distance correlation with different computational strategies involving submatrices. | Memory consumption | Time(sec.)
---|---|---
N | Direct | $L=500$ | $L=1000$ | Direct | $L=500$ | $L=1000$
1000 | 7.63MB | 7.63MB | 7.63MB | 0.32 | 0.94 | 0.48
2500 | 47.68MB | 7.63MB | 30.52MB | 1.58 | 4.40 | 5.40
5000 | 190.7MB | 7.63MB | 30.52MB | 6.57 | 13.08 | 15.47
10000 | 762.9MB | 7.63MB | 30.52MB | 26.83 | 45.91 | 53.00
25000 | CRASH | 7.63MB | 30.52MB | – | 290.13 | 309.91

##  5 Conclusions

In this paper, we have introduced an algorithm that automatically selects the
variates for a regression or classification model. The procedure operates in a
forward way adding a variable to the model at each iteration. The key of the
whole procedure is the extensive use of the distance correlation that presents
two important advantages: the choice of the variate is made without
considering a model and it is possible to compute this quantity for variates
of different nature (functional, multivariate, circular, …) as it is only
computed from distances. The simplicity of the latter is also its main
drawback because the number of operations (and memory consumption) is of a
quadratic order respect to the sample size. But fortunately, it is possible to
compute the distance correlation when the sample size is huge while the
consumption of resources is maintained under certain limitations. Our proposal
is presented is a very general way although in the applications in this paper,
we have restricted ourselves to additive models that offer a balanced
compromise between predictive ability and simplicity. The obtained results are
quite promising in scenarios where no competitors are available because no
other procedure can deal with variates of different nature in a homogeneous
way. As a final comment, the procedure was applied to a real problem related
with the Iberian Energy Market (Price and Demand) where the number of possible
covariates is really big. The algorithm was able to find synthetic regression
models offering interesting insights about the relationship among the response
and the covariates. The final selected models mix functional, scalar and
categorical information.

## Acknowledgement

The authors acknowledge financial support from Ministerio de Economía y
Competitividad grant MTM2016-76969-P and European Regional Development Fund
(ERDF). The authors are also grateful to the Guest Editor, Associate Editor
and two anonymous referees for their helpful and valuable comments.

## References

  * Akaike [1973] H. Akaike.  Maximum likelihood identification of Gaussian autoregressive moving average models.  _Biometrika_ , 60(2):255–265, 1973.
  * Bühlmann and Yu [2003] P. Bühlmann and B. Yu.  Boosting with the ${L}_{2}$ loss: regression and classification.  _J Am Stat Assoc_ , 98(462):324–339, 2003.
  * Du et al. [2012] P. Du, G. Cheng, and H. Liang.  Semiparametric regression models with additive nonparametric components and high dimensional parametric components.  _Comput Stat Data Anal_ , 56(6):2006–2017, 2012\.
  * Efron et al. [2002] B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani.  Least Angle Regression. Technical report, Department of Statistics, 2002.
  * Fan and Li [2001] J. Fan and R. Li.  Variable selection via nonconcave penalized likelihood and its oracle properties.  _J Am Stat Assoc_ , 96(456):1348–1360, 2001.
  * Febrero-Bande and González-Manteiga [2013] M. Febrero-Bande and W. González-Manteiga.  Generalized additive models for functional data.  _TEST_ , 22(2):278–292, 2013.
  * Ferraty and Vieu [2009] F. Ferraty and P. Vieu.  Additive prediction and boosting for functional data.  _Comput Stat Data Anal_ , 53(4):1400–1413, 2009\.
  * Hastie et al. [2015] T. Hastie, R. Tibshirani, and M. Wainwright.  _Statistical learning with sparsity: the LASSO and generalizations_.  CRC Press, 2015.
  * Lin and Zhang [2006] Y. Lin and H. H. Zhang.  Component selection and smoothing in smoothing spline analysis of variance models.  _Ann Stat_ , 34(5):2272–2297, 2006.
  * Liu et al. [2011] X. Liu, L. Wang, and H. Liang.  Estimation and variable selection for semiparametric additive partial linear models.  _Stat Sinica_ , 21(3):1225–1248, 2011.
  * Lyons [2013] R. Lyons.  Distance covariance in metric spaces.  _Ann Probab_ , 41(5):3284–3305, 2013.
  * Mallows [1973] C. L. Mallows.  Some comments on $C_{p}$.  _Technometrics_ , 15(4):661–675, 1973.
  * Müller and Yao [2008] H.-G. Müller and F. Yao.  Functional additive models.  _J Am Stat Assoc_ , 103(484):1534–1544, 2008\.
  * Peng et al. [2005] H. Peng, F. Long, and C. Ding.  Feature selection based on mutual information criteria of max-dependency, max-relevance, and min-redundancy.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 27(8):1226–1238, 2005.
  * Radchenko and James [2011] P. Radchenko and G. M. James.  Improved variable selection with Forward-Lasso adaptive shrinkage.  _Ann Appl Stat_ , 5(1):427–448, 2011.
  * Rizzo and Székely [2014] M. L. Rizzo and G. J. Székely.  energy: E-statistics. R package version 1.6.2, 2014.
  * Schwarz [1978] G. Schwarz.  Estimating the dimension of a model.  _Ann Stat_ , 6(2):461–464, 1978.
  * Stone [1979] M. Stone.  Comments on model selection criteria of Akaike and Schwarz.  _J Roy Stat Soc B Met_ , 41(2):276–278, 1979\.
  * Székely and Rizzo [2013] G. J. Székely and M. L. Rizzo.  The distance correlation t-test of independence in high dimension.  _J Multivariate Anal_ , 117:193–213, 2013.
  * Székely and Rizzo [2017] G. J. Székely and M. L. Rizzo.  The Energy of Data.  _Annu Rev Stat Appl_ , 4(1):447–479, 2017.
  * Székely et al. [2007] G. J. Székely, M. L. Rizzo, and N. K. Bakirov.  Measuring and testing dependence by correlation of distances.  _Ann Stat_ , 35(6):2769–2794, 2007.
  * Tibshirani [1996] R. Tibshirani.  Regression shrinkage and selection via the LASSO.  _J Roy Stat Soc B Met_ , 58(1):267–288, 1996\.
  * Wang et al. [2011] L. Wang, X. Liu, H. Liang, and R. J. Carroll.  Estimation and variable selection for generalized additive partial linear models.  _Ann Stat_ , 39(4):1827–1851, 2011.
  * Xue [2009] L. Xue.  Consistent variable selection in additive models.  _Stat Sinica_ , 19(3):1281–1296, 2009.
  * Yenigün and Rizzo [2015] C. D. Yenigün and M. L. Rizzo.  Variable selection in regression using maximal correlation and distance correlation.  _J. Stat. Comput. Sim._ , 85(8):1692–1705, 2015\.
  * Zou [2006] H. Zou.  The adaptive LASSO and its oracle properties.  _J Am Stat Assoc_ , 101(476):1418–1429, 2006\.