Equivalence relations and $L^{p}$ distances between time series with application to the Black Summer Australian bushfires

We begin with the following general question: can one define equivalence relations between time series, whereby two time series are declared equivalent if they differ only up to noise? This is our guiding heuristic; in this paper, we give one concrete way to understand this.

We begin with motivation: in set theory and throughout mathematics, an equivalence relation $\sim$ on a set $\mathcal{X}$ can simplify the set. Indeed, an equivalence relation partitions the original set into equivalence classes - subsets of $\mathcal{X}$ in which everything is mutually equivalent. By considering the set of equivalence classes $\mathcal{X}/\sim$, equivalent elements of the original $\mathcal{X}$ have been collapsed or divided down to the same element of $\mathcal{X}/\sim$; for this reason, the latter is called the quotient set.

If we wish to understand time series as equivalent if they differ only up to noise, and then wish to define distance measures that “cancel out the noise”, we would like to define equivalence relations $\sim$ such that the quotient set $\mathcal{X}/\sim$ is metrizable, namely can be imbued with an appropriate metric. Then, two equivalence classes of time series $A$ and $B$ would have an unambiguous distance $d$ apart. That is, any time series $X_{t}$ in the $A$ class and any time series $Y_{t}$ in the $B$ class would have a noise-invariant distance $d$ apart.

To concretely define these equivalence relations, we take the posture that time series consist of locally stationary segments separated by change points (or structural breaks) at which the statistical properties of the time series change, a frequent assumption in the statistical literature [1]. If two time series have the exact same structural breaks and appropriately chosen statistical properties (up to completely identical distributions) in each of the segments, then the difference between them is random noise. This is of course just one way to build a framework envisioned by the paragraphs above.

The following theorem proves that the above goals are possible. Formally, it combines Theorems 2.4 and 2.5:

The proof and construction of such equivalence relations build on several insights: the change point algorithms developed by [2, 3] and [4], the distances between change points of [5], and the construction of Lebesgue $L^{p}$ spaces as metrizable quotients by [6].

The most significant equivalence relation between functions in measure theory is that of almost everywhere equivalence: define $f\sim g$ if $f=g$ a.e., that is, $\{x:f(x)\neq g(x)\}$ is measure zero. We rely on this crucial equivalence relation for both the motivation and mathematics of our theorem and construction.

We recall the construction of the $L^{p}$ spaces. Let $\mathcal{L}^{p}([0,H])$ be the vector space of all Lebesgue-measurable functions $f$ on $[0,H]$ such that $|f|^{p}$ has finite integral. We have a quotient mapping, defined and outlined in [7],

We will use this notation in what follows. One then equips this quotient space with the $L^{p}$ norms, as defined in [7] or in Section 2 for our specific case.

Let $(X_{t})_{t\in I}$ be a real-valued time series over some time interval $I=\{0,1,...,H\}$, and $\mathfrak{X}$ the space of such time series. In this paper, we only record each data series once, so treat each $X_{t_{i}}$ as a single observed real number. Let $P_{f}\mathbb{R}$ be the set of finite subsets of $\mathbb{R}$. Throughout the paper, $p\geq 1$ is a chosen real number.

First, one selects a statistical attribute, such as mean, variance, or higher-order moment. For simplicity of exposition, we describe the procedure for the mean. We use the two-step algorithm of [16] to determine a set of change points $S=\{c_{1},...,c_{m}\}$ with respect to changes in the mean. This is an inductive deterministic procedure that uses iterated hypothesis testing. At each step, the null hypothesis of no change point existing is tested, and subject to preset parameters, the null hypothesis may be ruled invalid; this determines a change point. In addition, we set $c_{0}=0,c_{m+1}=H$. Properly interpreted, this gives a function

Different parameters yield different algorithms and hence different such functions. This procedure has algorithmically determined that $(X_{t})$ has a constant mean over each interval $[c_{i},c_{i+1}],i=0,...,m$. By using the change point algorithm or simply averaging $(X_{t})$ over these intervals, we compute and record the empirically observed means $\mu_{i}$. Now we define a function

as a weighted indicator function. This is a piecewise constant function that codifies the change points of the time series and its changing means. That is,

