11institutetext: University of Illinois at Urbana-Champaign, Urbana, IL
61801-2302, USA

# Attacks on Visualization-Based Malware Detection: Balancing Effectiveness
and Executability

Hadjer Benkraouda 11 0000-0001-5511-3182 Jingyu Qian 11 0000-0002-3953-5382
Hung Quoc Tran Benkraouda, Qian, and Tran share co-first authorship.11
0000-0001-7767-3180 Berkay Kaplan 11 0000-0002-4365-7606

###### Abstract

With the rapid development of machine learning for image classification,
researchers have found new applications of visualization techniques in malware
detection. By converting binary code into images, researchers have shown
satisfactory results in applying machine learning to extract features that are
difficult to discover manually. Such visualization-based malware detection
methods can capture malware patterns from many different malware families and
improve malware detection speed. On the other hand, recent research has also
shown adversarial attacks against such visualization-based malware detection.
Attackers can generate adversarial examples by perturbing the malware binary
in non-reachable regions, such as padding at the end of the binary.
Alternatively, attackers can perturb the malware image embedding and then
verify the executability of the malware post-transformation. One major
limitation of the first attack scenario is that a simple pre-processing step
can remove the perturbations before classification. For the second attack
scenario, it is hard to maintain the original malware’s executability and
functionality. In this work, we provide literature review on existing malware
visualization techniques and attacks against them. We summarize the limitation
of the previous work, and design a new adversarial example attack against
visualization-based malware detection that can evade pre-processing filtering
and maintain the original malware functionality. We test our attack on a
public malware dataset and achieve a 98% success rate.

###### Keywords:

Malware visualization Adversarial machine learning Binary rewriting.

##  1 Introduction

With the proliferation of connectivity and smart devices in all aspects of
human life, these devices have become increasingly targeted by malicious
actors. One of the most common forms of attack is through malicious software
(i.e., malware). According to AVTEST, one of the leading independent research
institute for IT security, to date, there have been more than 100 million new
malware applications in 2020 alone [20]. Inspired by the success of machine
learning in other fields, researchers have proposed using machine learning for
many security applications. With the rapid growth and evolving nature of new
malware applications, machine learning-based solutions are a natural fit for
malware detection and classification due to their robustness.

Several papers have designed malware detection systems using machine learning
(e.g., [35, 24]). The proposed solutions tackle malware detection from
different perspectives. The main differences are in the representation used
and the subsequent machine learning model selected for effective
classification. These representations include raw bytes, embeddings,
representative features, and binary visualization. Visualization methods, in
particular, have shown high accuracy in detecting malware compared to
conventional methods. These data reduction and visualization techniques for
detection have shown improvements in both speed and memory efficiency [29,
17]. Additionally, visualization-based techniques have achieved higher
detection accuracy, mainly attributed to the applicability of deep learning
techniques in detecting malware patterns [1]. We, therefore, focus our work on
visualization-based malware detection models.

Nevertheless, machine learning models are susceptible to adversarial example
attacks, which add imperceptible non-random perturbations to test samples,
causing machine learning models to misclassify them. Successful adversarial
examples have been seen to fool systems into misclassifying people [37], cause
systems to recognize street stop signs as speed limit signs [10], or cause
voice-controllable systems to misinterpret commands or perform arbitrary
commands [41].

Recent work has shown that machine learning-based techniques for malware
detection are also susceptible to adversarial examples. In these systems, the
attacks alter the malware binaries to cause the target model to classify the
malware sample as benign or vise versa. However, adversarial examples in this
domain are more challenging to produce. In addition to the constraint of
imperceptibility and minimal changes that conventional adversarial examples
must comply with, adversarial examples for malware detection must maintain the
original malware functionality. This means that the attacker cannot change the
bytes arbitrarily. Instead, the attacker has to understand the functionality
of the malware and perform careful changes.

There have been previous attempts to create adversarial examples against
visualization-based malware classifiers [22, 28]. These attacks produce
adversarial examples either by using conventional image-based techniques such
as the Fast Gradient Sign Method [15] or Carlini and Wagner method [4] or by
injecting byte values to unreachable regions within the binary. These attacks
are simplistic and can be detected and removed easily with minimal
countermeasures [25].

In this paper, we propose a new adversarial example attack that combines
binary rewriting and adversarial attacks in image classification. We target a
convolutional neural network (i.e., CNN) model for malware detection. Because
there is no open-sourced code for visualization-based malware detection, our
first phase of the project includes constructing the malware detection model
(Figure 1 left). We apply a similar CNN structure as previous work for
visualization-based malware detection and achieves an overall accuracy of 99%.
In the second phase of the project (Figure 1 right), we design an adversarial
example attack against this malware detection model. Our attack performs
additive changes to the original malware and ensures that the added
instructions are semantic NOPs, i.e., they do not change values of any
register or manipulate the program state. Our attack achieves a 98% success
rate on a public malware dataset. The success of the proposed attack reveals
that it is necessary for visualization-based malware detection to perform more
advanced and robust protection against adversarial examples other than simply
filtering the padding or the non-reachable header section.

Figure 1: The two project phases: constructing the malware detection model and
designing the adversarial example attack.

The rest of the paper is organized as follows. In Section 2, we introduce
background and related work on visualization-based malware detection and
adversarial machine learning. In Section 3, we provide our detailed design of
the attack against visualization-based malware detection and illustrate how we
solve the challenges of creating successful adversarial examples while
maintaining the original malware functionality. In Section 4, we discuss our
experiment setup and measure our attack success rate. In Section 5, we discuss
limitations of our attack and potential future work.

##  2 Background and Related Work

In this section, we introduce the background and related work on
visualization-based malware detection. We then discuss some traditional
methods to camouflage malware, adversarial machine learning and how attacks
against image classification and malware detection work. Finally, we include a
table for SoK of the papers that we mentioned, and point out their
limitations.

###  2.1 Malware Visualization

With the development of image processing technology, visualization-based
techniques are also proposed for malware detection and analysis. These
techniques can be applied directly to the binary without complicated
disassembly and execution process. Researchers have proposed approaches to
visualize malware as gray-scale or RGB-colored images. From these images,
machine learning techniques, such as CNN, can classify whether the tested
software is benign or malicious.

Figure 2 illustrates a typical approach to visualize the malware as an image.
The malware binary is grouped by 8-bit vectors. Each vector represents a value
from 0 to 255, which can be mapped to a gray-scale pixel value. The shape of
the final image depends on the width of the image, which is usually a tunable
parameter, and the size of the malware binary in bytes. This methodology can
be adapted to visualize the malware as an RGB-colored image, which considers
different feature types and represents them in different color channels [12].

Figure 2: Typical approach to visualize the binary in gray-scale [31]

Here we introduce a few projects from the literature that used several
different malware visualization techniques. Han et al. [18] proposed a malware
visualization method that converts the binary to a gray-scale bitmap image and
then generates the entropy graph for the entire malware binary. They then used
a histogram similarity measuring method to group malware within the same
malware family. Nataraj et al. [31] also visualized the malware binary as a
gray-scale image but extracted texture features to characterize and analyze
the malware. They used GIST to capture texture features, and apply the
K-nearest neighbors algorithm with Euclidean distance to classify the malware
into different malware families. Xiaofang et al. [40] mapped malware binaries
to gray-scale images, extracted a 64-dimension feature vector from the image,
and performed fingerprint matching to identify similar malware.

Unlike the other work, which converts the malware binary to a gray-scale
image, Fu et al. [12] took a different approach to visualize the malware
(Figure 3). Their approach extracted both local and global features to
generate an RGB-colored image for the malware. Specifically, they extracted
three types of features, including section entropy, raw byte values, and
relative size of each section to the whole file. For raw byte values, they use
the same approach to visualize malware as gray-scale images (Figure 2). Each
types of features occupies a single color channel (i.e., either red, green, or
blue). For the final classification process, they compared different machine
learning techniques, including random forest, K-nearest neighbors, and support
vector machine.

Han et al. [19] proposed a novel method to classify malware and malware
families. They extracted opcode instruction sequences from the malware binary
as features and generated the RGB-colored image matrix from these features.
The image matrices are compared with each other using selective area matching
[19] to group malware into malware families.

Another work interestingly extends the field by visualizing the behavior of
the malware instead of the malware binary itself, and suggests that any
feature of a malware can be visualized to determine its classification [36].
This work further indicates that the possibility of malware visualization is
limitless as a program has a lot of useful features ranging from its behavior
to its metadata. But, for this work, it specifically focuses on the malware
behavior by running an API call monitoring utility through a Virtual Machine
(i.e., VM) to examine the APIs used while the program is executed in user mode
[36]. While there are several other techniques to capture malware behavior,
such as capturing the network activity of the malware, or the changes in the
operating system’s resources, API monitoring has been chosen in this study due
to its conciseness and the shortcomings of other techniques that have been
discussed in detail [36]. Afterwards, the calls are mapped to hot colors,
specifically such as red, or orange, for classification. Finally, the
classification is mostly done through a similarity ratio to the tested
software against the known malware’s color mapping [36].

Figure 3: Visualize the malware as an RGB image, considering both local and
global features [12]

###  2.2 Traditional Malware Camouflage

Although there are several methods of malware detection based on
visualization, attackers can still employ various methods such as encryption
and obfuscation to hide their malware in the targeted software’s code and
counter static malware detection methods [7, 27, 34].

Malware encryption intends to encrypt the malware body to hide its intentions
and avoid static analysis detection so that a direct signature matching
defense cannot detect the malware [7, 27, 21]. It relies on a decryption loop
(a.k.a., decryptor) to decrypt the malicious payload and execute the malware.
Therefore, if the defense can find out the decryption loops in the malware, it
can decrypt the malicious code first and then perform a simple signature
matching to detect the malware. In addition, the original defense can be
easily augmented with a signature checking component to identify suspicious
decryption loops within the malware. Visualization-based defense can be even
better at detecting malware encryption because it can extract locality
features specific to suspicious decryption loops of the malware. Even a more
complicated malware encryption technique that picks different decryptors for
different malware (i.e., oligomorphism) only prolongs the detection time [34].

Figure 4: Polymorphism virus structure [34].

