The paper

Methods

The plan, code and parameters were registered before any full-pipeline null outcome was computed (prereg:f71f0995…). plan/ is byte-identical to the registration. code/analyze.py, code/fr.py and code/tmax.py are identical to the registered copies. deviations.json lists the one addition, code/summarize.py, which rounds the registered outputs for display.

Data. We used Foster and Rahmstorf's own inputs and outputs, fetched from Zenodo and checked by SHA-256 (data/external.json, CC BY 4.0):

• the code and data archive, whose ObservedData.csv holds monthly global mean surface temperature for 1880–2024 from NASA GISTEMP v4, NOAA GlobalTemp v6, HadCRUT5, Berkeley Earth and ERA5 (from 1940), together with stratospheric aerosol optical depth, the NINO3.4 index and the sunspot number;
• their adjusted series.

Adjustment. code/fr.py ports their R function Xogen.r to numpy:

• a lowess smooth spanning 20 years is subtracted from the temperatures;
• the residuals are regressed on the volcanic, NINO3.4 and sunspot series, lagged 0–24, 0–12 and 0–12 months, choosing the lag combination with the lowest AIC among 4,225, fitted over 1950–2024;
• this is followed by four rounds of backfitting, in which the lowess is refitted to the adjusted series.

Before registration, we ran their R code and the port on all five records. The port chose the same lags and gave the same adjusted values to the 4 decimals Xogen.r rounds to. The bundle checks this again against the published adjusted series (C1).

Changepoint test. code/tmax.py follows Beaulieu et al. (2024), whose test Foster and Rahmstorf apply:

• Data: annual means for 1970–2024 (N = 55).
• Model: for each changepoint k from 6 to 49 (10% trimming), a continuous piecewise-linear trend with Gaussian AR(1) errors, fitted by exact maximum likelihood profiled over the AR coefficient φ. Then T_k is the slope change divided by its generalized-least-squares standard error, and T_max = max |T_k|.
• Agreement with R: T_k agrees with R's arima(..., method = "ML") to within 0.02 on the series checked. It gives T_43 = 1.575 for HadCRUT 1970–2023, where Beaulieu et al. report 1.528 from their own copy of the record.
• Critical values: these come from AR(1) Monte Carlo at the φ fitted to a straight line plus AR(1). For the observed series that means 100,000 simulations at the fitted φ. For simulated series it means a table of 10,000 simulations at each φ from −0.40 to 0.80 in steps of 0.02, interpolated linearly.

Quadratic test. This is ordinary least squares of a quadratic in time, with the standard error of the quadratic coefficient multiplied by √ν, where ν = 1 + 2r₁/(1 − r₂/r₁) is computed from the residuals' lag-1 and lag-2 autocorrelations: the ARMA(1,1) correction of Foster and Rahmstorf (2011). It is applied to monthly values for 1970–2024 and separately to annual means. p-values are two-sided normal.

Null world. For each record, simulated monthly temperatures are g₀(t) + E(t) + ε(t):

• g₀(t) is a straight line fitted to the adjusted series for 1970–2024, joined before 1970 by that series' 20-year lowess, shifted to meet the line;
• E(t) is the volcanic, ENSO and solar terms with the coefficients and lags the procedure fits to the real record, applied to the real predictor series;
• ε(t) is stationary Gaussian noise, from one of three models fitted to the real adjusted series:
  • primary, null_trend: ARMA(1,1) by maximum likelihood, fitted to the adjusted series minus g₀ over 1970–2024. This is the noise as estimated under the null hypothesis, as a parametric test requires;
  • lowess: ARMA(1,1) fitted to the adjusted series minus its 20-year lowess over 1950–2024. This has less low-frequency variance, which favours the original finding;
  • long_memory: ARFIMA(1, d, 0) fitted to the same residuals as the primary model. The fractional order d comes from the log-periodogram regression of Geweke and Porter-Hudak (1983) over the lowest ⌊√n⌋ frequencies, clipped to [0, 0.45], with an AR(1) fitted to the fractionally differenced series.

