# Detecting Distributional Differences in Labeled Sequence Data with
Application to Tropical Cyclone Satellite Imagery

Trey McNeely  Galen Vincent  Kimberly M Wood  Rafael Izbicki  Ann B Lee
Department of Statistics and Data Science, Carnegie Mellon University
Department of Geosciences, Mississippi State University  Department of
Statistics, Federal University of São Carlos

###### Abstract

Our goal is to quantify whether, and if so how, spatio-temporal patterns in
tropical cyclone (TC) satellite imagery signal an upcoming rapid intensity
change event. To address this question, we propose a new nonparametric test of
association between a time series of images and a series of binary event
labels. We ask whether there is a difference in distribution between
(dependent but identically distributed) 24-h sequences of images preceding an
event versus a non-event. By rewriting the statistical test as a regression
problem, we leverage neural networks to infer modes of structural evolution of
TC convection that are representative of the lead-up to rapid intensity change
events. Dependencies between nearby sequences are handled by a bootstrap
procedure that estimates the marginal distribution of the label series. We
prove that type I error control is guaranteed as long as the distribution of
the label series is well-estimated, which is made easier by the extensive
historical data for binary TC event labels. We show empirical evidence that
our proposed method identifies archetypes of infrared imagery associated with
elevated rapid intensification risk, typically marked by deep or deepening
core convection over time. Such results provide a foundation for improved
forecasts of rapid intensification.

two-sample testing,

high-dimensional time series,

association studies,

remote sensing,

functional data,

weather forecasting,

###### keywords:

, , , , and

##  1 Introduction

A broad array of problems in the physical, environmental and biological
sciences feature high-dimensional time series $\\{\mathbf{X}_{t}\\}_{t\geq
0}$, associated with binary “labels” $\\{Y_{t}\\}_{t\geq 0}$ indicating an
event of interest. Examples include sequences of satellite or other remote
sensing data paired with natural events like the occurrence of an earthquake
or the rapid intensification of a hurricane; or multivariate
electroencephalographic (EEG) and magnetoencephalographic (MEG) data showing
brain activity paired with physiological events like the occurrence of a
stroke. Most research on this front concerns prediction of events (Luo et al.,
2014), measurement of event impact (Scharwächter and Müller, 2020a, b), or
detection of change-points (Aminikhanghahi and Cook, 2017; Evans and G’Sell,
2020) after events occur. Furthermore, joint analyses of time and event series
often assume that the time series is univariate, or model the relationship
between multiple scalar quantities as they change over time. There is a lack
of theoretically and computationally sound methods for (non-parametric)
association studies and statistical tests for dependent and high-dimensional
sequence data.

This work is motivated by the need to identify spatio-temporal patterns in the
convective evolution of tropical cyclone satellite imagery prior to a rapid
intensity change; see “Motivating Application.” Our immediate goal is not
operational forecasting or prediction per se, but rather gaining scientific
insight into the spatio-temporal evolution
$\mathbf{S}_{<t}=\\{\mathbf{X}_{t-T},\mathbf{X}_{t-T+1},\ldots,\mathbf{X}_{t}\\}$
of convective structure or satellite imagery $\mathbf{X}_{t}$, leading up to a
rapid intensity change event ($Y_{t}=1$), for some lead time $T,$ and
identifying whether it differs in distribution from sequences
$\mathbf{S}_{<t}$ that precede a non-event ($Y_{t}=0$).

From a statistical methodology standpoint, this problem amounts to a
challenging two-sample testing problem for high-dimensional dependent but
identically distributed (DID) data. From the observed time and event series,
we extract labeled sequence data $\\{(\mathbf{S}_{<t},Y_{t})\\}_{t\geq 0}$,
which we assume is a stationary process. That is, both $\mathbf{S}_{<t}$ and
$Y_{t}$ are auto-correlated and dependent for different instances of time $t$.
By the assumption of stationarity, the data $(\mathbf{S}_{<t},Y_{t})$ are
identically distributed. Given historical data, we test whether the
distributions of $\mathbf{S}_{<t}|Y_{t}=1$ and $\mathbf{S}_{<t}|Y_{t}=0$ are
the same. The challenge is to construct an efficient test of (1) that is valid
(controls type I error) for DID data, and that applies to different types of
sequence data (images, functions, and sequentially observed data from multiple
physical probes) with a minimum of assumptions. Finally, for many problems of
applied interest, scientists want to know not just whether sequences preceding
an event versus non-event are significantly different in distribution, but if
so, also how the two distributions are different. That is, if the null
hypothesis that the distributions of (stationary) sequences $\mathbf{S}|Y=1$
and $\mathbf{S}|Y=0$ are the same is rejected, the question is how to identify
the patterns in the state space $\mathcal{S}$ of $\mathbf{S}$ that contributed
to the rejection.111The stationarity assumption makes the state space
$\mathcal{S}$ well-defined, and allows us to drop the subindex $<t$ in our
notation. These patterns correspond to sequences $\mathbf{s}\in\mathcal{S}$
that are more or less likely to be associated with an event ($Y=1$) than by
chance.

Motivating Application: Tropical Cyclones. Tropical cyclones (TCs) are highly
structured storms which rank among the deadliest and costliest natural
disasters in the United States (Klotzbach et al., 2018). Cases of rapid
intensification (RI) and rapid weakening (RW) of such storms — defined for
this work as a change in maximum wind speed of at least 25 knots within 24
hours, $Y=1$ — are notoriously difficult to forecast (Kaplan and DeMaria,
2003; Kaplan, DeMaria and Knaff, 2010; Kaplan et al., 2015; Wood and Ritchie,
2015). RI prediction has thus been the “highest-priority forecast challenge”
identified by the National Hurricane Center (NHC) in the last decade (Gall et
al., 2013). RW events, while a lower priority than RI, are also associated
with above-average forecast errors and are of great interest to
meteorologists.

Models such as SHIPS-RII (Kaplan et al., 2015) have made great progress on
skillfully forecasting RI events, but these approaches rely on scalar
predictors (e.g., area-averaged vertical wind shear) at fixed points in time
and thus neglect the evolving spatial structure of the TC, and these
structural changes often influence such events. To address this gap,
meteorologists and forecasters seek interpretable patterns in the
spatiotemporal structure of physically-relevant 2D fields that could indicate
an elevated risk of RI. The first step in this search is to find interpretable
temporally-evolving sequences of spatial structure $\mathbf{S}_{<t}$ that
differ in distribution depending on whether a TC is undergoing a rapid
intensity change event ($Y_{t}=1$) or not ($Y_{t}=0$).

Our application examines one such 2D field: deep convection within the storm.
Convection, or deep thunderstorm-like clouds, is a key component of the
mechanism through which TCs extract energy from the ocean, meaning that
convection, in both strength and distribution, should be closely related to
storm intensity. We quantify convective structure using cloud-top temperature
as measured by infrared (IR) imagery from the Geostationary Operational
Environmental Satellites (GOES). As convection strengthens, the cloud tops are
pushed higher into the atmosphere where temperatures are lower. The
temperature of the cloud top can therefore be used as a proxy for the strength
of convection in the storm. We ask: “Do 24-hour sequences of convective
structure, $\\{\mathbf{S}_{<t}\\}_{t\geq 0}$, contain information about
upcoming intensity change, $\\{Y_{t}\\}_{t\geq 0}$? If so, how do the spatio-
temporal patterns of rapidly changing TCs differ from that of not rapidly
changing TCs?”

Both $\\{\mathbf{S}_{<t}\\}_{t\geq 0}$ and $\\{Y_{t}\\}_{t\geq 0}$ are highly
dependent (auto-correlated) time series in this application: environmental
fields such as convection change slowly, and rapid intensity change events are
by definition extended periods of change (typically 12-48 hours), meaning that
successive measurements of these variables at short time intervals (e.g., 3
hours) will be highly dependent. It is particularly challenging to perform
traditional two-sample tests in this DID setting due to the combination of low
sample size (671 TCs between 2000-2020), high-dimensional image data, and
strong temporal dependence.

Contribution and Relevance. Our contribution is two-fold: (i) On the
methodology side, we present a new statistical framework for detecting
arbitrary distributional differences in a high-dimensional setting with
labeled sequence data. The proposed two-sample test is valid in a DID setting
and provides local diagnostics as to the type of sequences
$\mathbf{s}\in\mathcal{S}$ that contribute to rejection of $H_{0}$ in Equation
2. (ii) On the applied side, we utilize our proposed framework to identify and
describe patterns of convective evolution in TCs prior to the onset of rapid
intensity change.

Figure 1: Leveraging prediction tools for statistical inference and scientific
insight. In our application, we seek relationships between high-dimensional TC
observations and TC behavior. We pose this as a two-sample testing problem,
but the high dimensionality and sequential correlation in our data make this
hypothesis difficult to test. By rewriting the test in terms of a prediction
problem, we are able to leverage powerful prediction methods to infer modes of
structural evolution in TCs which indicate rapid intensity change.

Figure 1 shows a flow chart of our overall approach. As indicated by the
vertical blue arrow to the left, we seek quantitative conclusions regarding
how TC behavior relates to convective structure as observed by GOES imagery
and extracted functions. We cast the two-sample test as a prediction problem
(Figure 1, right, for “Statistical Methods”). This allows us to leverage
powerful prediction techniques, such as convolutional neural nets (CNNs), to
gain scientific insight from high-dimensional functional or video data without
a prior dimension reduction (arrow back to “TC behavior” and “Intensity
Guidance”). To account for dependence in the labeled sequence
$\\{(\mathbf{S}_{<t},Y_{t})\\}_{t\geq 0}$, we develop a bootstrap regression
test (Algorithm 1), which yields a valid p-value accompanied by interpretable
diagnostics.

Our bootstrap test is ideal for TC studies: the low number of unique TCs for
which high-resolution satellite imagery is available prohibits efficient
inference via classical blocking schemes, whereas our method can take
advantage of the extensive historical record of event (label) sequences.
Theorem 1 shows that the bootstrap test is valid as long as we can estimate
the distribution of the label sequences well; this is made easier by the fact
that labels are binary.

Section 7 provides evidence that deep and/or deepening convection is a
necessary precursor to RI, but that other factors (e.g., low vertical shear)
must be present. This elevated risk of RI due to deepening convection is often
present _prior_ to the onset of intensification, demonstrating value to
forecasting pipelines. RW, meanwhile, is a more variable process and did not
return significant results.

Relation to Other Work. Event impact and causal inference for time series. Our
problem set-up is closest in flavor to association and causal inference
studies for testing the relationship between a time and event series
(Scharwächter and Müller, 2020b; Luo et al., 2014). The vast majority of these
works assume univariate time series, or test for each dimension in a
multivariate time series separately (Candes et al., 2018). The recent paper by
Scharwächter and Müller (2020b) leverages a two-sample testing approach for
high-dimensional data as in this work, albeit to study how a discrete binary
event history impacts a multivariate time series, rather than how the
evolution of a complex time series is associated with a later event. A key
methodological difference is that their work handles the lack of independence
by sampling data points that are far from each other, while our method allows
for the use of all available data. This distinction is key for applications
with limited data. Their proposed multiple testing procedure also does not
control the false positive error rate exactly, but heuristically.

Two-sample tests in high dimension. Recently, there has been a growing
interest in nonparametric two-sample tests in high dimension. Popular machine
learning-based approaches include classification accuracy tests (Kim et al.,
2021), kernel-based tests (Gretton et al., 2012), and divergence-based density
ratio tests (Moon and Hero, 2014; Kandasamy et al., 2015). We use the same
regression test statistic (Equation 4) as Kim, Lee and Lei (2019) to allow for
interpretable local diagnostics. However, Kim, Lee and Lei (2019) and the
above-mentioned two-sample testing papers only handle the standard
independent, identically distributed (IID) data setting (Figure 2(a)), whereas
the methods in this work apply to DID sequence data (Figure 2(c)).

Also related to our paper are modern tests of conditional independence between
two random vectors $Y$ and $Z$ given a third random vector $X$ (see Discussion
in Section 8.2) and tests of the conditional mean and quantile dependence of
$Y$ on $X$. Most high-dimensional research on this front (including model-X
knockoffs; Candes et al. 2018; Sesia, Sabatti and Candès 2019; Berrett et al.
2020) assumes IID data $(X_{i1},\ldots,X_{ip},Y_{i})\sim F_{X,Y}$ for
$i=1,\ldots,n$, and also assumes that the distribution of $Y$ given $X$
depends on only a small fraction of the $p$ covariates. The latter sparsity
assumption is reasonable for, e.g., genome-wide association studies, but not
for remote sensing applications with image and functional data. Another key
difference between our testing approach and so-called model-X methodologies,
which also use machine learning algorithms to approximate the distribution of
$Y$ given $X$ (Katsevich and Ramdas, 2020), is that model-X approaches
estimate or make assumptions regarding the distribution of $X$ (or $X$ given
$Z$), whereas we instead estimate the distribution of the response $Y$.

Bootstrap for time series. There exist many different types of bootstrap
methods for dependent data; see, e.g, Bühlmann (2002); Horowitz (2003); Kreiss
and Paparoditis (2011) for a review. The goal is often to model the data
distribution for parameter estimation, rather than as here to test for an
association between a label series and a high-dimensional time series. Our
framework also does not bootstrap the entire distribution of
$\\{(\mathbf{S}_{<t},Y_{t})\\}_{t\geq 0}$ but only the distribution of labels
$\\{Y_{t}\\}_{t\geq 0}$. For binary labels, this is a much easier estimation
problem than, for example, block-bootstrap of high-dimensional time series.