Polymorphism is a more sophisticated method to hide the malicious payload
based on malware encryption (Figure 4). It uses several types of
transformations on the decryptor, such as changing the order of instructions
with additional jump instructions to maintain the original semantics and
permuting the register allocation to deceive anti-virus software [7]. It also
typically injects junk or dead codes to further mutate the decryptor so that
it is hard to recognize it [34]. However, after enough emulation and a simple
string matching algorithm, the underlying encrypted sections of the malware
can still be revealed [34].

Figure 5: Metamorphism structure: a metamorphic engineer is responsible for
changing the malware instructions to equivalent ones probabilitically [34].

A major drawback of either malware encryption or polymorphism is that it
relies on an explicit decryptor to decrypt the malicious payload and execute
the malware. This leaves a significant mark on the malware that can be
relatively easily detected. On the other hand, metamorphism is a more advanced
technique to camouflage the malware without using any encrypted parts. Malware
metamorphism is a technique to mutate the malware binary using different
obfuscations by a metamorphic engine (Figure 5). In this way, the attacker
changes the syntax of the original malware but keeps the original malware
behavior. In particular, metamorphism allows the malware to change its opcode
with each execution of the infected program. Alam et al. [2] group some
typical obfuscations used in metamorphism into three categories. The first
category is the opcode level obfuscation, which includes instruction
reordering, dead code insertion, and register renaming [2]. The second
category is control flow level obfuscation, which includes changing the order
of instructions, and applying branch functions, opaque predicates, jump
tables, and exception tables [2]. The last category is obfuscation by self-
modifying code, which intends to change instructions during runtime in order
to hide malicious payload to avoid reverse engineering and detection by anti-
malware software [2].

It has been shown that malware metamorphism can easily defeat the signature-
based malware detection [6] because signature-based detection is unable to
capture the changes of the malware due to dynamic code obfuscation. However,
metamorphic malware is usually initiated from known malware, and with the
initial knowledge of existing malware, it is still possible to detect malware
metamorphism. Zhang et al. [42] proposed a defense to characterize the
semantics of the program and perform code pattern matching, based on static
analysis of control and data flow of call traces. Alam et al. [2] proposed a
metamorphic malware analysis framework that builds the behavioral signatures
to detect metamorphic malware in real-time. Chouchane et al. [5] proposed
another method to detect the existence of a metamorphic engineer by checking
the likelihood of a piece of code generated by some known metamorphic engines.
In addition, the metamorphic malware typically has a significant proportion of
binary related to the metamorphic engine, which can be recognized by a more
advanced detector, such as a visualization-based malware detector.

Go et al. discussed the importance of developing new approaches against
polymorphism and metamorphism specifically [13]. Their method converts the
binary to a grey-scale image and uses the ResNeXt CNN model to build
resiliency against such malware camouflage techniques [13]. Their paper does
not explicitly discuss their method’s effectiveness against obfuscation
attacks [13]. The researchers used the Malimg dataset but did not mention that
their dataset contained examples of obfuscation [13]. Islam et al. focused
more on obfuscation detection by attempting to integrate static analysis with
dynamic [21]. The paper acknowledges the ease of bypassing static analysis
with obfuscation but proposes integrating dynamic analysis using information
vectors derived from FLF, PSI and API calls [21]. Since obfuscating the
features of the code would result in outlier vectors, their approach can
detect such attacks [21].

###  2.3 Adversarial Machine Learning

In this section, we introduce adversarial machine learning. In general,
adversarial machine learning aims to fool the machine learning model by
carefully generating adversarial examples through evasion attacks or polluting
the training phase through poisoning attacks. We focus our discussion on
evasion attacks leveraged against image classification and malware detection
due to their many close similarities to visualization-based malware detection.

####  2.3.1 Attacking Image Classification.

Previously, image classification has been used in many applications not
related to binary classification. In these fields, multiple attacks have been
produced to cause errors in detection. Goodfellow et al. [15] illustrated the
fast gradient sign method (i.e., FGSM) to generate adversarial examples.
Figure 6 shows how FGSM is applied to cause the image classifier to mistakenly
classify a panda to a gibbon by adding carefully crafted perturbation. In
Carlini et al.[4], it was shown that the addition of random noise to images
can significantly reduce the accuracy of classifiers while being imperceptible
to the human eye. These attacks can be mitigated to some degree by denoising
techniques, as proposed in Liao et al.[26]. Nevertheless, such mitigation
efforts do not fully reverse the perturbations and may introduce more noise
accidentally. Furthermore, as shown in Eykholt et al.[11], classification can
also be interrupted by altering small sections of an image, where the image
would still be readable by a human eye.

Figure 6: Adversarial example generated using FGSM [15].

A notable difference in image classification for binary visualization is that
the validator for images lies in code execution instead of human recognition.
As a result, perturbations intended for adversarial attacks on binary
visualization must continue to function as an unchanged code sample. A code
sample that maintains functionality after modification can be said to maintain
executability, as the code will execute as intended.

####  2.3.2 Attacking Malware Detection.

Because malware detection from raw bytes relies heavily on the performance of
the machine learning model for classifying the image embedding of the malware
(e.g., RGB-colored image), it is also vulnerable to similar attacks against
image classification. However, there are more difficulties in generating valid
adversarial examples in the malware detection domain than in the image
classification domain. The main challenge is to solve the inverse feature-
mapping problem [33]. Pierazzi et al. [33] proposed a novel formalization of
problem-space attacks and a novel problem-space attack in the Android malware
domain. They used conditional statements that are never executed during
runtime to wrap the malicious code payload. In particular, they used opaque
predicates to ensure the obfuscated conditions always resolve to false but
look legitimate. They showed that their attack against Android is successful
against the Drebin classifier [3] and several SVM-based malware detectors.

Liu et al. proposed introducing perturbations to the visualized binary to
lower the success rate of ML-based malware detectors [28]. The paper
introduced a method that leverages gradient descent and L-norm optimization
methods. However, as it changes the image in potentially unexpected ways, it
cannot guarantee the executability of the perturbed malware.

Khormali et al. [22] showed simple adversarial examples to bypass
visualization-based malware detection. Their attack only attempts rudimentary
methods to generate adversarial examples, namely padding and injection [22].
It preserves the executability and functionality of the original malware, but
it is easy to detect and not scalable. Kolosnjaji et al. [23] proposed a
gradient-based attack to evade ML-based malware detection. Their attack also
injects padding bytes to the original malware but does not consider preserving
the functionality and does not target visualization-based malware detection.
Demetrio et al. [9] proposed a general framework to perform white-box and
black-box adversarial attacks on learning-based malware detection by injecting
the malicious payload to the DOS header. However, the authors also claimed
that their header attacks could be easily patched through filtering process
before classification. Grosse et al. [16] demonstrated the attack against a
deep neural network approach for Android malware detection. They crafted
adversarial examples by iteratively adding small gradient-guided perturbation
to the malware on application level instead of directly perturbing the binary.
They restrict the perturbation to a discrete set of operations that do not
interfere with the malware functionality. In addition, they discussed some
remedies against their attacks, including feature reduction, distillation, and
adversarial training (i.e., re-training the model with the addition of
adversarial examples). Their work focused more on attacking application-level
malware detection instead of visualization-based malware detection. Sharif et
al. [38] proposed an optimization-guided attack to mislead deep neural
networks for malware detection. Their attack is more invasive in that it
changes the reachable code section of the malware. Nevertheless, their attack
considers a limited set of transformation types of malware functions and does
not target visualization-based malware detection [38]. Overall, there is no
work proposing robust attacks against the visualization-based malware
detection that preserve both executability and functionality of the original
malware and are hard to detect through pre-processing. We seek to fill this
room in the research space.

###  2.4 SoK of Existing Literatures

Finally, we provide Table 1 as our SoK of existing malware visualization
techniques and attacks against them. For each reviewed work, we summarize its
methodology and list its limitations.

Table 1: SoK of reviewed papers: malware visualization and adversarial attacks against malware visualization software. Category |  Paper |  Methodology |  Limitations
---|---|---|---
Malware visualization on binary: Gray-Scale |  Han et al. [18] |  Converts binary to the bitmap image and generates the entropy graph from visualized malware |  Hard to classify packed malware binaries
Nataraj et al. [31] |  Extracts image texture features from visualized malware |  Relying on global texture features can be beat by attackers
Xiaofang et al. [40] |  Extracts a 64-dimensional feature vector and performs fingerprint matching to identify similar malware |  Relying on global image features can be beat by attackers
Malware visualization on binary: RGB-Colored |  Fu et al. [12] |  Combines entropy, relative section size, and raw bytes to generate an RGB-colored image |  Limited to PE format
Han et al. [19] |  Extracts opcode instruction sequences |  Classification is not yet automated
Malware visualization on behavioral features |  Shaid et al. [36] |  API call monitoring |  Does not consider network behavior and does not work directly on malware binary
Attacking malware detection |  Pierazzi et al. [33] |  General problem-space attack for inverse feature-mapping |  Attacks focus on Android malware and are not against neural network-based detector
Liu et al. [28] |  Gradient descent to perturb binary |  Not guarantees executability
Kormali et al. [22] |  Padding and injection |  Easy to detect and not scalable
Kolosnjaji et al. [23] |  Padding and injection |  Does not preserve functionality and does not target visualization-based malware detection
Demetrio et al. [9] |  Injects the malicious payload to DOS header |  Easy to patch through filtering process before classification
Grosse et al. [16] |  Iteratively adds small gradient-guided perturbations |  Only targets Android malware and does not attack on binary level
Sharif et al. [38] |  Manipulates code section guided by an optimization function |  Considers only a limited set of manipulation types

##  3 Robust Adversarial Example Attack against Visualization-Based Malware
Detection

In this section, we focus on the workflow of generating adversarial examples,
and we leave the construction of the malware detector to Section 4.1. We have
two main goals for our adversarial example attack. Firstly, we aim to find an
adversarial example generated from a single malware such that the malware
detector will misclassify it as benign software. Secondly, an adversarial
example generated in this way must maintain the functionality of the original
malware. An overview of the full workflow of our adversarial example
generation algorithm is shown in Figure 7. At a high level, the attack starts
with using a mask generator to add empty spaces to instruction boundaries
where perturbations are allowed. Then, the adversarial example generator
(i.e., AE generator) will generate the optimal adversarial example in the
image space. To ensure that the original malware functionality is not changed,
we use a NOP generator to produce a semantic NOP list and update the optimal
adversarial example to the closest matching viable one that preserves malware
functionality. If this processed adversarial example is still misclassified as
benign, then our attack succeeded. Otherwise, we relaunch the AE generator,
starting from the failed adversarial example, creating a new optimal AE, and
starting a new iteration. We iterate until we produce a successful adversarial
example or we reach a pre-set threshold of iterations. In the following sub-
sections, we discuss the details of each component in the workflow.

