The prediction error floor for millisecond pulsar timing residuals is bounded below by stochastic processes (spin noise, gravitational wave background, interstellar medium turbulence) that preclude deterministic prediction to better than 50 nanoseconds at horizons exceeding 1 hour, establishing a hard physical limit on pre-epoch key derivation in pulsar-anchored time-lock encryption.

For millisecond pulsars (MSPs) with spin periods P < 10 ms and timing baselines > 1 year, the root-mean-square (RMS) timing residual prediction error σ_pred satisfies σ_pred ≥ 50 ns for all prediction horizons Δt > 3600 s (1 hour), due to irreducible stochastic contributions from: (1) spin noise with power spectral density S_spin(f) ∝ f^{-α}, α ∈ [2,6]; (2) gravitational wave background (GWB) strain h_c ~ 10^{-15} at f = 1/yr; and (3) interstellar medium (ISM) dispersion measure (DM) variations with ΔDMRMS ~ 10^{-4} to 10^{-3} pc cm^{-3} yr^{-1}. This 50 ns floor constitutes a hard physical limit preventing deterministic pre-epoch key derivation in any pulsar-anchored time-lock cryptographic scheme with security parameter λ ≥ 128 bits.

PHASE 1 — Noise Floor Characterization (Days 1–45): Analyze existing PTA datasets (NANOGrav 15yr, PPTA DR3, EPTA DR2) to empirically measure prediction residuals as a function of horizon Δt for all MSPs with ≥ 5-year baselines. Fit Bayesian noise models (white + red + DM) and compute posterior predictive distributions at horizons Δt = {1h, 6h, 24h, 1wk, 1mo}.

PHASE 2 — Stochastic Component Isolation (Days 46–90): Decompose residuals into spin noise, GWB, and ISM contributions using: (a) chromatic index analysis across frequency bands; (b) Hellings-Downs spatial correlation for GWB; (c) structure function analysis for DM variations. Quantify each component's contribution to σ_pred(Δt).

PHASE 3 — Predictive Model Stress Test (Days 91–135): Train deterministic + ML hybrid predictors on 80% of each MSP's timeline; evaluate on 20% held-out data. Test predictors: (a) polynomial ephemeris extrapolation; (b) Gaussian process regression; (c) LSTM/Transformer sequence models. Record σ_pred for each Δt bin.

PHASE 4 — Cryptographic Bound Formalization (Days 136–160): Translate σ_pred distributions into information-theoretic bounds on key entropy. Compute H(K | TOA_predicted) for key lengths 128, 256 bits as a function of σ_pred. Determine minimum σ_pred required to break 128-bit security.

# Pulsar Timing Prediction Floor Validation — Core Architecture
# Dependencies: enterprise, enterprise_extensions, TEMPO2, PTMCMCSampler,
#               numpy, scipy, torch, gpytorch, optuna, astropy

import numpy as np
from enterprise.pulsar import Pulsar
from enterprise.signals import signal_base, white_signals, gp_signals
from enterprise_extensions import models, sampler
import torch
import torch.nn as nn

# ============================================================
# MODULE 1: DATA INGESTION & MSP SELECTION
# ============================================================
class MSPDataLoader:
    def __init__(self, par_files, tim_files, min_baseline_yr=5.0,
                 max_DM=100.0, max_period_ms=10.0):
        self.selection_criteria = {
            'baseline': min_baseline_yr,  # years
            'DM_max': max_DM,             # pc/cm^3
            'P_max': max_period_ms * 1e-3 # seconds
        }

    def load_and_filter(self):
        pulsars = []
        for par, tim in zip(self.par_files, self.tim_files):
            psr = Pulsar(par, tim, ephem='DE440', clk='TT(BIPM2021)')
            baseline = (psr.toas.max() - psr.toas.min()) / 365.25  # days->yr
            if (psr.period < self.selection_criteria['P_max'] and
                psr.DM < self.selection_criteria['DM_max'] and
                baseline >= self.selection_criteria['baseline']):
                pulsars.append(psr)
        return pulsars  # Expected: ~45 MSPs

# ============================================================
# MODULE 2: BAYESIAN NOISE MODEL FITTING
# ============================================================
class NoiseModelFitter:
    def __init__(self, pulsar, n_samples=1_000_000):
        self.psr = pulsar
        self.n_samples = n_samples

    def build_model(self):
        # White noise: EFAC + EQUAD per backend
        ef = white_signals.MeasurementNoise(efac=True)
        eq = white_signals.EquadNoise(log10_equad=True)
        # Red spin noise: power law
        rn = gp_signals.RedNoise(
            spectrum=models.powerlaw,
            components=30,
            prior_log10_A=(-16, -12),
            prior_gamma=(0, 7)
        )
        # DM noise: chromatic index kappa=2
        dm = gp_signals.DMNoise(
            spectrum=models.powerlaw,
            components=30
        )
        return ef + eq + rn + dm

    def run_mcmc(self):
        pta = signal_base.PTA([self.build_model()(self.psr)])
        samp = sampler.setup_sampler(pta, outdir=f'chains/{self.psr.name}')
        samp.sample(self.n_samples)
        return samp  # Returns posterior samples