TC analyses. Many TC analysis tools incorporate information from 2D fields via
scalar values such as area averages; for example, the operationally-used SHIPS
and SHIPS-RII forecast schemes include scalars for the fraction of pixels with
IR temperatures below -30∘C within the 50-200-km annulus (DeMaria and Kaplan,
1999; Kaplan et al., 2015). Such approaches discard complex, time-evolving
structure in the 2D fields. More recent analyses take spatial information into
account by applying dimension reduction techniques like functional principal
component analysis (PCA) to the field, such as for TC eye formation forecasts
in (Knaff and DeMaria, 2017). However, dimension reduction adds an extra layer
of abstraction between TC structure and subsequent TC behavior and can reduce
meteorologists’ ability to interpret the information. The ORB framework
(Organization, Radial structure, Bulk morphology) was introduced in McNeely et
al. (2020) to summarize convective structure into a dictionary of functional
features. The key objective was to utilize entire functions to quantify
structure rather than thresholded feature statistics, thereby enabling richer
descriptions of spatial structure while remaining interpretable. In this
paper, we leverage CNNs to model the relationship between TC intensity change
and the temporal evolution of one such continuous ORB function — the radial
profile of cloud-top temperatures as observed by GOES-IR imagery (Sanabia,
Barrett and Fine, 2014; McNeely et al., 2019, 2020). Our bootstrap test then
provides a powerful tool for directly assessing whether there is an
association between RI or RW events and TC structure or the TC environment in
terms of structural summaries (1D ORB functions) that are easily digestible to
meteorologists.

Outline. We begin by defining the problem set-up in Section 2, paying special
attention to different dependence structures in
$\\{(\mathbf{S}_{<t},Y_{t})\\}_{t\geq 0}$. In Section 3 we describe the TC
data. In Section 4 we lay out the details of our bootstrap test for
distributional differences in dependent sequence data. In Section 5 we provide
theoretical justification for validity of the bootstrap test. In Section 6 we
introduce a simulated toy example to empirically demonstrate the advantage of
a Markov Chain-based bootstrap test over traditional permutation testing. In
Section 7 we apply our method to study the evolution of convective structure
in TCs prior to a rapid intensity change. Finally, in Section 8 we discuss
potential extensions of our method.

##  2 Set-Up

Our goal is to detect distributional differences in labeled sequence data
$\\{(\mathbf{S}_{<t},Y_{t})\\}_{t\geq 0}$, where the “labels”
$Y_{t}\in\\{0,1\\}$ are binary, and the covariates
$\mathbf{S}_{<t}\in\mathcal{S}$ can represent high-dimensional quantities. We
formalize this question in the hypothesis

$\displaystyle H_{0}:p(\mathbf{s}_{<t}|Y_{t}=1)$
$\displaystyle=p(\mathbf{s}_{<t}|Y_{t}=0),\ \text{for all}\ t\ \text{and}\
\mathbf{s}_{<t},\ \text{versus}$ (1) $\displaystyle\
H_{1}:p(\mathbf{s}_{<t}|Y_{t}=1)$ $\displaystyle\neq
p(\mathbf{s}_{<t}|Y_{t}=0),\ \text{for some}\ t\ \text{and}\ \mathbf{s}_{<t}.$

We assume that the sequence $\\{(\mathbf{S}_{<t},Y_{t})\\}_{t\geq 0}$ is
stationary (Assumption 1), which allows us to rewrite the hypothesis as

$\displaystyle H_{0}:p(\mathbf{s}|Y=1)$ $\displaystyle=p(\mathbf{s}|Y=0),\
\text{for all}\ \mathbf{s}\in\mathcal{S},\ \text{versus}$ (2) $\displaystyle\
H_{1}:p(\mathbf{s}|Y=1)$ $\displaystyle\neq p(\mathbf{s}|Y=0),\ \text{for
some}\ \mathbf{s}\in\mathcal{S}.$

We will study three different settings (see Figure 2) which admit increasingly
complex dependence in $\\{\mathbf{S}_{<t}\\}_{t\geq 0}$ and
$\\{Y_{t}\\}_{t\geq 0}$. In Setting A (Figure 2(a)) there is no temporal
dependence, meaning the data $\\{(\mathbf{S}_{<t},Y_{t})\\}$ are IID. Testing
Equation 1 is still challenging in this setting when $\mathbf{S}_{<t}$ is
high-dimensional, but various methods and associated theory have been
developed to handle these challenges; see e.g. Kim, Lee and Lei (2019); Kim et
al. (2021). In Setting B (Figure 2(b)) there is temporal dependence in
$\\{\mathbf{S}_{<t}\\}$, but $\\{Y_{t}\\}$ are conditionally independent given
the associated value in $\\{\mathbf{S}_{<t}\\}$. In Setting C (Figure 2(c))
there is temporal dependence in both $\\{\mathbf{S}_{<t}\\}$ and
$\\{Y_{t}\\}$, regardless of the association between the two variables. We
expect the TC data to exhibit the structure of Setting C because intensity
change labels ($Y_{t}$) are not conditionally independent solely given the
convective structure of the storm ($\mathbf{S}_{<t}$). The effect of
convective activity on TC intensity change generally manifests within 24
hours, so a 24-hour history in $\mathbf{S}_{<t}$ should be sufficient (Rogers,
2010).

(a) Setting A: $\\{(\mathbf{S}_{<t},Y_{t})\\}_{t\geq 0}$ with no temporal
dependence between pairs $(\mathbf{S}_{<t},Y_{t})$ for different $t$.

(b) Setting B: $Y_{t}$ conditionally independent of $Y_{t-1}$ given
$\mathbf{S}_{<t}$; $\mathbf{S}_{<t}$ is autocorrelated.

(c) Setting C: $Y_{t}$ conditionally dependent on $Y_{t-1}$ given
$\mathbf{S}_{<t}$; $\mathbf{S}_{<t}$ and $Y_{t}$ are each autocorrelated.

Figure 2: Dependence settings. Directed acyclic graphs (DAGs) illustrating the
three dependence structures we explore. Note that each variable
$\mathbf{S}_{<t}$ can itself represent a temporal sequence of high-dimensional
functions or images, as in Figure 3.

##  3 Sequence Data from Tropical Cyclone Satellite Imagery

Figure 3: Evolution of TC convection as structural trajectories. The raw data for each trajectory $\mathbf{S}_{<t}$ is a sequence of a sequence of TC-centered cloud-top temperature images from GOES ($T_{b}$). We convert each GOES image into a radial profile ($X_{t}$). The 24-hour sequence of consecutive radial profiles, sampled every 30 minutes, defines a _structural trajectory_ or Hovmöller diagram ($\mathbf{S}_{<t}$). These trajectories serve as high-dimensional inputs to $\widehat{m}_{\text{post}}(\mathbf{s}_{<t})$. |  | NAL | ENP | Total | Year Range | Years
---|---|---|---|---|---|---
(i) | Training Data |  |  |  |  |
| All Sequences | 47,502 | 31,549 | 79,051 |  |
| RI Sequences | 7,015 | 6,742 | 13,757 |  |
| RW Sequences | 5,878 | 7,298 | 13,176 |  |
| Unique TCs | 209 | 185 | 394 | 2000-2012 | 13
(ii) | Test Data |  |  |  |  |
| All Sequences | 28,368 | 32,817 | 61,185 |  |
| RI Sequences | 3,965 | 6,386 | 10,351 |  |
| RW Sequences | 3,167 | 7,182 | 10,349 |  |
| Unique TCs | 125 | 152 | 277 | 2013-2020 | 8
(iii) | Synoptic Labels |  |  |  |  |
| All Labels | 14,683 | 15,274 | 29,957 |  |
| RI Labels | 1,850 | 2,462 | 4,312 |  |
| RW Labels | 1,643 | 2,534 | 4,177 |  |
| Unique TCs | 532 | 589 | 1,121 | 1979-2012 | 34
Table 1: Sample sizes: Data set summary for each category, (i) labeled
Sequences $(\mathbf{S}_{<t},Y_{t})$ used in training, (ii) unlabeled test
sequences $\mathbf{S}_{<t}$, and (iii) synoptic labels $Y_{t}$ used when
complete trajectories are not needed.

Analysis of TC convective structure relies on two types of observations:
sequences of longwave infrared imagery captured by GOES imagers and records of
TC intensity and location recorded in NOAA’s HURDAT2 database (Landsea and
Franklin, 2013).

Longwave infrared (IR) imagery ($\sim$10.3$\mu m$ wavelength) serves as a
proxy for convective strength: where IR-estimated cloud-top temperatures are
low, convection is strong. GOES longwave IR imagery is available through
NOAA’s MERGIR database (Janowiak, Joyce and Xie, 2020) at
30-minute$\times$4-km resolution over both the North Atlantic (NAL) and
Eastern North Pacific (ENP) basins from 2000–present. Every thirty minutes
during the lifetime of a storm, we download a $\sim$2,000 km$\times$2,000 km
“stamp” of IR imagery surrounding the TC location. Figure 3 (left) shows two
such stamps after an 800-km radius mask is applied.

TC location and intensity are given by the NHC’s HURDAT2 best track database,
which utilizes all available data on each TC (including data not available in
real time) to estimate critical characteristics over the lifetime of each TC.
HURDAT2 best tracks are provided at 6-hour (or synoptic) time resolution; we
linearly interpolate TC latitude and longitude between HURDAT2 data points to
estimate TC location at non-synoptic times. Since we are interested in the
behavior of mature TCs (as opposed to e.g. early development), we consider TC
genesis to be the first synoptic time at which intensity surpasses 35 kt and
lysis to be the last synoptic time at which intensity is at least 35 kt.

Structural trajectories via ORB. Direct utilization of IR imagery is
challenging in a meteorological setting, where interpretability is at a
premium. In previous works, we defined the ORB framework (Organization, Radial
structure, Bulk morphology) to summarize convective structure into a
dictionary of continuous functions instead of using traditional area-averaged
features common to TC analyses. This framework enables richer descriptions of
spatial structure without sacrificing interpretability (McNeely et al., 2019).
We have leveraged the ORB framework to analyze the evolution of TC convective
structure and demonstrated how projections of ORB functions onto a PCA basis
can be used to identify rapid intensification events (McNeely et al., 2020),
but we have not yet directly utilized temporally evolving, continuous
functions.

This work studies the temporal evolution of an entire ORB function; in this
case, the radial profile, defined as
$\overline{T}(r)=\frac{1}{2\pi}\int_{0}^{2\pi}T_{b}(r,\theta)d\theta$. The
radial profile $\overline{T}(r)$ captures the structure of cloud-top
temperatures $T_{b}$ as a function of radius $r$ from the TC center and serves
as an easily interpretable description of the depth and location of convection
near the TC core (Sanabia, Barrett and Fine, 2014; McNeely et al., 2020). The
radial profiles are computed at 5-km resolution from 0-400km ($d=80$) (Figure
3, center); we denote these summaries of convective structure at each time $t$
by $\mathbf{X}_{t}$. Finally, at each time $t$ we stack the preceding 24 hours
(48 profiles) into a _structural trajectory_
$\mathbf{S}_{<t}=\\{\mathbf{X}_{t},\mathbf{X}_{t-1},\dots,\mathbf{X}_{t-48}\\}$.
We visualize these trajectories with Hovmöller diagrams (Hovmöller 1949; see
Figure 3, right).

Labeling TC sequence data. HURDAT2 contains estimated TC intensities only at
synoptic times (0000 UTC, 0600 UTC, 1200 UTC, and 1800 UTC). We thus begin by
labeling these points $Y_{t}\in\\{0,1\\}$ based on whether the TC was
undergoing RI (or RW, for those analyses) at time $t$, where $Y=1$ indicates
occurrence of a rapid intensity change event. We then interpolate to non-
synoptic times by assigning label $Y_{t}=1$ if the an observation falls
between two consecutive synoptic $Y_{t}=1$ observations, and $Y_{t}=0$
otherwise. See McNeely et al. (2022a) (Section B) for further details on this
procedure.

There are three sample sizes of interest in this application: (i) the number
of labeled training sequences $\mathbf{S}_{<t}$ (further divided into 60%
train/40% validation), (ii) the number of test sequences $\mathbf{S}_{<t}$,
and (iii) the number of synoptic best track entries used when only labels
$Y_{t}$ are required (e.g., $\widehat{m}_{\text{seq}}$ in Algorithm 1(3)).
These sample sizes and associated years are summarized in Table 1.

##  4 Methods

Our TC problem set-up is difficult because of: (i) the complexity of the data
themselves, with one observation representing an entire sequence
$\mathbf{S}_{<t}$ of functions; (ii) dependence between labels $Y_{t}$ (and
sequences $\mathbf{S}_{<t}$) at nearby time points $t$; and finally (iii) the
need for scientific interpretability, or more precisely, statistical findings
which are easily digestible by TC scientists and forecasters.

We aim to test hypothesis 1 for DID sequence data
$\\{(\mathbf{S}_{<t},Y_{t})\\}_{t\geq 0}$ that satisfy the DAG of Setting C,
Figure 2(c). Our approach builds on Kim, Lee and Lei (2019), where the authors
present a regression approach for detecting differences in high-dimensional
IID data $\\{(\mathbf{S}_{i},Y_{i})\\}_{i=1}^{n}$, where
$\mathbf{S}_{i}\in\mathcal{S}:=\\{\mathbf{s}\in\mathbb{R}^{D}:p(\mathbf{s})>0\\}$,
for $Y_{i}\in\\{0,1\\}$. Their set-up is equivalent to Setting A, Figure 2(a).
The main idea is to rewrite the two-sample test in the equivalent formulation

$\displaystyle H_{0}:\mathbb{P}(Y=1|\mathbf{S}=\mathbf{s})$
$\displaystyle=\mathbb{P}(Y=1),\text{ for all }\mathbf{s}\in\mathcal{S},\text{
and }$ (3) $\displaystyle H_{1}:\mathbb{P}(Y=1|\mathbf{S}=\mathbf{s})$
$\displaystyle\neq\mathbb{P}(Y=1),\text{ for some }\mathbf{s}\in\mathcal{S}.$