Figure 7: The overview of the workflow of the adversarial example generation
algorithm.

###  3.1 Mask Generator

The first step of our attack workflow aims at controlling and locating where
the perturbations can be added. This step is provided to ensure both
executability and robustness to simple pre-processing defenses while
maintaining the original semantic operation. The central intuition of this
step is to allow additional instructions to be embedded within the code
section so that they are not easily distinguishable from the rest of the code
section and, in the meantime, ensure that these added instructions do not
introduce any changes to the original malware instructions. These additional
instructions will serve as empty spaces to allow perturbing the malware in the
image representation.

Figure 8: Mask generator adds empty space to allow perturbation

To achieve this, we create a mask generator (Figure 8). The algorithm of the
mask generator is as follows. First, we extract the code section of the
malware sample and identify the instruction boundaries. Next, we need to
decide the size of each perturbation block, its location, and frequency. The
attacker can set these parameters to satisfy certain size limitations of the
perturbation relative to the original malware size to make it harder to
detect. On the other hand, the attacker can also increase the frequency of the
perturbation blocks to make the attack easier. The perturbations can be
initialized to random inputs or naive NOPs. After that, these perturbations
are inserted into the expected instruction boundaries. In this way, we make
sure that the original malware instructions are not changed at all because we
never add extra bytes to the middle of the binary of any instruction. With
this malware augmented with the perturbation sequences initialized to naive
NOPs, we use a binary-to-image converter to represent the malware in the image
space. The binary-to-image converter treats the consecutive 8 bits of the
binary as a single pixel to build a PNG file. Figure 8 illustrates how our
mask generator adds empty space to allow perturbation. In this example, we add
a single NOP every two instructions. The image representation of the binary is
expanded by two pixels (i.e., the pixel marked with ‘*’) due to the two added
nops.

Besides the malware image, the mask generator produces a mask in the form of
an array of the same dimension as the augmented malware image. The mask flags
the locations where perturbations are allowed with ones, while the rest of the
array is filled with zeros. We name this mask the perturbation mask.

###  3.2 AE Generator

Once the perturbation mask and the augmented malware image are generated, we
launch a modified version of the CW attack [4] to generate the optimal
adversarial example (i.e., optimal AE) in the image space, which is
misclassified as benign by the malware detector. The only difference is the
application of the perturbation mask to further restrict the positions of
perturbations. The objective function is given as

$\min||M\delta||_{2}+C\cdot f(x+M\delta)\quad s.t.\quad x+M\delta\in[-1,1].$
(1)

Here $M$ is the perturbation mask, and $x$ is the augmented malware image
produced by the mask generator. The optimal AE is unlikely to correspond to
the perturbed malware that maintains the original malware functionality
because the CW attack only intends to attack the malware detector, which is
modeled as an image classifier. Our objective function does not place
additional restrictions to ensure that the perturbations will be either
executable or semantic NOPs after being reversed back to the binary, which is
required to preserve the malware functionality. One possible approach is to
apply the similar idea of ensuring printability for adversarial examples in a
real-world attack against image classification. However, the number of
semantic NOPs is much larger than the number of printable colors, which will
significantly increase the computation overhead of the objective function. On
the other hand, the optimal AE is a good starting point to guide us to
generate a viable adversarial example that keeps the malware functionality. We
also avoid complicating the objective function to reduce the runtime overhead
to generate the optimal adversarial examples.

###  3.3 NOP Generator

To generate the viable adversarial example from the optimal AE, we need to
replace the perturbations introduced by the CW attack with binaries that keep
the malware functionality (i.e., semantic NOPs). The replacement algorithm
will be discussed in Section 3.4. This section shows how we create the
semantic NOPs, which are semantically equivalent to a naive NOP, which does
not change the values stored in the registers and the function state, such as
the stack. We do not use the same strategy in [33] that used unreachable
conditional statements because that requires careful crafting of obfuscated
conditions to avoid being removed by the compiler optimization or the
filtering process, which makes the perturbation generation process more
complicated.

Figure 9: Semantic NOP seeds used to construct longer semantic NOPs.

We start by creating short semantic NOPs. We call these short semantic NOPs
semantic NOP seeds. Semantic NOP seeds fall into four categories (Figure 9):

  1. 1.

Arithmetic sequences: Some neutral sequence of arithmetic operators which can
be performed on any register. An example is adding zero to the register.

  2. 2.

Movement sequences: Instructions to move the register value back and forth or
jump to the defined locations to skip the regions that are not expected to
execute. Examples are moving the register value to itself or jumping to the
immediate next instruction.

  3. 3.

Logical sequences: Some neutral sequence of logical operators which can be
performed on any register. An example is ANDing 1 to the register.

  4. 4.

Other miscellaneous sequences: A sequence of simple NOPs or a sequence to
change and recover the flags.

Because the perturbation space is not pre-determined, we do not generate
semantic NOPs for any arbitrary size. Instead, we build the NOP generator to
combine semantic NOP seeds to make longer semantic NOPs that can fit larger
perturbation spaces. For instance, if the NOP generator is asked to compute
3-byte semantic NOPs and the naive NOP (i.e., byte 90 in heximal) is one of
the semantic NOP seeds, then it can combine three NOPs. Given the expected
size of each perturbation block, the NOP generator produces a list of semantic
NOPs of the given size.

It is necessary to keep the byte length of semantic NOP seeds as small as
possible so as to improve the ability for the AE optimizer to produce a viable
adversarial example that maintains the malware functionality. However, too
much restriction on the byte length also limits our ability to generate enough
semantic NOP seeds. In our design, we pick the minimum size of a semantic NOP
seed as one byte (i.e., the byte length for a naive nop) and the maximum size
of a semantic NOP seed as eight bytes. We provide our example semantic NOP
seeds and the corresponding byte size in Table 2. We admit that this is a
simple design choice and is far from the optimal selection. We also do not
comprehensively list all the possible operations given the operation type,
which can make the set of semantic NOP seed too large. However, our evaluation
results reveal that this selection is already enough for us to generate
adversarial examples for malware that maintains functionality with high
success rate.

Table 2: Example of semantic NOP seeds and their byte length. Operation Type | Example in Hex | Byte Length
---|---|---
Nop | 90 | 1
Move register to itself | 89c0 | 2
Jump to next instruction | 7700 | 2
Push and pop | 5058 | 2
Not and not | f7d0f7d0 | 4
Add/subtract 0 | 9c83c0009d | 5
Logical AND with 1 | 9c83e0ff9d | 5
Logical OR with 0 | 9c83c8009d | 5
Logical XOR with 0 | 9c83f0009d | 5
Negate and negate | 9cf7d8f7d89d | 6
Increment and decrement | 9cffc0ffc89d | 6
Add x and subtract x | 9c83c00183e8019d | 8

###  3.4 AE Optimizer

In this step, we build a module, the AE optimizer, which produces a viable
adversarial example that maintains the original malware functionality. The AE
optimizer takes in the perturbation mask, the optimal AE generated from the CW
attack, and the list of semantic NOPs produced by the NOP generator. Next, the
AE optimizer locates the allowed positions for perturbations in the optimal AE
using the perturbation mask. Subsequently, the Euclidean distance between the
instruction sequences in the allowed perturbation spaces and the semantic NOPs
is calculated using the following equation:

$d\left(p,q\right)=\sqrt{\sum_{i=1}^{n}\left(q_{i}-p_{i}\right)^{2}}$ (2)

Here $p$ and $q$ are the generated instruction sequence and the semantic NOPs.
This process identifies the semantic NOPs closest to each of the sequences in
the allowed perturbation space of the optimal AE. In the current
implementation, this process is done sequentially for each perturbation block;
however, it can be easily parallelized to improve the runtime performance.
After that, the semantic NOPs with the minimum distance are used to replace
the perturbation blocks in the optimal AE. The new adversarial example is
called the optimal viable AE.

Finally, we pass the optimal viable AE to our malware detector for
classification. If it is classified as benign, we stop the process because we
have already produced a successful adversarial example. If it is correctly
classified as malware, it will be used as the starting point for another
iteration of the CW attack, and the process is repeated. We expect that
starting from a failed optimal viable AE can direct us better to the
successful optimal viable AE. However, it is possible that the AE can be stuck
in the local optimum. Another more general approach is to start over the whole
process again from the visualized malware augmented with randomly initialized
semantic NOPs.

##  4 Evaluation

###  4.1 Experiment Setup

####  4.1.1 Malware Detection Model.

We use a CNN-based malware detector as our attack target. Because there is no
open-source code for the visualization-based malware detector, we build our
own CNN by following a similar structure from previous work. Specifically, we
re-built the structure described in [22, 23]. Our CNN is composed of an
adaptive average pooling layer to handle inputs of different dimensions, two
consecutive convolutional layers with max-pooling, and three fully connected
layers. We use ReLU as the activation function (Figure 10).

Figure 10: Malware detection neural network model structure.

####  4.1.2 Dataset.

We build our dataset by combining one public malware dataset and one public
benign software dataset. The malware dataset is from UCSB’s public malimg
dataset [30], which has a collection of around 9500 malware formatted as PNGs.
The benign dataset is a collection of around 150 benign software from the
Architecture Object Code Dataset [8]. The software was chosen from the AOCD
with various functionalities to ensure that the classifier did not learn
latent features representative of the type of benign code.

####  4.1.3 Model Training.

A relatively small subset of the malimg dataset is used to acquire a $50\%$
split between malware and benign images in the training data. This subset was
retrieved randomly from the malimg dataset, and examples were chosen without
considering their corresponding malware family. In doing this, we prevent the
classifier from learning any latent features from any particular malware
family and improve its accuracy for any malware class in the validation set.
Validating this model’s accuracy, we are able to confirm the classification
accuracy of state-of-the-art models, with a $100\%$ accuracy on identifying
benign software and a $99\%$ accuracy on correctly identifying malware
samples.

####  4.1.4 Attacking the Malware Detector.