Let $\mathcal{F}(X_{t}):=[f]\in L^{p}([0,H])$ be the image of $f$ under the almost everywhere equivalence relation. This step means that the values of $f$ at the change points $c_{i}$ do not matter, for a finite set of points is measure zero. Let $V$ be the subspace of $L^{p}([0,H])$ consisting of (the images of) piecewise constant functions. Properly interpreted, our procedure is the function

This defines an equivalence relation $\sim$ on $\mathfrak{X}$ and an embedding

We note that $V$ is a vector subspace of $L^{p}([0,H])$ and thus inherits its $L^{p}$ norm $\|.\|_{p}$, defined by

normalized so that $\|\mathbbm{1}_{[0,H]}\|=1$. When $p=2$, there is an inner product associated to the norm:

At this point we must make some careful definitions and observations. While $\mathfrak{X}$ has a vector space structure, the change point algorithms do not respect addition of time series, and hence $\mathcal{F}$ is not linear. The absence of linearity forces us to differentiate two separate definitions. First, the metric space structure on $V$ induces a metric on $\mathfrak{X}/\mathord{\sim}$ and a pseudo-metric on $\mathfrak{X}$ itself. Concretely, if $(X_{t}),(Y_{t})$ map to piecewise constant functions $f,g$ respectively, define $d_{p}(X_{t},Y_{t})=\|f-g\|_{p}$.

Secondly and separately, we can pullback the norm structure to $\mathfrak{X}$ as a measure of the overall magnitude of the time series, defining $\operatorname{mag}(X_{t})=\|\mathcal{F}(X_{t})\|_{p}=\|f\|_{p}$. We use the symbol mag, rather than $\|.\|$, as this is neither a norm nor semi-norm on $\mathfrak{X}$. Concretely, $\operatorname{mag}(X_{t})=d_{p}(X_{t},0)$ but $d_{p}(X_{t},Y_{t})\neq\operatorname{mag}(X_{t}-Y_{t}).$ That is, mag does not induce the metric $d_{p}.$ Only after passage to $V$, via the non-linear $\mathcal{F}$, does the norm induce the metric.

Now we prove that $d_{p}$ has all properties of a metric on $\mathfrak{X}/\mathord{\sim}$.

1. 1.

   It is immediate that $d_{p}(X_{t},Y_{t})\geq 0$ always.

2. 2.

   It is immediate that $d_{p}(X_{t},Y_{t})=d_{p}(Y_{t},X_{t})$.

3. 3.

   The triangle inequality $d_{p}(X_{t},Z_{t})\leq d_{p}(X_{t},Y_{t})+d_{p}(Y_{t},Z_{t})$ follows from Lemma 2.1.

4. 4.

   If $d_{p}(X_{t},Y_{t})=0$, then $(X_{t}),(Y_{t})$ are equal elements of $\mathfrak{X}/\mathord{\sim}$ by Lemma 2.2.

In particular, equivalence does not depend on the parameter $p$. Figure 1 depicts some examples of equivalent time series.

Combining all these results we have proved the following:

By an identical method, using change point algorithms for changes in variance, one can construct deterministic functions $\mathcal{F}_{var}$ that prove the following:

Other statistical properties induce other equivalence relations and maps $\mathcal{F}$, all producing variants of Theorem 1.1. There exist change point algorithms that detect changes in all stochastic properties simultaneously; these induce finer equivalence relations.

In this section, we investigate some desirable properties of our metrics $d_{p}$ relative to existing options, by both theoretical propositions and simulation. First, we briefly compare our new distance to the simple Euclidean distance between time series.

Heuristically, $X_{t},X^{\prime}_{t}$ and $Y_{t},Y^{\prime}_{t}$ differ only up to noise. Our distance $d_{p}$ filters out the noise and returns the same distance between $X_{t}$ and $Y_{t}$ as $X^{\prime}_{t}$ and $Y^{\prime}_{t}$. However, if the noise between $X_{t}$ and $X^{\prime}_{t}$ is “opposite” that of the noise between $Y_{t}$ and $Y^{\prime}_{t}$, the Euclidean metric may be sensitive to this and add these errors rather than remove them.

Next, we examine the performance of our metrics $d_{p}$ relative to previously proposed such distances between time series based on structural breaks. Our previous work [5] uses metrics and semi-metrics between the sets of change points to measure distance between time series; that work does so without using data of the mean or variance between these breaks. The following simulation study and proposition provide justification for why the distance $d_{p}$ may more appropriately measure distance between time series.