These hypotheses involve a regression function for the “class posterior”
$m_{\text{post}}(\mathbf{s}):=\mathbb{P}(Y=1|\mathbf{S}=\mathbf{s})$, and a
“class prior” $m_{\text{prior}}:=\mathbb{P}(Y=1)$. One then tests $H_{0}$
against $H_{1}$ using the test statistic

$\lambda=\sum_{\mathbf{s}\in\mathcal{V}}\left(\widehat{m}_{\text{post}}(\mathbf{s})-\widehat{m}_{\text{prior}}\right)^{2},$
(4)

where $\widehat{m}_{\text{post}}(\mathbf{s})$ is an estimate of
$m_{\text{post}}(\mathbf{s})$,
$\widehat{m}_{\text{prior}}=\frac{1}{n}\sum_{i=1}^{n}I(Y_{i}=1)$ is the class
proportion of the training sample, and $\mathcal{V}\subset\mathcal{S}$ is a
fixed finite set of evaluation data points. Depending on the choice of
regression method, the regression test based on $\lambda$ can adapt to
challenging non-standard data, like images and sequences of images or
functions.

Because the null distribution of $\lambda$ is typically unknown, Kim, Lee and
Lei (2019) compute p-values based on $\lambda$ by using a permutation
procedure. The procedure relies on the exchangeability of the labels $Y$ under
$H_{0}$ in Equation 3. However, there are at least two types of dependence in
our data $\\{(\mathbf{S}_{<t},Y_{t})\\}_{t\geq 0}$ which violate the
assumption of exchangeability: (i) autocorrelation in $\\{Y_{t}\\}_{t\geq 0}$
which is inherent or governed by unobserved quantities as in Setting C, Figure
2(c), and (ii) the presence of an observed, correlated confounding sequence
$\\{\mathbf{Z}_{t}\\}_{t\geq 0}$ (discussed in Section 8.2). In either case,
the theoretical guarantees of a valid test in Kim, Lee and Lei (2019) no
longer hold.

###  4.1 Accounting for Dependence in $Y_{t}|\mathbf{S}_{<t}$

How do we handle dependence in the labeled sequence data
$\\{(\mathbf{S}_{<t},Y_{t})\\}_{t\geq 0}$? Permutation tests essentially model
the distribution of $\\{Y_{t}\\}_{t\geq 0}$ assuming IID labels. One way to
admit dependence in the relabeling procedure is to instead assume a Markov
property of order $k$ on $\\{Y_{t}\\}_{t\geq 0}$; that is, to assume that the
random variable $Y_{t}$ depends only on the previous $k$ variables
$Y_{t-1},\dots,Y_{t-k}$. To estimate the null distribution of the test
statistic $\lambda$ (3), we draw new labels from a Markov autoregressive
model:

$\widetilde{Y}_{t}\sim\text{Bernoulli}\left(\widehat{\mathbb{P}}(Y_{t}=1|Y_{t-1}=\widetilde{Y}_{t-1},\dots,Y_{t-k}=\widetilde{Y}_{t-k})\right).$
(5)

The marginal distribution of the labels, denoted by $m_{\text{seq}}$, can be
estimated from a holdout sample of observed data $\\{Y_{t}\\}_{t\geq 0}$ with
a variety of methods, including binary Markov chains and random forests. As we
shall see, as long as the marginal estimate of $\\{Y_{t}\\}_{t\geq 0}$
converges in distribution to the true data-generating process as the size of
the holdout sample increases, then the bootstrap test detailed in Algorithm 1
will be asymptotically valid. The result holds even if
$m_{\text{post}}(\mathbf{s})$ and $m_{\text{prior}}$ are not well-estimated.
This is good news for TC analysis, as one usually has ample access to label
series data $Y_{t}$, whereas sample sizes for the sequences $\mathbf{S}_{<t}$
derived from high-resolution satellite images are smaller.

Proof of the validity of our bootstrap test for Equation 3 is given in Section
5, Theorem 1. Section 6 includes empirical results on the power of the test
for synthetic data with the DAG structures in Figure 2.

type of test (BST=TRUE for bootstrap test; BST=FALSE for permutation test);
train data $\\{(\mathbf{S}_{<t},Y_{t})\\}_{t\in\mathcal{T}_{1}}$ and
regression method for estimating
$m_{\text{post}}(\mathbf{s}):=\mathbb{P}(Y_{t}=1|\mathbf{S}_{<t}=\mathbf{s})$;
for bootstrap test: train data $\\{Y_{t}\\}_{t\in\mathcal{T}_{2}}$ and
regression method for estimating
$m_{\text{seq}}(Y_{t-1},\dots,Y_{t-k}):=\mathbb{P}(Y_{t}=1|Y_{t-1},\dots,Y_{t-k})$;
number of repetitions $B$; evaluation points $\mathcal{V}$.

p-value for testing $H_{0}:p(\mathbf{s}|Y=1)=p(\mathbf{s}|Y=0)\ \text{for
all}\ s\in\mathcal{S}$; local posterior differences
$\\{\lambda(\mathbf{s})\\}_{\mathbf{s}\in\mathcal{V}}$ at the evaluation
points $\mathbf{s}\in\mathcal{S}$.

// Estimate underlying probability distributions

(1) Estimate $m_{\text{prior}}:=\mathbb{P}(Y_{t}=1)$ with class proportion
$\widehat{m}_{\text{prior}}=\frac{1}{\lvert\mathcal{T}_{1}\rvert}\sum_{t\in\mathcal{T}_{1}}Y_{t}$.

(2) Regress $Y_{t}$ on $\mathbf{S}_{<t}$ using $\mathcal{T}_{1}$ to compute
$\widehat{m}_{\text{post}}$.

(3) if  _BST_ then

Regress $Y_{t}$ on $Y_{t-1},\dots,Y_{t-k}$ using $\mathcal{T}_{2}$ to compute
$\widehat{m}_{\text{seq}}$.

// Compute test statistic and estimate its null distribution

(4) Compute test statistic
$\mathbf{\lambda}=\sum_{\mathbf{s}\in\mathcal{V}}\lambda^{2}(\mathbf{s}),\text{
where
}\lambda(\mathbf{s})=\widehat{m}_{\text{post}}(\mathbf{s})-\widehat{m}_{\text{prior}}$.

(5) for  _$b\in\\{1,2,\dots,B\\}$_ do

\- Draw new train labels $\\{\widetilde{Y}_{t}\\}_{t\in\mathcal{T}_{1}}$ under
$H_{0}$:

if  _BST_ then

Draw an initial label sequence $\widetilde{Y}_{1},\dots,\widetilde{Y}_{k}$
from the empirical distribution.

Draw sequence of length 100$\times k$ from
$\widetilde{Y}_{t}\sim\text{Binom}(\widehat{m}_{\text{seq}}(\widetilde{Y}_{t-k},\dots,\widetilde{Y}_{t-1}))$
for burn-in.

Draw new labels
$\widetilde{Y}_{t}\sim\text{Binom}(\widehat{m}_{\text{seq}}(\widetilde{Y}_{t-k},\dots,\widetilde{Y}_{t-1}))$
for $t\in\mathcal{T}_{1}$.

else

Permute original labels $\\{Y_{t}\\}_{t\in\mathcal{T}_{1}}$.

\- Regress $\widetilde{Y}_{t}$ on $\mathbf{S}_{<t}$ using $\mathcal{T}_{1}$ to
compute $\widehat{m}_{\text{post}}^{(b)}$.

\- Recompute test statistic
$\widetilde{\lambda}^{(b)}=\sum_{\mathbf{s}\in\mathcal{V}}\left(\widetilde{\lambda}^{(b)}(\mathbf{s})\right)^{2},\text{
where
}\widetilde{\lambda}^{(b)}(\mathbf{s})=\widehat{m}^{(b)}_{\text{post}}(\mathbf{s})-\widehat{m}_{\text{prior}}$.

// Compute approximate p-value

(6) Compute $p$-value according to
$\widehat{p}=\frac{1}{B+1}\left(1+\sum_{b=1}^{B}\mathbb{I}\left(\widetilde{\lambda}^{(b)}>\lambda\right)\right).$

return $\widehat{p},\ \\{\lambda(\mathbf{s})\\}_{\mathbf{s}\in\mathcal{V}}$

Algorithm 1 Test for distributional differences in labeled sequence data.

###  4.2 Local Diagnostics

Suppose $H_{0}$ is rejected. That is, we detect that the two distributions
$p(\mathbf{s}|Y=0)$ and $p(\mathbf{s}|Y=1)$ are indeed different. How do we
then provide the scientist with interpretable diagnostics that explain how the
two distributions are different?

Classification accuracy tests (Kim et al., 2021) require a separate post-hoc
procedure to identify local distributional differences (Gretton et al., 2012;
Chakravarti et al., 2021). A key advantage of the regression test is that the
test statistic in Equation 4, by construction, is a sum of local posterior
differences. Indeed, for each evaluation point $\mathbf{s}\in\mathcal{V}$, we
compute the local posterior difference (LPD)

$\lambda(\mathbf{s})=\widehat{m}_{\text{post}}(\mathbf{s})-\widehat{m}_{\text{prior}}.$
(6)

A large value of $\lambda(\mathbf{s})$ indicates that the distributions
$p(\mathbf{s}|Y=0)$ and $p(\mathbf{s}|Y=1)$ are very different at
$\mathbf{s}\in\mathcal{V}$, which in turn contributes to a larger test
statistic $\lambda=\sum_{\mathbf{s}\in\mathcal{V}}\lambda^{2}(\mathbf{s})$ and
potential rejection of $H_{0}$.

######  Remark 1.

The posterior difference
$\lambda(\mathbf{s})=\mathbb{P}(Y=1|\mathbf{s})-\mathbb{P}(Y=1)$ can be viewed
as a scaled density difference:

$\lambda(\mathbf{s})=\frac{p(\mathbf{s}|Y=1)-p(\mathbf{s}|Y=0)}{w(\mathbf{s})}$
(7)

where ${w(\mathbf{s})=\frac{1}{1-\pi}\ p(\mathbf{s}|Y=1)+\frac{1}{\pi}\
p(\mathbf{s}|Y=0)}$ is a positive scaling function, and
${\pi:=\mathbb{P}(Y=1)}$ denotes the prior class probability or
${m}_{\text{prior}}$. This difference has several desired properties for
assessing local distributional differences: $\lambda(\mathbf{s})$ is always
bounded, unlike other popular discrepancy measures such as the density ratio,
$p(\mathbf{s}|Y=1)/p(\mathbf{s}|Y=0)$, and the density difference,
$p(\mathbf{s}|Y=1)-p(\mathbf{s}|Y=0)$, itself. Furthermore, the posterior
difference does not decay to zero as fast as $p(\mathbf{s})$, which leads to
high sensitivity to detect differences in low density regions; for example, in
the case of balanced classes, $\lambda(\mathbf{s})$ takes a value of
$+\frac{1}{2}$ when $p(\mathbf{s}|Y=1)\gg p(\mathbf{s}|Y=0)$, and a value of
$-\frac{1}{2}$ when $p(\mathbf{s}|Y=1)\ll p(\mathbf{s}|Y=0)$, regardless of
the actual magnitudes of $p(\mathbf{s})$, $p(\mathbf{s}|Y=1)$, and
$p(\mathbf{s}|Y=0)$.

In summary, our method tests for distributional differences in labeled
sequence data as follows:

  1. 1.

Decide on a suitable model for the marginal distribution of
$\\{Y_{t}\\}_{t\geq 0}$.

  2. 2.

Apply Algorithm 1 to compute the p-value for testing the hypotheses in
Equation 3.

  3. 3.

If $H_{0}$ is rejected, then examine local posterior differences (6) to
identify what patterns $\mathbf{s}$ in the state space $\mathcal{S}$ of
sequence data contributed the most to the rejection.

##  5 Theory

This section provides theoretical justification for our bootstrap procedure
for testing Equation 3. In particular, we show that Algorithm 1 controls the
type I error asymptotically.

######  Assumption 1 (Stationary sequence).

$\\{(\mathbf{S}_{<t},Y_{t})\\}_{t\geq 0}$ is a stationary sequence, where
$\mathbf{S}_{<t}\in\mathcal{S}$ and $Y_{t}\in\\{0,1\\}$

Assumption 1 is needed for the hypothesis in Equation 2 to be well-defined.

######  Assumption 2 (Conditional independence).

$\\{(\mathbf{S}_{<t},Y_{t})\\}_{t\geq 0}$ satisfies the DAG of Setting C
(Figure 2).

Assumption 2 encodes the conditional independences required for our method to
control type I error.

In this section, we denote the test statistic by

$\displaystyle\lambda(\mathcal{D})=\int(\widehat{m}_{\text{post}}(\mathbf{s})-\widehat{m}_{\text{prior}})^{2}dQ(\mathbf{s}),$
(8)

where $Q$ is any fixed measure over $\mathcal{S}$, and
$\widehat{m}_{\text{post}}$ and $\widehat{m}_{\text{prior}}$ are obtained
using a training set
$\mathcal{D}:=\\{(\mathbf{S}_{<t},Y_{t})\\}_{t\in\mathcal{T}_{1}}$. In
Algorithm 1, $Q$ is the distribution that assigns mass $1/|\mathcal{V}|$ for
each evaluation point $\mathbf{s}\in\mathcal{V}$, but the results we show here
apply to any $Q$. We also assume that the regression estimator we use is a
continuous function of the training set.

######  Assumption 3 (Continuous regression method).

$\widehat{m}_{\text{post}}$ is obtained by applying a regression estimator
that is a continuous function of $\mathcal{D}$.