We evaluate our attack using the malware from the malimg dataset [30]. For
each malware, we augment it with empty spaces per instruction initialized with
eight naive NOPs. If the AE generator fails to produce the optimal AE, we
consider it an attack failure. If the AE generator can produce the optimal AE,
we run the AE optimizer to replace perturbations with the closest semantic
NOPs. We set the iteration threshold to be ten. If the optimal viable AE is
generated within the ten iterations of the CW attack and AE optimization and
can be successfully misclassified as benign, we view it as a successful
attack.

###  4.2 Results

We first evaluate our attack on the malware family, “Dialplatform.B”, which
contains 174 malware samples. Malware in this family has moderately large code
section, which gives us more flexibility to generate adversarial examples in a
short time. The corresponding image of each malware is either $216\times 64$
or $432\times 64$ in pixels. The average number of instructions in the code
section for each malware is 622. All the malware is classified as malware by
our malware detector.

We successfully generate adversarial examples for 172 malware out of the total
174 malware (98.9%). For the two malware that we fail to generate an
adversarial example for, our algorithm workflow oscillates around a local
optimum, so our iteration of CW attack and AE optimization (steps 2, 4, 5 in
Figure 1). There are two potential methods to avoid the local minimum
oscillation issue. First, we can initialize the empty spaces between
instructions with random semantic NOPs. Randomized initialization has already
been applied frequently to solve similar local optimum oscillation problems
before. Second, in the AE optimizer step (Section 3.4), we could pick a sub-
optimal viable adversarial example if we detect local optimum oscillation. In
this way, our algorithm can break the local optimum while also searching for a
valid adversarial example.

For all of the adversarial examples, the functionality of the original malware
is preserved by construction since we only instrument the original binary with
semantic NOPs at instruction boundaries. All of the adversarial examples can
be generated within five iterations of the CW attack and AE optimization. The
running time to generate the adversarial example is given in Figure 11. On
average, it takes 77.8 seconds to generate the adversarial example for the
malware with around 600 instructions, and the time overhead is in AE
optimization step. The expansion rate of the original malware due to
augmentation is shown in Figure 12. On average, the perturbation size is
35.82% of the size of the original malware as seen from Figure 12. We argue
that though the expansion rate is high, the added perturbation is still hard
to filter out due to the flexibility of our semantic NOPs. The attackers can
even build a much larger semantic NOP set than our current version. They can
also consider using dead code branching. In this way, unless the users know
the semantic NOPs and the initial size of the malware, it is hard to filter
out our attacks.

We achieve similar results for other malware families with similar sizes of
the code section. On the other hand, our algorithm does not work well for
malware with very small code sections. We argue that the spaces to add
semantic NOPs are too few and too small to allow enough perturbation to cause
misclassification due to the small code section. We expect that enlarging the
mask size can improve the attack success rate, but this might defeat the
attacker’s purpose for distributing small and easy-to-transmit malware.
Another potential way is to combine perturbations in code sections with those
in other sections, such as data sections. However, we argue that keeping
executability can be hard and naive padding is relatively easy to filter out.
We leave attacks for malware with a small code section as future work.

Figure 11: The runtime evaluation for the adversarial example attack. Figure
12: The size of the added perturbation of the adversarial example attack.

To further test the end-to-end methodology for generating adversarial examples
for malware with a larger number of instructions, we also run the same
experiment with 11 malware from the “Lolyda.AA3” class. Each malware contains
at least 200,000 instructions. We achieve an attack success rate of 81.8%
(i.e., 9 out of 11). The two failures are also due to local optimum
oscillation issues that we face when generating adversarial examples for
malware in “Dialplatform.B” family. We expect random initialization and sub-
optimal AE optimizer can solve the problem and leave it to future work. On the
other hand, the major problem for generating adversarial examples for large
malware is the time overhead. In our experiments, the time overhead reaches
about six hours to generate a single adversarial example. As a proof-of-
concept attack and assuming the attacker can tolerate this time overhead, our
attack can still work properly, but the practicality of the attack in real
life is still under question considering the size overhead added to the
original malware as well. We expect parallelism can improve the running time.
Specifically, our current AE optimizer finds viable semantic NOP sequences for
each empty space defined by the mask sequentially. However, each empty space
is independent of the other if we perturb it with self-contained NOP
sequences. Therefore, parallelization can be easily applied in the AE
optimization step. We leave its implementation and evaluation as future work.

Our high attack success rate is a conservative representation of how
vulnerable visualization-based malware detection can be to the adversarial
examples. In our experiment, we set each perturbation block to be small (i.e.,
8 bytes). A more powerful attacker who would like to risk being detected more
can further increase the size of each perturbation block so that the attack
can be easier to launch.

##  5 Discussion

###  5.1 Limitations

While the proposed attack can successfully find adversarial examples with a
high success rate, due to the nature of the optimization algorithm, the time
necessary to find the optimal viable AE increases drastically with the size of
the original malware. In our evaluation, we performed adversarial example
attacks on the “DialPlatform.B” malware family, where each malware image is of
dimension $216\times 64$ or $432\times 64$ with an average of 622
instructions. Since the images are of reasonably small dimensions, the
potential perturbation spaces are fairly limited. As a result, an adversarial
example can be generated in less than two minutes. However, for larger
malware, as in the “Lolyda.AA3” family, producing a single optimal adversarial
example can take a few hours. As each perturbation block can be independently
replaced with the closest semantic NOPs, we expect parallelization to improve
the running time.

In addition, our attack only adds instruction to the malware’s code sections.
Therefore, when the code section for the original malware is small, it will be
hard to provide enough perturbations to cause the classifier to misclassify
the malware image. Similar issues occur when the data section is much larger
than the code section. One possible approach to find an easier way to generate
perturbation is to find out the hot zone for the machine learning detector on
the malware in the first place and then only adding enough semantic NOPs to
these locations. Another approach to solving this challenge is to design a
mechanism to perturb the data section without affecting the original
functionality of the malware.

There are some potential defenses against general attacks to malware
detection. Tong et al. [39] proposed a method to boost the robustness of
feature-space models by identifying the conserved features and constraining
them in adversarial training. Their defense is evaluated only on PDF malware,
and we plan to evaluate further the defense on more malware families and more
diverse model types.

###  5.2 Future Work

The proposed attack algorithm is a proof-of-concept and a work-in-progress,
and there are several research directions we would like to explore as future
work:

  1. 1.

In the current version of the attack, the size and frequency of the added
perturbations are chosen beforehand by the attacker. In the development of the
attack, we would like to explore the possibility of adding these two
parameters (i.e. size and frequency of the perturbations) to the optimization
process. We speculate that this can improve the performance of our AE
generation algorithm. First, it can lead to faster convergence into a viable
AE. Additionally, it can also lead to a smaller perturbation size that is
customized for each malware sample.

  2. 2.

Another avenue of improvement to our baseline algorithm is in speed
performance as the baseline algorithm does not leverage any speed
optimizations. Our AE generation algorithm can allow batching and
parallelization to enhance the speed performance. First, we plan to batch
instructions to be processed at once instead of processing individually.
Additionally, since Euclidean distance calculations do not have to be done
sequentially we plan to calculate all the distances in parallel.

  3. 3.

Another avenue that we would like to explore is the intrinsic properties of
NOPs. Specifically, we would like to study the difference between NOPs and
whether some NOPs are better than others. Additionally, we would like to draw
from the field of software engineering in creating semantically similar code
using modifications that do not change the semantics at a higher-level
language (e.g. adding empty loops or if statements that would never be true).

  4. 4.

In our current implementation, we restrict adding perturbation to the code
section. In the algorithm, we would like to explore the effects of adding
perturbations to other sections and understanding its effects on
executability, robustness to pre-processing, and maintaining semantic
operations.

  5. 5.

We plan to perform a comprehensive comparison between our proposed attack with
the state-of-the-art attacks with respect to speed, stealthiness, and
deviation from the original binary. Additionally, we want to evaluate our
attack’s success against the defense and detection mechanisms available (e.g.
VirusTotal).

  6. 6.

To have an end-to-end solution we would like to add a module that checks
executability of each of the binary (e.g., using _execve_ command).

  7. 7.

Our current work has already revealed that malware detection based on binary
visualization can be beaten by an adversarial example with a high success
rate. In the meantime, the attacker can also maintain the functionality of the
original malware. Therefore, our work motivates future directions to propose
defenses against our attack. Previous defenses against adversarial examples,
such as defensive distillation [32] and adversarial training [14], does not
usually focus on real-world tasks, such as defending malware detection models.
Whether these defenses are effective in this use case is worth exploring.

##  6 Conclusion

In this work, we provide a literature review on existing techniques and
attacks for malware detection. We summarize their limitations, and propose a
novel end-to-end method for generating adversarial examples against
visualization-based malware detection. We design our attack workflow to ensure
that the malware sample remains executable while maintaining its original
functionality. Additionally, we ensure that the added perturbation is robust
against pre-processing by inserting semantic NOPs in the reachable code
section. We construct a dataset that includes both malware and benign software
samples and use it to build an visualization-based ML malware detector that
achieves high accuracy. Next, we design a workflow that generates semantic NOP
sequences and use them to construct viable adversarial examples. Our results
show that it is possible to successfully generate adversarial examples that
can bypass a highly accurate visualization-based ML malware detector while
maintaining executability and without changing the code operation. Our work
motivates the design for more robust visualization-based malware detection
against carefully crafted adversarial examples.