Each simulated record is 1,200 months of burn-in plus the real record's length. It goes through the whole procedure, the adjustment re-estimated including the lag search and backfitting, and then both tests.

Outcomes and decision rules. We ran 2,000 simulated records for each record and noise model, all from seed 20261010, split into independent streams.

• Size: the fraction of null records whose nominal p falls below 0.05 or 0.01.
• Full-pipeline p: (1 + the number of null statistics at least as large as the observed one) / 2,001. The smallest possible value is 0.0005.
• Changepoint rule: "robust" if the full-pipeline p is below 0.01 in all five records, "not robust" if it is 0.05 or above in at least three, and otherwise "mixed".
• Quadratic rule: the same, with 0.02 in place of 0.01, matching the reported confidence of better than 98%.

The sensitivity noise models do not enter the decisions. sh code/run reproduces everything on python:3.12-slim with env/requirements.txt, offline, in about 15 minutes on one CPU core.

Results

Reproduction. The port reproduces the published adjusted series of all 5 records to within 0.0001 °C. Table 1 gives both tests on the adjusted and unadjusted series. Every adjusted record has a changepoint statistic beyond its nominal critical value at level 0.01. Among the unadjusted records, NOAA reaches nominal p = 0.041 and NASA 0.053, while the other three stay above 0.05.

Table 1. Both tests on the observed series, annual means for 1970–2024; p is the nominal autoregressive Monte Carlo value for the changepoint test, and the critical value is its nominal one at the stricter level (Methods); z is the corrected quadratic statistic.

Full-pipeline null, primary noise model. Table 2 gives the results from 2000 simulated records per dataset.

• Size: the changepoint test at nominal level 0.05 rejected 0.059 to 0.07 of null records, and at nominal level 0.01 it rejected 0.011 to 0.018.
• Critical values: the full-pipeline critical value of T_max at level 0.01 was 4.01 to 4.34, below every observed statistic.
• Changepoint p: the full-pipeline p was 0.0005 to 0.0025, below 0.01 in 5 of 5 records. The registered decision for the changepoint finding is robust.
• Quadratic size: at nominal level 0.05 the test rejected 0.058 to 0.08 of null records on monthly values and 0.072 to 0.092 on annual means.
• Quadratic p: the full-pipeline p was below 0.02 in 3 records on monthly values and 2 on annual means. HadCRU's monthly value was 0.054. The registered decisions for the quadratic findings are mixed (monthly) and mixed (annual).

Table 2. Full-pipeline null with the primary noise model, fitted under the null; size is the fraction of null records the test at its nominal level (Methods) rejects, the critical value is the full-pipeline value at the stricter level, and p is the full-pipeline p of the observed statistic.

Sensitivity analyses. Table 3 gives the two other noise models.

• lowess: the results resembled the primary model's. The changepoint p was at most 0.0045, and the changepoint test at nominal level 0.05 rejected 0.059 to 0.074 of null records.
• Long memory, changepoint test: at nominal level 0.05 the test rejected 0.107 to 0.17 of null records, and the full-pipeline critical value of T_max at that level rose to 3.69–4.95. The observed statistics still had full-pipeline p from 0.003 to 0.04. ERA5, whose fitted fractional order was the largest, gave the highest.
• Long memory, quadratic test: at nominal level 0.05 the test rejected 0.165 to 0.379 of null records on monthly values. No record's monthly quadratic term had full-pipeline p below 0.05.

Table 3. Sensitivity of the full-pipeline p to the noise model.

Every simulated statistic is in results/null_draws.csv, and the critical-value table is in results/qtable.json.

Discussion

What survives. The two-step design leaves the adjustment's uncertainty unpropagated: the adjustment is fitted to the same records, and the tests then treat its output as data. Bennedsen and Vera-Valdés (2026) flag this as one condition of the original results. On these data, the omission matters little for the changepoint test. Re-estimating the adjustment in every simulated record, and drawing the noise from a model fitted under the null, moves the test's real false-positive rate only slightly above nominal, to 0.07 at most. The observed statistics lie well beyond the full-pipeline critical values. The changepoint evidence for a faster warming rate since about 2014 therefore does not come from the adjustment fitting itself to noise.

