Methods

This page states what was done, in the order it was done, with the decisions that could have gone another way. The paper covers the same ground more briefly. Every number quoted on this site comes from results/canonical_results.json, which records where in the result tables each one was read from.

1. Data and units

Every raw file is SHA-256 hashed in data/raw/MANIFEST.json at download. Unit handling lives in one module (units.py) and is enforced, not assumed:

• An absolute temperature must lie between 150 and 350 K, or loading stops. A Celsius array mislabelled as Kelvin fails this immediately.
• An anomaly larger than ±60 K is refused.
• Berkeley absolute temperature is rebuilt as climatology[month] + anomaly and only then converted to Kelvin.
• Irradiance must lie in −1 to 1500 W m⁻²; NCEP's handful of −0.13 W m⁻² values are clipped to zero.
• Palaeo time is ka BP (thousands of years before 1950, positive into the past). La2004's negative kyr are negated on load, and the loader asserts the present-day elements.

Area weights are exact spherical band areas, sin(φ_north) − sin(φ_south), computed from cell edges. The same code serves the unequally spaced Gaussian grid and the regular grid. Check: area-weighted land fraction is 0.284 (NCEP) and 0.288 (Berkeley), against a true value of about 0.29.

2. Orbital geometry and insolation

Daily-mean top-of-atmosphere insolation follows Berger (1978): Q̄ = (S₀/π)(a/r)²[H₀ sin φ sin δ + cos φ cos δ sin H₀], with H₀ = arccos(−tan φ tan δ), including the polar-day and polar-night limits. S₀ = 1360.8 W m⁻² (Kopp & Lean 2011). The Berger convention of 1365 W m⁻² is used only for validation against published values.

The Sun's position is computed by two independent routes:

1. Meeus (1998) low-precision solar theory on real calendar dates, used for the instrumental era.
2. A Keplerian solution driven by La2004 elements, used for the palaeo era and for counterfactual orbits.

The two agree to 0.018° in solar longitude over a full year.

Validation targets:

Angle convention. La2004's ϖ is the heliocentric longitude of perihelion. The Sun's apparent longitude at perihelion is ϖ + 180°. A unit test pins this: perihelion must fall in early January, and an off-by-π error would move it to July.

Calendar. Monthly means are equal-time averages. For climatologies the March equinox is anchored at its 1981–2010 mean date (day 78.7 of a no-leap year). Averaging over equal steps in solar longitude would double-count Kepler's second law.

3. Regime A: the modern seasonal cycle

Target. The 1981–2010 monthly climatology of every grid cell, every zonal row (land and ocean separately) and 34 regions (8 latitude zones × {all, land, ocean} plus 10 named boxes). A regional month is kept only if at least 80% of the region's area has data.

Seasonal share. The fraction of monthly variance carried by the climatology is computed within the baseline years, where the month-dummy projection makes the seasonal and anomaly parts exactly orthogonal. It is the ceiling on what any purely seasonal model can explain. The share explained by the physical model is then exactly seasonal share × R²(model).

Physical model. A one-pole relaxation: a zero-dimensional energy balance with heat capacity C and feedback λ.

C dT/dt = k′Q(t) − λ(T − T₀), so T̂(ω) = k Q̂(ω)/(1 + iωτ), with τ = C/λ.

The daily insolation cycle is filtered in Fourier space and then averaged into calendar months, so monthly averaging is treated exactly. For fixed τ the model is linear in (T₀, k), which makes the fit an exact profile over 174 τ values between 0 and 365 days. Admissibility: the primary fit requires k ≥ 0 and τ ≤ 365 d. The unconstrained fit is reported as a diagnostic only. It wins mainly with k < 0 ("more sun, colder") or with decade-long τ, and those cases are excluded from every attribution statement.

Comparators with the same data: a one-harmonic (3-parameter) and two-harmonic (5-parameter) regression, and zero-lag insolation (2 parameters).

Eccentricity detectability (perfect-model test). A free gain and lag absorb most of the eccentricity term, which is close to an annual harmonic peaking at perihelion. So a p-value on ΔR² answers the wrong question. Instead:

1. Fit a circular-orbit (e = 0) model to the noise-free response of the real-forcing model.
2. The misfit of that fit, 1 − R², is the most R² eccentricity could ever add.
3. Divide it by the real-forcing model's misfit to the data.

Below 1, the eccentricity contribution is smaller than what the model already gets wrong. Uncertainty comes from resampling the 30 baseline years (B = 500).

Global mean. With e = 0 the global-mean insolation is constant through the year (the numerical standard deviation is 0.001 W m⁻²), so the global-mean insolation cycle is purely eccentricity-driven. It is compared with global-mean temperature. The annual harmonic is decomposed into area-weighted hemisphere × surface contributions.

Phase trends. The day of maximum of the annual harmonic is computed per year. Its trend is estimated by OLS with Newey–West HAC errors and measured relative to the calendar-anchored insolation phase. That comparison removes Gregorian leap-year drift of the equinox, about −0.074 d per decade. The precession-only prediction comes from La2004 elements 1 kyr apart. Windows: the full record, 1954–2007 and 1979 onward.

4. Regime B: orbital time scales

Gridding. Records are bin-averaged, not point-interpolated, onto a regular 2 ka grid (2.5 ka before the MPT). Averaging keeps high-frequency variance from aliasing into orbital bands. Empty bins fall back to interpolation and their fraction is reported; it is zero at 2 ka for every record.

Windows.

• Common 0–640 ka for all three records, so the circularity comparison uses the same length and the same estimator.
• 0–800 ka for LR04 and EPICA.
• 1250–2000 ka for LR04 before the mid-Pleistocene transition.