For this simulation study, we will make use of the following (semi)-metrics between finite sets $S$ and $T$. First, let

Then, let the Hausdorff metric be defined as

The modified Hausdorff semi-metric is defined as

The new family of semi-metrics introduced in [5] is defined by

These aforementioned (semi-)metrics between finite sets can readily be applied to measure distance between time series with respect to their change points. If $(X_{t})$ and $(Y_{t})$ are time series with change point sets $S,T$, respectively (under some change point algorithm), one can simply set $d(X_{t},Y_{t})$ to be any of $d_{H}(S,T),d^{1}_{MH}(S,T)$ or $d^{1}_{MJ}(S,T).$

We generate ten synthetic time series (labeled TS1 - TS10), which we display in Figure 2. These provide example where two different structural breaks (change points) may have markedly different changes in the mean. Black vertical lines indicate detected change points, with the orange line representing the piecewise constant functions $f$. When considering only structural breaks, the most similar time series are groups as follows: {TS1, TS2, TS3}, {TS4, TS5}, {TS6, TS7}, {TS8, TS9, TS10}. However, considering structural breaks and their mean, using our metric $d_{p}$ induced from $\mathcal{F}_{\mu}$, the most similar time series are {TS1, TS3} and {TS7, TS8}.

In Figure 3, we implement hierarchical clustering on the ten synthetic time series with our new $d_{p}$ metric, here with $p=1$, as well as existing distance measures in [5]. We see that only the $d_{p}$ metric, in Figure 3(d), appropriately clusters similar time series according to their changing values, while the other three figures identify erroneous similarities between the time series based solely on their structural breaks. That is, the other (semi-)metrics neglect the additional data that $d_{p}$ recognizes, demonstrating the superior performance of the latter.

Next, we prove a proposition that demonstrates that only our new metrics $d_{p}$ are continuous with respect to small deformations in the time series that may introduce new structural breaks. That is, the existing distance measures between change points may be excessively sensitive to new structural breaks, either caused by a deformation in the time series data or a faulty detection by a change point algorithm. We adopt the same suppositions of Lemma 2.3, illustrated in Figure 1.

Next, we prove a property that establishes a form of linearity for the metrics $d_{p}$ analogous to the Euclidean metric or other classical $L^{p}$ metrics.

The proof of this proposition is quite trivial, but understanding the primary condition is more interesting. The property that $\mathcal{F}(aX_{t}+b)=a\mathcal{F}(X_{t})+b\in V$ really has two components: first, that the change points of $aX_{t}+b$ are the same as $X_{t}$, assuming $a\neq 0$, and secondly, that the statistical quantity $\mathcal{F}$ is measuring scales commutes with scalar multiplication and the addition of constants. The property (33) has a near-identical interpretation. The former component is more subtle than it looks, relying on the fact that the change point algorithm use normalized test statistics. Fortunately, it is the case that essentially all change point algorithms do use normalized test statistics (and should!) so that means this condition holds for all procedures we will use. As for the latter component, the property (32) holds for the mean and hence $\mathcal{F}_{\mu}$ of Theorem 2.4. This is simply due to the fact that the mean of a random variable $Z$ satisfies $E(aZ+b)=aE(Z)+b$. On the other hand, the property (33) holds for the standard deviation and hence $\mathcal{F}_{\sigma}$ closely related to $\mathcal{F}_{var}$ of Theorem 2.5. This is also simply due to the fact that $\sigma(aZ+b)=|a|\sigma(Z)$.

The explanation above and the simple proof show that this proposition applies to the mean-mapping $\mathcal{F}_{\mu}$ and associated distances of Theorem 2.4 as well the standard deviation-mapping $\mathcal{F}_{\sigma}$, related to the variance-mapping $\mathcal{F}_{var}$ of Theorem 2.5. That is, $d_{p}$ respect linear operations between time series (scalar multiplication and the addition of a constant) for both mean and standard deviation.

We conclude this section by considering several invariances our $d_{p}$ metrics do and do not have, and some modifications one may make for a desired application. First, Proposition 3.1 immediately establishes noise-invariance, when properly defined, the guiding property of the entire framework. Next, as a result of Proposition 3.3, $d_{p}$ are not invariant to amplitude (caling $X_{t}\mapsto aX_{t}$) or offset in the $y$-direction (when $X_{t}$ is left invariant but $Y_{t}\mapsto Y_{t}+b$). However, both of these invariances can be attained by simply normalizing the time series prior to analysis (computing Z-scores) as required by the application. This is required in most domains [50, 42], and in fact we do perform such a normalization when computing normalized distances in Section 4.