# ============================================================
# MODULE 3: PREDICTION HORIZON ANALYSIS
# ============================================================
class PredictionHorizonAnalyzer:
    HORIZONS_SEC = [3600, 21600, 43200, 86400, 172800, 604800, 2592000]
    # 1h, 6h, 12h, 24h, 48h, 1wk, 1mo

    def leave_future_out_cv(self, toas, residuals, noise_params,
                             n_windows=50):
        results = {h: [] for h in self.HORIZONS_SEC}
        toa_range = toas.max() - toas.min()

        for _ in range(n_windows):
            # Random split point in middle 60% of data
            split_frac = np.random.uniform(0.3, 0.7)
            T_train = toas.min() + split_frac * toa_range

            train_mask = toas <= T_train
            for horizon in self.HORIZONS_SEC:
                # Select TOAs within [T_train, T_train + horizon]
                pred_mask = ((toas > T_train) &
                             (toas <= T_train + horizon / 86400.0))
                if pred_mask.sum() < 3:
                    continue

                # Predict using GP posterior mean
                pred_residuals = self._gp_predict(
                    toas[train_mask], residuals[train_mask],
                    toas[pred_mask], noise_params
                )
                errors = residuals[pred_mask] - pred_residuals
                results[horizon].append(np.sqrt(np.mean(errors**2)) * 1e9)  # ns

        return {h: np.array(v) for h, v in results.items()}

    def _gp_predict(self, t_train, r_train, t_pred, params):
        # Gaussian Process prediction using fitted noise model
        # Returns posterior mean at t_pred given training data
        from scipy.linalg import solve
        K_train = self._covariance_matrix(t_train, t_train, params)
        K_cross = self._covariance_matrix(t_pred, t_train, params)
        K_train_reg = K_train + params['white_var'] * np.eye(len(t_train))
        return K_cross @ solve(K_train_reg, r_train)

    def _covariance_matrix(self, t1, t2, params):
        # Power-law red noise covariance (Fourier sum approximation)
        tau = np.abs(t1[:, None] - t2[None, :])  # days
        A, gamma = params['red_A'], params['red_gamma']
        # Approximate: C(tau) = A^2 * (tau/yr)^{-(gamma-1)} / (gamma-1)
        yr = 365.25
        return A**2 / (gamma - 1) * (tau / yr + 1e-10)**(-(gamma - 1))

# ============================================================
# MODULE 4: TRANSFORMER PREDICTOR
# ============================================================
class PulsarTransformerPredictor(nn.Module):
    def __init__(self, d_model=256, nhead=8, num_layers=6,
                 seq_len=180, pred_len=1):
        super().__init__()
        self.embedding = nn.Linear(2, d_model)  # [TOA, residual]
        encoder_layer = nn.TransformerEncoderLayer(
            d_model=d_model, nhead=nhead, dim_feedforward=1024,
            dropout=0.1, batch_first=True
        )
        self.transformer = nn.TransformerEncoder(encoder_layer, num_layers)
        self.output_head = nn.Linear(d_model, pred_len)
        self.seq_len = seq_len

    def forward(self, x):
        # x: [batch, seq_len, 2] — normalized TOA gaps + residuals
        x_emb = self.embedding(x)
        encoded = self.transformer(x_emb)
        return self.output_head(encoded[:, -1, :])  # Predict next residual

    def predict_at_horizon(self, history_toas, history_res, target_toa):
        # Interpolate/extrapolate to target_toa
        # Returns predicted residual in nanoseconds
        self.eval()
        with torch.no_grad():
            x = self._prepare_input(history_toas, history_res, target_toa)
            return self.forward(x).item() * 1e9  # Convert to ns

# ============================================================
# MODULE 5: STOCHASTIC DECOMPOSITION
# ============================================================
class StochasticDecomposer:
    def structure_function(self, times, residuals, tau_bins):
        """Compute D(tau) = <[R(t+tau) - R(t)]^2>"""
        D = np.zeros(len(tau_bins))
        for i, tau in enumerate(tau_bins):
            pairs = []
            for j, t in enumerate(times):
                future = np.where(np.abs(times - (t + tau)) < tau * 0.1)[0]
                if len(future) > 0:
                    pairs.append((residuals[future[0]] - residuals[j])**2)
            D[i] = np.mean(pairs) if pairs else np.nan
        return D

    def fit_power_law(self, tau_bins, D_tau):
        """Fit D(tau) = A * tau^beta; return (A, beta, R^2)"""
        from scipy.optimize import curve_fit
        valid = ~np.isnan(D_tau) & (D_tau > 0)
        log_tau = np.log10(tau_bins[valid])
        log_D = np.log10(D_tau[valid])
        popt, _ = curve_fit(lambda x, a, b: a + b*x, log_tau, log_D)
        beta = popt[1]
        residuals = log_D - (popt[0] + beta * log_tau)
        ss_res = np.sum(residuals**2)
        ss_tot = np.sum((log_D - log_D.mean())**2)
        return 10**popt[0], beta, 1 - ss_res/ss_tot

# ============================================================
# MODULE 6: CRYPTOGRAPHIC BOUND CALCULATOR
# ============================================================
class CryptographicBoundCalculator:
    def __init__(self, key_bits=128, quantization_ns=10.0):
        self.key_bits = key_bits
        self.delta = quantization_ns  # ns

    def conditional_entropy(self, sigma_pred_ns):
        """
        H(K | TOA_predicted) for Gaussian prediction error.
        Key derived as K = hash(round(TOA / delta) * delta)
        Attacker knows predicted TOA with uncertainty sigma_pred.
        """
        import scipy.stats as stats
        # Number of TOA bins within ±3*sigma of prediction
        n_bins = int(6 * sigma_pred_ns / self.delta) + 1
        # Probability of each bin given Gaussian uncertainty
        bin_centers = np.arange(-n_bins//2, n_bins//2 + 1) * self.delta
        probs = stats.norm.pdf(bin_centers, 0, sigma_pred_ns)
        probs /= probs.sum()
        # Shannon entropy of bin distribution (attacker's uncertainty)
        H_bins = -np.sum(probs * np.log2(probs + 1e-300))
        # Total key entropy
        H_key = self.key_bits
        # Conditional entropy (remaining uncertainty after prediction)
        H_conditional = min(H_key, H_bins)
        return H_conditional

    def security_bits(self, sigma_pred_ns):