Spectra. Thomson multitaper. The multitaper half-bandwidth W = NW/(N·Δt) decides whether the 23- and 19-kyr lines can be separated. At NW = 4 over 800 kyr they cannot; the analysis uses NW = 3 (2.5 on the 640-kyr window) and states the resolution flag with every result. Band enrichment is the observed share of power in a band divided by the median share for 1000 AR(1) surrogates fitted to the record.

Coherence. Multitaper magnitude-squared coherence and phase. The null phase-randomises the record, not the forcing. The analytic K-taper level assumes the tested series is locally white. A record with its own spectral line at the forcing period breaks that assumption: two lines at one frequency are coherent whatever their phase. A unit test demonstrates that this makes the analytic level reject almost always, while record surrogates hold their size. Claims are restricted to physically meaningful pairs (obliquity at 41 kyr, precession at 23 and 19 kyr, eccentricity at 100 kyr, and the relevant insolation series). They are corrected by Benjamini–Hochberg within each record and window.

Tuning-artefact null. LR04's chronology is tuned to 65°N insolation (its own header says so), so its orbital coherence is partly circular by construction. To measure how much:

1. Take 200 phase-randomised surrogates of LR04 (its spectrum, no orbital phase).
2. Align each to a 65°N-insolation target (a one-pole response with τ = 15 ka) by slope-constrained dynamic time warping. Steps (1,1), (1,2) and (2,1) limit the implied change in sedimentation rate to between ×0.5 and ×2. A Sakoe–Chiba band caps the age adjustment at 0, 6, 10 or 20 ka.
3. The untuned Cheng record is re-tuned the same way.

DTW is more aggressive than LR04's tie-point tuning, so the results are upper bounds for a given age freedom.

Linear response. A one-pole relaxation of an insolation series, filtered in physical time on the long La2004 record so it is spun up inside the window. τ is chosen from 0–100 ka by blocked 5-fold cross-validation with a 10 ka embargo. The surrogate R² null re-fits every τ for each of 500 phase-randomised copies of the record. The statistical comparator is a multiple regression on obliquity, e·sin ϖ̃, e·cos ϖ̃ and e with a common τ. Standard errors are HAC, with the bandwidth set by the residual autocorrelation.

Wavelets. Morlet (ω₀ = 6) continuous wavelet transform after Torrence and Compo (1998), with the cone of influence and a χ² test against AR(1). It is used only where non-stationarity is the question: the mid-Pleistocene transition.

5. Model comparison

Three families, run on the same targets, folds and metrics:

• Physical: zero-lag insolation, and the one-pole model with and without a trend.
• Statistical: monthly climatology, harmonic regression with 1–2 harmonics, and Fourier regression with 1–6 harmonics chosen by inner CV. In Regime B: sinusoids at the La2004 periods and regression on orbital elements.
• Flexible: a Gaussian process on time (periodic × envelope + trend + noise; hyperparameters by L-BFGS-B capped at 60 iterations), and gradient boosting on orbital inputs alone or with time and month.

Evaluation schemes:

• Regime A: blocked 5-fold CV over contiguous years with a 12-month embargo; an extrapolation split (train ≤ 1985, test ≥ 1990); and transfer. Transfer predicts the seasonal cycle of latitude bands, or whole domains (another hemisphere, the ocean), that the model never saw. This is the one test here that probes mechanism rather than association.
• Regime B: blocked 5-fold CV with a 10 ka embargo, and extrapolation between the younger and older halves.

Skill is reported as R², RMSE, and MSE skill against a reference refitted on the same training fold (climatology in Regime A, the training mean in Regime B).

Neural forecasting is omitted on purpose. Three-parameter models already reach the seasonal ceiling, and orbital inputs carry no year-to-year information for a network to learn.

6. Sensitivity

2939 rows, one statistic per row, across these axes:

• temporal aggregation (monthly, bimonthly, seasonal);
• baseline period (1951–1980, 1981–2010, 1991–2020, full record);
• variable and dataset (NCEP T, NCEP surface SW, Berkeley T);
• orbital-parameter subset (full, obliquity-only, eccentricity-only);
• detrending;
• resampling step (1, 2, 4 ka) and multitaper NW (2–4);
• predictor subset (nine combinations);
• smoothing window (0–60 ka), with a naive and a smoothing-aware null;
• random age-model error (standard deviation 2–6 ka, correlation 20 ka).

A synthetic Slutsky–Yule experiment runs the smoothing tests on pure AR(1) noise.

7. Inference that survives dependence

• Serial correlation: correlations use the Bretherton et al. (1999) effective sample size. On independent AR(1) pairs with φ = 0.9 this cuts the false-positive rate from 49.7% (naive) to 4.3%.
• Regressions: Newey–West HAC standard errors.
• Resampling: block bootstrap.
• Spatial dependence: fields report their effective spatial degrees of freedom and use BH-FDR, never the raw number of cells.

8. Reproducibility

python scripts/00_fetch_data.py         # data, verified by checksum
python scripts/01_pilot_single_region.py
python scripts/02_seasonal_global.py ncep && python scripts/02_seasonal_global.py berkeley
python scripts/03_paleo.py
python scripts/04_models.py             # checkpointed; resumes after interruption
python scripts/05_sensitivity.py
python scripts/06_figures.py
python scripts/07_export_site.py
python scripts/verify_claims.py         # every tagged number in the paper against the registry
python -m pytest -q

All random draws use fixed seeds. Python 3.11; package versions are pinned in requirements.txt.