# Enumeration of Polyominoes & Polycubes Composed of Magnetic Cubes††thanks:
This work used HPE DSI/IT computational resources and was supported by the NSF
under Grant Nos. [IIS-1553063, IIS-1849303, CNS-1932572] .††thanks: 1Authors
are with the University of Houston, Houston, TX 77204 USA {ylu36, debiediger,
<EMAIL_ADDRESS>2Authors are with Southern Methodist University,
Dallas, TX 75275 USA {abhattacharjee<EMAIL_ADDRESS>
Yitong Lu1, Anuruddha Bhattacharjee2, Daniel Biediger1, Min Jun Kim2, Aaron T.
Becker1

(May 2020)

###### Abstract

This paper examines a family of designs for magnetic cubes and counts how many
configurations are possible for each design as a function of the number of
modules. _Magnetic modular cubes_ are cubes with magnets arranged on their
faces. The magnets are positioned so that each face has either magnetic south
or north pole outward. Moreover, we require that the net magnetic moment of
the cube passes through the center of opposing faces. These magnetic
arrangements enable coupling when cube faces with opposite polarity are
brought in close proximity and enable moving the cubes by controlling the
orientation of a global magnetic field. This paper investigates the 2D and 3D
shapes that can be constructed by magnetic modular cubes, and describes all
possible magnet arrangements that obey these rules. We select ten magnetic
arrangements and assign a “color” to each of them for ease of visualization
and reference. We provide a method to enumerate the number of unique
polyominoes and polycubes that can be constructed from a given set of colored
cubes. We use this method to enumerate all arrangements for up to 20 modules
in 2D and 16 modules in 3D. We provide a motion planner for 2D assembly and
through simulations compare which arrangements require fewer movements to
generate and which arrangements are more common. Hardware demonstrations
explore the self-assembly and disassembly of these modules in 2D and 3D.

##  I Introduction

Small-scale modular robots are often designed to be manipulated in large
groups by an external system. Magnetic forces are a popular choice for
manipulating these robots because the robots can be directly controlled simply
by incorporating magnets or ferromagnetic material in the robot body. In the
last few decades, many research projects have been conducted based on magnetic
actuation and coupling systems to explore reconfigurable modular robotics at
different length scales [1, 2, 3, 4]. Researchers have developed modules with
excellent locomotion and manipulation abilities, and have integrated advanced
control algorithms. Several reconfiguration strategies of modules, such as
autonomous and distributed stochastic self-assembly, disassembly, self-
reconfiguration by rotation in a planar workspace, locomotion through
reconfiguration on a cubic lattice, and universal and distributed
reconfiguration planning for square and hexagonal-lattice-based robots have
been investigated [5, 6].

Figure 1: Designs of 2D and 3D modular cubes. Underneath each cube the
arrangement is labeled as {left [L], front [F]/back [B], top [T]/under [U],
right [R]}.

Much related work has focused on the “Tilt model” where components move in a
straight line until they hit an obstacle, and combine when mating particles
are brought close together. Solutions often design workspace obstacles that
enable arbitrary reconfiguration of components pushed by a global force [7,
8]. We designed algorithms that enabled efficient construction of desired
shapes [9], and workspace obstacles that enabled sorting and classifying
constructed polyominoes [10].

This previous work did not consider the stability or robustness of the
structures generated. Recent work in [11, 12] examined the forces between
individual components and evaluated if a structure of module robots would be
stable — if it would topple or collapse under its own weight, especially
during the process of reconfiguration.

Our future goal is to design motion plans for magnetic modular cubes to enable
them to form shapes suitable for fast motion, yet with designed fracture
points so they can be easily disassembled at a goal and reconfigured into
tools or structures. As first steps toward that goal, this paper presents the
following contributions: (1) Design and demonstrations of magnetic modular
cubes (MMCs) which are simple in design, scalable, and highly reconfigurable
through controlled assembly-disassembly. (2) Analysis of all possible
magnetization profiles that match the specifications of our modular cubes. (3)
A method to count the number of unique configurations for a given set of
cubes. (4) Enumerate the unique configurations up to $n=20$ in 2D and $n=16$
in 3D. (5) Provide a motion planner for 2D assembly. (6) Monte Carlo
simulations to measure the frequency that certain shapes can be constructed
and the required number of movement steps to generate the shapes. (7) Conduct
hardware experiments that explore the self-assembly and disassembly of these
modules in 2D and 3D. For applications of these magnetic cubes, see [13].

##  II Design of magnetic modular cubes

We design the MMCs so they can be controlled by a global magnetic field. The
cubes are manufactured so they align with an applied magnetic field. The
global magnetic field provides a reference frame with a magnetic south to the
left and north to the right. We labeled each face of the cube as [L, F, B, T,
U, R] for _left_ , _front_ , _back_ , _top_ , _under_ , and _right_ faces. On
each face at most one axially magnetized permanent magnet is embedded. Each
cube must have a net magnetic orientation. The magnets are positioned with
poles facing outwards such that each face is either north or south and such
that the net magnetic moment of the cube from magnetic south to north passes
through the center of opposing faces. This net magnetic moment is required
because it enables “pivot walking” and “rolling” the cubes. Pivot walking is a
movement gait that alternately lifts the north or south pole of the cube so
the cube balances on the opposing edge, and then applies torque to spin the
cube on the balancing edge, as in [14, 15]. The rolling motion [16] can be
achieved by applying continuous rotational magnetic torque along an axis
perpendicular to the axis of net magnetization of a cube to roll it towards
the desired direction.

Because the magnetic field passes through the center of the cube face, the
magnetic field must be rotationally symmetric across the left-to-right axis.
Thus the magnetization of F==B and T==U. We can choose the L to be $(S$ or
$0)$, the F and B to be $(S,$ $0,$ or $N)$, the T and U to be $(S,$ $0,$ or
$N)$, and the R $(0$ or $N)$, for a total of $2\cdot 3\cdot 3\cdot 2=36$
arrangements. The arrangement is labeled as {L, F/B, T/U, R}. We then remove
the nine arrangements of the form $\\{0,*,*,0\\}$ because they have no net
magnetic moment. Therefore, there are 27 sets of suitable magnet arrangements.
2D designs of MMCs only allow connection along the $x$-axis and $y$-axis. No
magnets are embedded on the top [T] or under [U] faces. 3D designs of MMCs
enable connection along faces (see Fig. 1).

Due to the MMCs’ cubic design, magnetically connected structures of MMCs are
polyominoes in 2D and polycubes in 3D. Polyominoes and polycubes are a classic
topic in combinatorics. This paper draws on a rich literature and results on
enumerating polyominoes and polycubes [17, 18, 19, 20, 21] and on coloring
them [22]. Work on characterizing upper and lower bounds on the number of
polyominoes and polycubes is an active research area [23, 21], and will be
relevant to modular robotic construction.

Similarly, there has been a great deal of recent research on how to design
workspaces that exploit global control of magnetic particles to design
arbitrary shapes [8, 24, 25], and on the complexity of motion planning under
such constraints. These papers provide important solutions on how to design
workspaces that enable fast rearrangement of cubes, but they do not provide
insight on how to build those cubes. Moreover, these works tend to assume all
cubes will stick together if brought close enough in proximity. In contrast,
this paper focuses on the configurations made possible by the magnetic
profiles of the cubes.

##  III Generate all possible polyominoes in 2D

_Polyominoes_ are face-connected sets of unit cubes that lie on the square-
grid graph. Each cube is represented by an integer tuple ($x$, $y$). Cubes
($x_{1}$, $y_{1}$) and ($x_{2}$, $y_{2}$) are adjacent if
$|x_{1}-x_{2}|+|y_{1}-y_{2}|=1$. Two _free polyominoes_ are considered
distinct if they have different shapes, but they cannot be rotations of each
other. Two _fixed polyominoes_ are considered distinct if they have different
shapes or orientations. Because the global magnetic field provides an
orientation, we consider fixed polyominoes for our MMCs design. The number of
fixed polyominoes of size $n$ is denoted by $A_{n,2}\textsubscript{D}$. The
number of polyominoes grows exponentially [23]

$\lim_{n\to\infty}(A_{n,2}\textsubscript{D})^{\frac{1}{n}}=\lambda.$ (1)

Figure 2: Enumeration of colored cubes in 2D. The composition of cube types
determines the number of polyominoes that can be made. The white arrow
indicates the direction of south to north magnetic pole. Table I: Number of
fixed polyominoes $A_{n,2}\textsubscript{D}$, where ${n}$ is the number of
cubes. With MMC-m2D and MMC-c2D, ${C_{n,i}}$ is the number of colored
polyominoes with ${i}$ magenta cubes and $n$ total cubes. Row maximum in
orange.

We implemented a method to enumerate the fixed polyominoes and unique colored
polyominoes that can be constructed by a set of MMCs [26]. Based on the
hardware experiment environment (see Section V), we developed a 2D motion
planner using a bounded 11$\times$11 workspace that generates the set of valid
colored polyominoes reachable by global moves from an initial configuration.
The cubes are represented by axis-aligned, unit-length tiles in our motion
planner. Arrows on the tiles indicate the direction of the net magnetization
(from magnetic south pole to north pole). We select four magnet arrangements
that use at least two magnets and, for ease of visualization, assign them the
colors cyan, magenta, brown, and black (see Fig. 2 top right).

We enumerate the fixed polyominoes from a given number of cubes, then
calculate the set of valid colored polyominoes that can be generated from a
given set of MMCs. We determine the shortest movement sequences to generate
all reachable polyominoes from a given initial configuration. For each
polyomino, we analyze the fraction of starting configurations for which it is
reachable and the number of moves required to construct it.

###  III-A Enumeration of fixed polyominoes

To enumerate fixed polyominoes, we implemented the method used by Redelmeier
in [27]. We use a variant of this algorithm to enumerate colored polyominoes
constructed from MMCs. Because the magnets in these cubes define a coordinate
system, we enumerate all polyominoes with magnetic moments pointing in the
$+y$ direction. The algorithm works by recursively generating all cubes up to
a given size $n$. We start with the leftmost cube of the bottom row, placed at
the origin $(0,0)$. For a given configuration with $k$ cubes, the algorithm
generates a new configuration of size $k+1$ by adding a cube connected to the
polyomino. We maintain a list of the adjacent cubes, _list_adj_ , and a list
of which adjacent cubes have been selected, _poly_ , for each recursive call.
In the initial call, the list of adjacent cubes, _list_adj_ , consists of the
origin $(0,0)$, the cube above the origin $(0,1)$, and the cube to the right
of the origin $(1,0)$. We set _poly_[0] = 0. At each recursive call, if
$k\equiv n$, the procedure returns. Otherwise, a recursive call is made for
each cube in _list_adj_ greater than the last cube in _poly_. For each of
these calls, the cube is added to _poly_ , and any cubes adjacent to the new
cube are appended to _list_adj_ if they are not already in _list_adj_. To
avoid generating the same cube more than once, cubes are only added if they
are above the origin, or in the same row to the right of the origin:

$\left\\{(x,y)\leavevmode\nobreak\ |\leavevmode\nobreak\ (y>0)\textrm{ or
}(y\equiv 0\textrm{ and }x\geq 0)\right\\}$ (2)

###  III-B Enumeration of valid colored polyominoes from a set

To count the number of valid colored polyominoes from a set of colored cubes,
we modify the previous algorithm. We include the coloring information in
_poly_. Before each recursive call, we check if the new cube can be colored in
the given color. A check is successful if there are sufficient unused cubes of
the color, and if placing the cube does not violate any of the magnetic-
assembly rules. If the check is successful, a recursive call is performed.
Fig. 2 top left shows that for $n=3$ cubes, 6 fixed polyominoes can be
enumerated. Using a supply with one cyan and two magenta cubes produces more
polyominoes than one brown and two magenta cubes. However, using three colored
cubes enables constructing even more colored polyominoes.

Figure 3: Left: The self-assembly result for an arbitrary initial
configuration of four cubes with different designs. Right: Details of the
movements. Figure 4: Top Left: The number of polyominoes that can be made with
1 cyan and 3 magenta cubes from random starts. The red line is the maximum
number of possible configurations (26). This test was repeated for 1000
trials, and analyzed in the remaining plots. Bottom Left: Cumulative number of
unique reachable polyominoes. Red + signs show outliers. Middle Top: Total
number of steps to construct each polyomino (when the polyomino is reachable).
Middle Bottom: Frequency results showing how often each polyomino is
reachable. Right: Scatter plot of percentage reachable vs. number of steps.

Table I shows the number of fixed polyominoes that can be enumerated with up
to $n=20$ cubes. It also lists the number of valid polyominoes with magenta
and cyan cubes that can be constructed by coloring them with up to 20 cyan
cubes. With zero to two magenta cubes, the number of valid polyominoes is less
than the number of fixed polyominoes ($C_{n,i}$ for $i$ = 0, 1 or 2). This is
because two magenta or two cyan cubes can only be connected in series (with
arrows aligned tip-to-tail). However, a more extensive set of polyominoes with
different shapes and orientations can be constructed when the fraction of
magenta cubes is increased. The largest variety is possible when $i$ =
$\lfloor\frac{n}{2}\rfloor$. The top right corner of Table I shows examples of
valid colored polyominoes.

###  III-C Self-assembly algorithm

We provide a low-fidelity motion model that can compute reachable polyomino
configurations and the shortest movement sequences for a set of MMCs from an
initial configuration. We assume that all modules move at the same speed in
the same direction, unless they encounter a fixed obstacle. Without loss of
generality, we limit the movement directions to north, east, south, and west.

Fig. 3 shows an example of all the valid colored polyominoes can be generated
with 4 different cube designs. The initial configuration and the details of
the movements are shown on the right. Underneath each polyomino is listed the
shortest step sequence to construct it. At each step, we could move one unit
length in the four cardinal directions {n, e, s, w} or rotate {r}
counterclockwise (CCW) or clockwise (CW) a quarter-turn. We prune this search
by only allowing rotations if no translation moves have been made and by only
performing CW rotations. For example, in this shorthand, “2r, 4w, 6n” means
rotate 180∘ CW followed by four unit moves west, and then six unit moves
north. We used a breadth-first search (BFS) algorithm to discover these
reachable configurations. The root of the search is the initial configuration,
and we maintain a list of all relative configurations that have been reached.
A _relative configuration_ is represented by rotating the coordinate frame so
that the magnetic north points up and then translating all cubes. To further
prune the tree, we only switch direction if the relative configuration of the
cubes changes, and we terminate search branches if the present configuration
has occurred before. When the relative configuration changes, we call this an
_intermediate configuration_. Intermediate configurations occur when at least
one cube (but not all cubes) strikes an obstacle. When this occurs, the new
relative configuration is compared to the list and is only appended if it is
unique. After the polyomino is assembled, we compare the paths taken to save
the shortest step sequence to construct the structure.

Table II: Number of fixed polycubes $A_{n,3}\textsubscript{D}$, where ${n}$ is
the number of cubes. Row maximum in orange.
With MMC-r3D and MMC-b3D, ${C_{n,i}}$ is the number of colored polycubes with
${i}$ red cubes and $n$ total cubes.

The number of reachable polyominoes depends on the initial configuration. For
a given number of magenta and cyan cubes, a finite number of colored
polyominoes can be generated (see Table I), but some may be unreachable from
an arbitrary initial configuration. Different initial configurations usually
generate different numbers of reachable polyominoes. Fig. 4 top left shows
examples of 20 random initial configurations. For each configuration, the
number of unique reachable polyominoes was calculated. Three black boxes show
the valid reachable polyominoes that can be generated from three initial
configurations. However, from none of these 20 configurations could we
construct every one of the 26 valid colored polyominoes (see Table I). If the
desired polyomino is unreachable from a given initial configuration, one could
repeatedly generate new random initial configurations until the polyomino is
reachable. Rapidly rotating the magnetic field is one way to separate and
scramble the cubes into new configurations.

To test the feasibility of this, we ran 1000 trials using one cyan and three
magenta cubes. Fig. 4 bottom left shows the net unique polyominoes that can be
constructed after a given number of configurations. The median number of
reachable polyominoes from the first start is 14. With two restarts, it is
possible to cover on average 20.28 of the 26 possible polyominoes
configurations. In 999 out of 1000 trials, all 26 polyominoes were reachable
after 20 accumulated configurations. The middle top figure shows the number of
steps needed to construct each polyomino (if the polyomino is reachable) over
1000 trials. For these trials, we only counted translation steps and did not
count rotation steps because translations require multiple changes to the
magnetic field. The middle bottom figure shows the percentage reachable
results for each polyomino. Icons at the bottom of these figures show each
valid polyomino. I-shapes require on average the fewest steps to form.
L-shapes with a cyan cube in the corner are reachable most often (the highest
frequency is 786 out of 1000 trials). Z-shapes are harder to make because they
require assembly in a specific temporal order (moving all cubes to give them
the same $y$ coordinate, so to make a Z-shape, the motion must be stopped
before all cubes reach the boundary). Z-shapes also require on average the
most steps to form (see Fig. 4 right).

##  IV Generate all possible polycubes in 3D

A _polycube_ is a three-dimensional polyomino constructed by cubes attaching
face to face. Each cube is represented by an integer tuple $(x,y,z)$. Cubes 1
and 2 are adjacent if $|x_{1}-x_{2}|+|y_{1}-y_{2}|+|z_{1}-z_{2}|=1$. Polycubes
can be enumerated in two ways, depending on whether different orientations are
counted as one polycube or two. Two _free polycubes_ are considered distinct
if they have different shapes, but are not 3D rotations of each other. Two
_fixed polycubes_ are considered equivalent if one can be transformed into the
other by a translation. The number of fixed three-dimensional polycubes of
size $n$ is denoted by $A_{n,3}\textsubscript{D}$.

Figure 5: Enumeration of colored cubes in 3D. The composition of cube types
determines the number of polycubes that can be made. The white arrow indicates
the direction of south to north magnetic pole.  Figure 6: Left: Maximum number
of valid colored 3D polycubes with different sets of cube colors, where ${n}$
is the total number of cubes. Right: Enumeration results of valid colored
polyominoes and polycubes.

###  IV-A Enumeration of fixed polycubes

As with the enumeration of polyominoes in 2D, we start with the leftmost cube
of the bottom row and place it at the origin $(0,0,0)$. In the initial call,
the list of adjacent cubes, _list_adj_ , consists of the origin $(0,0,0)$, the
cube to the right of the origin $(1,0,0)$, the cube to the back of the origin
$(0,1,0)$, and the cube on the top of the origin $(0,0,1)$, we set _poly_[0] =
0. To avoid double-counting polycubes, cubes are only added at:

$\displaystyle\big{\\{}(x,y,z)\leavevmode\nobreak\ |\leavevmode\nobreak\ $
$\displaystyle(z>0)\textrm{ or }(z\equiv 0\textrm{ and }y>0)$
$\displaystyle\textrm{ or }(z\equiv 0\textrm{ and }y\equiv 0\textrm{ and
}x>0)\big{\\}}.$ (3)

###  IV-B Enumeration of valid colored polycubes from a set

We select six magnet arrangements that use at least four magnets and assign
them the colors in Fig. 5 top right. As with the enumeration of valid colored
polyominoes in 2D, to count the number of valid colored polycubes generated
from a set of MMCs, we color the cubes while enumerating fixed polycubes
during the recursive step. Fig. 5 top left shows that for $n=3$ cubes, 15
fixed polycubes can be enumerated. Middle left and bottom left figures show
examples of coloring with two colors. Using a supply with red and blue cubes
can produce more polycubes than red and green. With a supply of three colors,
a wider variety of valid colored polycubes can be constructed.

Table II shows the number of fixed polycubes that can be enumerated with up to
$n=16$ cubes and the number of valid polycubes that can be constructed by
coloring them with red and blue cubes. The top right corner of the table shows
examples of valid colored polycubes. Fig. 6 left shows the maximum number of
valid polycubes that can be constructed with five different sets of cubes in
3D, where ${n}$ is the total number of cubes. For each cube set, the left
column shows the maximum number of valid colored polycubes that can be
generated for a given total number of cubes. The right column in each set
(gray font) shows configurations in a shorthand, e.g. “1r4b” means 1 red and 4
blue cubes.

###  IV-C Rejecting non-magnetically connected components

Our enumeration of valid polyominoes (or polycubes) allows placing a cube face
without magnets next to a cube of any color, but there is no magnetic
connection. Our goal for the polyominoes (or polycubes) is to create
aggregates that can be manipulated in unison by a magnetic field. For the
purpose of this study, we require that every cube be magnetically connected to
form a single component. This requires connected-component analysis.
Connected-component labelling is normally used for distinguishing different
objects in a binary image [28]. A _magnetically-connected component_ is a set
of cubes, $C$, such that for any two cubes in $C$, there is a valid connected
path between them and this path is contained in $C$. A _valid connection_
between two MMCs means the two cubes share a face and on this face one cube
has a south pole outward and the other has a north pole outward. Before adding
a polycube of size $n$ to our enumeration, we first check that the polycube is
a single magnetically-connected component. We start with the $(0,0,0)$ cube,
and perform a breadth-first search on all valid-connected neighboring cubes.
If the total number of connected cubes is $n$, the polycube is counted,
otherwise it is rejected. All the polyominoes in Fig. 7 are the same shape and
are magnetically connected, but the lack of magnetic bonds introduces planar
cracks in the assembly, which could be exploited for controlled disassembly.
Future work could catalog these shapes during the enumeration process.

Figure 7: Examples of magnetically connected polyominoes with the same shape
but different bonds. Red lines show potential failure surfaces with no
magnetic connection.

###  IV-D Time complexity

As with the enumeration scheme by Redelmeier, our enumeration counts each
possible polyomino and polycube only one time, the enumeration complexity is
exponential, $O(c_{\text{2D}}^{n})$ and $O(c_{\text{3D}}^{n})$. Our
implementation is slightly less efficient if we allow cubes that do not have
magnets on each side, since we generate intermediate configurations that might
not form a magnetically-connected component. After constructions, the final
step is to check if the polycube is magnetically-connected, and reject those
that are not. The connected check requires $O(n)$ time. 3D construction is
more complex than 2D, so the constant $c_{\text{3D}}>c_{\text{2D}}$. Fig. 6
right shows the results of enumeration of valid colored polyominoes and
polycubes.

##  V Experimental Setup and Results

A large-scale nested triaxial Helmholtz coil system was used to conduct
experiments. This system was designed parametrically to create uniform
magnetic fields in $x$, $y$, and $z$. Each coil pair is connected to its own
programmable power supply which is controlled by a National Instruments data
acquisition (DAQ) board. The power supplies generate sinusoidal outputs to the
coils, which in turn creates a uniform rotating magnetic field with user-
specified magnitude and frequency. These signals are controlled using a
customized C++ program. The 3D magnetic field vector produced from the system
is represented by

$\textbf{{B}}=\begin{bmatrix}B_{x}\\\ B_{y}\\\ B_{z}\\\
\end{bmatrix}=\begin{bmatrix}A\cos(\alpha)\cos(\theta)\\\
A\cos(\alpha)\sin(\theta)\\\ A\sin(\alpha)\\\ \end{bmatrix},$ (4)

where B is the applied magnetic flux density, $B_{x}$, $B_{y}$, and $B_{z}$
are the three-dimensional components along $x$, $y$ and $z$-axes respectively,
$A$ is the amplitude, $\alpha$ is the pitch angle, and $\theta$ is the yaw
angle. Individual cubes that have a volume of 1cm3 were designed using a CAD
software (Onshape) and printed using polylactic acid (PLA) on a 3D printer
(Ultimaker 2 Extended+). A Matlab program collected magnetic field settings
throughout the experiment and recorded assembly behavior of modular subunits
with a digital camera.

Figure 8: 2D self-assembly and disassembly behavior of modular subunits with
eight cubes.

A magnetic flux of 10mT was enough to actuate the modular cubes and achieve
successful assembly. Pivot walking motion was used to actuate individual
modular cubes with precise motion, following predetermined paths to create a
target assembly. When the arrows on top faces of the modular cubes are aligned
side-by-side, this is defined as a parallel assembly. When the arrows are
aligned tip-to-tail, this is defined as a serial assembly. The self-assembly
is caused by the magnetic attraction force generated from embedded permanent
magnets in the modular cubes when they are moved close to each other. Fig. 8
shows 2D self-assembly and disassembly behavior of modular subunits for eight
cubes. Starting with two groups of modular subunits, with an equal number of
modular cubes in each group, serial and parallel assembly were explored to
form square- and  L -shapes (at 44s). However, the self-assembly configuration
depends on the motion path taken (i.e., parallel assembly or serial assembly)
and the choice of subunits. From 44s to 55s, the assembled structures could be
moved downwards by pivot walking, like an individual modular cube. From 55s to
56s, the sub-assemblies were rotated $90^{\circ}$ CW by changing the direction
of the static magnetic field. The sub-assemblies could then be moved right or
left to bring them into close proximity. The sub-assemblies were then joined
to create a large 9-shaped structure (at 62s). Then, it was disassembled into
two different shapes (at 82s), with a rapid $90^{\circ}$ CCW rotational
magnetic torque. The disassembly occurred at the weakest magnetic joint in the
structure. The modular structures can also be moved in the workspace using the
rolling motion [16] which allows flipping the structures (at 91s).

Figure 9: 3D self-assembly and disassembly behavior of modular subunits with
eight cubes.

Fig. 9 shows a representative 3D construction sequence using global magnetic
control. A 3D self-assembly was performed following a 2D self-assembly to
create a pair of sub-assemblies. Starting with a similar configuration as Fig.
8, but choosing a different motion path in the 2D self-assembly process formed
an  L -shape and a square-shape with different cube orderings at 36s. The sub-
assemblies were rotated CW (at 30s) and CCW (at 47s) for movement in different
directions. Once 2D structures were formed with pivot walking motion followed
by 2D self-assembly, the sub-units can be further joined in 3D by using
rolling motions. The rolling motion was achieved by applying continuous
rotational magnetic torque on the modular structures about a fixed axis until
the target location was reached (from 57s to 66s). While the sub-assemblies
were rolled to bring them in close proximity, a blue modular cube in the  L
-shaped structure was disassembled (at 65s). The blue cube reattached to the
square-shaped structure (at 66s) which resulted in an unplanned
reconfiguration. This spontaneous disassembly and reconfiguration could have
benefits, but could also be avoided with improved control. Once the
reconfigured sub-units move into close proximity, they self-assembled to
create a two layered 3D structure (at 79s). This 3D structure was also
actuated in the workspace with the pivot walking motion (from 79s to 90s).

##  VI Conclusion

This paper describes the 2D and 3D configurations of modular cubes with
magnetic faces controllable by a global magnetic field. Experimental results
demonstrate the self-assembly and disassembly behavior of magnetic modular
cubes in 2D and 3D. Computational modeling enumerates the set of possible
polyominoes and polycubes from a set of modular cubes. The 2D motion planner
computes the shortest movement sequences to generate all reachable polyomino
configurations.

Future work will focus on 3D planning, in particular making a high-fidelity
simulator for motion planning. It would be interesting to calculate the number
of valid colorings of a given polycube as a function of the supply of colored
cubes. During the enumeration process one could also catalog useful
configurations, ranking them on their mechanical stability, fracture modes, or
magnetic profile in addition to their shape.

## References

  * [1] J. W. Romanishin, K. Gilpin, and D. Rus, “M-blocks: Momentum-driven, magnetic modular robots,” in _IROS_. IEEE, 2013, pp. 4288–4295.
  * [2] S. Tasoglu, C. Yu, H. Gungordu, S. Guven, T. Vural, and U. Demirci, “Guided and magnetic self-assembly of tunable magnetoceptive gels,” _Nature Communications_ , vol. 5, no. 1, pp. 1–11, 2014.
  * [3] H. Xie, M. Sun, X. Fan, Z. Lin, W. Chen, L. Wang, L. Dong, and Q. He, “Reconfigurable magnetic microrobot swarm: Multimode transformation, locomotion, and manipulation,” _Science Robotics_ , vol. 4, no. 28, 2019\.
  * [4] W. Saab, P. Racioppo, and P. Ben-Tzvi, “A review of coupling mechanism designs for modular reconfigurable robots.” _Robotica_ , vol. 37, no. 2, pp. 378–403, 2019.
  * [5] R. Pfeifer, M. Lungarella, and F. Iida, “Self-organization, embodiment, and biologically inspired robotics,” _Science_ , vol. 318, no. 5853, pp. 1088–1093, 2007.
  * [6] K. Tomita, S. Murata, H. Kurokawa, E. Yoshida, and S. Kokaji, “Self-assembly and self-repair method for a distributed mechanical system,” _IEEE Transactions on Robotics and Automation_ , vol. 15, no. 6, pp. 1035–1045, 1999\.
  * [7] J. Balanza-Martinez, A. Luchsinger, D. Caballero, R. Reyes, A. A. Cantu, R. Schweller, L. A. Garcia, and T. Wylie, “Full tilt: Universal constructors for general shapes with uniform external forces,” in _Proceedings of the Thirtieth Annual ACM-SIAM Symposium on Discrete Algorithms_. SIAM, 2019, pp. 2689–2708.
  * [8] J. Balanza-Martinez, T. Gomez, D. Caballero, A. Luchsinger, A. A. Cantu, R. Reyes, M. Flores, R. Schweller, and T. Wylie, “Hierarchical shape construction and complexity for slidable polyominoes under uniform external forces,” in _Symposium on Discrete Algorithms_. SIAM, 2020, pp. 2625–2641.
  * [9] A. T. Becker, S. P. Fekete, P. Keldenich, D. Krupke, C. Rieck, C. Scheffer, and A. Schmidt, “Tilt assembly: Algorithms for micro-factories that build objects with uniform external forces,” _Algorithmica_ , vol. 82, no. 2, pp. 165–187, 2020.
  * [10] P. Keldenich, S. Manzoor, L. Huang, D. Krupke, A. Schmidt, S. P. Fekete, and A. T. Becker, “On designing 2d discrete workspaces to sort or classify polynminoes,” in _IROS_. IEEE, 2018, pp. 1–9.
  * [11] P. Hołobut and J. Lengiewicz, “Distributed computation of forces in modular-robotic ensembles as part of reconfiguration planning,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2017, pp. 2103–2109.
  * [12] B. Piranda, P. Chodkiewicz, P. Hołobut, S. Bordas, J. Bourgeois, and J. Lengiewicz, “Distributed prediction of unsafe reconfiguration scenarios of modular-robotic programmable matter,” _arXiv preprint arXiv:2006.11071_ , 2020.
  * [13] A. Bhattacharjee, Y. Lu, A. T. Becker, and M. J. Kim, “Magnetically-controlled modular cubes with reconfigurable self-assembly and disassembly,” _IEEE Transactions on Robotics_ , 2021.
  * [14] E. Al Khatib, A. Bhattacharjee, P. Razzaghi, L. W. Rogowski, M. J. Kim, and Y. Hurmuzlu, “Magnetically actuated simple millirobots for complex navigation and modular assembly,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, pp. 2958–2965, 2020.
  * [15] L. Rogowski, A. Bhattacharjee, X. Zhang, G. Kararsiz, H. Fu, and M. J. Kim, “Manipulation planning of magnetically actuated programmable cuboids for collaborative assembly,” in _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems, Las Vegas, NV, USA_ , 2020.
  * [16] C. Bi, M. Guix, B. V. Johnson, W. Jing, and D. J. Cappelleri, “Design of microscale magnetic tumbling robots for locomotion in multiple environments and complex terrains,” _Micromachines_ , vol. 9, no. 2, p. 68, 2018.
  * [17] G. Aleksandrowicz and G. Barequet, “Counting d-dimensional polycubes and nonrectangular planar polyominoes,” in _International Computing and Combinatorics Conference_. Springer, 2006, pp. 418–427.
  * [18] S. Luther and S. Mertens, “Counting lattice animals in high dimensions,” _Journal of Statistical Mechanics: Theory and Experiment_ , vol. 2011, no. 09, p. P09026, 2011.
  * [19] G. Barequet and M. Shalah, “Polyominoes on twisted cylinders,” in _Proceedings of the Twenty-Ninth Annual Symposium on Computational Geometry_ , 2013, pp. 339–340.
  * [20] A. R. Conway, “The design of efficient algorithms for enumeration,” _arXiv preprint arXiv:1610.09806_ , 2016.
  * [21] G. Barequet and M. Shalah, “Improved upper bounds on the growth constants of polyominoes and polycubes,” _arXiv preprint arXiv:1906.11447_ , 2019.
  * [22] E. D. Demaine, M. L. Demaine, D. Eppstein, and J. O’Rourke, “Some polycubes have no edge-unzipping,” _arXiv preprint arXiv:1907.08433_ , 2019.
  * [23] G. Barequet, G. Rote, and M. Shalah, “$\lambda>$ 4: An improved lower bound on the growth constant of polyominoes,” _Commun. ACM_ , vol. 59, no. 7, pp. 88–95, 2016.
  * [24] D. Caballero, A. A. Cantu, T. Gomez, A. Luchsinger, R. Schweller, and T. Wylie, “Building patterned shapes in robot swarms with uniform control signals,” in _Proceedings of CCCG_ , 2020.
  * [25] U. K. Cheang, F. Meshkati, H. Kim, K. Lee, H. C. Fu, and M. J. Kim, “Versatile microrobotics using simple modular subunits,” _Scientific Reports_ , vol. 6, p. 30472, 2016.
  * [26] Y. Lu and A. T. Becker, March 2021. [Online]. Available: https://github.com/RoboticSwarmControl/EnumeratingPolycubes/
  * [27] D. H. Redelmeier, “Counting polyominoes: yet another attack,” _Discrete Math._ , vol. 36, no. 2, pp. 191–203, 1981.
  * [28] L. He, X. Ren, Q. Gao, X. Zhao, B. Yao, and Y. Chao, “The connected-component labeling problem: A review of state-of-the-art algorithms,” _Pattern Recognition_ , vol. 70, pp. 25–43, 2017.


# Financial Time Series Analysis and Forecasting
with HHT Feature Generation and Machine Learning

Tim Leung Applied Mathematics Department, University of Washington, Seattle WA
98195. Email<EMAIL_ADDRESS>Corresponding author. Theodore Zhao Applied
Mathematics Department, University of Washington, Seattle WA 98195. Email:
<EMAIL_ADDRESS>
###### Abstract

We present the method of complementary ensemble empirical mode decomposition
(CEEMD) and Hilbert-Huang transform (HHT) for analyzing nonstationary
financial time series. This noise-assisted approach decomposes any time series
into a number of intrinsic mode functions, along with the corresponding
instantaneous amplitudes and instantaneous frequencies. Different combinations
of modes allow us to reconstruct the time series using components of different
timescales. We then apply Hilbert spectral analysis to define and compute the
associated instantaneous energy-frequency spectrum to illustrate the
properties of various timescales embedded in the original time series. Using
HHT, we generate a collection of new features and integrate them into machine
learning models, such as regression tree ensemble, support vector machine
(SVM), and long short-term memory (LSTM) neural network. Using empirical
financial data, we compare several HHT-enhanced machine learning models in
terms of forecasting performance.

##  1 Introduction

Market observations and empirical studies have shown that financial time
series such as market indices and asset prices are often driven by multiscale
factors, ranging from long-term economic cycles to rapid fluctuations in the
short term. This suggests that financial time series are potentially embedded
with different timescales. Traditionally, Fourier analysis has been a powerful
tool in finding a spectral representation of time series. However, the Fourier
bases, by assumption, are stationary and linear. On the other hand,
nonstationary and behaviors and nonlinear dynamics are often observed in
financial time series. These characteristics can hardly be captured by linear
models and call for an adaptive and nonlinear approach for analysis. For
decades, methods based on short time Fourier transform and wavelets have been
developed and applied to nonstationary time series, but there are still
challenges in capturing nonlinear dynamics, and the often prescribed
assumptions make the methods not fully adaptive. This gives rise to the need
for an adaptive and nonlinear approach for analysis.

One alternative approach in adaptive time series analysis is the Hilbert-Huang
transform (HHT) [1]. The HHT method can decompose any time series into
oscillating components with nonstationary amplitudes and frequencies using the
empirical mode decomposition (EMD). This fully adaptive method provides a
multiscale decomposition for the original time series, which gives richer
information about the time series. The instantaneous frequency and
instantaneous amplitude of each component are later extracted using the
Hilbert transform. The decomposition onto different timescales also and allows
for reconstruction up to different resolutions, providing a smoothing and
filtering tool that is ideal for noisy financial time series. The method of
HHT and its variations have been applied in numerous fields, from engineering
to geophysics [2].

There have been several studies on the variations and alternatives to EMD,
including optimization based methods [3, 4, 5, 6] and noise-assisted
approaches [7, 8]. For financial time series with high level of intrinsic
noise, we apply the complementary ensemble empirical mode decomposition
(CEEMD) [8]. Like EMD, CEEMD decomposes any time series – stationary or not –
into a number of intrinsic mode functions representing the local
characteristics of the time series at different timescales, but the timescale
separation is improved by resolving mode mixing in EMD [9]. The noise-assisted
approach is also more robust to intrinsic noise in the data. In previous
studies, CEEMD have been found to be useful for forecasting [10, 11] and
signal processing [12]. More recent application of EMD and the noise-assisted
variations includes [13, 14, 15, 16].

Applications of EMD and HHT to finance date back to the work by Huang and co-
authors on modeling mortgage rate data [17]. In previous studies, EMD has been
used for financial time series forecasting [18, 19] and for examining the
correlation between financial time series [20]. The studies that adopted the
full HHT procedure for machine learning appeared in the electrical engineering
field [21, 22, 23]. However, in most of the previous studies using EMD or HHT
as forecasting features, the decomposition was pre-implemented on the whole
time span, including training and testing data, which is not practical when
doing real-time forecasting. Moreover, both EMD and HHT lend information from
the past and the future to compute the current values, which can cause
information leakage in the forecasting process. On the other hand,
extrapolating prediction as in [18] is subject to the end effect of EMD, which
makes the decomposition error quite high at the end of the time span. It is
unknown from previous studies of how to overcome the challenge of end effect .

In this paper, using CEEMD and Hilbert spectral analysis, we derive the
instantaneous energy-frequency spectrum associated with financial time series
to examine the properties of various timescales embedded in the original time
series. Important interpretable multiscale outputs from the HHT procedure are
recognized as a new set of features, called HHT features. In other words, our
method drastically expands the feature set available for machine learning.
Compared to previous studies, our main contribution is a novel approach to
incorporates HHT features with end effect correction into machine learning
methods for extrapolating prediction on financial time series. We also studied
the feature selection problem in using HHT for machine learning.

We experiment the methods using daily data for the S&P 500 index, CBOE
volatility (VIX) index, SPDR gold exchange-traded fund (GLD), and 10-year
treasury yield (TNX). We show that while these time series have contrasting
path behaviors (upward trending vs mean-reverting), their instantaneous
energy-frequency spectra reveal that they share very similar average
frequencies on both long and short timescales. In feature selection part, our
comparison analysis identifies the most useful HHT features for prediction. We
highlight the improved prediction performance of our method when benchmarked
to the prediction using only the original time series. We then discussed the
end effect of EMD and proposed a novel framework of extrapolating forecasting
with end effect correction. The prediction was implemented on the above-
mentioned dataset. Our results are potentially useful not only for prediction
but also portfolio construction [24].

The rest of the paper is structured as follows. First, we review the methods
of EMD and CEEMD and provide the algorithms for IMF extraction in Sect. 2. The
implementation of the methods is illustrated with real financial data in Sect.
2.5. We also discuss the end effect of the decomposition in Sect. 2.6 . Next
in Sect. 3, we introduce HHT features for financial machine learning. The
training and evaluation framework for extrapolating prediction is proposed. We
also discuss the incorporation of HHT features into several machine learning
models, including regression tree ensemble (RTE), support vector machine
(SVM), and especially the structure design of long short-term memory (LSTM)
neural network. In Sect. 4 we discuss the forecasting performance of the HHT-
enhanced machine learning (ML) models and feature selection of the HHT
features and IMF components. Section 4.2 discusses the end effect problem and
extrapolating forecasting with end effect correction.

##  2 Methodology

In this section, we present our methodology for processing a financial time
series and generating new features for machine learning. The major components
are the complementary ensemble empirical mode decomposition (CEEMD) and
Hilbert spectral analysis. We also discuss a procedure for time series
filtering and reconstruction from the multiscale decomposition.

###  2.1 Empirical Mode Decomposition

EMD is the first step of our multistage procedure. For any given time series
$x(t)$ observed over a period of time $[0,T]$, we decompose it in an iterative
way into a finite sequence of oscillating components
$c_{1}(t),\cdots,c_{n}(t)$, plus a nonoscillatory trend called the residue
term $r_{n}(t)$. Precisely, we have

$x(t)=\sum_{j=1}^{n}c_{j}(t)+r_{n}(t).$ (1)

To ensure that $c_{j}(t)$ have the proper oscillatory properties, the concept
of intrinsic mode function (IMF) is applied [1]. The intrinsic mode functions
are real functions in time that admit well-behaved and physically meaningful
Hilbert transform. Specifically, each IMF is defined by the following two
criteria:

  * •

No local oscillation: the number of extrema and the number of zero crossings
must be equal or at most differ by one.

  * •

Symmetric: the maxima of the function defined by the upper envelope and the
minima defined by the lower envelope must sum up to zero at any time
$t\in[0,T]$.

The definition of an IMF guarantees pure oscillation while allowing time-
dependent frequency and amplitude. Mathematically, we can express am IMF
$c(t)$ as

$c(t)=a(t)\cos(\theta(t)),$ (2)

where $a(t)\geq 0$ is the instantaneous amplitude, and $\theta(t)$ is the
phase function with $\theta^{\prime}(t)\geq 0$. Such property makes it ready
for time-frequency analysis under the Hilbert spectral analysis, which will be
discussed in Sect. 2.3.

As is standard [1, 25], we consider a sifting process that decomposes any time
series into a finite set of IMFs that oscillate on different timescales, plus
a nonoscillatory residue term. The key idea of the method is as follows: look
for the finest oscillation by finding all the local maxima and minima, and
then subtract the remaining trend, until the oscillation satisfies the IMF
conditions. Each IMF discovered is removed sequentially from the time series
until a nonoscillatory residue remains. The residue is a constant or monotonic
function, or has at most one maximum or minimum. The algorithm is summarized
as follows:

  * •

Initialize the residue term as $r_{0}(t)=x(t)$ and set $j=1$.

  * •

While $r_{j-1}(t)$ does not satisfies the nonoscillatory condition, i.e.
$r_{j-1}(t)$ has more than one maximum or minimum, do the following sifting
process to extract the oscillation as an IMF denoted by $c_{j}(t)$.

  * •

Initialize the component as $c_{j}(t)=r_{j-1}(t)$.

    * –

Interpolate the maxima of $c_{j}(t)$ using a cubic spline as the upper
envelope $u(t)$, and interpolate the minima of $c_{j}(t)$ using a cubic spline
as the lower envelope $l(t)$. Compute the mean of the upper and lower
envelopes $m(t)=\frac{1}{2}\left(u(t)+l(t)\right)$.

    * –

Iterate $c_{j}(t)\leftarrow c_{j}(t)-m(t)$.

    * –

Stop when $c_{j}(t)$ satisfies the criteria of an IMF. Let $c_{j}(t)$ be the
$j$-th component, and iterate $r_{j}(t)=r_{j-1}(t)-c_{j}(t)$ and $j\leftarrow
j+1$.

  * •

Return the IMFs $c_{1}(t),\cdots,c_{n}(t)$, and residue $r_{n}(t)$.

The sifting process described above is purely empirical and adaptive with very
little assumption on the temporal change, making it ideal for nonstationary
time series.

###  2.2 Complementary Ensemble Empirical Mode Decomposition

The EMD algorithm is designed to decompose a time series into multiple
components from high frequency to low frequency. However, the phenomenon of
mode mixing, first illustrated in [9], could occur in some cases. Mode mixing
is defined as either one IMF consisting of widely disparate scales, or signals
of similar scales reside in several IMF components [7]. Such a phenomenon
poses potential challenges on the interpretation of the IMFs. This problem is
particularly relevant for financial time series when high degree of
nonstationarity and noise are observed.

To resolve the mode mixing issue, we consider the ensemble empirical mode
decomposition (EEMD) proposed by [7] for our analysis. EEMD is a noise-
assisted signal processing technique that extracts each mode from an ensemble
mean computed based on $N$ trials. In each trial $i$, an i.i.d. white noise
$w_{i}$ with a zero mean and finite variance $\sigma$ is added to the original
time series $x(t)$, and $x(t)+w_{i}(t)$ is referred as the _signal_ in this
trial. The original EMD algorithm is then applied to the signal, outputting
the IMFs $c_{ij}(t),\ j=1,\cdots,n$, and the residual term $r_{in}(t)$.
Finally, the ensemble mean of the IMFs and residual terms across all the $N$
trails is regarded as the true mode extraction. The resulting ensemble mean of
the IMF components are given by

$c_{j}(t)=\frac{1}{N}\sum_{i=1}^{N}c_{ij}(t),\quad
r_{n}=\frac{1}{N}\sum_{i=1}^{N}r_{in}(t).$ (3)

Using (1), we can write

$\sum_{j=1}^{n}c_{ij}(t)+r_{in}(t)=x(t)+w_{i}(t).$ (4)

Then, the ensemble mean of the IMF components sums up to

$\displaystyle\sum_{j=1}^{n}c_{j}(t)+r_{n}(t)$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}(x(t)+w_{i}(t))=x(t)+\frac{1}{N}\sum_{i=1}^{N}w_{i}(t)\,.$

Therefore, the components from the ensemble mean provide a decomposition which
converges to the original time series $x(t)$ almost surely at the rate of
$\mathcal{O}(\frac{1}{\sqrt{N}})$. As such, a large ensemble size is typically
desired if possible. Implementing EMD with a large number of trials can be
potentially prohibitive in terms of computational cost and speed. To address
this issue, an alternative called the complementary ensemble empirical mode
decomposition (CEEMD) method is introduced [8]. CEEMD expands the ensemble by
adding a complementary negative noise $-w_{i}(t)$ for each trial, expanding
the total ensemble size to $2N$. Again, the components from the ensemble mean
sum up to equal the original time series:

$\sum_{j=1}^{n}c_{j}(t)+r_{n}(t)=x(t)+\frac{1}{2N}\sum_{i=1}^{N}(w_{i}(t)-w_{i}(t))=x(t).$
(5)

This holds regardless of the choice of $N$, thus reducing the need to have a
very large ensemble size. Besides reducing mode mixing, both EEMD and CEEMD
are also more robust to intrinsic noise in the original data, as these methods
automatically average out extra independent noises in the process.

Due to the high non-stationarity and noise level in financial time series, we
consider using CEEMD as the noise-assisted decomposition to resolve the mode
mixing problem. For the rest of the paper, we use CEEMD to extract the IMF
components from time series. Further HHT features such as instantaneous
frequency and amplitude are also derived from the CEEMD result.

###  2.3 Hilbert Spectrum

By definition, an IMF lends itself to the second stage of HHT, Hilbert
spectral analysis. An oscillating real-valued function can be viewed as the
projection of an orbit on the complex plane onto the real axis. For any
function in time $X(t)$, the Hilbert transform is given by

$Y(t)=\mathcal{H}[X](t):=\frac{1}{\pi}\int_{-\infty}^{+\infty}\frac{X(s)}{t-s}\,ds,$
(6)

where the improper integral is defined as the Cauchy principle value. The
transform exists for any function in $L^{p}$ [26]. As a result, $Y(t)$
provides the complementary imaginary part of $X(t)$ to form an analytic
function in the upper half-plane defined by

$Z(t)=X(t)+iY(t)=a(t)e^{i\theta(t)},$ (7)

where

$a(t)=\|Z(t)\|=\sqrt{X^{2}(t)+Y^{2}(t)}$ (8)

and

$\theta(t)=\arg Z(t)=\arctan(\frac{Y(t)}{X(t)}).$ (9)

As is standard [27, 28], for a function of the form in(2), if the amplitude
$a(t)$ and the frequency $\theta^{\prime}(t)$ are slow modulations, the
Hilbert transform will give a $\pi/2$ shift to the phase $\theta(t)$.
Therefore, the $a(t)$ given by (8) is exactly the instantaneous amplitude, and
the $\theta(t)$ given by (9) is exactly the instantaneous phase function. The
instantaneous frequency is then defined as the $2\pi$-standardized rate of
change of the phase function, that is,

$f(t)=\frac{1}{2\pi}\dot{\theta}(t)=\frac{1}{2\pi}\frac{d}{dt}\left(\arctan(\frac{Y(t)}{X(t)})\right).$
(10)

By the definition of IMF, it has a well-defined Hilbert transform, which means
that $a(t)$ and $f(t)$ computed reflect the amplitude and frequency of the
oscillation, thus the Hilbert transform is physically meaningful. Applying
Hilbert transform to each of the IMF components individually yields a sequence
of analytic signals [1]:

$c_{j}(t)+i\mathcal{H}[c_{j}](t)=c_{j}(t)+i\hat{c}_{j}(t)=a_{j}(t)e^{i\theta_{j}(t)},$
(11)

for $j=1,\cdots,n$. We refer to the $c_{j}(t)+i\hat{c}_{j}(t),\ j=1,\cdots,n$
as the complex IMFs, where the IMF components are the real parts, and the
corresponding Hilbert transforms serve as the imaginary parts. In turn, the
original time series can be represented as

$x(t)=\mathfrak{Re}\ \sum_{j=1}^{n}a_{j}(t)e^{i\int^{t}2\pi
f_{j}(s)ds}+r_{n}(t).$ (12)

This decomposition can be seen as a sparse spectral representation of the time
series with time-varying amplitude and frequency. In other words, each IMF
represents a generalized Fourier expansion that are suitable for nonlinear and
nonstationary financial time series. In summary, the procedure generates a
series of complex functions that are analytic in time, along with their
associated instantaneous amplitudes and instantaneous frequencies. These
components capture different time scales and resolutions embedded in the time
series and are used for time series filtering and reconstruction.

Lastly, the Hilbert spectrum is defined by

$H(f,t)=\sum_{j=1}^{n}H_{j}(f,t),\quad\mbox{where}\ {\displaystyle
H_{j}(f,t)={\begin{cases}a_{j}(t),&f=f_{j}(t),\\\
0,&{\text{otherwise.}}\end{cases}}}$ (13)

The instantaneous energy of the $j$-th component is defined as

$E_{j}(t)=|a_{j}(t)|^{2}.$ (14)

We examine the behavior of the Hilbert spectrum through the pair
$(f_{j}(t),E_{j}(t))$ (see Fig. 7). The Hilbert spectrum will also be
incorporated into machine learning models for time series forecasting.

###  2.4 Time Series Filtering and Reconstruction

Financial time series often exhibits characteristics on different timescales,
from long-term economic cycles to daily rapid fluctuations. Moreover,
financial data can be noisy, so we used a noise-added approach to make the
decomposition more robust. CEEMD decomposes the time series onto different
timescales, which sum up to the original time series.

By construction, the iterative sifting process identifies the finest
structures, and then extract longer and longer scales. As a result, the first
few components have higher frequencies which are more noisy, and the last few
components have lower frequencies representing long term structures.
Considering the noisy and complex nature of financial time series, such scale
separation is valuable. In fact, EMD can be used as a filter for time series
[17]. In the reconstruction of the original time series using the IMF
components, we can choose a subset of modes as a filter for desired
information. By removing the first few high frequency components, we create a
low-pass filter; that is,

$x_{L}^{(m)}(t)=\sum_{j=m}^{n}c_{j}(t)+r_{n}(t).$ (15)

This reconstruction using only the last few components can serve as a
smoothing of the time series. Similarly, we can also build a high-pass filter
with

$x_{H}^{(m)}(t)=\sum_{j=1}^{m}c_{j}(t),$ (16)

which captures the high-frequency local behaviors, and can also be used to
estimate the noise level or volatility.

###  2.5 Implementation with Financial Data

In this section, we use four financial time series of different type and in
different sectors for implementation, namely the S&P 500 index, the CBOE
volatility index VIX, the gold price exchange-traded fund (ETF) GLD, and the
10-year treasury yield (TNX), over a 10-year period from 4/1/2010 to
3/31/2020. For S&P 500 and GLD, let $s(t)$ be the value of the time series and
$x(t):=\log(s(t))$ be its log price. For VIX and 10-year treasury rate, let
$x(t)$ be the time series itself. In each case, CEEMD is applied to $x(t)$ to
extract the IMFs and residual from the time series, the Hilbert spectrum of
which is later computed with HHT. To reduce the error in instantaneous
frequency computation, we have used the robust locally weighted scatterplot
smoothing (robust LOWESS) [29]. The energy-frequency spectrum is then shown to
display the instantaneous values of $(f_{j}(t),E_{j}(t))$.

In Fig. 1, we implement CEEMD using the 10-year historical data for the
S&P500, GLD, VIX, and 10-year treasury yield. The top row in each plot shows
the original time series $x(t)$, followed by the IMFs with decreasing
frequencies, ended with the residual term showing the overall trend. The modes
on different timescales can indicate different behavior of the time series.
Taking the recent COVID-19 event for example, when S&P 500 began the sharp
decline from 2/21/2020 on, modes 3, 4, 5 and 6 also show notable peak around
the time period, but on earlier and earlier dates (2/25/2020, 1/21/2020,
12/26/2019, and 11/15/2019). This indicates that a mode on a longer timescale
signals trends earlier than other modes. The opposite behavior is observed for
VIX – the lower-frequency modes (modes 3 to 6) exhibit a steep increase
towards the day when VIX spikes.

Figure 1: IMFs and residual terms extracted from CEEMD of S&P 500 (top left),
VIX (top right), GLD (bottom left), and 10-year treasury yield (bottom right),
from 4/1/2010 to 3/31/2020. The top row in each plot shows the original time
series. The second to last but one rows show the IMF modes of the
corresponding time series. The bottom row of each plot shows the residual term
of the time series.

In Fig. 2, we illustrate the low-pass filtered reconstruction of the four time
series. This involves applying (15) using different collections of components.
Specifically, we have used the last 3 components and last 5 components
including the residual terms. We can see that, with more components included
in the reconstruction, the resulting time series resembles the original time
series on a finer timescale. Compared to some other time series smoothing
techniques, such as moving average, the CEEMD-based filter achieves smoothing
without lags. Furthermore, CEEMD is capable of reconstructing time series with
both trending or mean-reverting dynamics as shown in Fig. 2.

Figure 2: Reconstruction of financial time series: S&P 500, VIX, GLD and
10-year treasury yield. For S&P 500 and GLD, the decomposition was implemented
on the log price of the indices. The low-pass filtered time series by (15)
were then taken exponential to show the reconstruction.

Using formula (6), we implement Hilbert transform to the IMFs extracted from
CEEMD. Following (11), we take the IMF components as the real parts and the
corresponding Hilbert transform as the imaginary parts, and plot them together
in Fig. 3 as the complex IMFs which are analytic signals. We see that each
real part (IMF component) and the corresponding imaginary part (Hilbert
transform) are oscillations coupled together, where negative (positive)
Hilbert transform corresponds to a increasing (decreasing) IMF component. The
coupling behavior indicates the complex IMFs can be useful for time series
forecasting.

Figure 3: Complex IMFs derived from formula (6), S&P 500 (top left), VIX (top
right), GLD (bottom left), and 10-year treasury yield (bottom right). The real
parts (IMF components) are plotted together with their corresponding imaginary
parts (Hilbert transforms of IMFs).

Figure 4: The instantaneous frequency associated with each IMF component
derived from formula (10) of S&P 500, VIX, GLD, and TNX. From top to bottom,
the instantaneous frequency correspond to modes with lower frequencies.

The instantaneous frequency $f_{j}(t)$ from (10) associated with each IMF
component $c_{j}(t)$ of the four financial time series are presented in Fig.
4. From the 1st to 6th mode (top to bottom), the fluctuations of instantaneous
frequency tend to decrease, and the overall levels of instantaneous frequency
are reduced. Fig. 5 shows the instantaneous amplitude $a_{j}(t)$ from (8)
associated with each IMF component $c_{j}(t)$ of the four financial time
series. From the 1st to 6th mode (top to bottom), the fluctuations of
instantaneous amplitude is reduced significantly. As such, the higher-order
modes reflects a clearer long-term trend compared to lower-order modes. The
fluctuations of the instantaneous amplitude and frequency are also less rapid
than the IMF component itself, satisfying the slow modulation condition.

Figure 5: The instantaneous amplitude associated with each IMF component
derived from formula (8) of S&P 500, VIX, GLD, and TNX. From top to bottom,
the instantaneous amplitude correspond to modes with lower frequencies.

###  2.6 End Effect of EMD

EMD-based decomposition methods and HHT are adaptive for nonstationary
systems. However, the algorithms are also known for the issue of end effect,
whereby the decomposition and Hilbert transform have higher error when getting
closer to the ends of the time series. This is because the sifting algorithm
in Section 2.1 uses interpolation in each iteration, requiring an anchor point
beyond the time span. The end effect can be relieved when the time series has
enough length, and we are only interested in the properties in the interior.
However, the error at the end can be a potential problem and may an additional
challenge for forecasting.

To characterize the end effect, we need to consider the decomposition with the
position in the time frame. Let’s define the end effect factor as the position
of point $t$ on the decomposition interval $[T_{1},T_{2}]$:

$\lambda(t)=\frac{2t-(T_{1}+T_{2})}{T_{2}-T_{1}}.$ (17)

Note that $\lambda=0$ at the middle point of the time frame, and $\lambda=\pm
1$ at the ends. Suppose the time series has the true IMF mode decomposition
$x(t)=\sum_{j=1}^{n}c^{*}_{j}(t)+r^{*}_{j}(t)$, then the decomposition found
by the CEEMD will be:

$(c_{1}(t),\cdots,c_{n}(t),r_{n}(t))\sim\mathcal{D}(c^{*}_{1}(t),\cdots,c^{*}_{n}(t),r^{*}_{n}(t),\lambda(t)),$
(18)

which is a draw from the distribution due to the random nature of the noise-
assisted algorithm. Denote the shorthand notice
$\mathbf{c}(t)=(c_{1}(t),\cdots,c_{n}(t),r_{n}(t))$. The error vector of
$\mathbf{c}(t)$ is:

$\mathcal{E}(\mathbf{c}^{*}(t),\lambda(t))=\mathcal{D}(\mathbf{c}^{*}(t),\lambda(t))-\mathbf{c}^{*}(t).$
(19)

Note that $\mathbf{1}^{T}\mathcal{E}(\mathbf{c}^{*}(t),\lambda(t))=0$ at any
$t$, due to the exact decomposition in (5), and
$Var(\mathcal{E}_{j}(\mathbf{c}^{*},\lambda))$ increase with $|\lambda|$ due
to the end effect.

###  2.7 Empirical Energy-Frequency Spectrum

For each financial time series, there is an instantaneous frequency profile.
This allows us to compare the frequencies between any pair of indices, such as
S&P 500 vs VIX in Fig. 6. Each point on the plot is the pair of instantaneous
frequencies of S&P 500 and VIX recorded at the same time. For each mode, there
is one cluster of points (instantaneous frequencies) and the black cross marks
the mean value of each cluster. These marks are located extremely close to the
solid straight line of unit slope, which is a reference line for identical
frequencies. This means that, while S&P 500 and VIX have drastically different
path behaviors, their dynamics share very similar frequencies on average
across both long and short timescales. Similar observations are found for GLD
when compared against VIX.

Figure 6: Instantaneous frequency: S&P 500 vs VIX (left) and GLD vs VIX
(right).

Next, we show in Fig. 7 the instantaneous energy-frequency spectrum of the S&P
500 (top left), VIX (top right), GLD (bottom left), and 10-year treasury yield
(bottom right) across all IMF components. Each point on the plots is a pair of
$(f_{j}(t),E_{j}(t))$, for mode $j=1,\cdots,n$, and time $t=1,\cdots,T$. For
each mode, the black cross marks the mean of each cluster of points. The
clusters from the right to left correspond to higher modes.

For S&P 500, from IMF 1 to 6, the average instantaneous frequency decreases
from $2.122\times 10^{-1}$ to $1.287\times 10^{-3}$ respectively, indicating
cycles on timescales from 5 trading days (weekly behavior) to over 750 trading
days (about 3 years). Since the decomposition is implemented over a 10-year
period, the residual term reflects the long-term trend, possibly driven by
economic cycles. The solid line obtained from linear regression has a negative
slope for all the four financial time series. This means that the
instantaneous energy $E(t)$ tends to decrease as the instantaneous frequency
$f(t)$ increases.

Comparing the four indices, S&P 500, GLD and 10-year treasury rate exhibit a
more rapid dissipation of energy with respect to frequency than VIX. This
reveals that the rapid fluctuations in S&P 500, GLD and 10-year treasury rate
tend to have much smaller magnitude than those for VIX. Moreover, S&P 500 has
a significantly wider range of instantaneous energy than VIX. This is
consistent with the fact that S&P500 has risen from 1100 to 3000 while VIX has
stayed between 10 and 40 over the 10-year period. In terms of power spectrum,
S&P 500, GLD and 10-year treasury rate are closer to the famous $1/f$ power
spectrum, while VIX shows significantly different behavior. More
interestingly, the residue with respect to the linear regression on the log-
log plot of VIX shows significant concavity. These difference all indicate the
volatility process follows a very different mechanism than the price
processes.

For higher modes (lower frequencies), the instantaneous energy and frequency
have smoother paths over time, giving rise to more continuous trajectories on
the left-side of the plot. This is consistent with the fact that the
fluctuation of instantaneous amplitude and instantaneous frequency are reduced
from lower to higher modes, as shown in Fig. 5 and 4. Using the language of
dynamical systems, we can view the energy-frequency spectra as movements on
the phase plane $(f,E)$, which serves as an alternative representation of the
dynamics.

Figure 7: Instantaneous energy-frequency spectrum for S&P 500 (top left), VIX
(top right), GLD (bottom left), and 10-year treasury yield (bottom right). The
black cross marks the mean of the cluster of points for each mode.

###  2.8 Comparison to Other Methods

We end this section by highlighting the differences among the HHT-based
method, Fourier transform, and wavelet transform. As summarized in Table 1,
HHT allows for nonlinearity, nonstationarity. The HHT bases are determined in
adaptive form with respect to the time series, rather than being prespecified
as in the Fourier and wavelet transforms. In addition, HHT produces a sparse
spectrum with finitely many of modes, which is a useful property for our study
of financial time series. The adaptiveness and sparsity also allow for
practical feature generation for time series forecasting, which is the main
reason for our use of HHT in this study.

| Nonstationarity | Nonlinearity | Basis | Spectrum
---|---|---|---|---
Fourier | No | No | A priori | Global; Dense
Wavelet | Yes | No | A priori | Regional; Dense
HHT | Yes | Yes | Adaptive | Local; Sparse
Table 1: Comparison between different transforms for time series.

##  3 Machine Learning Using HHT Features

The implementation of HHT-enhanced ML predictions involves two main steps. In
the first step, HHT is applied to the original time series $x(t)$. CEEMD
decomposes $x(t)$ into $n$ IMFs $c_{1}(t),\cdots,c_{n}(t)$ plus the residual
term $r_{n}(t)$, and then Hilbert transform on the IMFs yields the
corresponding imaginary counterparts $\hat{c}_{1}(t),\cdots,\hat{c}_{n}(t)$,
along with the instantaneous amplitude $a_{1}(t),\cdots,a_{n}(t)$, and
instantaneous frequency $f_{1}(t),\cdots,f_{n}(t)$. These outputs represent
the _HHT feature set_ derived from the time series $x(t)$ (see Fig. 8). The
following sections discuss the second step of training and testing machine
learing models on HHT features.

Figure 8: The collection of HHT features derived from the original time series
$x(t)$.

###  3.1 Training and Testing

After the extracting of the HHT features, a machine learning model is chosen
to analyze the times series and make predictions for future time periods.
Specifically, the goal is to predict the change in the next period $\Delta
x(t+1):=x(t+1)-x(t)$, going forward in time $t$. Traditional time series
prediction takes the values of $x$ in the past $\tau$ periods
$x^{(t,\tau)}:=(x(t-\tau+1),\cdots,x(t))$ as input. In contrast, our approach
uses only the HHT features as predictors, rather than the original time
series. To this end, we use the following shorthand notations:

$\displaystyle c^{(t,\tau)}_{j}$
$\displaystyle:=(c_{j}(t-\tau+1),\cdots,c_{j}(t)),$ (20) $\displaystyle
c^{(t,\tau)}_{(j,l)}$
$\displaystyle:=(c^{(t,\tau)}_{j},\ldots,c^{(t,\tau)}_{k}),$ (21)

and likewise for other HHT features $\hat{c},a$ and $f$.

The objective is to minimize the expected $l_{2}$-loss of $\Delta x(t+1)$
prediction, given information up to time $t$. This leads to the optimization
problem:

$\min_{g\in\mathcal{G}}\mathbb{E}[(\Delta
x(t+1)-g(c^{(t,\tau)}_{(1,n)},\hat{c}^{(t,\tau)}_{(1,n)},a^{(t,\tau)}_{(1,n)},f^{(t,\tau)}_{(1,n)}))^{2}|\mathcal{F}_{t}],$

where $\mathcal{G}$ is the functional class specified by the ML model so that
each $g$ is a prediction function determined within the training period. The
filtration of information up to time $t$ is denoted by $\mathcal{F}_{t}$. For
implementation with real data, we do not know the distribution of input
features and the target output, so the expectation is estimated by the sample
average. Specifically, we split the time series into a training section from
$t=1$ to $t=T_{1}$, and a testing section after that from $t=T_{1}+1$ to
$t=T_{1}+T_{2}$. The ML models learn on the training set by solving

$\min_{g\in\mathcal{G}}\frac{1}{T_{1}}\sum_{t=1}^{T_{1}}(\Delta
x(t+1)-g(c^{(t,\tau)}_{(1,n)},\hat{c}^{(t,\tau)}_{(1,n)},a^{(t,\tau)}_{(1,n)},f^{(t,\tau)}_{(1,n)}))^{2}.$

The optimal predictive function found, $g_{*}$, is then tested on the unseen
data. The performance is measured in terms of mean squared error (MSE):

$\frac{1}{T_{2}}\sum_{t=T_{1}+1}^{T_{1}+T_{2}}(\Delta
x(t+1)-g_{*}(c^{(t,\tau)}_{(1,n)},\hat{c}^{(t,\tau)}_{(1,n)},a^{(t,\tau)}_{(1,n)},f^{(t,\tau)}_{(1,n)}))^{2}.$

###  3.2 Extrapolating Prediction

The training and testing framework discussed above is useful in validating the
predictive power of HHT and perform feature selection. However, in real-time
forecasting, we are only able to do the decomposition up to the current time
$T$, and the testing is the one-shot extrapolation into the future $T+1$ where
no decomposition is available. The above-mentioned procedure and most of the
previous studies suffered from information leakage from the future that
enhanced the prediction at the current time, due to the interpolation property
described in Section 2.6.

As discussed above, the traditional train/test split does not work for time
series prediction using HHT features. To test the extrapolation power of HHT
features, we propose the following procedure of training and testing:

  * •

For $t=T_{1}+1,\cdots,T_{1}+T_{2}$:

    * –

Implement CEEMD and HHT on $x(s)$ for $s\in[t-T,t-1]$.

    * –

Train machine learning model on $[t-T,t-1]$.

$\min_{g\in\mathcal{G}}\frac{1}{T}\sum_{s=t-T}^{t-1}(\Delta
x(s)-g(c^{(s-1,\tau)}_{(1,n)},\hat{c}^{(s-1,\tau)}_{(1,n)},a^{(s-1,\tau)}_{(1,n)},f^{(s-1,\tau)}_{(1,n)}))^{2}.$

    * –

Test one-shot prediction of
$\hat{x}(t)=g_{*}(c^{(t-1,\tau)}_{(1,n)},\hat{c}^{(t-1,\tau)}_{(1,n)},a^{(t-1,\tau)}_{(1,n)},f^{(t-1,\tau)}_{(1,n)})$.

  * •

Evaluate testing mean square error
$\frac{1}{T_{2}}\sum_{t=T_{1}+1}^{T_{1}+T_{2}}(x(t)-\hat{x}(t))^{2}$.

###  3.3 Machine Learning Models

Our approach is generally applicable to all ML models that take
multidimensional features as inputs. To illustrate our methodology, we
implement several ML models for time series predictions. Specifically, we
tested three different types of nonlinear regression models, namely the
regression tree ensemble (RTE), the support vector machine (SVM) regression,
and the long short-term memory (LSTM) neural network.

The algorithm framework for RTE and SVM are generically applied to
classification or regression tasks. The models are described in [30, 31] for
RTE with least square boosting and in [32, 33] for SVM for regression. The
LSTM model, on the other hand, requires design of the network structure. In
this section, we describe the network structure we used for this study, which
is visualized in Fig. 9

As a type of recurrent neural network (RNN), LSTM learns the long-term
dependency in sequenced data [34]. At each time step $t$, the LSTM block takes
in new observations $c_{j}(t),\hat{c}_{j}(t),\omega_{j}(t),a_{j}(t),\
j=1,\cdots n$. The new inputs together with the network state information from
the last step generate the hidden units $h_{1}(t),\cdots,h_{M}(t)$. The LSTM
block and the hidden units then propagate to the next LSTM block. For deeper
networks, we apply a drop-out layer to the hidden layers and then feed into
another layer of LSTM.

Our neural network structure consists of two layers of LSTM each with a drop-
out layer. The network then links to two fully-connected layers, each with a
ReLU layer for activation. The number of hidden units in the layers are as
follows: 100 - 200 for the first LSTM layer, 100 - 200 for the second LSTM
layer, 50 - 100 for the first fully-connected layer, and 50 - 100 for the
second fully-connected layer. For each data set, we tune the number of hidden
units via cross-validation.

For the regression task to predict $y(t)=\Delta x(t+1)$, the hidden units
$h(t)$ are connected to the output $y(t)$ through a multi-layer neural
network. In our study, we used multiple fully-connected layers with the
activation function of rectified linear unit (ReLU). The reason to use ReLU is
to create nonlinearity in the prediction and induce sparsity to choose a
subset of useful signals. To train the neural network, we use the Adam solver
[35] and train multiple times with different randomization. The model with the
lowest root mean square error (RMSE) on the training set is used for testing.

Figure 9: The structure of LSTM with HHT features.

##  4 Forecasting Experiments

In this section, we implement the machine learning forecasting experiments on
the historical prices of S&P 500 index, CBOE volatility index (VIX), SPDR gold
ETF (GLD), and the 10-year treasury yield (TNX). In the first part, we
experiment feature selection for the HHT features and IMF components. The
prediction follows the standard decomposition and train/test split framework.
In Section 4.2, we discuss the extrapolation prediction and the associated end
effect. The extrapolating prediction using LSTM is implemented.

###  4.1 Features and IMFs Selection

We implement the ML models using HHT features during the 10-year period
4/1/2010–3/31/2020. The data are split into a training set 4/1/2010–3/31/2018
and an unseen test set 4/1/2018–3/31/2020. The goal is short-time prediction
for the next day change $\Delta x(t+1)$, based on the features up to time $t$.
The comparison is between using the original time series feature
$x^{(t,\tau)}$, and with the HHT feature information
$c^{(t,\tau)}_{(1,n)},\hat{c}^{(t,\tau)}_{(1,n)},a^{(t,\tau)}_{(1,n)},f^{(t,\tau)}_{(1,n)}$.
We set the lookback window $\tau=5$ for regression tree ensemble and SVM,
while LSTM takes all the past features into the recurrent neural network, and
each LSTM block controls the past information flow and memory level. We train
the three ML models using four separate sets of features:

  1. 1.

original time series: $x(t)$,

  2. 2.

IMF components: $c_{1}(t),\cdots,c_{n}(t)$,

  3. 3.

complex IMF: $c_{1}(t),\hat{c}_{1}(t),\cdots,c_{n}(t),\hat{c}_{n}(t)$,

  4. 4.

full HHT features:
$c_{1}(t),\hat{c}_{1}(t),f_{1}(t),a_{1}(t),\cdots,c_{n}(t),\hat{c}_{n}(t),f_{n}(t),a_{n}(t)$,

For each time series, HHT produces a large set of interpretable features
corresponding to different time scales. The feature size is increased by $4n$
times, where $n$ is the number of IMF components. Through our experiments, we
determine the most useful set of features for predictions using different ML
methods, namely regression tree ensemble (RTE), support vector machine (SVM),
and long short-term memory (LSTM) neural network. For performance comparison,
we measure the mean squared errors (MSE) of the predictions over the whole
horizon, based on the four financial time series S&P 500, VIX, GLD, and
10-year treasury rate.

Tables 2 – 5 summarize the experimental results from using RTE, SVM, and LSTM
based on S&P 500 (log price), VIX index, GLD (log price) and 10-year treasury
rate respectively. As benchmarks, we show the MSEs from predictions with the
original time series $x$ as input feature, using the corresponding machine
learning model. In all the rest parts of tables, the predictions are made
using only HHT features and do not use the original time series $x$ directly.
Each row corresponds to using only a subset of HHT features, and each column
corresponds to the set of IMFs included. The first and last few components are
also referred to as the high-pass (high-frequency) and low-pass (low-
frequency) filters respectively.

Tables 2–5 show that generally the largest feature set (bottom row) lead to
much reduced errors. This suggests that HHT is capable of identifying the most
important features from a time series. Surprisingly, the feature pair
$(c,\hat{c})$, which consists of only the IMF components and the complementary
imaginary parts derived from Hilbert transform in (6), yields similar or
better performance than a much larger feature sets, e.g. ($c,\hat{c},f,a$),
which indicates that the most useful information in HHT is included in the
complex IMFs. Considering the much smaller feature set, the complex IMFs are
more preferable in terms of information criterion. Both feature sets
outperform the feature set with only real IMF components, which has been
explored in many previous studies of machine learning prediction with EMD.

In terms of IMF selection, all three ML models perform better with all or the
first few (high-frequency) IMFs than the last (low-frequency) IMFs. In many
cases, the machine learning model performs best when the last IMFs are
excluded (first 5). This is intuitive since low-frequency components may not
be as useful as high-frequency components for predicting short-term movements,
and a machine learning model with high complexity can potentially overfit on
the training set with redundant information.

Comparing different machine learning models, we see that LSTM does not perform
as good as the other more generic models in some scenarios. This is due to the
fact that the prediction framework is under the train/test split within the
same decomposition. The generic models lose predictive power compared to LSTM,
as we get into extrapolating prediction in the next section.

Model: RTE
---
Features | All IMFs | First 5 | First 4 | Last 5 | Last 4
$x$ (original time series) | 2.3036
$c$ | 1.2614 | 1.2185 | 1.3640 | 2.1418 | 2.2442
$c,\hat{c}$ | 0.8203 | 1.0508 | 0.9540 | 2.0625 | 2.2048
$c,f$ | 1.2845 | 1.1128 | 1.0747 | 2.0878 | 2.2539
$c,a$ | 1.2990 | 1.1378 | 1.0581 | 2.1720 | 2.2481
$c,f,a$ | 1.0137 | 1.3420 | 1.2366 | 2.0850 | 2.2892
$c,\hat{c},f,a$ | 0.8331 | 0.8713 | 1.0297 | 2.0815 | 2.2242
Model: SVM
Features | All IMFs | First 5 | First 4 | Last 5 | Last 4
$x$ (original time series) | 2.3313
$c$ | 2.0456 | 2.0545 | 2.0646 | 2.2685 | 2.3087
$c,\hat{c}$ | 1.9943 | 2.0214 | 1.9987 | 2.2421 | 2.2826
$c,f$ | 2.1161 | 1.9825 | 2.0030 | 2.2913 | 2.3100
$c,a$ | 2.1867 | 2.1490 | 2.1105 | 2.2889 | 2.3006
$c,f,a$ | 2.1842 | 2.0817 | 2.1089 | 2.2953 | 2.3064
$c,\hat{c},f,a$ | 2.1446 | 2.0589 | 1.9897 | 2.2769 | 2.3023
Model: LSTM
Features | All IMFs | First 5 | First 4 | Last 5 | Last 4
$x$ (original time series) | 2.3285
$c$ | 1.1400 | 1.2738 | 1.1537 | 2.0841 | 2.1510
$c,\hat{c}$ | 0.6348 | 0.5938 | 0.4876 | 1.9315 | 2.1239
$c,f$ | 1.2031 | 1.4041 | 1.0702 | 2.0743 | 2.1610
$c,a$ | 1.4561 | 1.1516 | 1.3415 | 2.0796 | 2.2130
$c,f,a$ | 1.2626 | 1.4305 | 1.1992 | 2.0707 | 2.1788
$c,\hat{c},f,a$ | 0.6834 | 0.9725 | 1.0058 | 1.9217 | 2.1233
Table 2: S&P 500 prediction error (MSE, $\times 10^{-4}$) with different
feature sets and components for three machine learning models (RTE (top), SVM
(middle), LSTM (bottom)). The benchmark (top row for each model) is the
prediction error from using the plain time series $x$ only. Model: RTE
---
Features | All IMFs | First 5 | First 4 | Last 5 | Last 4
$x$ (original time series) | 6.2796
$c$ | 2.8104 | 3.0372 | 3.9181 | 5.7601 | 6.2365
$c,\hat{c}$ | 1.9488 | 1.9982 | 1.8169 | 5.6598 | 6.1334
$c,f$ | 2.4096 | 2.4606 | 2.4603 | 5.7786 | 6.2631
$c,a$ | 3.7021 | 3.0241 | 3.9229 | 5.8065 | 6.2387
$c,f,a$ | 2.4866 | 2.7127 | 2.1297 | 5.7865 | 6.2145
$c,\hat{c},f,a$ | 2.1461 | 2.4642 | 1.5964 | 5.5183 | 6.1631
Model: SVM
Features | All IMFs | First 5 | First 4 | Last 5 | Last 4
$x$ (original time series) | 6.4127
$c$ | 5.9377 | 6.0057 | 5.6676 | 6.2427 | 6.2916
$c,\hat{c}$ | 5.9248 | 5.8020 | 5.6557 | 6.1872 | 6.3019
$c,f$ | 5.8231 | 5.7919 | 5.5276 | 6.2204 | 6.2988
$c,a$ | 6.2414 | 6.0873 | 5.9045 | 6.2802 | 6.2950
$c,f,a$ | 5.9618 | 6.0003 | 5.6598 | 6.2615 | 6.2985
$c,\hat{c},f,a$ | 5.8611 | 5.8534 | 5.3969 | 6.1940 | 6.3017
Model: LSTM
Features | All IMFs | First 5 | First 4 | Last 5 | Last 4
$x$ (original time series) | 6.3131
$c$ | 3.7682 | 3.9556 | 3.4594 | 5.5977 | 6.0955
$c,\hat{c}$ | 2.7004 | 2.5687 | 2.5492 | 5.3324 | 6.1218
$c,f$ | 5.4949 | 5.2656 | 5.2457 | 6.4108 | 6.2548
$c,a$ | 4.3976 | 4.3648 | 4.4501 | 6.2045 | 6.2070
$c,f,a$ | 7.0323 | 5.7723 | 4.9104 | 6.2129 | 6.4870
$c,\hat{c},f,a$ | 4.4516 | 4.6962 | 3.9247 | 5.8894 | 6.2638
Table 3: VIX index prediction error (MSE) with different feature sets and
components for three machine learning models (RTE (top), SVM (middle), LSTM
(bottom)). The benchmark (top row for each model) is the prediction error from
using the plain time series $x$ only. Model: RTE
---
Features | All IMFs | First 5 | First 4 | Last 5 | Last 4
$x$ (original time series) | 6.9816
$c$ | 2.7139 | 2.7363 | 2.5476 | 5.5666 | 6.7445
$c,\hat{c}$ | 1.0644 | 1.1821 | 1.0336 | 5.2963 | 6.4202
$c,f$ | 1.7681 | 1.9592 | 1.9780 | 6.0007 | 6.8471
$c,a$ | 3.1483 | 2.7417 | 2.5566 | 5.8702 | 6.7719
$c,f,a$ | 2.1982 | 2.0045 | 2.2453 | 5.7184 | 6.6546
$c,\hat{c},f,a$ | 0.9626 | 0.9145 | 0.8760 | 5.5237 | 6.5540
Model: SVM
Features | All IMFs | First 5 | First 4 | Last 5 | Last 4
$x$ (original time series) | 7.2381
$c$ | 3.8314 | 3.8783 | 3.6813 | 6.2574 | 7.3537
$c,\hat{c}$ | 2.6886 | 2.6979 | 2.7458 | 5.8944 | 6.8402
$c,f$ | 3.4269 | 3.4790 | 3.3913 | 6.0189 | 6.8042
$c,a$ | 4.5545 | 4.4103 | 4.3258 | 6.4397 | 7.0054
$c,f,a$ | 4.3807 | 4.1316 | 3.8366 | 6.5134 | 6.8918
$c,\hat{c},f,a$ | 3.0045 | 2.6902 | 2.6795 | 6.0785 | 6.8850
Model: LSTM
Features | All IMFs | First 5 | First 4 | Last 5 | Last 4
$x$ (original time series) | 6.9792
$c$ | 4.8152 | 3.4164 | 3.1468 | 5.7257 | 6.8414
$c,\hat{c}$ | 1.4325 | 1.3451 | 1.4513 | 5.1071 | 6.5997
$c,f$ | 2.8661 | 3.2673 | 2.8033 | 6.2854 | 6.8597
$c,a$ | 3.2262 | 3.1938 | 3.5799 | 5.7604 | 6.8049
$c,f,a$ | 3.2504 | 3.4774 | 3.6283 | 5.7533 | 6.7521
$c,\hat{c},f,a$ | 1.3953 | 1.5691 | 1.3248 | 5.1094 | 6.6098
Table 4: GLD prediction error (MSE, $\times 10^{-5}$) with different feature
sets and components for three machine learning models (RTE (top), SVM
(middle), LSTM (bottom)). The benchmark (top row for each model) is the
prediction error from using the plain time series $x$ only. Model: RTE
---
Features | All IMFs | First 5 | First 4 | Last 5 | Last 4
$x$ (original time series) | 2.6011
$c$ | 1.4615 | 1.2206 | 1.1800 | 2.2609 | 2.5347
$c,\hat{c}$ | 0.6244 | 0.6473 | 0.6545 | 2.1302 | 2.4472
$c,f$ | 1.2036 | 1.1138 | 1.0990 | 2.2331 | 2.5236
$c,a$ | 1.3131 | 1.3719 | 1.3811 | 2.2872 | 2.5832
$c,f,a$ | 1.0673 | 1.0318 | 1.0657 | 2.2663 | 2.5696
$c,\hat{c},f,a$ | 0.6299 | 0.6798 | 0.8444 | 2.1907 | 2.4539
Model: SVM
Features | All IMFs | First 5 | First 4 | Last 5 | Last 4
$x$ (original time series) | 2.7631
$c$ | 1.6064 | 1.5610 | 1.5968 | 2.4099 | 2.6059
$c,\hat{c}$ | 1.3893 | 1.3023 | 1.3737 | 2.2918 | 2.5383
$c,f$ | 1.9820 | 1.5531 | 1.4357 | 2.4251 | 2.5881
$c,a$ | 1.9156 | 1.7213 | 1.6954 | 2.4124 | 2.5908
$c,f,a$ | 2.0474 | 1.8175 | 1.6666 | 2.4550 | 2.5877
$c,\hat{c},f,a$ | 1.7609 | 1.4761 | 1.3948 | 2.3599 | 2.5489
Model: LSTM
Features | All IMFs | First 5 | First 4 | Last 5 | Last 4
$x$ (original time series) | 2.6071
$c$ | 1.1932 | 1.1545 | 1.1304 | 2.2402 | 2.4243
$c,\hat{c}$ | 0.4184 | 0.4348 | 0.5390 | 1.9995 | 2.4374
$c,f$ | 1.2942 | 1.1846 | 1.1972 | 2.2446 | 2.4780
$c,a$ | 1.1745 | 1.1426 | 1.0941 | 2.2503 | 2.4284
$c,f,a$ | 1.2464 | 1.1732 | 1.1668 | 2.2991 | 2.5692
$c,\hat{c},f,a$ | 0.4633 | 0.3816 | 0.4190 | 2.0137 | 2.4758
Table 5: 10-year treasury yield prediction error (MSE, $\times 10^{-3}$) with
different feature sets and components for three machine learning models (RTE
(top), SVM (middle), LSTM (bottom)). The benchmark (top row for each model) is
the prediction error from using the plain time series $x$ only.

###  4.2 Extrapolation and End Effect

As discussed in Section 2.6, the decomposition of $x(t)$ depends on the
position of $t$ in the time frame, characterized by the end effect factor
$\lambda(t)$ define in (17). Since the distribution of the decomposed
component is according to (18), depending on $\lambda$, the prediction model
then becomes

$x_{t+1}\sim
g(c^{(t,\tau)}_{(1,n)},\hat{c}^{(t,\tau)}_{(1,n)},a^{(t,\tau)}_{(1,n)},f^{(t,\tau)}_{(1,n)},\lambda(t)),$
(22)

which means we need to incorporate $\lambda$ into the predictors.

Furthermore, we notice that the distribution of the input and output in the
time series forecasting task is not homogeneous, which requires an evolving
prediction model like LSTM. The generic machine learning models RTE and SVM
have poor performance on the testing error as the distribution is different
from the training data. For this reason, we only show the results on LSTM with
the end effect correction factor $\lambda$ as the additional feature.

Table 6 shows the extrapolating prediction results with the end effect
correction. We see that the different HHT feature sets all improved the
prediction accuracy from using the plain time series and from the benchmark
naive prediction $x_{T+1}=x(T)$. Furthermore, the complex IMF components and
the full HHT feature set not only achieved lower mean square error, but also
required less hidden units in the network and less training time. This shows
the HHT features are better at explaining the dynamics and presenting the
data.

S&P 500 |  Benchmark: 5.5770 |  |
---|---|---|---
Features | MSEs | Avg. # hidden units | Avg. train time (s)
$x$ (original time series) | 5.5794 | 500.0 | 6.6166
$c$ | 5.0029 | 553.3 | 7.9682
$c,\hat{c}$ | 5.1114 | 373.1 | 5.3871
$c,\hat{c},f,a$ | 5.1444 | 322.6 | 4.1803
GLD |  Benchmark: 1.5698 |  |
Features | MSEs | Avg. # hidden units | Avg. train time (s)
$x$ (original time series) | 1.5781 | 550.0 | 7.7274
$c$ | 1.3497 | 552.0 | 7.9372
$c,\hat{c}$ | 1.3304 | 348.1 | 5.4421
$c,\hat{c},f,a$ | 1.3872 | 327.4 | 4.1814
VIX |  Benchmark: 87.141 |  |
Features | MSEs | Avg. # hidden units | Avg. train time (s)
$x$ (original time series) | 86.828 | 400.0 | 6.8851
$c$ | 79.581 | 401.0 | 5.6356
$c,\hat{c}$ | 79.265 | 399.5 | 5.6866
$c,\hat{c},f,a$ | 75.450 | 515.7 | 8.3581
TNX |  Benchmark: 37.299 |  |
Features | MSEs | Avg. # hidden units | Avg. train time (s)
$x$ (original time series) | 37.691 | 400.0 | 7.3313
$c$ | 35.749 | 504.8 | 8.8135
$c,\hat{c}$ | 34.044 | 380.7 | 6.2190
$c,\hat{c},f,a$ | 31.954 | 431.4 | 7.1112
Table 6: Out-of-sample prediction error (MSE, $\times 10^{-4}$) and
computational cost of LSTM with different features. The benchmark is the error
of naive prediction $x_{T+1}=x_{T}$. The first rows are the prediction with
the plain time series.

Fig. 10 shows the 1-month rolling MSE of LSTM prediction using different
feature sets. The experiment is implemented on S&P 500 GLD, VIX and TNX from
1/1/2020 - 10/29/2020. The HHT feature sets consistently outperform the
benchmark prediction with plain time series as feature (dashed lines).

Figure 10: 1-month rolling MSE of LSTM extrapolating prediction using HHT
features. The benchmark (dashed line) is the prediction with plain time series
as the feature.

##  5 Conclusion

We have presented the CEEMD method for multiscale analysis of nonstationary
financial time series. The key outputs of this method are the series of
intrinsic mode functions, along with the time-varying instantaneous amplitudes
and instantaneous frequencies. Different combinations of modes allow us to
reconstruct the time series using components of different timescales. Using
Hilbert spectral analysis, we compute the associated instantaneous energy-
frequency spectrum to illustrate the properties of various timescales embedded
in the original time series. Another use of the CEEMD method is to generate a
collection of unique features that can be integrated into machine learning
models, such as regression tree ensemble, support vector machine (SVM), and
long short-term memory (LSTM) neural network. Through a series of examples
with empirical financial data, we show how HHT features can enhance the
machine learning models in terms of forecasting performance.

## References

  * [1] N. E. Huang, Z. Shen, S. R. Long, M. C. Wu, H. H. Shih, Q. Zheng, N.-C. Yen, C. C. Tung, and H. H. Liu, “The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis,” Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, vol. 454, no. 1971, pp. 903–995, 1998.
  * [2] N. E. Huang, Hilbert-Huang Transform and Its Applications, vol. 16.  World Scientific, 2014.
  * [3] T. Y. Hou and Z. Shi, “Adaptive data analysis via sparse time-frequency representation,” Advances in Adaptive Data Analysis, vol. 3, no. 01n02, pp. 1–28, 2011.
  * [4] T. Y. Hou and Z. Shi, “Data-driven time-frequency analysis,” Appl. Comput. Harmon. Anal, vol. 35, no. 2, pp. 284–308, 2013.
  * [5] T. Y. Hou, M. P. Yan, and Z. Wu, “A variant of the EMD method for multi-scale data,” Advances in Adaptive Data Analysis, vol. 1, no. 04, pp. 483–516, 2009.
  * [6] B. Huang and A. Kunoth, “An optimization based empirical mode decomposition scheme,” Journal of Computational and Applied Mathematics, vol. 240, pp. 174–183, 2013.
  * [7] Z. Wu and N. E. Huang, “Ensemble empirical mode decomposition: a noise-assisted data analysis method,” Advances in Adaptive Data Analysis, vol. 1, no. 01, pp. 1–41, 2009.
  * [8] J.-R. Yeh, J.-S. Shieh, and N. E. Huang, “Complementary ensemble empirical mode decomposition: A novel noise enhanced data analysis method,” Advances in Adaptive Data Analysis, vol. 2, no. 02, pp. 135–156, 2010\.
  * [9] N. E. Huang, Z. Shen, and S. R. Long, “A new view of nonlinear water waves: the Hilbert spectrum,” Annual Review of Fluid Mechanics, vol. 31, no. 1, pp. 417–457, 1999.
  * [10] M. Niu, Y. Wang, S. Sun, and Y. Li, “A novel hybrid decomposition-and-ensemble model based on CEEMD and GWO for short-term PM2.5 concentration forecasting,” Atmospheric Environment, vol. 134, pp. 168–180, 2016.
  * [11] L. Tang, W. Dai, L. Yu, and S. Wang, “A novel CEEMD-based EELM ensemble learning paradigm for crude oil price forecasting,” International Journal of Information Technology & Decision Making, vol. 14, no. 01, pp. 141–169, 2015.
  * [12] J. Li, C. Liu, Z. Zeng, and L. Chen, “GPR signal denoising and target extraction with the CEEMD method,” IEEE Geoscience and Remote Sensing Letters, vol. 12, no. 8, pp. 1615–1619, 2015.
  * [13] M. K. Saxena, S. J. Raju, R. Arya, R. B. Pachori, S. Ravindranath, S. Kher, and S. M. Oak, “Empirical mode decomposition-based detection of bend-induced error and its correction in a raman optical fiber distributed temperature sensor,” IEEE Sensors Journal, vol. 16, no. 5, pp. 1243–1252, 2015.
  * [14] S. Nalband, R. Sreekrishna, and A. A. Prince, “Analysis of knee joint vibration signals using ensemble empirical mode decomposition,” Procedia Computer Science, vol. 89, pp. 820–827, 2016.
  * [15] S. Nalband, A. Prince, and A. Agrawal, “Entropy-based feature extraction and classification of vibroarthographic signal using complete ensemble empirical mode decomposition with adaptive noise,” IET Science, Measurement & Technology, vol. 12, no. 3, pp. 350–359, 2017.
  * [16] S. Nalband, C. Valliappan, A. A. Prince, and A. Agrawal, “Time-frequency based feature extraction for the analysis of vibroarthographic signals,” Computers & Electrical Engineering, vol. 69, pp. 720–731, 2018.
  * [17] N. E. Huang, M.-L. Wu, W. Qu, S. R. Long, and S. S. Shen, “Applications of Hilbert–Huang transform to non-stationary financial time series analysis,” Applied Stochastic Models in Business and Industry, vol. 19, no. 3, pp. 245–268, 2003.
  * [18] N. Nava, T. Di Matteo, and T. Aste, “Financial time series forecasting using empirical mode decomposition and support vector regression,” Risks, vol. 6, no. 1, p. 7, 2018.
  * [19] J. Wang and J. Wang, “Forecasting stochastic neural network based on financial empirical mode decomposition,” Neural Networks, vol. 90, pp. 8–20, 2017\.
  * [20] N. Nava, T. Di Matteo, and T. Aste, “Dynamic correlations at different time-scales with empirical mode decomposition,” Physica A: Statistical Mechanics and its Applications, vol. 502, pp. 534–544, 2018.
  * [21] V. Kurbatsky, D. Sidorov, V. Spiryaev, and N. Tomin, “Using the Hilbert-Huang transform for ANN prediction of nonstationary processes in complex power systems,” in 8 th World Energy System Conference “WESC, pp. 106–110, 2010.
  * [22] V. G. Kurbatskii, D. N. Sidorov, V. A. Spiryaev, and N. V. Tomin, “On the neural network approach for forecasting of nonstationary time series on the basis of the Hilbert-Huang transform,” Automation and Remote Control, vol. 72, no. 7, pp. 1405–1414, 2011.
  * [23] V. G. Kurbatsky, D. N. Sidorov, V. A. Spiryaev, and N. V. Tomin, “Forecasting nonstationary time series based on Hilbert-Huang transform and machine learning,” Automation and Remote Control, vol. 75, no. 5, pp. 922–934, 2014.
  * [24] J. Zhang, T. Leung, and A. Aravkin, “Sparse mean-reverting portfolios via penalized likelihood optimization,” Automatica, vol. 111, p. 108651, 2020\.
  * [25] G. Rilling, P. Flandrin, P. Goncalves, et al., “On empirical mode decomposition and its algorithms,” in IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing, vol. 3, pp. 8–11, 2003.
  * [26] E. C. Titchmarsh, Introduction to the theory of Fourier integrals.  Clarendon Press, 1948.
  * [27] E. Bedrosian, “A product theorem for Hilbert transforms,” Proceedings of the IEEE, vol. 51, no. 5, pp. 868–869, 1963.
  * [28] A. Nuttall and E. Bedrosian, “On the quadrature approximation to the Hilbert transform of modulated signals,” Proceedings of the IEEE, vol. 54, no. 10, pp. 1458–1459, 1966.
  * [29] W. S. Cleveland, “Robust locally weighted regression and smoothing scatterplots,” Journal of the American Statistical Association, vol. 74, no. 368, pp. 829–836, 1979.
  * [30] L. Breiman, Classification and Regression Trees.  Routledge, 2017.
  * [31] T. Hastie, R. Tibshirani, and J. Friedman, The elements of statistical learning: data mining, inference, and prediction.  Springer Science & Business Media, 2009.
  * [32] H. Drucker, C. J. Burges, L. Kaufman, A. J. Smola, and V. Vapnik, “Support vector regression machines,” in Advances in Neural Information Processing Systems, pp. 155–161, 1997.
  * [33] V. Vapnik, The nature of statistical learning theory.  Springer Science & Business Media, 2013.
  * [34] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” Neural Computation, vol. 9, no. 8, pp. 1735–1780, 1997.
  * [35] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” ArXiv Preprint arXiv:1412.6980, 2014.