Moreover, let
$\mathcal{D}^{t_{2}}_{0}:=\\{(\mathbf{S}_{<t},Y^{0}_{t})\\}_{t\in\mathcal{T}_{\text{1}}}$
denote a random draw from the data set used to estimate the null distribution
of $\lambda$, where
$\\{Y^{0}_{t}\\}_{t\in\mathcal{T}_{\text{1}}}\sim\text{G}_{\widehat{\mathbf{p}}_{t_{2}}}$
and G is a distribution over $\\{0,1\\}^{|\mathcal{T}_{\text{1}}|}$ indexed by
the parameter $\widehat{\mathbf{p}}_{t_{2}}$, which is estimated using a
holdout set $\mathcal{D}^{\prime}=\\{Y_{t}\\}_{t\in\mathcal{T}_{\text{2}}}$,
with $t_{2}=|\mathcal{T}_{\text{2}}|$. In the method described in Section 4.1,
$\text{G}_{\widehat{\mathbf{p}}_{t_{2}}}$ is the Markov autoregressive model
$\widehat{m}_{\text{seq}}$, but this model can be more general. We require it
to converge to the true distribution of the marginal process
$\\{Y_{t}\\}_{t\in\mathcal{T}_{1}}$ when the null hypothesis holds:

######  Assumption 4 (Consistency of the marginal distribution estimator).

The estimator $\widehat{\mathbf{p}}_{t_{2}}$ is such that, if the null
hypothesis is true,

$\text{G}_{\widehat{\mathbf{p}}_{t_{2}}}\xrightarrow[t_{2}\longrightarrow\infty]{\enskip\mbox{Dist}\enskip}\text{G}^{*},$

where $\text{G}^{*}$ is the true generating process of
$\\{Y_{t}\\}_{t\in\mathcal{T}_{1}}$.

In the following, we show two examples where Assumption 4 holds.

######  Example 1.

Under Settings A and B (Figure 2), $Y_{t}$’s are IID under the null
hypothesis. Thus, $G^{*}$ is necessarily a product of IID Bernoulli random
variables with some parameter $p$. Now, let $G_{\widehat{\mathbf{p}}_{t_{2}}}$
be the product of IID Bernoulli random variables with parameter given by
$p_{t_{2}}:=(t_{2})^{-1}\sum_{t\in\mathcal{T}_{2}}Y_{t}$. The law of large
numbers implies that
$p_{t_{2}}\xrightarrow[t_{2}\longrightarrow\infty]{\enskip\mbox{a.s.}\enskip}p$.
Thus, the cumulative distribution function of $Y_{t}^{0}$, given by

$F_{Y_{t}^{0}}(y_{t})=\begin{cases}0,&\text{if $y_{t}<0$}\\\
1-p_{t_{2}},&\text{if $0\leq y_{t}<1$}\\\ 1&\text{otherwise}\end{cases}$

is such that
$F_{Y_{t}^{0}}(y_{t})\xrightarrow[t_{2}\longrightarrow\infty]{\enskip\enskip}F_{Y_{t}}(y_{t})$.
It follows that

$\displaystyle\mathbb{P}\left(Y^{0}_{t}\leq y_{t},\ \forall
t\in\mathcal{T}_{1}\right)$
$\displaystyle=\Pi_{t\in\mathcal{T}_{1}}F_{Y_{t}^{0}}(y_{t})$
$\displaystyle\xrightarrow[t_{2}\longrightarrow\infty]{\enskip\enskip}\Pi_{t\in\mathcal{T}_{1}}F_{Y_{t}}(y_{t})=\mathbb{P}\left(Y_{t}\leq
y_{t},\ \forall t\in\mathcal{T}_{1}\right)$

and therefore Assumption 4 holds.

######  Example 2 (Markov Chain).

If (under $H_{0}$) the process $\\{Y_{t}\\}_{t\geq 0}$ is an irreducible and
ergodic stationary $k$-order Markov Chain, then the maximum likelihood
estimators of the transition probabilities converge almost surely to the true
transition probabilities (Grimmett and Stirzaker, 2020). The same reasoning of
Example 1 then implies that Assumption 4 holds for such estimator under
Setting C.

The following theorem shows that, under $H_{0}$, the test statistic has
approximately the same distribution as the test statistic evaluated at the
generated data $\mathcal{D}^{t_{2}}_{0}$.

######  Theorem 1.

Assume 1, 2, 3 and 4. Under the null hypothesis,

$\lambda(\mathcal{D}^{t_{2}}_{0})\xrightarrow[t_{2}\longrightarrow\infty]{\enskip\mbox{Dist}\enskip}\lambda(\mathcal{D})$

It follows from Theorem 1 that type I error is controlled asymptotically.

######  Corollary 1 (Type I error control).

Let

$\widehat{p}^{\
t_{2}}_{B}(\mathcal{D}):=\frac{1}{B+1}\left(1+\sum_{b=1}^{B}\mathbb{I}\left(\lambda(\mathcal{D}^{(b)})>\lambda(\mathcal{D})\right)\right)$

be the Monte Carlo p-value for $H_{0}$, where
$\mathcal{D}^{(1)},\ldots,\mathcal{D}^{(B)}\stackrel{{\scriptstyle\text{{IID}}}}{{\sim}}\mathcal{D}^{t_{2}}_{0}$.
Assume that Assumptions 1, 2, 3 and 4 hold. Then, under the null hypothesis,
for any $0<\alpha<1$,

$\lim_{t_{2}\longrightarrow\infty}\lim_{B\longrightarrow\infty}\mathbb{P}\left(\widehat{p}^{\
t_{2}}_{B}(\mathcal{D})\leq\alpha\right)=\alpha.$

See McNeely et al. (2022a) (Section C) for proofs.

##  6 Performance of Tests on Synthetic Data

In this section, we use synthetic data to examine the performance (validity,
power and diagnostics) of our Markov chain bootstrap test for the data
dependence settings in Figure 2(c).

###  6.1 Synthetic Univariate Sequence Data

For simplicity, we first consider a scalar covariate $S_{t}$ of interest. We
then create dependent sequences $\\{S_{t},Y_{t}\\}_{t\geq 0}$ with a logistic
generative model:

$\displaystyle Y_{t}|S_{t}\sim$ $\displaystyle\,\text{Bernoulli}(p_{t})$
$\displaystyle p_{t}=\,$ $\displaystyle\text{logistic}(\gamma
H_{\delta}(S_{t})+U_{t})$ $\displaystyle H_{\delta}(S)=$
$\displaystyle\begin{cases}0&\lvert S\rvert<\delta\\\ S&\lvert
S\rvert\geq\delta\end{cases}$ (9) $\displaystyle S_{t}=$ $\displaystyle
U_{t}^{\prime},\ \ U_{t}^{\prime}\sim\text{AR}_{\phi^{\prime}}(1)$
$\displaystyle U_{t}\sim$ $\displaystyle\text{AR}_{\phi}(1).$

The variables $U_{t}$ and $U^{\prime}_{t}$ are spurious variables (not
included in the DAGs), which induce autocorrelation in the binary response
variable $Y_{t}$ and the covariate $S_{t}$, respectively. In our toy example,
we assume that both $U_{t}$ and $U^{\prime}_{t}$ are given by autoregressive
models of order $1$; more specifically, by AR(1) models of the form
$U_{t}=\phi U_{t-1}+\sqrt{1-\phi^{2}}\epsilon_{t}$ and
$U^{\prime}_{t}=\phi^{\prime}U^{\prime}_{t-1}+\sqrt{1-\phi^{\prime
2}}\epsilon^{\prime}_{t}$, where
$\epsilon_{t},\epsilon^{\prime}_{t}\overset{iid}{\sim}N(0,1)$, and
$\phi,\phi^{\prime}{\in[0,1]}$ are parameters for the 1-lag autocorrelation in
$U_{t}$ and $U^{\prime}_{t}$, respectively. Increasing $\phi^{\prime}$ thus
increases the autocorrelation (but not the variance) of the variable of
interest $S_{t}$, while increasing $\phi^{\prime}$ increases the
autocorrelation (but not the variance) of the spurious variable $U_{t}$.

The parameter $\gamma\geq 0$ determines the signal strength, or the strength
of the dependence of $Y_{t}$ on $S_{t}$. Testing $H_{0}$ in Equation 1 is
equivalent to testing $H_{0}:\gamma=0$. Ideally, our method should also be
able to identify local regions in the sample space $\mathcal{S}$ where the two
distributions are different or the same. To assess our method’s local
performance, we hence include a hard thresholding operator $H_{\delta}(\cdot)$
in Equation 9, which, regardless of the signal strength $\gamma$, creates a
region in $\mathcal{S}=\mathbb{R}$ where $\mathbb{P}(Y=1|S=s)=\mathbb{P}(Y=1)$
for $s\in(-\delta,\delta)$.

The two parameters $\phi$ and $\phi^{\prime}$ allow us to create synthetic
data with the dependence structures in Figure 2. More specifically,
Setting A: $\phi^{\prime}=\phi=0$
Setting B: $\phi^{\prime}>0$ induces autocorrelation in $\\{S_{t}\\}_{t\geq
0}$ via $U^{\prime}_{t}$; $\phi=0$
Setting C: $\phi^{\prime}>0$ induces autocorrelation in $\\{S_{t}\\}_{t\geq
0}$ via $U^{\prime}_{t}$, while $\phi>0$ induces autocorrelation in
$\\{Y_{t}\\}_{t\geq 0}$ via $U_{t}$.

######  Remark 2.

In this toy example, $S_{t}=U^{\prime}_{t}$ with $U^{\prime}_{t}$ observed.
More generally, $\mathbf{S}_{<t}$ can depend on unmeasured variables
$\mathbf{U}^{\prime}_{t}$, as well as confounding variables $\mathbf{Z}_{t}$,
which are connected to both $\mathbf{S}_{<t}$ and $Y_{t}$. The latter setting
is discussed in Section 8.2. We also note that the spurious variables
$\mathbf{U}_{t}$ and $\mathbf{U}^{\prime}_{t}$ can have more complex temporal
dependence (than in the example), as indicated by the fully connected
sequences in Figure 2(c).

###  6.2 Test Results for Synthetic Data

The logistic generative model in Equation 9 provides a variety of controls
over the dependence structure of $\\{(S_{t},Y_{t})\\}_{t\geq 0}$. For our
synthetic experiments, we implement Algorithm 1 with either the MC bootstrap
or the permutation test for $\lvert\mathcal{V}\rvert=250$ evaluation points.
We estimate the regression function
$m_{\text{post}}(s)=\mathbb{P}(Y_{t}=1|{S}_{t}={s})$ using train data
$\\{(S_{t},Y_{t})\\}_{t\in\mathcal{T}_{1}}$ and a Nadaraya-Watson (NW) kernel
estimator with an Epanechnikov kernel and the bandwidth chosen as the sample
standard deviation of $s$ divided by $\lvert\mathcal{T}_{1}\rvert^{1/5}$ (Li
and Racine, 2007). For the bootstrap test, we estimate the label distribution
$m_{\text{seq}}$ using an order $k=4$ Markov chain and labels
$\\{Y_{t}\\}_{t\in\mathcal{T}_{2}}$.

Validity: Testing $H_{0}$ in Equation 1 is equivalent to testing
$H_{0}:\gamma=0$ (no signal strength). To examine validity, we set $\gamma=0$
and simulate 500 independent data sets for each experiment, or combination of
“setting” (A, B, C) and “test method” (permutation or MC bootstrap test with
$\lvert\mathcal{T}_{1}\rvert=\lvert\mathcal{T}_{2}\rvert=\lvert\mathcal{V}\rvert=250$).
That is, each experiment returns 500 p-values. If the test controls type I
error, we expect these p-values to be approximately uniformly distributed.
Figure 4 assesses validity by plotting the difference between the empirical
and (uniform) theoretical quantiles against the theoretical quantiles; this is
equivalent to a standard quantile-quantile plot with the diagonal subtracted.
As a baseline, we provide a 95% confidence interval of this difference based
on $10,000$ Monte Carlo simulations of 500 uniformly distributed random
variables.

The permutation test (left panel) is valid under Settings A and B, where
$Y_{t}$ and $Y_{t-1}$ are independent after conditioning on $S_{t}$. However,
under Setting C, the p-values tend to have lower values than a uniform
distribution, corresponding to higher-than-nominal type I errors at most
significance levels. The MC bootstrap test (right panel) controls the type I
error at all significance levels for all three dependence settings, indicating
that our adjustment to account for the dependence in $Y_{t}\mid S_{t}$
achieved the desired result.

Figure 4: Synthetic sequence data: Validity of permutation test (left) versus
Markov Chain bootstrap test (right). Under $H_{0}:\gamma=0$, a valid test is
expected to return uniformly distributed p-values. Each curve corresponds to a
different experimental setting (A, B, or C), and shows the difference between
empirical and uniform theoretical quantiles for 500 repetitions; see text for
details. The gray region represents a 95% pointwise confidence interval
derived from Monte Carlo samples of 500 uniform deviates. Under Setting C
(labels dependent even after conditioning on predictors; purple curve), the
permutation test (left panel) does not control the type I error, but the
Markov Chain bootstrap test (right panel) does. (Setting A:
$\phi=\phi^{\prime}=0$. Setting B: $\phi=0,\ \phi^{\prime}=0.8$. Setting C:
$\phi=\phi^{\prime}=0.8$. Sample sizes
$\lvert\mathcal{T}_{1}\rvert=\lvert\mathcal{T}_{2}\rvert=\lvert\mathcal{V}\rvert=250$.)

Power: We next examine how the power of the test $H_{0}:\gamma=0$ versus
$H_{1}:\gamma\neq 0$ depends on

  1. i.

signal strength $\gamma$ (Figure 5, top),

  2. ii.

train sample size $\lvert\mathcal{T}_{1}\rvert$ (Figure 5, bottom),

  3. iii.

autocorrelation $\phi$ in labels $\\{Y_{t}\\}_{t\geq 0}$, or equivalently,
correlation in $\\{Y_{t}\rvert S_{t}\\}_{t\geq 0}$ (Figure 6, left), and

  4. iv.

autocorrelation $\phi^{\prime}$ in predictors $\\{S_{t}\\}_{t\geq 0}$ (Figure
6, right).