Other invariances must be established by modifying our $d_{p}$ distances, which should be possible by combining with existing approaches. For example, [42] elegantly multiply the Euclidean metric by an ordered ratio of the complexities of two time series to establish a complexity-invariant distance. That approach should work analogously for our $d_{p}$ metrics. Phase invariance is a complex problem, for which typically the best approach is testing all possible alignments [45, 115]. A similar problem exists for (global) warping, namely (uniform) scaling in the $t$-direction: one must test all possibilities within a given range [43]. Such approaches could be incorporated for our $d_{p}$ distances.

In this section, we describe a general procedure to use the function $\mathcal{F}$ and the induced distances $d_{p}$ to compute and analyze distances in a collection of time series. We will subsequently apply this procedure to analyze the Australian bushfires in the next section. Suppose $(X_{t}^{(1)}),...,(X_{t}^{(n)})$ are time series over some time period $[0,H]$. We suppress the time $t$ in our notation.

First, we choose an appropriate change point algorithm for the data, including threshold parameters, and a desired statistical quantity, such as mean or variance. This specifies a precise instance of $\mathcal{F}$. Let the piecewise constant functions associated to the time series $X^{(i)}$ be $f_{i}$, so that $\mathcal{F}(X^{(i)})=[f_{i}],i=1,...,n$. Also, we consider the normalized piecewise constant functions $\hat{f_{i}}=\frac{f_{i}}{\|f_{i}\|_{p}}$. We form three different matrices:

Given a distance matrix $D$ one may produce an affinity matrix $A$ by:

$D^{us}$ is a distance matrix, thus named as it consists of unscaled distances. $D^{\text{norm}}$ is also a distance matrix, consisting of distances between normalized vectors, again with $p\geq 1$ suppressed from the notation. We refer to $\Omega$ as an alignment matrix as it measures angles between the $f_{i}$ as vectors in an inner product space. We note that all entries of $\Omega$ lie in $[-1,1]$ with all diagonal elements equal to $1.$ To the distance matrices $D^{us}$ and $D^{\text{norm}}$ we associate affinity matrices $A^{us}$ and $A^{\text{norm}}$, respectively.

In this section, we analyze AQI time series from $n=52$ locations across NSW [107]. Each time series consists of hourly data from October 20, 2019 - January 20, 2020, to yield time series of length $H=2210$. In the rare case of missing data points, we take the value from the most recent prior hour. We record latitudes and longitudes of the measurement stations, from which we compute the matrix of geographical distances $G$ using the Haversine formula [116].

Mathematically, this data is a suitable opportunity to apply our $d_{p}$ metrics derived from $L^{p}$ distances between piecewise constant functions generated from change point algorithms. The change point procedure helps stabilize the data more appropriately, measuring average air quality over intervals of days rather than the huge peaks seen in the raw AQI time series. Since we expect air quality to remain somewhat locally stationary, these change point algorithms are more appropriate than more flexible but less targeted non-parametric models.

Within this multivariate context, we perform two separate analyses. First, our new $d_{p}$ distance matrices, computed according to the algorithmic framework of Section 2, enables the identification of similarity clusters and anomalies. Second, calculating the geographical distances between measuring sites gives a contextual distance matrix $G$. As described in Section 4, we may form associated consistency matrices. With each of these matrices, we apply hierarchical and spectral clustering, thereby identifying anomalies in the spread of poor air relative to space and time.

Henceforth, we set $p=1$.

We have proposed a new theoretical and practical framework for understanding equivalence and computing distances between time series. Mathematically, we aim to identify and quotient by the relation that two time series differ only up to noise, and computationally, we want to remove noise from calculations of distances $d_{p}$ and magnitudes $\operatorname{mag}$. Our procedure builds on earlier work to measure distances between time series via sets of structural breaks, and our new method produces more reliable and continuous measures. We have applied the function $\mathcal{F}_{\mu}$ and these metrics to analyze data with extraordinary peaks, determining piecewise constant functions $f$ designed to understand the true movement of the mean (or variance) throughout the time series. This procedure is highly flexible: we could swap variance for mean, vary our change point algorithm parameters, or use completely different change point algorithms altogether.