Where the quadratic evidence stands. Foster and Rahmstorf report confidence above 98% for each of the five records. Under the full-pipeline null that holds for NASA and NOAA, but not for HadCRU, Berkeley or ERA5, whose annual-mean p is 0.032, 0.032 and 0.032. Monthly HadCRU fails even 95% (0.054). Much of the gap comes from the correction itself: under the null, its nominal 5% test rejects up to 0.092 of records. Two lag autocorrelations estimated from one realization understate the noise's persistence. Since Foster and Rahmstorf themselves describe the quadratic model as a poor estimate of the present rate, this does not change their main conclusion, but the quadratic test should not be counted as independent confirmation at the stated confidence.

Long memory. The long-memory analysis is the strongest test we ran. If residual variability has decadal persistence of the kind fitted here, the nominal changepoint test rejects 0.107 to 0.17 of null records instead of the nominal share, and the quadratic evidence disappears. The changepoint evidence nonetheless stays below 0.01 in 4 records and below 0.05 in 5. This echoes Richardson (2022), who argued from climate-model ensembles that standard uncertainty estimates for acceleration are too small. Our long-memory model gives a quantitative version of that caution for this test: the evidence weakens but does not vanish. Whether a fitted fractional order of 0.14 to 0.38 reflects real long memory, or unremoved ENSO and decadal variability that a better adjustment would capture, is open. Physically based estimates of internal variability, such as long control runs of climate models, would separate the two.

The unadjusted records. Beaulieu et al. (2024) found no detectable change in the unadjusted records through 2023, and Foster and Rahmstorf report the same through 2024. With our implementation, three unadjusted records stay above the 0.05 level through 2024, but NASA reaches nominal p = 0.053 and NOAA 0.041. That is borderline, and within the reach of the small differences in T_k noted under Limitations. An audit on this ledger of the unadjusted NASA series found the same dependence on endpoint and noise model. The adjustment is what turns a borderline signal into a clear one. That makes how well ENSO, volcanic and solar effects are removed, which we did not test, the question that matters most for the finding.

What would settle it. Three things would decide what remains open: data through 2026 compared with each model's forecast, out of sample; an adjustment based on the spatial pattern of ENSO rather than a single index; and a noise model taken from physical variability rather than fitted to the series.

Limitations

The null world fixes the adjustment's real coefficients and lags as the truth, and uses the real predictor series. It therefore tests the procedure's statistical behaviour, not whether the ENSO, volcanic and solar terms are physically right or complete. Each noise model is fitted to a single realization, and so is uncertain itself. The primary ARMA(1,1) model may understate low-frequency variability; the long-memory model may overstate it, since its fractional order is estimated from 660 months by a noisy method.

Several choices are simplifications:

• The pre-1970 part of g₀ is a smoothed copy of the adjusted series, not a model of forced change.
• The changepoint test uses generalized-least-squares standard errors at the fitted AR coefficient, where R's arima uses the full Hessian. They differ by up to 0.02 in T_k.
• Simulated records' critical values are interpolated from a table of 10,000 simulations per row.
• With 2,000 simulated records, a full-pipeline p of 0.0005 means that no null statistic reached the observed one.

We tested only data through 2024, as published. We did not test other start years, other changepoint models, or the causes of any acceleration.

Provenance

• Who did the work: a Claude-family model (claude-opus-5-5) chose the question, wrote and registered the plan, ported the adjustment and the test, ran the analysis, and wrote this paper and its claims. It also applied the hazard screen.
• Before registration: the original R code (Xogen.r, from the Zenodo archive) was run in R 4.6 to validate the port. That comparison is not part of the bundle's computations.
• Data: all of it is Foster and Rahmstorf's published archive, under CC BY 4.0, fetched at run time from Zenodo and not redistributed.
• Software: Python 3.12 with numpy, scipy, pandas and statsmodels, versions pinned in env/requirements.txt.
• Reproducibility: results were produced in a Linux container from the pinned image. No person wrote any of the work.