At each fixed value of $\gamma$, we perform 1000 simulations. Power is then
estimated as the fraction of rejected null hypotheses at the $\alpha=0.05$
level. To ensure validity, we choose $\lvert\mathcal{T}_{1}\rvert\geq 250$ as
before for the MC bootstrap test. (The permutation test is valid by
construction under Settings A and B, but not Setting C.)

Figure 5:  Synthetic sequence data: Power as a function of signal strength and
training sample size. Top: The power of all tests increases with the signal
strength $\gamma$, regardless of dependence setting. The MC bootstrap test has
similar power as the permutation test, but the former test can be applied to
the more challenging Setting C with dependent labels $Y_{t}|S_{t}$. Sample
sizes $|\mathcal{T}_{1}|=|\mathcal{T}_{2}|=|\mathcal{V}|=250$. Bottom: The
power of all tests, at the alternative $\gamma=0.5$, increases with the train
sample size $|\mathcal{T}_{1}|$, regardless of dependence setting. Sample
sizes $|\mathcal{T}_{2}|=|\mathcal{V}|=250$. The filled regions represent 95%
pointwise confidence intervals for binomial proportions. (Setting A:
$\phi=\phi^{\prime}=0$. Setting B: $\phi=0,\ \phi^{\prime}=0.8$. Setting C:
$\phi=\phi^{\prime}=0.8$.) Setting C is not shown for the permutation test, as
it is not valid. Figure 6: Synthetic sequence data: Power of bootstrap test as
function of autocorrelation of $Y$ (left) and $S$ (right). The red diamond-
shaped markers correspond to Setting A with IID data. The power to detect the
signal $\gamma=0.5$ is independent of the correlation in $Y_{t}|S_{t}$ (value
of $\phi$), but decreases with correlation in $S_{t}$ (larger values of
$\phi^{\prime}$). The filled regions represent 95% pointwise confidence
intervals for binomial proportions. (Sample sizes
$|\mathcal{T}_{1}|=|\mathcal{T}_{2}|=|\mathcal{V}|=250$.)

As expected, the power increases with the signal strength $\gamma$ for all
tests and dependence settings (Figure 5, top). When both tests are valid, the
MC test has the same power as the permutation test. The practical implication
is that, even if one _thinks_ Setting B is a good approximation to the problem
at hand, there are benefits to applying the MC bootstrap test: one can achieve
similar power with the advantage of having robustness in the event that the
labels are dependent after conditioning on predictors.

Figure 5 (bottom) indicates that the power of the tests may be determined by
the quality of the regression estimator
$\widehat{m}_{\text{post}}(\mathbf{s})$: indeed, Figure 5 (bottom) shows that
the power at a fixed alternative ($\gamma=0.5$) increases with the train
sample size $\lvert\mathcal{T}_{1}\lvert$. The latter result is consistent
with theorem 3.3 of Kim, Lee and Lei (2019), which states for a regression
permutation test under Setting A, that if the chosen regression method
$\widehat{m}_{\text{post}}(\mathbf{s})$ has a small mean integrated squared
error, then the power of testing (2) is large over a wide region of
alternative hypotheses.

Figure 7: Synthetic sequence data: Local posterior difference in null region.
The gray curves in each panel/setting show the true LPD, which is zero in the
local null (no signal) region $s\in(-0.25,0.25)$. The solid and dashed colored
curves represent the mean and one standard deviation estimates of the LPD over
200 simulated data sets. These estimates are on average close to the true LPD,
with the dispersion decreasing for a local nonparametric estimator (like the
NW kernel estimator) as the train sample size increases. The gray bar at the
bottom marks the standard deviation of $s$; the large variance of the
estimated LPDs far from $s=0$ is partially due to the concentration of data
near $s=0$. (Setting A: $\phi=\phi^{\prime}=0$. Setting B: $\phi=0,\
\phi^{\prime}=0.8$. Setting C: $\phi=\phi^{\prime}=0.8$. Sample sizes
$\lvert\mathcal{T}_{1}\rvert$ vary,
$\lvert\mathcal{T}_{2}\rvert=\lvert\mathcal{V}\rvert=250$.)

Figure 6 brings insight on how dependence in $\\{(S_{t},Y_{t})\\}_{t\geq 0}$
affect the power of the MC bootstrap test. The red diamond-shaped markers
represent Setting A with IID sequence data ($\phi=\phi^{\prime}=0$).
Increasing correlation in the labels $Y_{t}|S_{t}$ (larger values of $\phi$)
has no effect on power, while the test remains valid as long as
$\widehat{m}_{\text{seq}}$ is accurate; this further emphasizes the previous
result that the bootstrap test is robust to correlation in $Y_{t}|S_{t}$
without sacrificing power. Meanwhile, increasing correlation in $S_{t}$
(larger values of $\phi^{\prime}$) reduces power; this follows from a reduced
effective sample size which in turn reduces the quality of
$\widehat{m}_{\text{post}}$.

Local Posterior Differences: For our synthetic example, the hard thresholding
operator $H_{\delta}(\cdot)$ induces a region $s\in(-\delta,+\delta)$ where
$p(s|Y=1)=p(s|Y=0)$. If the null hypothesis in Equation 2 is rejected, then
the estimated LPDs can identify the regions of large versus small
distributional differences, as long as the regression estimator
$\widehat{m}_{\text{post}}$ is consistent and the train sample size $|T_{1}|$
is sufficiently large. Figure 7 shows the average and one standard deviation
estimates of the LPD for a NW kernel estimator over 200 simulations.

##  7 Relating Evolution of TC Convection to Rapid Intensity Change

In our TC study, each observation consists of (i) a 24-hour sequence
$\mathbf{S}_{<t}=\\{\mathbf{X}_{t},\mathbf{X}_{t-1},\dots,\mathbf{X}_{t-48}\\}$
of one-dimensional radial profile functions
$\mathbf{X}_{t}=\frac{1}{2\pi}\int_{0}^{2\pi}T_{b,t}(r,\theta)d\theta$,
sampled every 30 minutes for a total of 48 profiles (see Figure 3), and (ii) a
binary label $Y_{t}\in\\{0,1\\}$ for the entire sequence222We check that the
dominant principal components of $\mathbf{X}_{t}$ and the continuous
intensities used to derive $Y_{t}$ are stationary via augmented Dickey-Fuller
tests; the p-value of each test (including tests of the first 3 ORB
coefficients in McNeely et al. 2020) are all $<10^{-20}$. We conclude that
Assumption 1 is reasonable for these data.. Since all individual sample points
that are part of a rapid intensity change event are labeled as $Y=1$, a
24-hour sequence $\mathbf{S}_{<t}$ with sequence label $Y_{t}=1$ could either
be part of an ongoing RI/RW event (if $Y_{t}=1$ falls near the end of an
event), or be part of the lead-up to RI/RW (if $Y_{t}=1$ falls near the
beginning of an event). Analyses of the latter case — such as approaches that
can identify archetypal modes of structural evolution preceding the onset of
RI/RW — are particularly valuable to forecasting of RI/RW events.

We divide our TC study into three parts:

  1. 1.

Analysis by event type (RI versus not RI, or RW versus not RW) within each
basin, North Atlantic (NAL) or eastern North Pacific (ENP), for a total of
four different two-sample tests.

  2. 2.

Case studies of three tropical cyclones (Hurricanes Nicole, Jose, and Rosa).

  3. 3.

Analysis of a subset of our data that consists of sequences immediately
preceding RI onset. (We ask whether our regression two-sample test can find
archetypical evolutionary modes preceding RI onset, and if so, what the “lead
time” between typical patterns and the RI onset would be.)

As in the synthetic example, we estimate $m_{\text{seq}}$ in Algorithm 1 with
a Markov chain of order $k=8$. The sequence data $\mathbf{S}_{<t}$ are however
much more complex than in our synthetic example. This is where we benefit from
a more complex regression method for estimating $m_{\text{post}}$; here we fit
a convolutional neural network to the 24-hour sequence data. Further details
on how we estimate $m_{\text{seq}}$ and $m_{\text{post}}$ can be found in
McNeely et al. (2022a) (Section A).

###  7.1 Analysis by Event Type and Basin

Figure 8: Analysis by event type and basin. The MC bootstrap test rejects
$H_{0}$ (that is, it detects significant differences in convection) for the
RI-NAL, RI-ENP, and RW-ENP models; the test is not rejected for the RW-NAL
model. Analyses of local posterior difference (LPD) are shown for the first
three models (see panels i-iii). Left column: Two-dimensional PCA map of
sequence data. One point in the map represents a 24-hour structural trajectory
(sequence of radial profiles) $\mathbf{S}_{<t}$ with the color coding for the
estimated LPD. Right column: Six 24-hour structural trajectories at locations
A \-  F in the PCA map, shown as Hovmöller diagrams (recall Figure 3). The
examples for each study are selected at random from each combination of LPD
sign (positive LPD:  A \-  C ; negative LPD:  D \-  F ) and TC intensity
(Tropical Storm:  A ,  D ; Category 1-2:  B ,  E ; Category 3-5:  C ,  F ).
See text for a discussion of the results.

Significance test. We start by testing
$H_{0}:{p(\mathbf{s}_{t}|Y_{t}=1)}={p(\mathbf{s}_{t}|Y_{t}=0)}$ by event type
and basin. For rapid intensification (RI), the MC bootstrap test rejects
$H_{0}$ at level $0.05$ for both the NAL and ENP basins, meaning that we
indeed detect a significant difference (p$<$0.01) in 24-hour sequences
$\mathbf{S}_{<t}$ of convective structure leading up to RI versus not-RI
events. For rapid weakening (RW), the MC bootstrap test rejects $H_{0}$ in the
ENP basin, but not in the NAL basin.

Our results are consistent with scientists’ understanding of TCs; for rapid
intensification, a TC exhibits a narrow range of convective patterns “primed”
to efficiently convert heat energy to mechanical energy across the storm,
hence the structural difference in convection for RI versus not-RI events.
Rapid weakening, on the other hand, is a more complex process driven by
several factors external to the TC, such as vertical wind shear, which may not
be fully captured by convective structure. In addition, RW is expected to be
more difficult to detect in the NAL basin due to the broader range of possible
environmental configurations and the increased rarity of over-water RW in the
basin.

Local posterior difference. Next we investigate what kind of structural
patterns lead to the rejection of $H_{0}$ for the RI-NAL, RI-ENP, and RW-ENP
models. Figure 8 (left) shows a two-dimensional embedding of the sequence data
via principal component analysis (PCA; computed separately for each basin).
Each point represents a 24-hour sequence $\mathbf{S}_{<t}$ colored by its
local posterior different (LPD). Note that PCA is only used for purposes of
visualization; the test itself is performed on the entire sequence of radial
profiles without a prior dimension reduction step. Figure 8 (right) shows
examples of Hovmöller diagrams for six 24-hour sequences $\mathbf{S}_{<t}$
sorted by LPD and TC intensity.

TCs are known to have different distributions of $\mathbf{S}_{<t}$ for
different basins. Nevertheless, we identify the same type of evolutionary
patterns of convective structure for RI-NAL (panel i) and RI-ENP (panel ii):
Positive LPD or “high chance of RI” (see diagrams  A \-  C for i and ii) tends
to occur for cold cloud tops near the core (dark blue at smaller radius),
growing in coverage and depth of convection with time (dark blue region
extending to larger radii when going from -24 to 0 hours). Meanwhile, negative
LPD or “low chance of RI” (see diagrams  D \-  F for i and ii) tends to occur
when TCs already possess a well-defined eye (narrow yellow region near the
center) or exhibit decaying core convection (dark blue region decreasing in
size when going from -24 to 0 hours). While such patterns can be directly
quantified and studied in future works, this work remains focused on
exploration of entire radial profiles.

Finally, RW-ENP results (panel iii) are not exactly opposite of the RI-ENP
results (panel ii), meaning that “high chance of RW” patterns might not mirror
“low chance of RI” patterns, and vice versa. In particular, TCs commonly form
strong convective cores without eyes (iii-  E ) prior to intensification, form
an eye (end of i-  C ) during RI, then rapidly weaken by dissipating entirely
(ii-  D , iii-  A , iii-  B ) without reforming a cold, eyeless core because
the reduction of intensity is accompanied by a collapse of convection
throughout the TC. Thus, RI/RW-ENP results are not symmetric. Unfortunately,
the RW-ENP model also predominantly captures the trivial result that
currently-intense TCs are more likely to weaken; see Section 8.2 for a
discussion of potential corrections.

###  7.2 Hurricane Case Studies

Figure 9: Case studies of Hurricanes Nicole [NAL 2016, i], Jose [NAL 2017,
ii], and Rosa [ENP 2018, iii]. Left: Structural trajectories of each storm
through its entire lifetime, with the sidebar to the right showing the LPD at
each 30-minute observation; each LPD value is evaluated using the preceding
24-hour sequence for an RI model (trained on the train sample for that basin),
with positive values in red and negative values in blue. Right: Storm
intensity over time, where filled versus empty circles mark RI and RW events,
respectively. The physical track of the storm across the North Atlantic is
shown inset, with the same time points labeled.
_Panel (i)_ : For Hurricane Nicole, the LPDs capture rapidly intensifying
periods (  A ,  C ,  D ), collapsing convection  B , and the decay of the TC
eye  E .
_Panel (ii)_ : Hurricane Jose was subjected to high vertical wind shear near
Sept. 9, which our model does not account for; while the core convection of
the TC remained poised for RI, high shear instead caused rapid weakening.
_Panel (iii)_ : Hurricane Rosa exhibited two interesting phenomena captured by
the LPDs. First, it experienced a pause in its rapid intensification  A ; the
LPDs indicate an ongoing RI threat at this time, consistent with the
resumption of intensification 18 hours later. Second, beginning on Sept. 28,
the TC underwent an eyewall replacement cycle, associated with weakening. The
LPDs mark this shift at the TC’s peak intensity, as well as the brief period
(following  C ) where eyewall replacement is completed prior to landfall.