Several future implementation opportunities exist in this equivalence framework. One could develop and explore notions of weak equivalence if two time series are equivalent under $\mathcal{F}_{\mu}$ or $\mathcal{F}_{var}$, and strong equivalence if they are equivalent under both, or more strongly, if they have identical distributions within their change point segments. Moving beyond piecewise constant functions, two autoregressive progresses with the same parameters $\phi_{t}=A\phi_{t-1}+w_{t}$ and $\psi_{t}=A\psi_{t-1}+w^{\prime}_{t}$ could be declared equivalent, if $w_{t}$ and $w^{\prime}_{t}$ are each “white noise” variables, namely $N(0,\epsilon^{2})$ Gaussian distributions for sufficiently small $\epsilon$. Algorithms could be developed to empirically detect such similarity in the distributions from an observed sample set, define equivalence relations, and hence produce reliable distance measures on the induced quotient space. Developing further theoretical and algorithmic ways of determining quotient spaces equipped with metrics may have many applications in determining appropriately adjusted measures between time series. Within a statistical or nonlinear dynamics context, algorithmic frameworks designed to learn appropriate notions of equivalence “up to noise” would also be of great interest to the community. Finally, our distance measures could be incorporated with other advances in distance analysis between time series, such as the search for various different invariances (amplitude, phase, complexity, and others) as discussed at the end of Section 3.

In our analysis of the NSW bushfires, we have determined anomalies not just in the air quality data itself, but also those sites that stand out in the otherwise strong consistency between air quality and geographical distances. By using several different discrepancy matrices, we are able to uncover structure that one alone cannot reveal. Katoomba is a location of particular interest in our analysis. It was initially revealed as highly anomalous in our unscaled distance matrix $D^{us}$ but less so in the normalized matrix $D^{\text{norm}}$. The initial logical explanation for this is that the feature most anomalous with Katoomba was its considerable scale, reflected in the location’s overall magnitude, the highest among all measuring stations, which is then nullified when normalizing distances are computed in $D^{\text{norm}}$. However, curiously, Katoomba again appears as an outlier in $\text{Con}^{\text{norm}}$, which uses normalized distances to measure consistency with geographical distance. Examining Katoomba’s time series again (Figure 4(e)) we note that Katoomba experienced its greatest severity in air quality later on in our period of analysis, from mid-December onward, whereas the Sydney stations (which Katoomba is relatively close to) generally experienced severity prior to this.

There is also utility in examining the results of a range of clustering methods. For example, only spectral clustering identifies Oakdale as an anomalous element in its own cluster when considering unscaled distances’ consistency. We have also identified Narrabri as an interesting case. Unlike Katoomba, this location was never identified as a significant anomaly with respect to geography or AQI, considering distances or alignment. And yet, it is anomalous in the consistency between geographical and AQI distances. It is one of very few locations where these two measures are not closely aligned relative to other sites. Indeed, in Figure 4(c) we can see it closely resembles the profile of the Sydney sites despite being quite distant, to the north of the state. In addition, the holistic analysis of various distance and consistency matrices yields several high-level insights. First, there is broad consistency between AQI and geographical affinity matrices, which measure the relationship between sites with respect to these two aspects, as evidenced with the low consistency matrix norms in all implementations. Second, there is clearly sensitivity in the anomalies detected based on the distance measures that are used.

The mathematical analysis of relationships and anomalies between air quality measuring stations could prompt other researchers to conduct closer meteorological or geological examinations of different locations. Such research could shed light on the anomalous nature of Katoomba, Narrabri and other sites discussed in this manuscript. The fact that Sydney’s air quality was much more moderate after mid-December while nearby Katoomba’s became much more severe are too persistent with time to be explained by coincidences such as wind direction, and further exploration of this reasons could provide further insight on how bushfires and air pollution spread. This and other mathematical analysis could encourage researchers to notice hereto-unknown features of the land and its susceptibility or resistance to extreme bushfires, allowing better allocation of resources and understanding of vulnerability to bushfires. As such extreme events increase, we urge more interdisciplinary work between nonlinear dynamics and environmental researchers into such questions.

The data analyzed in this study are publicly available at [107]. A cached copy is conveniently available at https://github.com/MaxMenzies/BushfireData.