## References

  * [1] Abuhamad, M., Abuhmed, T., Mohaisen, A., Nyang, D.: Large-scale and language-oblivious code authorship identification. Proceedings of the 2018 ACM SIGSAC Conference on Computer and Communications Security (2018)
  * [2] Alam, S., Horspool, R.N., Traore, I., Sogukpinar, I.: A framework for metamorphic malware analysis and real-time detection. computers & security 48, 212–233 (2015)
  * [3] Arp, D., Spreitzenbarth, M., Hubner, M., Gascon, H., Rieck, K., Siemens, C.: Drebin: Effective and explainable detection of android malware in your pocket. In: Ndss. vol. 14, pp. 23–26 (2014)
  * [4] Carlini, N., Wagner, D.: Towards evaluating the robustness of neural networks. In: 2017 ieee symposium on security and privacy (sp). pp. 39–57. IEEE (2017)
  * [5] Chouchane, M.R., Lakhotia, A.: Using engine signature to detect metamorphic malware. In: Proceedings of the 4th ACM workshop on Recurring malcode. pp. 73–78 (2006)
  * [6] Christodorescu, M., Jha, S.: Testing malware detectors. ACM SIGSOFT Software Engineering Notes 29(4), 34–44 (2004)
  * [7] Christodorescu, M., Jha, S., Seshia, S.A., Song, D., Bryant, R.E.: Semantics-aware malware detection. In: 2005 IEEE Symposium on Security and Privacy (S&P’05). pp. 32–46. IEEE (2005)
  * [8] Clemens, J.: Automatic classification of object code using machine learning. Digital Investigation 14, S156–S162 (08 2015). https://doi.org/10.1016/j.diin.2015.05.007
  * [9] Demetrio, L., Coull, S.E., Biggio, B., Lagorio, G., Armando, A., Roli, F.: Adversarial exemples: A survey and experimental evaluation of practical attacks on machine learning for windows malware detection. arXiv preprint arXiv:2008.07125 (2020)
  * [10] Evtimov, I., Eykholt, K., Fernandes, E., Kohno, T., Li, B., Prakash, A., Rahmati, A., Song, D.: Robust physical-world attacks on machine learning models. ArXiv abs/1707.08945 (2017)
  * [11] Eykholt, K., Evtimov, I., Fernandes, E., Li, B., Rahmati, A., Xiao, C., Prakash, A., Kohno, T., Song, D.: Robust physical-world attacks on deep learning visual classification. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR) (June 2018)
  * [12] Fu, J., Xue, J., Wang, Y., Liu, Z., Shan, C.: Malware visualization for fine-grained classification. IEEE Access 6, 14510–14523 (2018)
  * [13] Go, J.H., Jan, T., Mohanty, M., Patel, O.P., Puthal, D., Prasad, M.: Visualization approach for malware classification with resnext. In: 2020 IEEE Congress on Evolutionary Computation (CEC). pp. 1–7. IEEE (2020)
  * [14] Goodfellow, I.J., Shlens, J., Szegedy, C.: Explaining and harnessing adversarial examples. arXiv preprint arXiv:1412.6572 (2014)
  * [15] Goodfellow, I.J., Shlens, J., Szegedy, C.: Explaining and harnessing adversarial examples. CoRR abs/1412.6572 (2015)
  * [16] Grosse, K., Papernot, N., Manoharan, P., Backes, M., McDaniel, P.: Adversarial perturbations against deep neural networks for malware classification. arXiv preprint arXiv:1606.04435 (2016)
  * [17] Han, K.S., Lim, J.H., Kang, B., Im, E.G.: Malware analysis using visualized images and entropy graphs. International Journal of Information Security 14, 1–14 (2014)
  * [18] Han, K.S., Lim, J.H., Kang, B., Im, E.G.: Malware analysis using visualized images and entropy graphs. International Journal of Information Security 14(1), 1–14 (2015)
  * [19] Han, K., Lim, J.H., Im, E.G.: Malware analysis method using visualization of binary files. In: Proceedings of the 2013 Research in Adaptive and Convergent Systems, pp. 317–321 (2013)
  * [20] Institute, A.T.: Malware, https://www.av-test.org/en/statistics/malware/
  * [21] Islam, R., Tian, R., Batten, L.M., Versteeg, S.: Classification of malware based on integrated static and dynamic features. Journal of Network and Computer Applications 36(2), 646–656 (2013)
  * [22] Khormali, A., Abusnaina, A., Chen, S., Nyang, D., Mohaisen, A.: Copycat: practical adversarial attacks on visualization-based malware detection. arXiv preprint arXiv:1909.09735 (2019)
  * [23] Kolosnjaji, B., Demontis, A., Biggio, B., Maiorca, D., Giacinto, G., Eckert, C., Roli, F.: Adversarial malware binaries: Evading deep learning for malware detection in executables. In: 2018 26th European Signal Processing Conference (EUSIPCO). pp. 533–537. IEEE (2018)
  * [24] Krcál, M., Svec, O., Bálek, M., Jasek, O.: Deep convolutional malware classifiers can learn from raw executables and labels only. In: ICLR (2018)
  * [25] Krügel, C., Robertson, W.K., Valeur, F., Vigna, G.: Static disassembly of obfuscated binaries. In: USENIX Security Symposium (2004)
  * [26] Liao, F., Liang, M., Dong, Y., Pang, T., Hu, X., Zhu, J.: Defense against adversarial attacks using high-level representation guided denoiser. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR) (June 2018)
  * [27] Liu, L., Wang, B.s., Yu, B., Zhong, Q.x.: Automatic malware classification and new malware detection using machine learning. Frontiers of Information Technology & Electronic Engineering 18(9), 1336–1347 (2017)
  * [28] Liu, X., Zhang, J., Lin, Y., Li, H.: Atmpa: Attacking machine learning-based malware visualization detection methods via adversarial examples. In: 2019 IEEE/ACM 27th International Symposium on Quality of Service (IWQoS). pp. 1–10. IEEE (2019)
  * [29] Makandar, A., Patrot, A.: Malware class recognition using image processing techniques. 2017 International Conference on Data Management, Analytics and Innovation (ICDMAI) pp. 76–80 (2017)
  * [30] Nataraj, L., Karthikeyan, S., Jacob, G., Manjunath, B.S.: Malware images: Visualization and automatic classification. In: Proceedings of the 8th International Symposium on Visualization for Cyber Security. VizSec ’11, Association for Computing Machinery, New York, NY, USA (2011). https://doi.org/10.1145/2016904.2016908, https://doi.org/10.1145/2016904.2016908
  * [31] Nataraj, L., Karthikeyan, S., Jacob, G., Manjunath, B.S.: Malware images: visualization and automatic classification. In: Proceedings of the 8th international symposium on visualization for cyber security. pp. 1–7 (2011)
  * [32] Papernot, N., McDaniel, P., Wu, X., Jha, S., Swami, A.: Distillation as a defense to adversarial perturbations against deep neural networks. In: 2016 IEEE symposium on security and privacy (SP). pp. 582–597. IEEE (2016)
  * [33] Pierazzi, F., Pendlebury, F., Cortellazzi, J., Cavallaro, L.: Intriguing properties of adversarial ml attacks in the problem space. In: 2020 IEEE Symposium on Security and Privacy (SP). pp. 1332–1349. IEEE (2020)
  * [34] Rad, B.B., Masrom, M., Ibrahim, S.: Camouflage in malware: from encryption to metamorphism. International Journal of Computer Science and Network Security 12(8), 74–83 (2012)
  * [35] Raff, E., Barker, J., Sylvester, J., Brandon, R., Catanzaro, B., Nicholas, C.: Malware detection by eating a whole exe. In: AAAI Workshops (2018)
  * [36] Shaid, S.Z.M., Maarof, M.A.: Malware behavior image for malware variant identification. In: 2014 International Symposium on Biometrics and Security Technologies (ISBAST). pp. 238–243. IEEE (2014)
  * [37] Sharif, M., Bhagavatula, S., Bauer, L., Reiter, M.: Accessorize to a crime: Real and stealthy attacks on state-of-the-art face recognition. Proceedings of the 2016 ACM SIGSAC Conference on Computer and Communications Security (2016)
  * [38] Sharif, M., Lucas, K., Bauer, L., Reiter, M.K., Shintre, S.: Optimization-guided binary diversification to mislead neural networks for malware detection. arXiv preprint arXiv:1912.09064 (2019)
  * [39] Tong, L., Li, B., Hajaj, C., Xiao, C., Zhang, N., Vorobeychik, Y.: Improving robustness of $\\{$ML$\\}$ classifiers against realizable evasion attacks using conserved features. In: 28th $\\{$USENIX$\\}$ Security Symposium ($\\{$USENIX$\\}$ Security 19). pp. 285–302 (2019)
  * [40] Xiaofang, B., Li, C., Weihua, H., Qu, W.: Malware variant detection using similarity search over content fingerprint. In: The 26th Chinese Control and Decision Conference (2014 CCDC). pp. 5334–5339. IEEE (2014)
  * [41] Zhang, G., Yan, C., Ji, X., Zhang, T., Zhang, T., Xu, W.: Dolphinattack: Inaudible voice commands. Proceedings of the 2017 ACM SIGSAC Conference on Computer and Communications Security (2017)
  * [42] Zhang, Q., Reeves, D.S.: Metaaware: Identifying metamorphic malware. In: Twenty-Third Annual Computer Security Applications Conference (ACSAC 2007). pp. 411–420. IEEE (2007)


WGMwgm QEqe EPep PMSpms BECbec DEde

[type=editor, auid=000,bioid=1,orcid=0000-0001-6634-1340] Conceptualization,
Methodology, Software, Writing- Original draft preparation.

Investigation, Writing- Reviewing and Editing.

# Simple El Niño prediction scheme using the signature of climate time series

Nozomi Sugiura  Shinya Kouketsu  Research Institute for Global Change, Japan
Agency for Marine-Earth Science and Technology, Yokosuka, Japan

###### Abstract

El Niño is a typical example of a coupled atmosphere–ocean phenomenon, but it
is unclear whether it can be described quantitatively by a correlation between
relevant climate events. To provide clarity on this issue, we developed a
machine learning-based El Niño prediction model that uses the time series of
climate indices. By transforming the multidimensional time series into the
path signature, the model is able to properly evaluate the order and
nonlinearity of climate events, which allowed us to achieve good forecasting
skill (mean square error = 0.596 for 6-month prediction). In addition, it is
possible to provide information about the sequence of climate events that tend
to change the future NINO3.4 sea surface temperatures. In forecasting
experiments conducted, changes in the North Pacific Index and several NINO
indices were found to be important precursors. The results suggest that El
Niño is predictable to some extent based on the correlation of climate events.

###### keywords:

signature El Niño time series analysis machine learning

##  1 Introduction

El Niño is an important climate phenomenon that has an immense socio-economic
impact. Consequently, its onset/offset mechanism has been garnering intense
scientific interest (Neelin et al., 1998; Wallace et al., 1998; Timmermann et
al., 2018), and the relationships with the other climate modes and events has
been investigated (Bjerknes, 1969; Alexander et al., 2002; White et al., 2014)
for many years. The prediction of El Niño events is still under investigation
from various perspectives, including statistical inference from past time
series of climate records and results from climate models, initialization from
climate models, and data assimilation using oceanic or coupled
atmospheric–oceanic models.