Thus far, we have analyzed the collection of 24-hour sequences in the
2013-2020 test sample as a whole; however, forecasters monitor individual
storms in real time for signals of RI. Here, we take an in-depth look at our
results for three individual TCs in the test set: Hurricanes Nicole, Jose, and
Rosa. Each of these storms display distinct evolutionary modes. We track each
TC through its lifetime to investigate the relationship between the evolution
of convective structure, the LPDs, and intensity change.

Figure 9, left, shows the evolution of radial profiles for each storm.
Appended to the right of these profiles is a sidebar showing the LPDs from the
RI model for the associated basin. The right panels display the evolution of
each storm’s intensity, where filled and hollow markers indicate RI and RW
events, respectively. The storm’s physical track over the ocean is shown as an
inset. These three hurricane case studies respectively highlight: (i) signals
which lead RI, (ii) the effect of vertical wind shear, and (iii) the
appearance of an eyewall replacement cycle (a process by which a second
eyewall forms, robbing the TC of energy as it shrinks to replace the original
eyewall).

Panel (i) depicts Hurricane Nicole [2016]. The TC underwent RI  A on Oct. 6
before its intensity stalled while influenced by the outflow from Hurricane
Matthew [2016], which induced vertical wind shear  B . Several days later,
Nicole re-intensified  C ,  D before again weakening and then transitioning
into an extratropical system  E . The LPDs in Figure 9 appear to align with
intensity changes and the emergence of deep convection  A ,  C and eye
formation  D . The structural trajectories immediately preceding  A and  C are
particularly interesting: the TC has not yet begun to intensify rapidly, but
24-h sequences including and prior to  A and  C have strongly positive LPDs.
These results indicate that structural trajectories of radial profiles in the
North Atlantic may contain signals of RI prior to onset.

Panel (ii) shows Hurricane Jose [2017] and highlights the importance of
vertical wind shear. This TC exhibited deep convection in the core for nearly
two weeks, remaining at elevated RI risk according to the RI-NAL posterior
differences. However, after an initial period of RI, high vertical wind shear
disrupted the TC structure around Sept. 9. The TC decayed from 135 kt to about
70 kt and never appreciably intensified again despite several periods of
elevated LPDs prior to  D . The underlying regression
$\widehat{m}_{\text{post}}$ does not account for the vertical wind shear which
prevented the TC’s intensification. See Section 8.2 for a discussion about how
to potentially account for external factors such as vertical wind shear.

Finally, panel (iii) depicts the short-lived eastern North Pacific TC
Hurricane Rosa [2018], which underwent an extended period of RI before
beginning an eyewall replacement cycle. This evolution included two
interesting phenomena. First, the TC experienced a pause in its rapid
intensification  A ; the LPDs indicate an ongoing RI threat at this time,
consistent with the resumption of intensification 12-18 hours later. Second,
beginning on Sept. 28, the TC underwent an eyewall replacement cycle, which
manifests as a expansion of the eye accompanied by an evening out of
convection across the storm. Such cycles typically result in a decrease of the
TC’s maximum sustained winds. The LPDs mark this shift at the TC’s peak
intensity, as well as the brief period  C where eyewall replacement is
completed prior to landfall.

###  7.3 Structural Trajectories Preceding Rapid Intensification

Figure 10: Detecting RI prior to onset. _Left_ : PCA map of a subset of the
RI-NAL test sample, representing 24-hour sequences prior to RI onset. We show
36 events with at least 24 hours of non-missing radial profiles $X_{t}$ (that
is, at least one sequence $\mathbf{S}_{<t}$). _Right_ : Examples of 48-hour
structural trajectories preceding RI onset. The local posterior differences,
shown in the bar to the left of each map, largely signal TCs with deep core
convection in the last 24 hours as “primed for RI”.

In the previous section, case studies indicated that positive LPDs can lead RI
events. That is, they may signal convective structure primed for RI prior to
RI onset. That core convection can predict RI is known; at least one RII
forecast model used by the NHC includes the fraction of GOES pixels between 50
and 200 km with temperatures below $-30^{\circ}$C as a predictor (Kaplan et
al., 2015). We would thus expect our LPDs to contain some signal leading RI.

In this section, we hone in on a subset of the RI-NAL test sample, which
represents complete 24-hour sequences $\mathbf{S}_{<t}$ leading up to RI onset
points (such as points i-  A and i-  C in Figure 9). This leaves us with a
total of 36 sequences, visualized in the PCA map of Figure 10, left. We then
investigate whether we are able to detect RI within the 48 hours prior to
onset, and if so, what the lead time might be. As before, a positive posterior
difference indicates $p(\mathbf{s}_{t}|Y_{t}=1)>p(\mathbf{s}_{t}|Y_{t}=0)$.
Because $\mathbf{S}_{<t}$ here encodes the recent evolution of convective
structure, we interpret LPD $\lambda(\mathbf{s}_{t})\gg 0$ prior to RI as an
indication that the TC at time $t$ exhibits convective structures primed for
RI onset. That is, $\lambda(\mathbf{s}_{t})\gg 0$ in these cases indicates an
elevated risk of RI based on convective evolution alone.

In 30 (83%) out of 36 onsets, the RI-NAL LPDs indicate an above-average RI
threat at onset. In 28 (78%) and 24 (66%) onsets, LPDs indicate an above-
average RI threat 6 or 12 hours prior to onset, respectively. Some example
trajectories are shown in Figure 10, right. Prior to RI onset, TCs tend to
exhibit deep core convection, evidence of strong outflow (indicated by the
downward slant of $\mathbf{S}_{<t}$), and often pronounced diurnal cycles
(oscillations in $\mathbf{S}_{<t}$ over time). These results indicate that
structural trajectories or LPDs from our test could serve as valuable inputs
to RI forecast models such as the SHIPS-RII (Kaplan et al., 2015).

##  8 Limitations and Potential Extensions of Methods

We end our paper by discussing limitations of the current analysis, and
suggest potential extensions of our work to answer different questions about
the behavior of TCs in the same binary (RI versus not-RI) setting.

###  8.1 Adding Multiple Data Sources, Environmental Variables, and Other
Functional Features

Deep convection as revealed by IR imagery is only one ingredient required for
RI onset. Our current analysis does not include environmental variables, such
as vertical wind shear, which are known to affect RI. Our analysis also
neglects asymmetric patterns in radial profiles, and key ORB functional
features. These limitations however only pertain to the presented TC analysis,
rather than being limitations of the proposed methodology per se. Indeed, our
testing framework for distributional differences for a binary response
variable extend to settings with multiple data sources and different data
types:

For a joint analysis of several variables, one directly adds predictors to the
regression model $\widehat{m}_{\text{post}}$ so that the sequence
$\mathbf{S}_{<t}$ consists of, for example, both structural features (e.g.,
ORB functions) and environmental features (e.g., vertical wind shear or sea
surface temperatures). Furthermore, the CNN model can easily handle multiple
inputs; including other observation bands (e.g., the water vapor-sensitive
$6.9\mu$m GOES band or microwave imagery from polar-orbiting satellites), and
ORB functions other than radial profiles (e.g., the size of temperature level
sets as a function of the temperature threshold; see McNeely et al. 2019,
2020). These inputs are then included as channels, which are combined at the
dense layer (with separate feature extraction layers), or added together at
the logistic layer (a neural additive model with separate feature and dense
layers for each input).

In particular, the current TC analysis of GOES IR-imagery can be extended to
capture asymmetry in convective structure by separating the radial profiles
into four quadrants. Our ongoing work on structural forecasting illustrates
the promise of such an approach (McNeely et al., 2022b).

###  8.2 Adjusting for Other Variables

Our presented framework detects potential associations between high-
dimensional time series and binary labels; any detected association might of
course be due to variables that confound the relationship between the label
and the covariates of interest. Our general methodology can however be
extended to account for confounding variables with the test of independence
then generalizing to a test of conditional independence, and the corresponding
regression problem including additional covariates:

Suppose that we want to detect distributional differences in sequence data
$\\{\mathbf{S}_{<t}\\}_{t\geq 0}$ preceding an event $Y_{t}=1$ versus a non-
event $Y_{t}=0$, after adjusting for the effect of other variables with
sequence data $\\{\mathbf{Z}_{<t}\\}_{t\geq 0}$. For example, wind shear might
confound the relationship between convective structures and RI or RW. Assuming
a stationary process $\\{(\mathbf{S}_{<t},\mathbf{Z}_{<t},Y_{t})\\}_{t\geq
0}$, hence omitting $t$, we test conditional independence

$\displaystyle H_{0}:p(\mathbf{s}|Y=1,\mathbf{Z}=\mathbf{z})$
$\displaystyle=p(\mathbf{s}|Y=0,\mathbf{Z}=\mathbf{z}),\ \text{for all}\
\mathbf{s}\in\mathcal{S}\text{ and all }\mathbf{z},\ \text{versus}$ (10)
$\displaystyle\ H_{1}:p(\mathbf{s}|Y=1,\mathbf{Z}=\mathbf{z})$
$\displaystyle\neq p(\mathbf{s}|Y=0,\mathbf{Z}=\mathbf{z}),\ \text{for some}\
\mathbf{s}\in\mathcal{S}\text{ or }\mathbf{z}.$

These hypotheses are equivalent to

$\displaystyle
H_{0}:\mathbb{P}(Y=1|\mathbf{S}=\mathbf{s},\mathbf{Z}=\mathbf{z})=\mathbb{P}(Y=1|\mathbf{Z}=\mathbf{z}),$
$\displaystyle\text{ for all}\ \mathbf{s}\in\mathcal{S}\text{ and all
}\mathbf{z},$ (11) $\displaystyle
H_{1}:\mathbb{P}(Y=1|\mathbf{S}=\mathbf{s},\mathbf{Z}=\mathbf{z})\neq\mathbb{P}(Y=1|\mathbf{Z}=\mathbf{z}),$
$\displaystyle\text{ for some }\mathbf{s}\in\mathcal{S}\text{ or }\mathbf{z}.$

Analogous to Equation 4, we can define a regression test statistic $\lambda$
based on the difference between an estimate of
$\mathbb{P}(Y=1|\mathbf{S}=\mathbf{s},\mathbf{Z}=\mathbf{z})$ and
$\mathbb{P}(Y=1|\mathbf{Z}=\mathbf{z})$. Note that if $\mathbf{Z}$ is
associated with both $\mathbf{S}$ and $Y$, then a permutation test is not
valid even for Setting A with IID data, because of lack of exchangeability
under $H_{0}$. One solution for (IID as well as DID) sequence data
$\\{(\mathbf{S}_{<t},\mathbf{Z}_{<t},Y_{t})\\}_{t\geq 0}$ is to extend our
bootstrap test to a procedure where one estimates the distribution of the
label series $\\{Y_{t}\\}_{t\geq 0}$ conditional on $\mathbf{Z}$.

Regarding TC rapid intensity change, the admission of confounders would
improve the interpretability of both the RI and RW tests. In the case of RI,
wind shear $\mathbf{Z}$ is a powerful environmental predictor and can inhibit
intensification $Y$ of a TC with otherwise favorable structure $\mathbf{S}$
(e.g., Hurricane Jose [2017], Figure 9, ii-C). Meanwhile, the results for RW
(Figure 8(iii)) appear to weakly capture the obvious relationship: stronger
storms are more likely to rapidly weaken. By accounting for the effect of
current intensity ($\mathbf{Z}$), we could better assess the relationship
between structural evolution ($\mathbf{S}$) and intensity change ($Y$).

###  8.3 Local P-Values

In this work, we refer to the local posterior difference (Equation 6) as a
local diagnostic, rather than as a local p-value, because empirical results
show that we do not control the type I error of a point-wise test
$H_{0}(\mathbf{s}):\mathbb{P}(Y=1|\mathbf{S}=\mathbf{s})=\mathbb{P}(Y=1)$ at
current sample sizes. The LPD value can, however, in principle be used to test
the local null

$\displaystyle
H_{0}^{\epsilon}(\mathbf{s}):\mathbb{P}(Y=1|\mathbf{S}=\mathbf{s}^{\prime})=\mathbb{P}(Y=1)\text{
for all }\mathbf{s}^{\prime}\in B(\mathbf{s};\epsilon),$ (12)

where $B$ is a ball of radius $\epsilon$ centered at $\mathbf{s}$. In McNeely
et al. (2022a), we show that, under DAG B, the local p-values (13) are valid
if the regression estimator $\widehat{m}_{\text{post}}$ only uses the
observations in $\mathcal{D}$ such that $\mathbf{S}\in
B(\mathbf{s};\epsilon)$. The latter assumption holds for regression estimators
that are based on partitions, such as tree-based estimators (random forests,
boosting methods), as well as smoothing kernel estimators with finite support.

###  8.4 Bootstrapping the Label Series

To estimate the null distribution of the test statistic $\lambda$ (3), we
currently assume a $k$-step Markov chain and draw new labels from the Markov
autoregressive model (5). There are other ways one can bootstrap the label
distribution, including adopting sampling schemes that model long-range
dependence in the label sequences. For a review of bootstrap methods for
dependent data, see for example Bühlmann (2002), Horowitz (2003), and Kreiss
and Paparoditis (2011).

##  9 Conclusions

We describe a statistical framework for analyzing the relationship between
complex high-dimensional data $\\{\mathbf{X}_{t}\\}_{t\geq 0}$ and labels
$\\{Y_{t}\\}_{t\geq 0}$. For DID sequence data
$\\{(\mathbf{S}_{t},Y_{t})\\}_{t\geq 0}$, where
$\mathbf{S}_{<t}=\\{\mathbf{X}_{t-T},\mathbf{X}_{t-T+1},\ldots,\mathbf{X}_{t}\\}$
and $T>0$, we propose a two-sample test (Equation 2 and Algorithm 1) with
minimal assumptions beyond stationarity. The test relies on two simple key
ideas: (i) a test statistic based on the posterior difference
$\mathbb{P}(Y=1|\mathbf{S})-\mathbb{P}(\mathbf{S})$, which we estimate using a
machine learning algorithm suitable for the data at hand (empirical results
indicate that the test power depends on the quality of the regression
estimate; Section 6 and Figure 5); and (ii) a bootstrap test, where we
estimate the marginal distribution of $\\{Y_{t}\\}_{t\geq 0}$ (consistency
guarantees asymptotic validity; Theorem 1). Our framework provides
interpretable diagnostics in local posterior differences (Section 7) and can
be extended to include longer-range dependence structures (Section 8.4),
multiple data sources (Section 8.1), and potential confounding variables
(Section 8.2).

