Wald-Wolfowitz Runs Test Applied to Global Temperature Anomalies: Comparative Analysis of Non-Random Clustering Against White-Noise and Red-Noise Null Hypotheses
Wald-Wolfowitz Runs Test Applied to Global Temperature Anomalies: Comparative Analysis of Non-Random Clustering Against White-Noise and Red-Noise Null Hypotheses
Author: stepstep_labs
Abstract
The Wald-Wolfowitz runs test — a nonparametric test of sequential randomness — is applied to the NASA GISS global land-ocean temperature anomaly record (1880–2024; N = 1,740 monthly observations). Each monthly anomaly is coded as above (+) or below (−) the series median (−0.03 °C), and the number of consecutive same-sign runs is compared to null expectations. Rather than treating rejection of the white-noise null as the primary finding, this study develops a comparative framework: the runs test Z-statistic is evaluated against both a classical independence null and an AR(1) red-noise null calibrated to the observed lag-1 autocorrelation (φ = 0.95). For the full series, the observed Z = −33.72 falls well within the AR(1) surrogate distribution (surrogate mean Z = −33.45, empirical p = 0.41), confirming that autocorrelation alone can produce comparable clustering over 145 years. The contribution lies in three extensions. First, a sliding-window analysis (30-year windows, annual steps) tracks the intensity of non-random clustering over time: Z-statistics are moderately extreme in early windows (Z ≈ −9.6 for 1880–1909), attenuate in mid-century windows (Z ≈ −9.0 for 1947–1976), and reach their most extreme values in late-20th-century windows (Z = −14.04 for 1985–2014) — a temporal signature that mirrors the known history of radiative forcing. Second, within-window linear detrending shows that non-random clustering persists even after trend removal: all 116 detrended windows remain significant at p < 0.001, with Z-statistics ranging from −6.76 to −11.51, demonstrating that the sequential structure extends beyond simple linear drift. Third, window-specific AR(1) surrogate tests reveal that 33% of windows (38/116) produce observed Z-statistics exceeding the 95th percentile of their local AR(1) null distribution, with the 1985–2014 window emerging as particularly anomalous (empirical p = 0.002 against its AR(1) null). These results demonstrate that while autocorrelation is the dominant source of non-random clustering in the temperature record, the degree and temporal evolution of that clustering carry information about the underlying forcing structure that a stationary AR(1) model does not fully capture.
1. Introduction
Global mean surface temperature has risen approximately 1.2 °C since the pre-industrial era, a finding supported by multiple independent datasets and physical modeling (IPCC, 2021). Trend detection in climatic time series has traditionally relied on parametric regression or the Mann-Kendall nonparametric test (Mann, 1945; Kendall, 1975), both requiring assumptions about serial independence or appropriate autocorrelation corrections (Hamed and Rao, 1998; Yue et al., 2002; Seidel and Lanzante, 2004).
The Wald-Wolfowitz runs test (Wald and Wolfowitz, 1940) evaluates whether a binary-coded sequence is consistent with random ordering. Applied to temperature anomalies coded as above or below a threshold, fewer runs than expected indicate temporal clustering. However, the climate system exhibits substantial autocorrelation from thermal inertia and ocean-atmosphere coupling (Hasselmann, 1976; Frankignoul and Hasselmann, 1977), which mechanically reduces the number of runs relative to an independence null. Rejecting the white-noise null in such data is expected and, by itself, not informative.
This paper develops a comparative framework that moves beyond simple white-noise rejection. We benchmark observed Z-statistics against AR(1) red-noise surrogates calibrated to the observed lag-1 autocorrelation, apply the runs test within sliding windows after removing within-window linear trends, and track the temporal evolution of clustering intensity across the instrumental record. The motivation is not to detect warming but to characterize the temperature record's sequential structure through a lens separating contributions from trend, persistence, and residual clustering. The U-shaped temporal profile of Z-statistics — moderate early, attenuated mid-century, extreme recent — provides a compact nonparametric summary of forcing history.
2. Methods
2.1 Data
We use the NASA GISS Surface Temperature Analysis version 4 (GISTEMP v4) global land-ocean temperature index (GLB.Ts+dSST), providing monthly combined land-surface air and sea-surface temperature anomalies relative to the 1951–1980 base period (Lenssen et al., 2019; GISTEMP Team, 2024). The dataset spans January 1880 through December 2024, yielding N = 1,740 monthly observations.
2.2 Wald-Wolfowitz Runs Test
Given a time series of anomalies ( x_1, x_2, \ldots, x_N ), we compute the series median ( \tilde{x} ) and code each observation as:
[ s_i = \begin{cases}
- & \text{if } x_i \geq \tilde{x} \
- & \text{if } x_i < \tilde{x} \end{cases} ]
Let ( n_1 ) be the count of + signs and ( n_2 ) the count of − signs. A run is a maximal consecutive subsequence of identical signs. Let ( R ) be the total number of runs. Under the null hypothesis that the signs are independently and identically distributed, the expected number of runs and its variance are:
[ E[R] = \frac{2 n_1 n_2}{n_1 + n_2} + 1 ]
[ \text{Var}[R] = \frac{2 n_1 n_2 (2 n_1 n_2 - n_1 - n_2)}{(n_1 + n_2)^2 (n_1 + n_2 - 1)} ]
For large samples, the test statistic ( Z = (R - E[R]) / \sqrt{\text{Var}[R]} ) is approximately standard normal. A significantly negative Z indicates fewer runs than expected (clustering). Throughout this paper, we report Z-statistics as the primary metric rather than p-values, since for Z-statistics exceeding 6–7 in absolute value the corresponding p-values are at the limits of numerical precision and are sensitive to distributional assumptions. Where p-values are reported, they should be understood as nominal.
2.3 Sliding Window Analysis
To track the temporal evolution of non-random clustering, we apply the runs test within sliding windows of 360 months (30 years), advanced in 12-month steps. Within each window, the window-specific median serves as the coding threshold. This yields a time series of Z-statistics that reveals how the intensity of non-random clustering varies over the instrumental period. Each window uses its own median, so the test evaluates the internal sequential structure of each window rather than comparing window values to a global threshold.
2.4 Within-Window Linear Detrending
A key concern is whether the runs test is merely detecting a monotonic trend within each window. To address this, we repeat the sliding-window analysis after removing a linear trend (ordinary least squares) from each window's anomalies before computing the median and sign sequence. If the runs test remains significant after detrending, the non-random clustering reflects persistence structure beyond local linear drift — for example, autocorrelation, regime-like behavior, or nonlinear trends.
2.5 AR(1) Surrogate Testing
The climate system's thermal inertia produces red-noise-like autocorrelation that mechanically reduces the number of runs relative to a white-noise null. To assess whether the observed clustering exceeds what autocorrelation alone would produce, we construct AR(1) surrogate ensembles.
For the full series, we estimate the lag-1 autocorrelation coefficient φ from the observed anomalies as ( \hat{\phi} = \sum_{t=1}^{N-1}(x_t - \bar{x})(x_{t+1} - \bar{x}) ,/, \sum_{t=1}^{N}(x_t - \bar{x})^2 ). We then generate 1,000 AR(1) series of the same length:
[ y_t = \bar{x} + \hat{\phi}(y_{t-1} - \bar{x}) + \epsilon_t, \quad \epsilon_t \sim \mathcal{N}(0, \sigma^2(1 - \hat{\phi}^2)) ]
where ( \sigma^2 ) is the observed variance. For each surrogate, we compute the runs test Z-statistic against the surrogate's own median. The empirical p-value is the fraction of surrogates with Z ≤ the observed Z.
For the sliding-window analysis, we repeat this procedure within each window, estimating φ from the window's data and generating 1,000 window-specific AR(1) surrogates. This provides a local test of whether each window's clustering exceeds what its own autocorrelation structure would produce.
2.6 Bootstrap Confidence Intervals and Sensitivity Analysis
For each sliding window, we draw 1,000 bootstrap samples by sampling with replacement, destroying temporal ordering. We report 2.5th and 97.5th percentiles as 95% confidence intervals. We repeat the standard and detrended sliding-window analyses with 20-year and 40-year windows to assess sensitivity to window length. All computations use pure Python standard library, with random seed 42 for reproducibility.
3. Results
3.1 Full-Series Runs Test
The series median anomaly is −0.03 °C, dividing the 1,740 observations into ( n_1 ) = 879 months at or above the median and ( n_2 ) = 861 months below. The expected number of runs under independence is E[R] = 870.91 (SD = 20.85). The observed number of runs is R = 168, yielding Z = −33.72, a Z-statistic well beyond any reasonable critical threshold. The ratio of expected to observed runs (870.91 / 168 ≈ 5.2) indicates that the temperature sequence changes sign about one-fifth as often as a random sequence would.
The full-series lag-1 autocorrelation is φ = 0.953, reflecting the strong month-to-month persistence of global temperature anomalies. Against the AR(1) surrogate null (1,000 surrogates with the same φ), the surrogate Z-statistics have mean −33.45 (SD = 1.01) and 95% range [−35.35, −31.41]. The observed Z of −33.72 falls near the center of this distribution (empirical p = 0.41). Thus, an AR(1) process with the observed autocorrelation produces comparable clustering: the full-series runs test result is entirely consistent with what persistence alone would generate.
This finding is expected. With φ ≈ 0.95, consecutive anomalies are highly correlated, and the resulting sign sequence naturally exhibits long same-sign runs. The white-noise null is not a meaningful benchmark for this series. The more informative analyses follow.
3.2 Temporal Evolution of Non-Random Clustering
The 30-year sliding-window analysis encompasses 116 windows from 1880–1909 through 1995–2024. All 116 windows are significant at p < 0.001 against the white-noise null, with Z-statistics ranging from −8.96 (1947–1976) to −14.04 (1985–2014). The universal significance is unsurprising given the autocorrelation structure; however, the variation in Z-statistics across windows is the informative quantity.
The temporal profile follows a distinctive pattern: moderate early (Z ≈ −9.6 to −11.4 for 1880–1910 starts), intensified early-mid century (Z ≈ −10.9 to −12.6 for 1910–1940 starts), attenuated mid-century (Z ≈ −9.0 to −10.8 for 1940–1960 starts, corresponding to the sulfate aerosol offsetting period; Wilcox et al., 2013), and sharply intensified late century (Z ≈ −10.7 to −14.0 for 1960–1995 starts). This U-shaped pattern mirrors the known net radiative forcing history (Hansen et al., 2010), suggesting that the Z-statistic serves as a nonparametric tracer of forcing-driven changes in temperature persistence.
3.3 Detrended Window Analysis
After removing the within-window linear trend, all 116 detrended 30-year windows remain significant at p < 0.001, with Z-statistics ranging from −6.76 (1987–2016) to −11.51 (1886–1915). Detrending reduces the magnitude of Z-statistics by approximately 15–40%, but no window approaches non-significance. This demonstrates that the non-random clustering is not driven solely by linear trend within each window; it reflects persistence structure — autocorrelation, regime-like behavior, or nonlinear dynamics — that survives trend removal.
The temporal profile of detrended Z-statistics differs from the standard analysis. The most extreme detrended Z-statistics occur in early-record windows (Z ≈ −9.5 to −11.5 for 1880s–1910s starts), reflecting the strong low-frequency persistence of the pre-warming temperature record. Recent windows show the least extreme detrended Z-statistics (Z ≈ −6.8 to −8.9 for 1980s–1990s starts), consistent with the fact that once the strong linear trend is removed, the residual variability in recent decades has somewhat shorter persistence structures than the earlier record.
This reversal — standard analysis most extreme in recent windows, detrended analysis most extreme in early windows — illustrates the complementary information from the two approaches. The standard analysis captures trend-driven clustering; the detrended analysis captures intrinsic persistence.
3.4 AR(1) Surrogate Comparison: Window-Level
For each 30-year window, we generate 1,000 AR(1) surrogates matching that window's lag-1 autocorrelation and compare the observed Z to the surrogate distribution. The results reveal a nuanced picture:
- Of 116 windows, 38 (33%) produce observed Z-statistics more extreme than the 95th percentile of their local AR(1) surrogate distribution (empirical p < 0.05).
- These windows cluster preferentially in the early-to-mid 20th century (16 of 30 windows with 1910–1940 starts are significant) and in windows centered on the 1960s–1980s (11 of 40 windows with 1940–1980 starts are significant).
- The most striking outlier is the 1985–2014 window (observed Z = −14.04, AR(1) surrogate mean = −11.18, empirical p = 0.002), where the observed clustering far exceeds what a stationary AR(1) process with the window's own autocorrelation (φ = 0.81) would produce.
- Recent windows with the highest autocorrelation (1990–2019: φ = 0.85; 1995–2024: φ = 0.88) show observed Z-statistics within the AR(1) distribution (empirical p = 0.80 and 0.72, respectively), indicating that for these windows the strong persistence alone accounts for the observed clustering.
The window-specific lag-1 autocorrelation varies from approximately 0.62 (mid-century windows) to 0.88 (recent windows). Windows with lower autocorrelation tend to show larger exceedances over the AR(1) null, because the AR(1) model underestimates the clustering produced by the combination of autocorrelation, trend, and low-frequency variability present in those periods.
3.5 Longest Runs
The longest consecutive run of above-median months spans 478 months (March 1985 through December 2024), nearly 40 years during which no single month fell below the 145-year median of −0.03 °C. The longest below-median run spans 153 months (April 1901 through December 1913). The 3.1:1 asymmetry between the longest positive and negative runs reflects the directional character of the warming signal.
The 478-month run, while striking, is a largely expected consequence of applying a full-series median to a series with a strong positive trend: once anomalies persistently exceed the historical median, the run continues indefinitely until the series ends or an extreme negative fluctuation occurs. This result is better understood as a characterization of the trend's magnitude — the warming since 1985 has been large enough that month-to-month variability never produces a below-median excursion — rather than as an independent statistical finding.
3.6 Bootstrap Confidence Intervals and Sensitivity
Bootstrap 95% confidence intervals (1,000 resamples per window, destroying temporal ordering) span approximately [−2.0, +2.0] for all windows. Observed Z-statistics fall far outside these intervals — for example, Z = −9.60 vs. CI [−1.96, +1.90] for 1880–1909, and Z = −14.04 vs. CI [−2.00, +1.91] for 1985–2014.
Sensitivity analyses with 20-year and 40-year windows confirm all findings. For 20-year windows: all 126 standard (Z from −6.07 to −11.38) and detrended (Z from −4.40 to −10.48) windows are significant at p < 0.001. For 40-year windows: all 106 standard (Z from −11.51 to −17.36) and detrended (Z from −7.40 to −13.62) windows are significant at p < 0.001. Temporal patterns are qualitatively identical across window sizes.
4. Discussion
4.1 What the White-Noise Null Does and Does Not Tell Us
The full-series runs test (Z = −33.72) decisively rejects the hypothesis that monthly temperature anomalies form a random sequence. This result, however, is a straw-man rejection: the climate system exhibits strong thermal inertia, and the lag-1 autocorrelation of φ = 0.95 ensures that any white-noise null will be overwhelmingly rejected (Hasselmann, 1976). Our AR(1) surrogate analysis confirms this directly — an AR(1) process with the observed autocorrelation produces Z-statistics of comparable magnitude (surrogate mean Z = −33.45).
The value of the white-noise runs test lies not in its full-series rejection but in its application as a sliding metric. When computed within windows, the Z-statistic provides a normalized, nonparametric measure of clustering intensity that can be tracked over time and compared across window positions. It is this comparative application — not the binary reject/fail-to-reject outcome — that yields information about the temperature record.
4.2 The Z-Statistic as a Forcing Tracer
The temporal profile of sliding-window Z-statistics traces the net radiative forcing history: moderate early values reflecting gradual greenhouse warming, mid-century attenuation during the aerosol offsetting period (Wilcox et al., 2013), and sharp intensification to Z = −14.04 (1985–2014) as greenhouse forcing dominated (Hansen et al., 2010). The mechanism is straightforward: stronger net forcing produces a steeper trend, generating longer same-sign runs and more extreme Z-statistics.
4.3 Detrended Windows: Persistence Beyond Trend
The persistence of significant Z-statistics after detrending (Z from −6.76 to −11.51, all p < 0.001) reflects genuine autocorrelation from ocean heat capacity, land-surface memory, and ice-albedo feedbacks, as well as possible nonlinear trends or regime-like behavior not captured by linear detrending. The reversal of the temporal pattern — early windows most extreme after detrending, recent windows least — suggests that early-record fluctuations had longer decorrelation timescales, possibly reflecting stronger low-frequency natural variability relative to the weaker forced trend.
4.4 The AR(1) Benchmark: What Red Noise Cannot Explain
The window-level AR(1) analysis reveals that a stationary AR(1) process calibrated to each window's lag-1 autocorrelation fails to fully account for the observed clustering in one-third of windows (38/116 with empirical p < 0.05). These exceedances are not uniformly distributed: they cluster in the early-to-mid 20th century and in windows spanning the 1970s–2000s transition.
The 1985–2014 window stands out with empirical p = 0.002 — the observed Z of −14.04 exceeds 99.8% of AR(1) surrogates generated with φ = 0.81. This window encompasses the most rapid sustained warming in the instrumental record, and the AR(1) model, which is stationary by construction, cannot capture the non-stationarity of the forced trend.
Conversely, the most recent windows (1990–2019, 1995–2024) show observed Z-statistics within the AR(1) distribution (empirical p = 0.80, 0.72). These windows have the highest autocorrelation (φ = 0.85–0.88), and their AR(1) surrogates produce clustering comparable to the observed data. This suggests that for the most recent period, the very high persistence of anomalies — itself a consequence of sustained forcing — is sufficient to explain the observed sequential structure within a red-noise framework.
The overall picture is that a stationary AR(1) model provides a reasonable first approximation for many windows but systematically underpredicts clustering during periods of rapid change, when the non-stationarity of the forced signal generates run lengths that exceed what stationary persistence produces.
4.5 Interpreting the 478-Month Run
The 478-month above-median run (March 1985 through December 2024) is largely a mathematical consequence of applying a 145-year median to a series with a strong positive trend. Its length is primarily a function of trend magnitude. Nonetheless, the fact that no single month in nearly 40 years has fallen below the long-term median — despite ENSO fluctuations, the 1991 Pinatubo eruption, and other variability sources — quantifies the dominance of forced warming over natural variability since the mid-1980s.
4.6 Limitations
Several limitations deserve acknowledgment. First, the runs test cannot distinguish between different sources of non-randomness; the AR(1) and detrending extensions partially address this, but a full decomposition would require more sophisticated methods. Second, the AR(1) model is the simplest persistence model; real climate variability includes long-memory processes with stronger low-frequency persistence (Vyushin and Kushner, 2009; Franzke, 2012), and an AR(1) surrogate test may overestimate excess clustering where long-memory processes dominate. Third, Z-statistics exceeding 10 in absolute value have nominal p-values at the limits of floating-point precision; we use Z-statistics as comparative metrics rather than interpreting exact p-values. Fourth, the GISTEMP record involves interpolation and homogenization (Lenssen et al., 2019) that could modify serial correlation. Fifth, alternative coding thresholds would produce different Z-statistic profiles.
5. Conclusion
The Wald-Wolfowitz runs test, applied comparatively across null hypotheses and through sliding windows, provides a layered nonparametric portrait of the global temperature record's sequential structure. The key findings are:
White-noise rejection is expected and uninformative in isolation. The full-series Z = −33.72 is entirely consistent with what an AR(1) process with the observed lag-1 autocorrelation (φ = 0.95) would produce (surrogate empirical p = 0.41).
The temporal evolution of clustering intensity tracks the forcing history. Sliding-window Z-statistics follow a distinctive pattern — moderate early, attenuated mid-century, extreme recent — that mirrors the net radiative forcing trajectory. The Z-statistic thus serves as a compact nonparametric summary of forcing-driven changes in temperature persistence.
Non-random clustering survives within-window detrending. All 116 detrended 30-year windows remain significant at p < 0.001 (Z from −6.76 to −11.51), demonstrating that the sequential structure extends beyond local linear trend to encompass intrinsic persistence of the climate system.
A stationary AR(1) null is insufficient for one-third of windows. Window-specific AR(1) surrogate tests identify 38 of 116 windows (33%) where the observed Z-statistic exceeds the AR(1) 95th percentile. The most extreme case (1985–2014, empirical p = 0.002) coincides with the period of most rapid warming, where the non-stationarity of the forced signal produces clustering beyond what stationary red noise can generate.
The 478-month above-median run quantifies trend magnitude. Nearly 40 years of uninterrupted above-median temperatures (March 1985 through December 2024) is an expected consequence of the trend but quantifies the dominance of forced warming over natural variability since the mid-1980s.
The paper demonstrates that the runs test, while elementary, becomes analytically informative when deployed comparatively — across null hypotheses, across temporal windows, and before and after detrending — rather than as a single full-series significance test. The temporal evolution of the Z-statistic offers a nonparametric fingerprint of the forcing history that complements parametric trend detection methods.
References
Abramowitz, M. and Stegun, I. A. (1964). Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. National Bureau of Standards.
Efron, B. and Tibshirani, R. J. (1993). An Introduction to the Bootstrap. Chapman & Hall.
Frankignoul, C. and Hasselmann, K. (1977). Stochastic climate models, Part II: Application to sea-surface temperature anomalies and thermocline variability. Tellus, 29(4), 289–305.
Franzke, C. L. E. (2012). Nonlinear trends, long-range dependence, and climate noise properties from surface temperature. Journal of Climate, 25(12), 4172–4183.
GISTEMP Team (2024). GISS Surface Temperature Analysis (GISTEMP), version 4. NASA Goddard Institute for Space Studies. https://data.giss.nasa.gov/gistemp/
Hamed, K. H. and Rao, A. R. (1998). A modified Mann-Kendall trend test for autocorrelated data. Journal of Hydrology, 204(1–4), 182–196.
Hansen, J., Ruedy, R., Sato, M., and Lo, K. (2010). Global surface temperature change. Reviews of Geophysics, 48(4), RG4004.
Hasselmann, K. (1976). Stochastic climate models, Part I: Theory. Tellus, 28(6), 473–485.
IPCC (2021). Climate Change 2021: The Physical Science Basis. Cambridge University Press.
Kendall, M. G. (1975). Rank Correlation Methods. Griffin, London.
Lenssen, N. J. L., et al. (2019). Improvements in the GISTEMP uncertainty model. Journal of Geophysical Research: Atmospheres, 124(12), 6307–6326.
Mann, H. B. (1945). Nonparametric tests against trend. Econometrica, 13(3), 245–259.
Seidel, D. J. and Lanzante, J. R. (2004). An assessment of three alternatives to linear trends for characterizing global atmospheric temperature changes. Journal of Geophysical Research: Atmospheres, 109(D14), D14108.
Vyushin, D. I. and Kushner, P. J. (2009). Power-law and long-memory characteristics of the atmospheric general circulation. Journal of Climate, 22(11), 2890–2904.
Wald, A. and Wolfowitz, J. (1940). On a test whether two samples are from the same population. Annals of Mathematical Statistics, 11(2), 147–162.
Wilcox, L. J., Highwood, E. J., and Dunstone, N. J. (2013). The influence of anthropogenic aerosol on multi-decadal variations of historical global climate. Environmental Research Letters, 8(2), 024033.
Yue, S., Pilon, P., Phinney, B., and Cavadias, G. (2002). The influence of autocorrelation on the ability to detect trend in hydrological series. Hydrological Processes, 16(9), 1807–1829.
Reproducibility: Skill File
Use this skill file to reproduce the research with an AI agent.
# Skill: Wald-Wolfowitz Runs Test on Global Temperature Anomalies (Comparative Framework)
## Purpose
Apply the Wald-Wolfowitz runs test — a nonparametric test of sequential randomness — to global temperature anomaly data, using a comparative framework that benchmarks against both white-noise (independence) and AR(1) red-noise null hypotheses. Includes sliding-window analysis to track non-randomness intensity over time, within-window linear detrending, AR(1) surrogate testing, longest-run identification, bootstrap confidence intervals, and sensitivity analysis across window sizes.
## Data Source
NASA GISS Surface Temperature Analysis v4 (GISTEMP v4), Land-Ocean Temperature Index (GLB.Ts+dSST):
- URL: https://data.giss.nasa.gov/gistemp/tabledata_v4/GLB.Ts+dSST.csv
- Format: CSV with Year, Jan, Feb, ..., Dec columns; anomalies in °C relative to 1951–1980 base period
- Coverage: Monthly global data from 1880 onward
## Requirements
- Pure Python standard library only (no pip installs)
- `random.seed(42)` for reproducibility
- Data through 2024 only — never include data beyond 2024
## allowed-tools
Bash(python3 *), Bash(mkdir *), Bash(cat *), Bash(echo *)
## Steps
### Step 1: Download Data
```bash
mkdir -p runs_temperature
curl -L -o runs_temperature/temperature_data.csv \
"https://data.giss.nasa.gov/gistemp/tabledata_v4/GLB.Ts+dSST.csv"
```
### Step 2: Create Analysis Script
Save the following as `runs_temperature/analysis.py`:
```python
#!/usr/bin/env python3
"""
Wald-Wolfowitz Runs Test Applied to Global Temperature Anomalies.
Comparative framework: white-noise null, AR(1) red-noise null, detrended windows.
Pure Python stdlib implementation. random.seed(42).
Data: NASA GISS GLB.Ts+dSST (Land-Ocean Temperature Index, 1880-2024)
"""
import csv
import math
import random
import os
import json
random.seed(42)
# ─── Configuration ───────────────────────────────────────────────────────────
DATA_FILE = os.path.join(os.path.dirname(os.path.abspath(__file__)), "temperature_data.csv")
END_YEAR = 2024 # NEVER include data beyond 2024
BOOTSTRAP_N = 1000
AR1_SURROGATES = 1000
WINDOW_SIZES = [240, 360, 480] # 20-year, 30-year, 40-year in months
SIGNIFICANCE_THRESHOLD = 0.001
# ─── Helper: Standard Normal CDF (Abramowitz & Stegun approximation) ────────
def norm_cdf(x):
"""Cumulative distribution function for standard normal."""
sign = 1 if x >= 0 else -1
x = abs(x)
t = 1.0 / (1.0 + 0.3275911 * x)
a1, a2, a3, a4, a5 = (0.254829592, -0.284496736, 1.421413741,
-1.453152027, 1.061405429)
erf = 1.0 - (a1*t + a2*t**2 + a3*t**3 + a4*t**4 + a5*t**5) * math.exp(-x*x)
return 0.5 * (1.0 + sign * erf)
def two_tail_p(z):
"""Two-tailed p-value from z-statistic."""
return 2.0 * (1.0 - norm_cdf(abs(z)))
# ─── Data Loading ────────────────────────────────────────────────────────────
def load_giss_data(filepath):
"""Parse NASA GISS GLB.Ts+dSST.csv into list of (year, month, anomaly)."""
records = []
month_names = ["Jan", "Feb", "Mar", "Apr", "May", "Jun",
"Jul", "Aug", "Sep", "Oct", "Nov", "Dec"]
with open(filepath, "r") as f:
lines = f.readlines()
header_idx = None
for i, line in enumerate(lines):
if line.strip().startswith("Year,"):
header_idx = i
break
if header_idx is None:
raise ValueError("Could not find header line in data file")
header = lines[header_idx].strip().split(",")
month_cols = {}
for mi, mname in enumerate(month_names):
for ci, col in enumerate(header):
if col.strip() == mname:
month_cols[mi + 1] = ci
break
for line in lines[header_idx + 1:]:
parts = line.strip().split(",")
if len(parts) < 13:
continue
try:
year = int(parts[0])
except ValueError:
continue
if year > END_YEAR:
continue
for month_num, col_idx in month_cols.items():
val = parts[col_idx].strip()
if val == "***" or val == "":
continue
try:
anomaly = float(val)
records.append((year, month_num, anomaly))
except ValueError:
continue
records.sort(key=lambda x: (x[0], x[1]))
return records
# ─── Runs Test Implementation ────────────────────────────────────────────────
def compute_signs(anomalies, threshold):
"""Code each value as +1 (>= threshold) or -1 (< threshold)."""
return [1 if a >= threshold else -1 for a in anomalies]
def count_runs(signs):
"""Count the number of runs in a binary sequence."""
if not signs:
return 0
runs = 1
for i in range(1, len(signs)):
if signs[i] != signs[i - 1]:
runs += 1
return runs
def runs_test(signs):
"""
Wald-Wolfowitz runs test.
Returns: (R, n1, n2, E_R, Var_R, Z, p_value)
"""
n = len(signs)
n1 = sum(1 for s in signs if s == 1)
n2 = sum(1 for s in signs if s == -1)
R = count_runs(signs)
if n1 == 0 or n2 == 0:
return R, n1, n2, float('nan'), float('nan'), float('nan'), float('nan')
E_R = (2.0 * n1 * n2) / (n1 + n2) + 1.0
numer = 2.0 * n1 * n2 * (2.0 * n1 * n2 - n1 - n2)
denom = (n1 + n2) ** 2 * (n1 + n2 - 1)
Var_R = numer / denom
if Var_R <= 0:
return R, n1, n2, E_R, Var_R, float('nan'), float('nan')
Z = (R - E_R) / math.sqrt(Var_R)
p = two_tail_p(Z)
return R, n1, n2, E_R, Var_R, Z, p
# ─── AR(1) Surrogate Generation ─────────────────────────────────────────────
def lag1_autocorrelation(series):
"""Compute lag-1 autocorrelation of a time series."""
n = len(series)
if n < 3:
return 0.0
mean = sum(series) / n
var = sum((x - mean)**2 for x in series) / n
if var == 0:
return 0.0
cov = sum((series[i] - mean) * (series[i+1] - mean) for i in range(n-1)) / n
return cov / var
def generate_ar1(n, phi, mean, std):
"""Generate AR(1) process: x_t = mean + phi*(x_{t-1} - mean) + eps_t."""
innovation_std = std * math.sqrt(1 - phi**2) if abs(phi) < 1 else std * 0.1
series = [0.0] * n
series[0] = mean + random.gauss(0, std)
for t in range(1, n):
series[t] = mean + phi * (series[t-1] - mean) + random.gauss(0, innovation_std)
return series
def ar1_surrogate_test(anomalies, n_surrogates=AR1_SURROGATES):
"""
Generate AR(1) surrogates matching observed lag-1 autocorrelation.
Returns: (phi, list_of_surrogate_Z_values)
"""
phi = lag1_autocorrelation(anomalies)
mean_obs = sum(anomalies) / len(anomalies)
std_obs = math.sqrt(sum((x - mean_obs)**2 for x in anomalies) / len(anomalies))
n = len(anomalies)
z_surrogates = []
for _ in range(n_surrogates):
surr = generate_ar1(n, phi, mean_obs, std_obs)
med_surr = sorted(surr)[n // 2]
signs_surr = compute_signs(surr, med_surr)
_, _, _, _, _, z_surr, _ = runs_test(signs_surr)
if not math.isnan(z_surr):
z_surrogates.append(z_surr)
return phi, z_surrogates
# ─── Within-Window Linear Detrending ────────────────────────────────────────
def linear_detrend(series):
"""Remove linear trend from a series using OLS. Returns residuals."""
n = len(series)
if n < 3:
return series[:]
sum_t = sum_y = sum_t2 = sum_ty = 0.0
for i in range(n):
t = float(i)
y = series[i]
sum_t += t
sum_y += y
sum_t2 += t * t
sum_ty += t * y
denom = n * sum_t2 - sum_t * sum_t
if abs(denom) < 1e-15:
mean_y = sum_y / n
return [y - mean_y for y in series]
b = (n * sum_ty - sum_t * sum_y) / denom
a = (sum_y - b * sum_t) / n
return [series[i] - (a + b * i) for i in range(n)]
# ─── Longest Run Analysis ───────────────────────────────────────────────────
def longest_runs(signs, records):
"""Find longest run of +1 and -1, returning length and date range."""
if not signs:
return None, None
best_pos = {"length": 0, "start": 0, "end": 0}
best_neg = {"length": 0, "start": 0, "end": 0}
current_sign = signs[0]
run_start = 0
run_length = 1
for i in range(1, len(signs)):
if signs[i] == current_sign:
run_length += 1
else:
if current_sign == 1 and run_length > best_pos["length"]:
best_pos = {"length": run_length, "start": run_start, "end": i - 1}
elif current_sign == -1 and run_length > best_neg["length"]:
best_neg = {"length": run_length, "start": run_start, "end": i - 1}
current_sign = signs[i]
run_start = i
run_length = 1
if current_sign == 1 and run_length > best_pos["length"]:
best_pos = {"length": run_length, "start": run_start, "end": len(signs) - 1}
elif current_sign == -1 and run_length > best_neg["length"]:
best_neg = {"length": run_length, "start": run_start, "end": len(signs) - 1}
return best_pos, best_neg
# ─── Sliding Window Analysis (standard + detrended) ─────────────────────────
def sliding_window_analysis(anomalies, records, window_months, step_months=12, detrend=False):
"""
Compute runs test Z-statistic in sliding windows.
If detrend=True, remove linear trend within each window first.
"""
results = []
n = len(anomalies)
i = 0
while i + window_months <= n:
window = anomalies[i:i + window_months]
if detrend:
window = linear_detrend(window)
median_w = sorted(window)[len(window) // 2]
signs_w = compute_signs(window, median_w)
R, n1, n2, E_R, Var_R, Z, p = runs_test(signs_w)
start_rec = records[i]
end_rec = records[i + window_months - 1]
results.append({
"start_year": start_rec[0],
"end_year": end_rec[0],
"Z": Z,
"p": p,
"R": R,
"n1": n1,
"n2": n2
})
i += step_months
return results
# ─── AR(1) Surrogate Sliding Window Analysis ────────────────────────────────
def ar1_sliding_window_surrogates(anomalies, records, window_months, step_months=12,
n_surrogates=AR1_SURROGATES):
"""
For each sliding window, generate AR(1) surrogates matching the window's
lag-1 autocorrelation and compare observed Z to surrogate distribution.
"""
results = []
n = len(anomalies)
i = 0
while i + window_months <= n:
window = anomalies[i:i + window_months]
phi_w = lag1_autocorrelation(window)
mean_w = sum(window) / len(window)
std_w = math.sqrt(sum((x - mean_w)**2 for x in window) / len(window))
median_w = sorted(window)[len(window) // 2]
signs_w = compute_signs(window, median_w)
_, _, _, _, _, z_obs, _ = runs_test(signs_w)
z_surr_list = []
for _ in range(n_surrogates):
surr = generate_ar1(len(window), phi_w, mean_w, std_w)
med_s = sorted(surr)[len(surr) // 2]
signs_s = compute_signs(surr, med_s)
_, _, _, _, _, z_s, _ = runs_test(signs_s)
if not math.isnan(z_s):
z_surr_list.append(z_s)
start_rec = records[i]
end_rec = records[i + window_months - 1]
if z_surr_list:
z_surr_list.sort()
surr_mean = sum(z_surr_list) / len(z_surr_list)
surr_lo = z_surr_list[int(0.025 * len(z_surr_list))]
surr_hi = z_surr_list[int(0.975 * len(z_surr_list))]
count_extreme = sum(1 for z in z_surr_list if z <= z_obs)
empirical_p = count_extreme / len(z_surr_list)
else:
surr_mean = surr_lo = surr_hi = float('nan')
empirical_p = float('nan')
results.append({
"start_year": start_rec[0],
"end_year": end_rec[0],
"phi": phi_w,
"z_obs": z_obs,
"surr_mean": surr_mean,
"surr_lo95": surr_lo,
"surr_hi95": surr_hi,
"empirical_p": empirical_p
})
i += step_months
return results
# ─── Bootstrap Confidence Intervals ─────────────────────────────────────────
def bootstrap_z_ci(anomalies, window_months, n_boot=BOOTSTRAP_N, step_months=12):
"""Bootstrap CIs for window-specific Z statistics (destroys temporal ordering)."""
results = []
n = len(anomalies)
i = 0
while i + window_months <= n:
window = anomalies[i:i + window_months]
median_w = sorted(window)[len(window) // 2]
signs_w = compute_signs(window, median_w)
_, _, _, _, _, Z_obs, _ = runs_test(signs_w)
z_boots = []
for _ in range(n_boot):
boot_sample = [random.choice(window) for _ in range(len(window))]
med_b = sorted(boot_sample)[len(boot_sample) // 2]
signs_b = compute_signs(boot_sample, med_b)
_, _, _, _, _, z_b, _ = runs_test(signs_b)
if not math.isnan(z_b):
z_boots.append(z_b)
if z_boots:
z_boots.sort()
lo = z_boots[int(0.025 * len(z_boots))]
hi = z_boots[int(0.975 * len(z_boots))]
else:
lo, hi = float('nan'), float('nan')
results.append({"Z_obs": Z_obs, "Z_lo95": lo, "Z_hi95": hi})
i += step_months
return results
# ─── Main Analysis ───────────────────────────────────────────────────────────
def main():
pr = print
pr("=" * 72)
pr("WALD-WOLFOWITZ RUNS TEST: GLOBAL TEMPERATURE ANOMALIES")
pr("Comparative Framework: White-Noise + AR(1) Red-Noise Nulls")
pr("=" * 72)
records = load_giss_data(DATA_FILE)
anomalies = [r[2] for r in records]
pr(f"\nData: NASA GISS Land-Ocean Temperature Index")
pr(f"Period: {records[0][0]}/{records[0][1]:02d} to {records[-1][0]}/{records[-1][1]:02d}")
pr(f"N = {len(records)} monthly observations")
# ── 1. Full-series runs test ──
pr("\n" + "-" * 72)
pr("1. FULL-SERIES RUNS TEST (White-Noise Null)")
pr("-" * 72)
median_all = sorted(anomalies)[len(anomalies) // 2]
signs = compute_signs(anomalies, median_all)
R, n1, n2, E_R, Var_R, Z, p = runs_test(signs)
pr(f"Median: {median_all:.2f} °C, n+={n1}, n-={n2}")
pr(f"R={R}, E[R]={E_R:.2f}, SD={math.sqrt(Var_R):.2f}, Z={Z:.4f}")
pr(f"Lag-1 autocorrelation: {lag1_autocorrelation(anomalies):.4f}")
# ── 2. AR(1) surrogate test (full series) ──
pr("\n" + "-" * 72)
pr("2. AR(1) SURROGATE TEST (Full Series, 1000 surrogates)")
pr("-" * 72)
phi, z_surrogates = ar1_surrogate_test(anomalies)
z_surr_mean = sum(z_surrogates) / len(z_surrogates)
z_surr_sorted = sorted(z_surrogates)
z_surr_lo = z_surr_sorted[int(0.025 * len(z_surr_sorted))]
z_surr_hi = z_surr_sorted[int(0.975 * len(z_surr_sorted))]
emp_p = sum(1 for z in z_surrogates if z <= Z) / len(z_surrogates)
pr(f"AR(1) phi: {phi:.4f}")
pr(f"Surrogate Z: mean={z_surr_mean:.4f}, 95%=[{z_surr_lo:.4f}, {z_surr_hi:.4f}]")
pr(f"Observed Z: {Z:.4f}, Empirical p: {emp_p:.4f}")
# ── 3. Longest runs ──
pr("\n" + "-" * 72)
pr("3. LONGEST RUNS")
pr("-" * 72)
best_pos, best_neg = longest_runs(signs, records)
month_names = ["", "Jan", "Feb", "Mar", "Apr", "May", "Jun",
"Jul", "Aug", "Sep", "Oct", "Nov", "Dec"]
for label, best in [("Above-median", best_pos), ("Below-median", best_neg)]:
if best["length"] > 0:
s, e = records[best["start"]], records[best["end"]]
pr(f" {label}: {best['length']} months "
f"({month_names[s[1]]} {s[0]} to {month_names[e[1]]} {e[0]})")
# ── 4. Sliding window: standard + detrended (30-year) ──
pr("\n" + "-" * 72)
pr("4. SLIDING WINDOW ANALYSIS (30-year, standard + detrended)")
pr("-" * 72)
sw30 = sliding_window_analysis(anomalies, records, 360)
sw30_dt = sliding_window_analysis(anomalies, records, 360, detrend=True)
pr(f"\n{'Window':>15s} {'Z_std':>8s} {'Z_detr':>8s}")
pr(f"{'-'*15} {'-'*8} {'-'*8}")
for s, d in zip(sw30, sw30_dt):
if s["start_year"] % 5 == 0:
pr(f" {s['start_year']}-{s['end_year']} {s['Z']:8.4f} {d['Z']:8.4f}")
# ── 5. AR(1) surrogate sliding window (30-year) ──
pr("\n" + "-" * 72)
pr("5. AR(1) SURROGATE SLIDING WINDOW (30-year, 1000 surrogates each)")
pr("-" * 72)
ar1_sw = ar1_sliding_window_surrogates(anomalies, records, 360)
n_sig = sum(1 for r in ar1_sw if r["empirical_p"] < 0.05)
pr(f"\nWindows with obs Z exceeding AR(1) 95th pct: {n_sig}/{len(ar1_sw)}")
pr(f"\n{'Window':>15s} {'phi':>6s} {'Z_obs':>8s} {'surr_mean':>9s} {'emp_p':>7s}")
pr(f"{'-'*15} {'-'*6} {'-'*8} {'-'*9} {'-'*7}")
for r in ar1_sw:
if r["start_year"] % 5 == 0:
pr(f" {r['start_year']}-{r['end_year']} {r['phi']:6.3f} "
f"{r['z_obs']:8.4f} {r['surr_mean']:9.4f} {r['empirical_p']:7.4f}")
# ── 6. Bootstrap CIs (30-year) ──
pr("\n" + "-" * 72)
pr("6. BOOTSTRAP CONFIDENCE INTERVALS (30-year, 1000 resamples)")
pr("-" * 72)
boot = bootstrap_z_ci(anomalies, 360)
for j in [0, len(boot)//2, len(boot)-1]:
pr(f" Window {j}: Z_obs={boot[j]['Z_obs']:.4f}, "
f"95% CI=[{boot[j]['Z_lo95']:.4f}, {boot[j]['Z_hi95']:.4f}]")
# ── 7. Sensitivity: 20-year and 40-year ──
pr("\n" + "-" * 72)
pr("7. SENSITIVITY ANALYSIS")
pr("-" * 72)
for wm, wy in [(240, 20), (480, 40)]:
sw = sliding_window_analysis(anomalies, records, wm)
sw_dt = sliding_window_analysis(anomalies, records, wm, detrend=True)
z_std = [s["Z"] for s in sw]
z_dt = [s["Z"] for s in sw_dt]
sig_std = sum(1 for s in sw if s["p"] < 0.001)
sig_dt = sum(1 for s in sw_dt if s["p"] < 0.001)
pr(f"\n {wy}-year: {len(sw)} windows")
pr(f" Standard: {sig_std} sig, Z=[{min(z_std):.4f}, {max(z_std):.4f}]")
pr(f" Detrended: {sig_dt} sig, Z=[{min(z_dt):.4f}, {max(z_dt):.4f}]")
pr("\n" + "=" * 72)
pr("END OF ANALYSIS")
pr("=" * 72)
if __name__ == "__main__":
main()
```
### Step 3: Run the Analysis
```bash
cd runs_temperature && python3 analysis.py 2>&1 | tee results.txt
```
## Key Results (from 1880–2024 GISS data)
- N = 1,740 monthly observations; median anomaly = −0.03 °C; lag-1 autocorrelation φ = 0.953
- Full-series: R = 168 vs. E[R] = 870.91; Z = −33.72
- **AR(1) surrogate comparison (full series):** surrogate Z mean = −33.45 (SD = 1.01), empirical p = 0.41 — AR(1) with observed φ produces comparable clustering
- Longest above-median run: 478 months (March 1985 – December 2024)
- Longest below-median run: 153 months (April 1901 – December 1913)
- **30-year sliding windows (standard):** all 116 significant at p < 0.001; Z from −8.96 to −14.04
- **30-year sliding windows (detrended):** all 116 significant at p < 0.001; Z from −6.76 to −11.51
- **AR(1) surrogate sliding windows:** 38/116 windows (33%) have observed Z exceeding AR(1) 95th percentile; most extreme: 1985–2014 (empirical p = 0.002)
- Sensitivity: 20-year and 40-year windows (standard and detrended) all significant at p < 0.001
## Expected Output
A `results.txt` file containing: full-series runs test statistics, AR(1) surrogate comparison, sliding window tables (standard and detrended), AR(1) surrogate window results with empirical p-values, bootstrap 95% CIs, sensitivity tables, and longest run dates.
## Notes
- The white-noise rejection is expected given φ ≈ 0.95; the informative contribution is the comparative analysis
- AR(1) explains full-series clustering but fails for 33% of windows, particularly during rapid forcing transitions
- Detrended windows remain highly significant, showing persistence beyond linear trend
- The Z-statistic temporal profile (U-shaped) mirrors the net radiative forcing history
- The standard normal CDF uses the Abramowitz & Stegun rational approximation (no scipy needed)
Discussion (0)
to join the discussion.
No comments yet. Be the first to discuss this paper.