Although previous studies suggest that predictions with expensive climate
models outperform purely statistical predictions (Wu et al., 2021),
statistical predictions still appear to have value because of their simplicity
(Penland and Magorian, 1993). Recently, elaborate and quite skillful
predictions have been performed based on machine learning, or deep learning,
with the use of past oceanic sea surface temperatures (SSTs) and subsurface
information (Wang et al., 2020; Ham et al., 2019). Hu et al. (2021)
extensively evaluated the climate network method that utilizes the
relationship between spatio-temporal points and concluded that it has some
prediction skill over one year. Dijkstra et al. (2019) also reported that
machine learning models can improve the prediction skill over one year.
Nonetheless, few practical prediction studies have employed only the series of
multidimensional climate indices as learning datasets. As a remarkable
exception, Yan et al. (2020) used the NINO3.4 and Southern Oscillation indices
exclusively and successfully performed a skillful prediction of the NINO3.4
index via a temporal convolutional network. In light of the well-known fact
that the climate events and variabilities in the extratropics can also
modulate El Niño and the Southern Oscillation (ENSO) variability (Vimont et
al., 2003; Nakamura et al., 2006), the use of the various climate indices,
which represent the typical patterns of the ocean surface variables (e.g., sea
surface temperature and sea level pressure) in the course of climate changes,
can simply clarify the relationships between ENSO and the other climate modes
as well as improve the prediction efficiency of ENSO. However, a possible
disadvantage of statistical predictions is that they typically provide little
information about the correlations among these climate events through their
evolution.

To rectify this issue, we propose a new statistical method that is simple,
considerably skillful, and provides process information to explain how climate
events evolved. This study was conducted to develop a practical machine
learning-based El Niño prediction scheme using the past time series of climate
indices. The key ingredient that enables the faithful interpretation of the
past time series, including their nonlinearities, is the signature of paths,
which is a central concept in rough path theory (Lyons et al., 2007).

Although several studies have been published on the methodology of time series
analysis using the signature method (Morrill et al., 2021), there appears to
be no application of the method to global-scale climate events, and thus this
paper opens a new field of research in geosciences.

##  2 Methodology

In this study, we apply supervised learning to a time series of past climate
indices. We utilize the fact that each segment of the time series is an
explanatory variable that is equipped with future values at that time, which
can be regarded as objective variables. The most significant aspect of the
proposed method is that each segment of the time series is transformed into a
signature. Therefore, our case study simply employs the simplest setting,
utilizing the signature method to concentrate on the proof-of-concept. In this
section, after explaining the theoretical basis of why the signatures are
relevant, we present the machine learning procedure based on that theory.
Then, we discuss how to interpret the results and, finally, we present the
parameters used.

###  2.1 Approximating a function on a set of paths

In prediction study based on the learning of time series, it is crucial to
properly construct a predictor that link past time series segments to future
values. A predictor is represented as a continuous function on a set of
multidimensional paths. To secure the performance of the predictor, it is
essential to choose an appropriate basis for the function of the path because
it determines the expressive power of the function. Note that our concern is
not the basis for a path but for a function of paths. In this sense, the most
mathematically justified candidate for the basis is the signature (Lyons et
al., 2007; Sugiura and Hosoda, 2020).

For a $d$-dimensional path $X=X_{[s,t]}:[s,t]\to\mathbb{R}^{d}$ that maps
$\tau$ to $X_{\tau}$, the $0$-th to $n$-th iterated integrals are defined
recursively, as follows (Lyons et al., 2007):

$\displaystyle\mathcal{S}^{()}(X_{[s,t]})$ $\displaystyle=1,$ (1)
$\displaystyle\mathcal{S}^{(i_{1}\cdots i_{n})}(X_{[s,t]})$
$\displaystyle=\int_{s}^{t}\mathcal{S}^{(i_{1}\cdots
i_{n-1})}(X_{[s,t_{n}]})dX^{(i_{n})}_{t_{n}},\quad
i_{1},\cdots,i_{n}=1,\cdots,d.$ (2)

The signature $\mathcal{S}(X)$ of path $X$ is the collection of all the
iterated integrals, and the operation $\mathcal{S}:X\mapsto\mathcal{S}(X)$ is
called the signature transform. In particular, its truncation up to the $n$-th
iterated integrals is called the step-$n$ signature,
$\mathcal{S}_{n}(X)\in\bigoplus_{k=0}^{n}(\mathbb{R}^{d})^{\otimes k}$, which
means that the multi-index $I$ in component $\mathcal{S}_{n}^{(I)}(X_{[s,t]})$
runs across

$\displaystyle I$
$\displaystyle\in\phi\cup\\{1,\cdots,d\\}\cup\cdots\cup\\{1,\cdots,d\\}^{n}.$
(3)

Now, let $C(K,\mathbb{R})$ be the space of a continuous function on a compact
set $K$ of paths with at least one monotonous coordinate. Subset $A\subset
C(K,\mathbb{R})$ is defined as

$\displaystyle A$
$\displaystyle=\left\\{g:X\mapsto\sum_{I}w^{(I)}\mathcal{S}^{(I)}(X)\middle|w\in\bigoplus_{k=0}^{n}(\mathbb{R}^{d})^{\otimes
k}\text{~{}for some $n\geq 0$}\right\\}.$ (4)

Then, $A$ satisfies the following conditions:

  1. 1.

Because the step-$n$ signature transform $K\ni X\mapsto\mathcal{S}_{n}(X)$ is
continuous for any $n>0$, $A\subset C(K,\mathbb{R})$.

  2. 2.

$g_{1},g_{2}\in A\text{ and
}\lambda_{1},\lambda_{2}\in\mathbb{R}\implies\lambda_{1}g_{1}+\lambda_{2}g_{2}\in
A$.

  3. 3.

Constant-valued function $\mathbf{1}\in A$.

  4. 4.

Based on the shuffle identity (Lyons et al., 2007), $g_{1},g_{2}\in A\implies
g_{1}g_{2}\in A$.

  5. 5.

Based on the uniqueness theorem (Levin et al., 2013), for all $X,Y\in A$ with
$X\neq Y$, there exists $g\in A$ that satisfies $g(X)\neq g(Y)$.

From these conditions, we can apply the Stone–Weierstrass theorem (Stone,
1937) to the subset $A$ and conclude that $A$ is dense in $C(K,\mathbb{R})$,
which means that any function $f\in C(K,\mathbb{R})$ is uniformly approximated
by a function $g\in A$ with arbitrary accuracy (Levin et al., 2013; Fermanian,
2021).

From the above reasoning, we can construct a nonlinear predictor as a linear
combination of iterated integrals for each segment of the multidimensional
time series.

###  2.2 Procedure for machine learning of time series

In the proposed approach, the predictor is constructed as follows. Suppose we
have a time series of several climate indices defined at each calendar month
of $t_{m},~{}m=1,2,\cdots,M$. Any segment of the time series over a period,
say six months, can be viewed as a multidimensional path, which can be
represented by the signature.

For this supervised learning, the object variable is the NINO3.4 index
$y_{m+m_{a}}$ at time $\tau=t_{m+m_{a}}$, while the explanatory variables
$x_{m}$ are the signature for the segment of time series $X$ in the period
$[t_{m-m_{b}+1},t_{m}]$. The approximation property described in the previous
section allows us to express the object variable as a linear combination of
the explanatory variables:

$\displaystyle y_{m+m_{a}}$
$\displaystyle=y_{m}+\left<w_{m},x_{m}\right>+\epsilon,$ (5) $\displaystyle
x_{m}$ $\displaystyle:=\mathcal{S}_{n}(X_{[t_{m-m_{b}+1},t_{m}]}),$ (6)

where $\mathcal{S}_{n}(X_{[t_{0},t_{1}]})$ denotes the order-$n$ signature for
the $d$-dimensional time series in the interval $[t_{0},t_{1}]$,
$\left<a,b\right>$ denotes the scalar product $\sum_{I}a^{(I)}b^{(I)}$,
$w_{m}=\\{w_{m}^{(I)}|I=\text{multi-index}\\}$ is the weight vector for the
predictor, $\epsilon$ is a random variable representing prediction error,
$t_{m}$ is the starting time of the prediction, $t_{m-m_{b}+1}$ is the
starting point of the path segment, and $t_{m+m_{a}}$ is the target time for
prediction.

Figure 1: Schematic view of training and prediction flow, assuming that
$m_{b}$=$m_{a}$=6. Hatched squares represent transforming into the signature.
Predicted value at time index $m$+6 will be compared with the validation data
if available.

Before converting into the signature, a zero vector is added at the beginning
of each series $X_{[t_{m-m_{b}+1},t_{m}]}$ to account for the magnitude of the
value at the starting point (Morrill et al., 2021). We computed the signature
by using the Python library ${\tt esig}$ (Kormilitzin, 2017).

In the control case, we instead used the time series as is without converting
it into the signature:

$\displaystyle x_{m}$
$\displaystyle:=X_{[t_{m-m_{b}+1},t_{m}]}=\left(X_{t_{m-m_{b}+1}},X_{t_{m-m_{b}+2}},\cdots,X_{t_{m}}\right).$
(7)

This corresponds to an auto-regressive (AR) model.

Using the training dataset available up to time $t_{m}$,

$\displaystyle D_{m}$
$\displaystyle=\left\\{\left(x_{\mu},y_{\mu+m_{a}}-y_{\mu}\right)\middle|~{}\mu\in[m_{b},m-m_{a}]\right\\},$
(8)

we first estimate the optimal weight $w=w_{m}$ that minimizes the cost
function with an $L_{1}$-penalty term:

$\displaystyle J_{m}(w)$
$\displaystyle=\frac{1}{2|D_{m}|}\sum_{\mu=m_{b}}^{m-m_{a}}\left(y_{\mu+m_{a}}-y_{\mu}-\left<w,x_{\mu}\right>\right)^{2}+\alpha\sum_{I}\left|w^{(I)}\right|,$
(9)

where $|D_{m}|=m-m_{a}-m_{b}+1$ is the number of samples in $D_{m}$, and $I$
is the multi-index. The optimization problem is solved by the Lasso model fit
with least angle regression (Pedregosa et al., 2011), which is suitable for
problems with many parameters. We then predict a future NINO3.4 index as
$\widehat{y}_{m+m_{a}}=y_{m}+\left<w_{m},x_{m}\right>$ and compare it to
$y_{m+m_{a}}$. In other words, a cross-validation is made against the
validation data:

$\displaystyle D^{\prime}_{m}$
$\displaystyle=\left\\{\left(x_{m},y_{m+m_{a}}-y_{m}\right)\right\\}.$ (10)

We repeat the above procedure after incrementing the time index $m$ by $1$.

Figures 1 shows the schematic view of training and prediction flow. In this
flow, the weight $w_{m}$ is obtained by using the training dataset $D_{m}$,
and then the prediction from time $t_{m}$ using the signature $x_{m}$ yields
the value $\widehat{y}_{m+m_{a}}$, which is subject to comparison with the
validation data $y_{m+m_{a}}$. Note that the size of the training data
$|D_{m}|=m-m_{a}-m_{b}+1$ depends on the starting time $t_{m}$. The prediction
error can be obtained from the statistics of
$\widehat{y}_{m+m_{a}}-y_{m+m_{a}}$ for various starting times. By taking this
approach, where training and forecasts are done progressively by moving the
starting time of the forecast hiding future at that moment, each forecast is
assured to be a fair cross-validation.

###  2.3 Diagnosing the dominant event sequences

One difficulty with regular machine learning is that it does not provide
sufficient reasoning for the results. However, the signature-based method
allows us to mathematically extract from the path those properties that are
important in the prediction.

To diagnose the dominant event sequences that contribute to the prediction, we
compute the standard partial regression coefficients (SPRCs) $r_{m}^{(I)}$,
which represent the sensitivity of normalized value $y_{\mu+m_{a}}$ in the
future to each component of the normalized signature $x_{\mu}^{(I)}$ in the
past, as

$\displaystyle r_{m}^{(I)}$
$\displaystyle=\frac{\sigma_{x_{m}^{(I)}}}{\sigma_{y_{m}}}w_{m}^{(I)},$ (11)

where $\sigma_{y_{m}}$ and $\sigma_{x_{m}^{(I)}}$ denote the standard
deviations of $y_{\mu+m_{a}}-y_{\mu}$ and $x_{\mu}^{(I)}$, respectively, among
the samples in $D_{m}$, which represents the learning data in the period from
time $t_{m_{b}}$ to time $t_{m-m_{a}}$.

###  2.4 Setting of experimental parameters

We used a climate time series composed of $d=12$ indices in Table 1 retrieved
from NOAA cite (NOAA, 2021). The time series starts at $t_{1}$ (January of
1900), and ends at $t_{1459}$ (July of 2021).

Table 1: Twelve climate indices and their abbreviations Abbrev. | Climate Indices | References
---|---|---
NINO34 | Nino 3.4 (5N-5S, 170W-120W) SST | Rayner et al. (2003)
NINO12 | Nino 1+2 (0-10S, 90W-80W) SST | Rayner et al. (2003)
NINO3 | Nino 3 (5N-5S, 150W-90W) SST | Rayner et al. (2003)
NINO4 | Nino 4 (5N-5S, 160E-150W) SST | Rayner et al. (2003)
DMI | Dipole Mode Index | Hameed and Yamagata (2003)
AMO | Atlantic Multidecadal Oscillation index | Enfield et al. (2001)
NPI | North Pacific Index | Trenberth and Hurrell (1994)
SOI | Southern Oscillation Index | Ropelewski and Jones (1987)
NAO | North Atlantic Oscillation (NAO) index | Jones et al. (1997)
TPI | Tripole Index | Henley et al. (2015)
AO | Arctic Oscillation index | Thompson and Wallace (1998)
MON | Date elapsed (mid-day in month divided by 365) | -

The standard lead time for prediction is $6$ months ($m_{a}=6$), whereas each
past segment is of length $6$ months ($m_{b}=6$). The experiment duration was
from the prediction starting at $t_{961}$ (January of 1980), to the one
starting at $t_{1453}$ (January of 2021).

We used iterated integrals up to level $n=3$, which means that the total
number of terms in the linear combination was $N=(d^{n+1}-1)/(d-1)=1885$, and
the intensity of the $L_{1}$ penalty term was tuned to $\alpha=2.0$.

##  3 Results

Figure 2 shows the result of $6$-month prediction, whereas fig 3 shows it as
anomalies from climatology, These figures indicate the overall superiority of
prediction using the signature. The prediction error for each target month is
shown in Fig. 4. It is obvious that the predictions for July to September were
much better than those in the control case; however, they were comparable in
the other months. The overall prediction skill was $0.596\mathrm{K}$ for the
signature case and $0.663\mathrm{K}$ for the control case.

Wang et al. (2020) proposed an operator-theoretic technique called kernel
analog forecasting (KAF), which has a rigorous connection with Koopman
operator theory for dynamical systems, yielding statistically optimal
predictions as conditional expectations. They also compared it to the linear
inverse model (LIM), which is a linear evolution operator for modes. Note that
both methods employ as explanatory variables the dominant modes in
spatiotemporal SST, which we did not use in our study. Here, we use KAF and
LIM for the comparison of forecasting skill. For the comparison with KAF and
LIM, the root-mean square (rms) errors for 6-month prediction in the period
from 1998 to 2017 were computed. The signature model, AR model, KAF model, and
LIM had rms values of $0.617$, $0.686$, $0.62$, and $0.75\mathrm{K}$,
respectively. This comparison result suggests that the forecasting skill of
the signature model is comparable to that of the KAF model.

The spring prediction barrier is defined in Lai et al. (2018) as follows: “…
models have problems in predicting Boreal winter tropical Pacific sea surface
temperature (SST) when forecasts start in Boreal spring (February–May). This
is called the spring predictability barrier. ” Similarly, Zheng and Zhu (2010)
pointed out that “… errors have the largest values and the fastest growth
rates initialized before and during the NH spring.” In light of these
definitions, the spring predictability barrier, i.e., poor prediction skill
when starting from February and March, seems to disappear as indicated by the
rms error values in the target months of August to September.

Table 2 shows the dominant event sequences among iterated integrals. The
events with the first to the third indices are shown in each row. If the same
index appears twice in a row, then the event is intense. The top sequence in
the period from 1900 to 2020 is an intense NPI change followed by a Niño 1+2
SST change. The key indices are NPI and various NINO indices. In particular,
NPI, an atmospheric process, is involved in all the dominant sequences, which
should be a manifestation that El Niño is a coupled atmospheric–oceanic
process. In addition, the comparison between statistics for two different
periods suggests that the Nino1+2 index, corresponding to the region of
coastal South America, is becoming more important as a precursor in the 21st
century. Summarizing the above, fig. 5 illustrates how the dominant climate
events occur that will lead to changes in the future NINO3.4.

Although the main indices in terms of the iterative integrals are related to
water temperature in the NINO regions, it appears that not only these but also
various climate indices contribute incrementally, and the predictor is built
on the balance of them. In fact, if we perform an experiment with the main 5
indices (NINO12, NINO3, NINO34, NINO4, and NPI), the rms error for the
prediction is $0.666\mathrm{K}$, with no improvement from the control case. In
this respect, the inference structure is considered to be different from other
SST-based predictions such as KAF and LIM.

Table 2: Top five dominant event sequences among iterated integrals. “1st” denotes the first index for the corresponding iterated integral: $\mathcal{S}^{(i_{1}i_{2}i_{3})}(X)=\int_{s}^{t}\int_{s}^{t_{3}}\int_{s}^{t_{2}}dX_{t_{1}}^{(i_{1})}dX_{t_{2}}^{(i_{2})}dX_{t_{3}}^{(i_{3})}$. Events happen from first to third: $t_{1}<t_{2}<t_{3}$. If the same index appears twice in a row, then the event is intense. “SPRC” represents the standard partial regression coefficients (Eq. 11).  | Learning data from Jan. 1900 to Dec. 1999
---|---
| SPRC | Indices in iterated integrals
No. | $r_{m}^{(i_{1}i_{2}i_{3})}$ | 1st ($i_{1}$) | 2nd ($i_{2}$) | 3rd ($i_{3}$)
1 | $5.55$ | NPI | NINO3 | NPI
2 | $-3.70$ | NINO3 | NPI | NPI
3 | $3.63$ | NPI | NPI | NINO12
4 | $-3.57$ | NINO34 | NINO34 | NINO34
5 | $3.34$ | NINO34 | NPI | NPI
| Learning data from Jan. 1900 to Dec. 2020
| SPRC | Indices in iterated integrals
No. | $r_{m}^{(i_{1}i_{2}i_{3})}$ | 1st ($i_{1}$) | 2nd ($i_{2}$) | 3rd ($i_{3}$)
1 | $4.53$ | NPI | NPI | NINO12
2 | $-4.17$ | NINO34 | NINO34 | NPI
3 | $-3.37$ | NPI | NINO34 | NINO12
4 | $-3.25$ | NINO34 | NPI | NINO12
5 | $2.88$ | NPI | NINO3 | NPI
Figure 2: Comparison of NINO3.4 for 6-month predictions. Red: signature case;
blue: control case. Horizontal axis is the target year and month, and vertical
axis is temperature in $\mathrm{{}^{\circ}C}$.  Figure 3: The same as Fig. 2
but shown as anomalies, which are defined as the difference from the past
30-yr mean of monthly values. Red: signature case; blue: control case.
Horizontal axis is the target year and month, and vertical axis is temperature
anomalies in $\mathrm{{}^{\circ}C}$.  Figure 4: Prediction error for each
target month. Red: signature case; blue: control case. Horizontal axis is the
target month (1 = January, 2 = February, $\cdots$, 12 = December), and
vertical axis is rms error in $\mathrm{K}$.  Figure 5: Typical climate event
flows for predicting future Nino3.4 index. Arrows indicate time order. Key
indices include NINO12 (Nino1+2 SST), NINO34 (Nino3.4 SST), NINO4 (Nino4 SST),
and NPI (North Pacific Index).

##  4 Conclusions

We developed a model that can statistically predict El Niño using only the
time series of past multidimensional climate indices. By converting the time
series into the signature, the accuracy of the machine learning algorithm is
improved and, thereby, the NINO3.4 SST can be predicted to some extent six
months in advance. An important byproduct of this approach is that the
correlation of climate events can be read from the dominant iterative
integral. For example, it was suggested that variations in the NPI, NINO12,
and other indices occur in a certain order, which leads to variations in the
NINO3.4 SST. It was also found that the signature method can learn the
nonlinear development of El Niño more accurately than the traditional AR model
and, thus, is less sensitive to the spring barrier of predictability. Future
research is required to improve the scheme by incorporating more detailed
oceanographic information, evaluating uncertainties, and considering other
factors.