###  9.1 TC Results

We detect a distributional difference between sequences leading up to RI
versus not-RI events in both the North Atlantic and eastern North Pacific
basins ($p<0.01$). Local posterior differences for RI-NAL and RI-ENP indicate
that specific types of convection — deep and deepening core convection — are
present both before and during RI (Figures 8, 9). Furthermore, we observe that
particular convective structures are necessary for RI (Figure 10) and thus
useful indicators of future RI, but they are not sufficient to trigger RI on
their own (Figure 9, ii), as the TC environment (e.g., vertical wind shear and
ocean heat content) must also support intensification. Thus, while our current
results have apparent value for RI forecasting, an analysis of structural
trajectories alongside environmental factors such as vertical wind shear
promises better understanding of RI and may improve analysis of RW events as
well.

When posing the same question regarding RW versus not-RW events, we do not
detect a difference in the NAL basin ($p=0.18$); this is expected, as RW is
more likely to be driven by a variety of internal and external factors not
well captured by convective distribution whereas RI generally _requires_ TC
convective structure to be capable of sustaining rapid energy uptake. However,
the ENP basin is characterized by a narrow spatial region of conditions
favorable to TC development such as warm ocean waters and moist air; this
homogeneity leads to a significant signal for RW in the ENP basin ($p<0.01$).

###  9.2 Broader Impact

While we apply our methods to meteorology, the proposed statistical framework
is applicable to any labeled DID sequence data. Labeled video or other
sequence data are common in automation, medical monitoring, and multiple
domains in the physical sciences; in many of these areas, the ability to
identify high-risk patterns in spatio-temporal data could prove
transformative. The flexibility of our framework in admitting a fusion of
multiple data sources and adjustment for other variables is also vital to
analyzing the complex systems in the environmental and physical sciences and
many other domains.

## Supplemental text

The online supplement contains a description of our CNN regressor in Section
7, our algorithm for labeling RI and RW events, and a proof of Theorem 1.

## Data and Code Availability

All data used are publicly available: HURDAT2 at
https://www.nhc.noaa.gov/data/#hurdat and MERGIR at
https://disc.gsfc.nasa.gov/datasets/GPM_MERGIR_1/summary, which can both be
accessed via provided code. All code used to produce the results in this paper
is available at https://github.com/ihmcneely/ORB2sample and as supplementary
material.

## References

  * Aminikhanghahi and Cook (2017) [author] Aminikhanghahi, SamanehS. and Cook, Diane JD. J. (2017). A survey of methods for time series change point detection. Knowledge and information systems 51 339–367.
  * Berrett et al. (2020) [author] Berrett, Thomas BT. B., Wang, YiY., Barber, Rina FoygelR. F. and Samworth, Richard JR. J. (2020). The conditional permutation test for independence while controlling for confounders. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 82 175–197.
  * Bühlmann (2002) [author] Bühlmann, PeterP. (2002). Bootstraps for time series. Statistical science 52–72.
  * Candes et al. (2018) [author] Candes, EmmanuelE., Fan, YingyingY., Janson, LucasL. and Lv, JinchiJ. (2018). Panning for gold: ‘model-X’ knockoffs for high dimensional controlled variable selection. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 80 551–577.
  * Chakravarti et al. (2021) [author] Chakravarti, PurvashaP., Kuusela, MikaelM., Lei, JingJ. and Wasserman, LarryL. (2021). Model-Independent Detection of New Physics Signals Using Interpretable Semi-Supervised Classifier Tests. arXiv preprint arXiv:2102.07679.
  * DeMaria and Kaplan (1999) [author] DeMaria, MarkM. and Kaplan, JohnJ. (1999). An Updated Statistical Hurricane Intensity Prediction Scheme (SHIPS) for the Atlantic and Eastern North Pacific Basins. Weather and Forecasting 14 326-337.
  * Evans and G’Sell (2020) [author] Evans, CiaranC. and G’Sell, MaxM. (2020). Sequential changepoint detection for label shift in classification. arXiv preprint arXiv:2009.08592.
  * Gall et al. (2013) [author] Gall, RobertR., Franklin, JamesJ., Marks, FrankF., Rappaport, Edward NE. N. and Toepfer, FrederickF. (2013). The hurricane forecast improvement project. Bulletin of the American Meteorological Society 94 329–343.
  * Gretton et al. (2012) [author] Gretton, ArthurA., Borgwardt, Karsten MK. M., Rasch, Malte JM. J., Schölkopf, BernhardB. and Smola, AlexanderA. (2012). A kernel two-sample test. The Journal of Machine Learning Research 13 723–773.
  * Grimmett and Stirzaker (2020) [author] Grimmett, GeoffreyG. and Stirzaker, DavidD. (2020). Probability and random processes. Oxford University Press.
  * Horowitz (2003) [author] Horowitz, Joel LJ. L. (2003). Bootstrap methods for Markov processes. Econometrica 71 1049–1082.
  * Hovmöller (1949) [author] Hovmöller, ErnestE. (1949). The trough-and-ridge diagram. Tellus 1 62–66.
  * Janowiak, Joyce and Xie (2020) [author] Janowiak, JohnJ., Joyce, BobB. and Xie, PingpingP. (2020). NCEP/CPC L3 Half Hourly 4km Global (60S - 60N) Merged IR V1. 10.5067/P4HZB9N27EKU
  * Kandasamy et al. (2015) [author] Kandasamy, KirthevasanK., Krishnamurthy, AkshayA., Poczos, BarnabasB., Wasserman, Larry AL. A. and Robins, James MJ. M. (2015). Nonparametric von Mises Estimators for Entropies, Divergences and Mutual Informations. In NIPS 15 397–405.
  * Kaplan and DeMaria (2003) [author] Kaplan, JohnJ. and DeMaria, MarkM. (2003). Large-scale characteristics of rapidly intensifying tropical cyclones in the North Atlantic basin. Weather and forecasting 18 1093–1108.
  * Kaplan, DeMaria and Knaff (2010) [author] Kaplan, JohnJ., DeMaria, MarkM. and Knaff, John A.J. A. (2010). A Revised Tropical Cyclone Rapid Intensification Index for the Atlantic and Eastern North Pacific Basins. Weather and Forecasting 25 220-241.
  * Kaplan et al. (2015) [author] Kaplan, JohnJ., Rozoff, Christopher MC. M., DeMaria, MarkM., Sampson, Charles RC. R., Kossin, James PJ. P., Velden, Christopher SC. S., Cione, Joseph JJ. J., Dunion, Jason PJ. P., Knaff, John AJ. A., Zhang, Jun AJ. A. et al. (2015). Evaluating environmental impacts on tropical cyclone rapid intensification predictability utilizing statistical models. Weather and Forecasting 30 1374–1396.
  * Katsevich and Ramdas (2020) [author] Katsevich, EugeneE. and Ramdas, AadityaA. (2020). A theoretical treatment of conditional independence testing under model-x. arXiv preprint arXiv:2005.05506 4\.
  * Kim, Lee and Lei (2019) [author] Kim, IlmunI., Lee, Ann BA. B. and Lei, JingJ. (2019). Global and local two-sample tests via regression. Electronic Journal of Statistics 13 5253–5305.
  * Kim et al. (2021) [author] Kim, IlmunI., Ramdas, AadityaA., Singh, AartiA. and Wasserman, LarryL. (2021). Classification accuracy as a proxy for two-sample testing. The Annals of Statistics 49 411–434.
  * Klotzbach et al. (2018) [author] Klotzbach, Philip J.P. J., Bowen, Steven G.S. G., Pielke, RogerR. and Bell, MichaelM. (2018). Continental U.S. Hurricane Landfall Frequency and Associated Damage: Observations and Future Risks. Bulletin of the American Meteorological Society 99 1359-1376. 10.1175/BAMS-D-17-0184.1
  * Knaff and DeMaria (2017) [author] Knaff, John AJ. A. and DeMaria, Robert TR. T. (2017). Forecasting tropical cyclone eye formation and dissipation in infrared imagery. Weather and Forecasting 32 2103–2116.
  * Kreiss and Paparoditis (2011) [author] Kreiss, Jens-PeterJ.-P. and Paparoditis, EfstathiosE. (2011). Bootstrap methods for dependent data: A review. Journal of the Korean Statistical Society 40 357–378.
  * Landsea and Franklin (2013) [author] Landsea, Christopher W.C. W. and Franklin, James L.J. L. (2013). Atlantic Hurricane Database Uncertainty and Presentation of a New Database Format. Monthly Weather Review 141 3576-3592. 10.1175/MWR-D-12-00254.1
  * Li and Racine (2007) [author] Li, QiQ. and Racine, Jeffrey ScottJ. S. (2007). Nonparametric econometrics: theory and practice. Princeton University Press.
  * Luo et al. (2014) [author] Luo, ChenC., Lou, Jian-GuangJ.-G., Lin, QingweiQ., Fu, QiangQ., Ding, RuiR., Zhang, DongmeiD. and Wang, ZheZ. (2014). Correlating events with time series for incident diagnosis. In Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining 1583–1592.
  * McNeely et al. (2019) [author] McNeely, TreyT., Lee, Ann BA. B., Hammerling, DoritD. and Wood, KimberlyK. (2019). Quantifying the Spatial Structure of Tropical Cyclone Imagery.
  * McNeely et al. (2020) [author] McNeely, TreyT., Lee, Ann BA. B., Wood, Kimberly MK. M. and Hammerling, DoritD. (2020). Unlocking GOES: A Statistical Framework for Quantifying the Evolution of Convective Structure in Tropical Cyclones. Journal of Applied Meteorology and Climatology 59 1671–1689.
  * McNeely et al. (2022a) [author] McNeely, TreyT., Vincent, GalenG., Izbicki, RafaelR., Wood, Kimberly MK. M. and Lee, Ann BA. B. (2022a). Supplemental Materials for Detecting Distributional Differences in Labeled Sequence Data with Application to Tropical Cyclone Satellite Imagery.
  * McNeely et al. (2022b) [author] McNeely, TreyT., Khokhlov, PavelP., Dalmasso, NiccoloN., Wood, Kimberly MK. M. and Lee, Ann BA. B. (2022b). Structural Forecasting for Short-term Tropical Cyclone Intensity Guidance. arXiv preprint arXiv:2206.00067.
  * Moon and Hero (2014) [author] Moon, KevinK. and Hero, AlfredA. (2014). Multivariate f-divergence estimation with confidence. Advances in neural information processing systems 27 2420–2428.
  * Rogers (2010) [author] Rogers, RobertR. (2010). Convective-scale structure and evolution during a high-resolution simulation of tropical cyclone rapid intensification. Journal of the Atmospheric Sciences 67 44–70.
  * Sanabia, Barrett and Fine (2014) [author] Sanabia, Elizabeth R.E. R., Barrett, Bradford S.B. S. and Fine, Caitlin M.C. M. (2014). Relationships between Tropical Cyclone Intensity and Eyewall Structure as Determined by Radial Profiles of Inner-Core Infrared Brightness Temperature. Monthly Weather Review 142 4581-4599. 10.1175/MWR-D-13-00336.1
  * Scharwächter and Müller (2020a) [author] Scharwächter, ErikE. and Müller, EmmanuelE. (2020a). Does terrorism trigger online hate speech? On the association of events and time series. The Annals of Applied Statistics 14 1285–1303.
  * Scharwächter and Müller (2020b) [author] Scharwächter, ErikE. and Müller, EmmanuelE. (2020b). Two-sample testing for event impacts in time series. In Proceedings of the 2020 SIAM International Conference on Data Mining 10–18. SIAM.
  * Sesia, Sabatti and Candès (2019) [author] Sesia, MatteoM., Sabatti, ChiaraC. and Candès, Emmanuel JE. J. (2019). Gene hunting with hidden Markov model knockoffs. Biometrika 106 1–18.
  * Wood and Ritchie (2015) [author] Wood, Kimberly MK. M. and Ritchie, Elizabeth AE. A. (2015). A definition for rapid weakening of North Atlantic and Eastern North Pacific tropical cyclones. Geophysical Research Letters 42 10–091.

##  Appendix A Estimating ${m}_{\text{seq}}$ and ${m}_{\text{post}}$

Estimating $m_{\text{seq}}$: RI/RW events extending beyond 48 consecutive
hours are extremely rare ($\approx 5\%$ of events). As such, we estimate
$m_{\text{seq}}$ with a Markov chain of length $k=8$ _at the 6-hour synoptic
resolution_. When drawing labels $\\{\widetilde{Y}\\}_{t\geq 0}$ for a storm,
we draw from the Markov chain at synoptic times (0000 UTC, 0600 UTC, 1200 UTC,
and 1800 UTC), then fill to the 30-minute resolution as described in Section 3
($\widetilde{Y}_{t}=1$ between _consecutive_ synoptic $\widetilde{Y}_{t}=1$).
This approximates the distribution of RI event lengths well up 48 hours with a
slightly heavy tail beyond 48 hours. We note that each TC requires a separate
sequence of draws from $m_{\text{seq}}$; we thus initialize and burn-in the MC
separately for each TC (see Algorithm 1 in main paper for details).

Estimating $m_{\text{post}}$: Each radial profile is sampled at a 5-km
resolution from 0-5km out to 395-400km (for 80 samples in the radius
dimension). Our model architecture is identical for all four combinations of
basin (NAL/ENP) and event type (RI/RW). Our CNN consists of 3 Convolution-
Maximum Pooling layers (8 channels, 8 channels, and 4 channels respectively;
convolution size $3\times 3$, maximum pooling size $2\times 2$) followed by a
fully connected layer and a sigmoid output layer. See Table 2 for the full
architecture.

When we initially estimate $m_{\text{post}}$, we train all parameters. For
each subsequent fitting of $\widehat{m}_{\text{post}}^{(b)}$, we freeze the
feature layers (i.e. layers 1, 4, and 7), only re-initializing the fully
connected and output layers (i.e. the “regression” layers, 11 and 13) for
retraining. For each of the four basin-event type combinations, we train 101
CNN models ($B=100$). $\widehat{m}_{\text{post}}$ achieves approximately
75%-80% balanced accuracy on held-out test data for RI in each basin, and
approximately 60%-65% for RW, while $\widehat{m}_{\text{post}}^{(b)}$ (trained
on random labels) only has $\approx 50\%$ balanced accuracy on held-out test
data in all cases.

Layer (Type-ID) | Output Shape | Parameters
---|---|---
Conv2d-1: | [N, 8, 50, 82] | 80
ReLU-2: | [N, 8, 50, 82] | 0
MaxPool2d-3: | [N, 8, 25, 41] | 0
Conv2d-4: | [N, 8, 27, 43] | 584
ReLU-5: | [N, 8, 27, 43] | 0
MaxPool2d-6: | [N, 8, 13, 21] | 0
Conv2d-7: | [N, 4, 15, 23] | 292
ReLU-8: | [N, 4, 15, 23] | 0
MaxPool2d-9: | [N, 4, 7, 11] | 0
Dropout-10: | [N, 308] | 0
Linear-11: | [N, 256] | 79,104
ReLU-12: | [N, 256] | 0
Linear-13: | [N, 2] | 514
Sigmoid-14: | [N, 2] | 0
Total params: |  | 80,574
Table 2: Parameters of the CNN used to estimate $m_{\text{post}}$. The model
consists of three convolutional blocks (layers 1-9). Dropout (layer 10) with
probability 50% is used to regularize the model.

##  Appendix B Rapid Intensity Change

Operationally, RI/RW events are defined by a threshold _and_ a time. For
example, the rapid intensification index (RII; Kaplan et al. (2015)) provides
the probability of intensification $c$ in a given time period; this includes
24-hour periods such as 25 kt in 24 hr as well as shorter (or longer) periods
such as 20 kt in 12 hr (or 55 kt in 48 hr). In lieu of selecting a
discretization of time in this manner, we present an algorithm for identifying
RI/RW events of any length with a daily rate of intensification of at least
$c$ per 24 hours (Algorithm 2). By admitting events both shorter and longer
than 24 hours, we prevent the overlap of RI and RW events. See Figure 11 for
an example. The operational (25 kt in 24 hr) definition for RI and RW results
in significant overlap between RI and RW events; for example, the transition
from RI to RW on July 15 for Hurricane Carlos [2009] sees an 18-hour overlap
between operationally defined RI and RW events. Our labeling procedure
restricts the set of $Y=1$ times for each event by trimming non-intensifying
(or non-weakening, in the case of RW) times from the beginning and end of each
event, while also admitting sufficiently rapid intensification (or weakening)
in shorter time periods than 24 hours.

In words, Algorithm 2 labels events as such:

  1. 1.

Identify 24-hour periods containing intensification of at least $c$.

  2. 2.

Trim the beginning and end of each identified 24-hour RI period such that they
begin and end with 6-hour intensification.

  3. 3.

Label $Y_{t}=1$ if any trimmed RI period contains $t$.

RW is identified symmetrically by reversing the input sequence (intensities)
in time prior to labelling. By virtue of step 2, RI and RW events can only
share endpoints, but cannot overlap.

Figure 11: Comparison of RI/RW labeling via 24-hour operational windows (left)
versus via Algorithm 2 (right). Algorithm 2 restricts event labels to only
intensifying/weakening times as appropriate, resulting in more logical, non-
overlapping labels for the endpoints (e.g., July 15 in the above figure).

Require: Sequence of contiguous intensity (i.e. maximum wind speed)
observations $\\{W_{i}\\}_{t=1}^{T}$ for a storm and an intensity change
threshold $c$.

Initialize: $\mathbf{Z}=\mathbf{0}_{T,T}$, $Y=\mathbf{0}_{T}$

for  _$t$ in $1:T-1$_ do

$\Delta_{t,1}=W_{t+1}-W_{t}$; // lead-1 intensity change

if  _$t\leq T-4$_ then

$\Delta_{t,4}=\max(W_{t},\dots,W_{t+4})-W_{t}$; // lead-4 intensity change

end if

end for

$\mathcal{A}=\\{t:\Delta_{t,4}\geq c\\}$; // 24-hour windows containing RI

$\mathcal{B}=\\{t:\Delta_{t,1}>0\\}$; // 6-hour windows containing
intensification

for  _$t$ in $1:T-4$_ do

if  _$t\in\mathcal{A}$_ then // label 24-hour RI window

$Z_{t,t},\dots,Z_{t,t+4}=1$;

$h=4$;

while  _$h\geq 1$ and $t+h-1\in\mathcal{B}$_ do

$Z_{t,t+h}=0$;// trim non-intensification from end of event

$h=h-1$;

end while

$h=1$;

while  _$h <4$ and $t+h\in\mathcal{B}$_ do

$Z_{t,t+h}=0$;// trim non-intensification from start of event

$h=h+1$;

end while

end if

end for

for  _$t$ in $1:T$_ do

$Y_{t}=\max(Z_{t,1},\dots,Z_{t,T})$; // points only need to be valid for one
start

end for

Output: $\\{Y_{i}\\}_{t=1}^{T}$

Note: The above algorithm identifies rapid intensification. To identify rapid
weakening instead, reverse the input sequence $\\{W_{i}\\}_{t=1}^{T}$ at
initialization, then reverse the output sequence $\\{Y_{i}\\}_{t=1}^{T}$ at
output.

Algorithm 2 Identifying observations that fall within rapid intensification
events.

##  Appendix C Proofs

###### Proof of Theorem 1.

Let $F_{G_{\widehat{\mathbf{p}}_{t_{2}}}}$ be the cumulative distribution
function of $\\{Y^{0}_{t}\\}_{t\in\mathcal{T}_{\text{1}}}$. Then, under the
null hypothesis,

$\displaystyle\mathbb{P}\left(\mathbf{S}_{<t}\leq\mathbf{s}_{<t},Y^{0}_{t}\leq
y_{t},\ \forall t\in\mathcal{T}_{1}\right)$
$\displaystyle=\mathbb{P}\left(\mathbf{S}_{<t}\leq\mathbf{s}_{<t},\ \forall
t\in\mathcal{T}_{1}\right)\mathbb{P}\left(Y^{0}_{t}\leq y_{t},\ \forall
t\in\mathcal{T}_{1}\right)$ By construction
$\displaystyle=\mathbb{P}\left(\mathbf{S}_{<t}\leq\mathbf{s}_{<t},\ \forall
t\in\mathcal{T}_{1}\right)F_{G_{\widehat{\mathbf{p}}_{t_{2}}}}\left(y_{t},\
\forall t\in\mathcal{T}_{1}\right)$ By definition
$\displaystyle\xrightarrow[t_{2}\longrightarrow\infty]{\enskip\enskip}\mathbb{P}\left(\mathbf{S}_{<t}\leq\mathbf{s}_{<t},\
\forall t\in\mathcal{T}_{1}\right)F_{G^{*}}\left(y_{t},\ \forall
t\in\mathcal{T}_{1}\right)$ Assumption 4
$\displaystyle=\mathbb{P}\left(\mathbf{S}_{<t}\leq\mathbf{s}_{<t},\ \forall
t\in\mathcal{T}_{1}\right)\mathbb{P}\left(Y_{t}\leq y_{t},\ \forall
t\in\mathcal{T}_{1}\right)$ By definition
$\displaystyle=\mathbb{P}\left(\mathbf{S}_{<t}\leq\mathbf{s}_{<t},Y_{t}\leq
y_{t},\ \forall t\in\mathcal{T}_{1}\right).$ Assumption 2

It follows that

$\mathcal{D}^{t_{2}}_{0}\xrightarrow[t_{2}\longrightarrow\infty]{\enskip\mbox{Dist}\enskip}\mathcal{D}.$

The conclusion of the theorem follows from the continuous mapping theorem and
Assumption 3. ∎

###### Proof of Corollary 1.

By the law of large numbers, for every realization $D$ of $\mathcal{D}$,

$\widehat{p}^{\
t_{2}}_{B}(D)\xrightarrow[B\longrightarrow\infty]{\enskip\mbox{a.s.}\enskip}\
\mathbb{P}\left(\lambda(\mathcal{D}^{t_{2}}_{0})>\lambda({D})\right).$

Now, Theorem 1 implies that

$\displaystyle\lim_{t_{2}\longrightarrow\infty}\mathbb{P}\left(\lambda(\mathcal{D}^{t_{2}}_{0})>\lambda({D})\right)$
$\displaystyle=1-\lim_{t_{2}\longrightarrow\infty}\mathbb{P}\left(\lambda(\mathcal{D}^{t_{2}}_{0})\leq\lambda({D})\right)$
$\displaystyle=1-\lim_{t_{2}\longrightarrow\infty}F_{\lambda(\mathcal{D}^{t_{2}}_{0})}\left(\lambda({D})\right)=1-F_{\lambda(\mathcal{D})}\left(\lambda({D})\right).$

Conclude that

$\lim_{t_{2}\longrightarrow\infty}\lim_{B\longrightarrow\infty}\widehat{p}^{\
t_{2}}_{B}(\mathcal{D})=1-F_{\lambda(\mathcal{D})}\left(\lambda(\mathcal{D})\right)$

almost surely, and thus this convergence also holds in distribution. The
conclusion follows from the fact that
$F_{\lambda(\mathcal{D})}\left(\lambda(\mathcal{D})\right)$ is a uniform
random variable, and therefore so is
$1-F_{\lambda(\mathcal{\mathcal{missing}}{D})}\left(\lambda({D})\right)$. ∎

##  Appendix D Local p-values

In this section we show that LPDs can be used to test the local null
hypothesis

$\displaystyle
H_{0}^{\epsilon}(\mathbf{s}):\mathbb{P}(Y=1|\mathbf{S}=\mathbf{s}^{\prime})=\mathbb{P}(Y=1)\text{
for all }\mathbf{s}^{\prime}\in B(\mathbf{s};\epsilon),$ (13)

where $B$ is a ball of radius $\epsilon$ centered at $\mathbf{s}$.

This local test needs further assumptions to asymptotically control type I
error probabilities, which we describe in what follows.

######  Assumption 5 (Local regression estimator).

$\widehat{m}_{\text{post}}$ only uses the observations in $\mathcal{D}$ such
that $\mathbf{S}_{<t}\in B(\mathbf{s};\epsilon)$.

Assumption 5 holds for regression estimators that are based on partitions,
such as tree-based estimators (e.g., random forests, boosting methods) or
smoothing kernel estimators with kernels with bounded support.

######  Assumption 6 (Conditional independence).

$\\{(\mathbf{S}_{<t},Y_{t})\\}_{t\geq 0}$ satisfies the DAG of Setting B
(Figure 2).

######  Theorem 2.

Assume that the measure $Q$ in the test statistic from Equation 8 places all
of its mass on $B(\mathbf{s};\epsilon)$. Also assume that Assumptions 1, 2, 5
and 6 hold. Then, under the local null hypothesis in Equation 13, if the
marginal distribution is estimated using the estimator of Example 1, then, for
any $0<\alpha<1$

$\lim_{t_{2}\longrightarrow\infty}\lim_{B\longrightarrow\infty}\mathbb{P}\left(\widehat{p}^{\
t_{2}}_{B}(\mathcal{D})\leq\alpha\right)=\alpha.$

###### Proof.

Let $\mathcal{I}=\\{t\in\mathcal{T}_{1}:\mathbf{S}_{<t}\in B(s;\epsilon)\\}$.
Under $H_{0}^{\epsilon}(\mathbf{s})$,

$\displaystyle\mathbb{P}(\mathbf{S}_{<t}\leq\mathbf{s}_{<t}\ $
$\displaystyle\forall t\in\mathcal{T}_{1}\mbox{ and }Y^{0}_{t}\leq y_{t}\
\forall t\in\mathcal{I}|\mathcal{I})$
$\displaystyle=\mathbb{P}\left(\mathbf{S}_{<t}\leq\mathbf{s}_{<t}\ \forall
t\in\mathcal{T}_{1}\right|\mathcal{I})\mathbb{P}\left(Y^{0}_{t}\leq y_{t}\
\forall t\in\mathcal{I}|\mathcal{I}\right)$ By construction
$\displaystyle\xrightarrow[t_{2}\longrightarrow\infty]{\enskip\enskip}\mathbb{P}\left(\mathbf{S}_{<t}\leq\mathbf{s}_{<t}\
\forall t\in\mathcal{T}_{1}|\mathcal{I}\right)\mathbb{P}\left(Y_{t}\leq y_{t}\
\forall t\in\mathcal{I}|\mathcal{I}\right)$ Example 1
$\displaystyle=\mathbb{P}(\mathbf{S}_{<t}\leq\mathbf{s}_{<t}\ \forall
t\in\mathcal{T}_{1}\mbox{ and }Y_{t}\leq y_{t}\ \forall
t\in\mathcal{I}|\mathcal{I}).$ Assumption 6

The result follows from noticing that, under Assumption 5,
$\lambda(\mathcal{D})$ depends on $\mathcal{D}$ only through
$\\{\mathbf{S}_{<t}\\}_{t\in\mathcal{T}_{1}}$ and
$\\{(\mathbf{S}_{<t},Y_{t})\\}_{t\in\mathcal{I}}$ and by repeating the
arguments of the proofs of Theorem 1 and Corollary 1.

∎