The predictions obtained by this method are not marked with error bars, but
because we know the prediction error for each month as shown in Fig. 4, we can
consider these values as the prediction error. However, as it is not possible
to give a forecast error for each forecast individually, an ensemble could be
created by bootstrapping or other methods to improve this point, which may
also lead to a factory for forecast accuracy.

The length of the path segment used for the 6-month forecast, 6 months, was
confirmed in preliminary experiments (not shown) to be appropriate, but the
length of the path segment required for forecasts with other lead times may be
different. It is also necessary to confirm whether the step-$3$ signature is
optimal.

The dominant iterated integral for prediction may change from time to time
depending on the period covered, as shown in Table 2. It needs to be carefully
considered how this relates to the decadal changes in the Niño mechanism.
These points remain as future work.

##  5 Acknowledgments

This study was funded by JST-PROJECT-20218919.

Code availability section

Name: enso_signature

Contact<EMAIL_ADDRESS>+81-46-867-9054

Hardware requirements: CPU

Program language: Python

Software required: Python libraries esig and scikit-learn

Program size: 188 lines

The source codes are available for downloading at the link:
https://github.com/nozomi-sugiura/enso_signature

## References

  * Alexander et al. (2002) Alexander, M.A., Bladé, I., Newman, M., Lanzante, J.R., Lau, N.C., Scott, J.D., 2002\.  The atmospheric bridge: The influence of ENSO teleconnections on air–sea interaction over the global oceans.  Journal of climate 15, 2205–2231.
  * Bjerknes (1969) Bjerknes, J., 1969.  Atmospheric teleconnections from the equatorial Pacific.  Monthly weather review 97, 163–172.
  * Dijkstra et al. (2019) Dijkstra, H.A., Petersik, P., Hernández-García, E., López, C., 2019.  The application of machine learning techniques to improve El Niño prediction skill.  Frontiers in Physics 7, 153\.
  * Enfield et al. (2001) Enfield, D.B., Mestas-Nuñez, A.M., Trimble, P.J., 2001.  The atlantic multidecadal oscillation and its relation to rainfall and river flows in the continental u.s.  Geophysical Research Letters 28, 2077–2080.  URL: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2000GL012745, doi:https://doi.org/10.1029/2000GL012745, arXiv:https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2000GL012745.
  * Fermanian (2021) Fermanian, A., 2021.  Embedding and learning with signatures.  Computational Statistics & Data Analysis 157, 107148.
  * Ham et al. (2019) Ham, Y.G., Kim, J.H., Luo, J.J., 2019.  Deep learning for multi-year enso forecasts.  Nature 573, 568--572.  doi:https://doi.org/10.1038/s41586-019-1559-7.
  * Hameed and Yamagata (2003) Hameed, S., Yamagata, T., 2003\.  Possible impacts of indian ocean dipole mode events on global climate.  Climate Research - CLIMATE RES 25, 151--169.  doi:10.3354/cr025151.
  * Henley et al. (2015) Henley, B.J., Gergis, J., Karoly, D.J., Power, S., Kennedy, J., Folland, C.K., 2015\.  A Tripole Index for the Interdecadal Pacific Oscillation.  Climate Dynamics 45, 3077--3090.  doi:10.1007/s00382-015-2525-1.
  * Hu et al. (2021) Hu, X., Eichner, J., Faust, E., Kantz, H., 2021\.  Benchmarking prediction skill in binary el niño forecasts.  Climate Dynamics , 1--15.
  * Jones et al. (1997) Jones, P.D., Jonsson, T., Wheeler, D., 1997.  Extension to the north atlantic oscillation using early instrumental pressure observations from gibraltar and south-west iceland.  International Journal of Climatology 17, 1433--1450.  doi:https://doi.org/10.1002/(SICI)1097-0088(19971115)17:13<1433::AID-JOC203>3.0.CO;2-P.
  * Kormilitzin (2017) Kormilitzin, A., 2017.  the-signature-method-in-machine-learning.  https://github.com/kormilitzin/.
  * Lai et al. (2018) Lai, A.W.C., Herzog, M., Graf, H.F., 2018.  Enso forecasts near the spring predictability barrier and possible reasons for the recently reduced predictability.  Journal of Climate 31, 815--838.
  * Levin et al. (2013) Levin, D., Lyons, T., Ni, H., 2013.  Learning from the past, predicting the statistics for the future, learning an evolving system.  ArXiv e-prints arXiv:1309.0260.
  * Lyons et al. (2007) Lyons, T.J., Caruana, M., Lévy, T., 2007.  Differential Equations Driven by Rough Paths. volume 1908 of Lecture Notes in Mathematics.  Springer.
  * Morrill et al. (2021) Morrill, J., Fermanian, A., Kidger, P., Lyons, T., 2021\.  A Generalised Signature Method for Multivariate Time Series Feature Extraction.  arXiv:2006.00873.
  * Nakamura et al. (2006) Nakamura, T., Tachibana, Y., Honda, M., Yamane, S., 2006\.  Influence of the Northern Hemisphere annular mode on ENSO by modulating westerly wind bursts.  Geophysical Research Letters 33.
  * Neelin et al. (1998) Neelin, J.D., Battisti, D.S., Hirst, A.C., Jin, F.F., Wakata, Y., Yamagata, T., Zebiak, S.E., 1998.  Enso theory.  J. Geophys. Res. 103, 14 261--14 290.
  * NOAA (2021) NOAA, 2021.  Climate Timeseries at PSL.  URL: https://psl.noaa.gov/gcos_wgsp/Timeseries/.
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., Duchesnay, E., 2011\.  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research 12, 2825--2830.
  * Penland and Magorian (1993) Penland, C., Magorian, T., 1993\.  Prediction of Niño 3 Sea Surface Temperatures Using Linear Inverse Modeling.  Journal of Climate 6, 1067 -- 1076.  URL: https://journals.ametsoc.org/view/journals/clim/6/6/1520-0442_1993_006_1067_ponsst_2_0_co_2.xml, doi:10.1175/1520-0442(1993)006<1067:PONSST>2.0.CO;2.
  * Rayner et al. (2003) Rayner, N.A., Parker, D.E., Horton, E.B., Folland, C.K., Alexander, L.V., Rowell, D.P., Kent, E.C., Kaplan, A., 2003\.  Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century.  Journal of Geophysical Research: Atmospheres 108\.  URL: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2002JD002670, doi:https://doi.org/10.1029/2002JD002670, arXiv:https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2002JD002670.
  * Ropelewski and Jones (1987) Ropelewski, C.F., Jones, P.D., 1987\.  An extension of the tahiti--darwin southern oscillation index.  Monthly weather review 115, 2161--2165.
  * Stone (1937) Stone, M., 1937.  Applications of the theory of boolean rings to general topology.  Transactions of the American Mathematical Society 41, 375--481.
  * Sugiura and Hosoda (2020) Sugiura, N., Hosoda, S., 2020\.  Machine learning technique using the signature method for automated quality control of argo profiles.  Earth and Space Science 7, e2019EA001019.  URL: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2019EA001019, doi:https://doi.org/10.1029/2019EA001019, arXiv:https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2019EA001019. e2019EA001019 10.1029/2019EA001019.
  * Thompson and Wallace (1998) Thompson, D.W.J., Wallace, J.M., 1998\.  The arctic oscillation signature in the wintertime geopotential height and temperature fields.  Geophysical Research Letters 25, 1297--1300.  URL: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/98GL00950, doi:https://doi.org/10.1029/98GL00950, arXiv:https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/98GL00950.
  * Timmermann et al. (2018) Timmermann, A., An, S.I., Kug, J.S., Jin, F.F., Cai, W., Capotondi, A., Cobb, K.M., Lengaigne, M., McPhaden, M.J., Stuecker, M.F., et al., 2018\.  El niño--southern oscillation complexity.  Nature 559, 535--545.
  * Trenberth and Hurrell (1994) Trenberth, K., Hurrell, J., 1994\.  Decadal atmosphere-ocean variations in the Pacific.  Climate Dynamics 9, 303--319.  doi:https://doi.org/10.1007/BF00204745.
  * Vimont et al. (2003) Vimont, D.J., Wallace, J.M., Battisti, D.S., 2003.  The seasonal footprinting mechanism in the Pacific: Implications for ENSO.  Journal of Climate 16, 2668--2675.
  * Wallace et al. (1998) Wallace, J.M., Rasmusson, E.M., Mitchell, T.P., Kousky, V.E., Sarachik, E.S., von Storch, H., 1998\.  On the structure and evolution of enso related climate variability in the tropical pacific: Lessons from toga.  J. Geophys. Res. 103, 14 241--14 259.
  * Wang et al. (2020) Wang, X., Slawinska, J., Giannakis, D., 2020.  Extended-range statistical ENSO prediction through operator-theoretic techniques for nonlinear dynamics.  Sci Rep. 10, 2636\.  doi:10.1038/s41598-020-59128-7.
  * White et al. (2014) White, C.J., Hudson, D., Alves, O., 2014.  ENSO, the IOD and the intraseasonal prediction of heat extremes across Australia using POAMA-2.  Climate dynamics 43, 1791--1810.
  * Wu et al. (2021) Wu, X., Okumura, Y.M., Deser, C., DiNezio, P.N., 2021\.  Two-Year Dynamical Predictions of ENSO Event Duration during 1954--2015.  Journal of Climate 34, 4069 -- 4087.  URL: https://journals.ametsoc.org/view/journals/clim/34/10/JCLI-D-20-0619.1.xml, doi:10.1175/JCLI-D-20-0619.1.
  * Yan et al. (2020) Yan, J., Mu, L., Wang, L., Ranjan, R., Zomaya, A.Y., 2020\.  Temporal convolutional networks for the advance prediction of ENSO.  Scientific reports 10, 1--15.
  * Zheng and Zhu (2010) Zheng, F., Zhu, J., 2010.  Spring predictability barrier of ENSO events from the perspective of an ensemble prediction system.  Global and Planetary Change 72, 108--117.