Quantifying Model Uncertainties

---

Lecture 31

April 15, 2026

Overview

• Model Background
• Correlated Discrepancies
• Bootstrapping Residuals
• Upcoming Schedule
• References

Model Background

Planetary Energy Balance

Representation of Planetary Energy Balance

Source: Reprinted from A Climate Modeling Primer, A. Henderson-Sellers and K. McGuffie, Wiley, pg. 58, (1987) via https://www.e-education.psu.edu/meteo469/node/137.

An Energy Balance Model (EBM)

\[\begin{align*} \overbrace{\frac{dH}{dt}}^{\text{change in heat}} &= \overbrace{F}^{\text{RF}} - \overbrace{\lambda T}^{\substack{\text{change in} \\ \text{temperature}}} \\[1em] \Rightarrow C\frac{dT}{dt} &= F - \lambda T - \gamma(T-T_D)\\ C_D\frac{dT_D}{dt} &= \gamma(T-T_D) \end{align*}\]

Two Layer EBM Schematic

Source: Palmer et al. (2018)

EBM Discretization

Use Euler discretization:

\[\begin{align*} T(t+1) &= T(t) + \frac{F(t) - \lambda T(t) - \gamma(T(t) - T_D(d))}{C} \Delta t \\[0.5em] T_D(t+1) &= T_D(t) + \frac{\gamma (T(t) - T_D(t))}{C_D} \Delta t \end{align*}\]

Equilibrium Climate Sensitivity (ECS)

Under steady-state conditions (constant \(F\) and \(dT/dt = 0\)), \[T = \frac{F}{\lambda}.\]

When we double atmospheric CO₂, we refer to the equilibrium temperature \(S\) as the equilibrium climate sensitivity:

\[S = \underbrace{F_{2\times \text{CO}_2}}_{\approx 4 \text{W/m}^2}/\lambda\]

Probability Models for Simulation Models

Most common setting (e.g. Brynjarsdóttir & O’Hagan (2014)):

\[\mathbf{y} = \underbrace{\eta(\mathbf{x}; \theta)}_{\text{model}} + \underbrace{\delta(x)}_{\text{discrepancy}} + \underbrace{\varepsilon}_{\text{error}}\]

Model-Data Discrepancy

Systematic disagreements between model predictions and system state (vs. random observation errors).

For example, \(\delta\) might capture:

• Bias (e.g.: model consistently over/underpredicts);
• Accumulations of error (e.g.: persistent model underestimates);
• Partial observations (e.g.: do you count every animal?)

Correlated Discrepancies

MLE for Gaussian Discrepancy

function ebm(rf_nonaerosol, rf_aerosol; p=(3.2, 1.0, 1.3, 100.0, 800.0, -0.1))
    # set up model parameters
    S, γ, α, d, D, T₀ = p # this unpacks the parameter tuple into variables
    F2xCO₂ = 4.0 # radiative forcing [W/m²] for a doubling of CO₂
    λ = F2xCO₂ / S

    c = 4.184e6 # heat capacity/area [J/K/m²]
    C = c*d # heat capacity of mixed layer (per area)
    CD = c*D # heat capacity of deep layer (per area)
    F = rf_nonaerosol + α*rf_aerosol # radiative forcing
    Δt = 31558152. # annual timestep [s]

    T = zero(F)
    T[1] = T₀
    TD = zero(F)
    for i in 1:length(F)-1
        T[i+1] = T[i] + (F[i] - λ*T[i] - γ*(T[i]-TD[i]))/C * Δt
        TD[i+1] = TD[i] + γ*(T[i]-TD[i])/CD * Δt
    end
    # return after normalizing to reference period
    return T
end

ebm_wrap(params) = ebm(forcing_non_aerosol_85[hind_idx], forcing_aerosol_85[hind_idx], p = params)

function gaussian_iid_homosked(params, temp_dat, temp_err, m)
    S, γ, α, d, D, T₀, σ = params 
    ebm_sim = m((S, γ, α, d, D, T₀))
    ll = sum(logpdf.(Normal.(ebm_sim, sqrt.(σ^2 .+ temp_err.^2)), temp_dat))
    return ll
end

lower = [1.0, 0.5, 0.0, 50.0, 200.0, temp_lo[1], 0.0]
upper = [5.0, 1.5, 2.0, 200.0, 1000.0, temp_hi[1], 10.0]
p0 = [3.0, 1.0, 1.0, 100.0, 800.0,temp_obs[1], 5.0]

result = Optim.optimize(params -> -gaussian_iid_homosked(params, temp_obs, temp_sd, ebm_wrap), lower, upper, p0)
θ_iid = result.minimizer
θ_iid_rd = round.(θ_iid; digits=1)
parnames = ["S", "γ", "α", "d", "D", "T₀", "σ"]
param_df = DataFrame(Parameters=parnames, MLE=θ_iid_rd)
pretty_table(param_df[1:4, :]; backend=:markdown)

Parameters String	MLE Float64
S	3.4
γ	1.5
α	1.0
d	77.8

Parameters String	MLE Float64
D	623.7
T₀	-0.1
σ	0.1

Gaussian Discrepancy Hindcast

n_samples = 10_000
temp_iid = ebm_wrap(θ_iid) # simulate IID best fit
# simulate projections with discrepancy and errors
temp_iid_proj = zeros(n_samples, length(temp_sd))
for i = 1:length(temp_sd)
    temp_iid_err = rand(Normal(0, sqrt.(θ_iid[end]^2 .+ temp_sd[i]^2)), n_samples)
    temp_iid_proj[:, i] = temp_iid[i] .+ temp_iid_err
end
# calculate quantiles
temp_iid_q = mapslices(col -> quantile(col, [0.05, 0.5, 0.95]), temp_iid_proj; dims=1)

p = scatter(time_obs, temp_obs, color=:black, label="Observations", ylabel="(°C)", xlabel="Year", title="Temperature Anomaly")
plot!(p, time_obs, temp_iid_q[2, :], lw=3, ribbon=(temp_iid_q[2, :] - temp_iid_q[1, :], temp_iid_q[3, :] - temp_iid_q[2, :]), color=:red, fillalpha=0.3, label="Gaussian Discrepancy")
plot!(size=(1200, 500))

Figure 1: MLE Fit for Gaussian discrepancy

Analyzing Residual Assumptions

resids_homogauss = temp_obs - temp_iid

p1 = qqnorm(resids_homogauss, xlabel="Theoretical Values", ylabel="Empirical Values", title="Normal Q-Q Plot", size=(600, 450))
pacf_homogauss = pacf(resids_homogauss, 1:5)
p2 = plot(1:5, pacf_homogauss, marker=:circle, line=:stem, linewidth=3, markersize=8, tickfontsize=16, guidefontsize=18, legend=false, ylabel="Partial Autocorrelation", xlabel="Time Lag", title="Partial Autocorrelation Plot", size=(600, 450))
hline!(p2, [0], color=:black, linestyle=:dash)

display(p1)
display(p2)

(a) Residual diagnostics

(b)

Figure 2

AR(1) Residuals

Autocorrelated

\[ \begin{align*} y_t &= \text{EBM}(\theta; F_t) + \delta_t + \varepsilon_t \\ \delta_t &= \rho \delta_{t-1} + \omega_t \\ \varepsilon_t &\sim N(0, \sigma^2_{\text{obs}, t}) \\ \omega_t &\sim N(0, \sigma^2) \end{align*} \]

Independent

\[ \begin{align*} y_t &= \text{EBM}(\theta; F_t) + \delta_t + \varepsilon_t \\ \delta_t &\sim N(0, \sigma^2) \\ \varepsilon_t &\sim N(0, \sigma^2_{\text{obs}, t}) \end{align*} \]

Likelihood Function

Without observation errors (\(\varepsilon\)), can whiten residuals:

\[ \begin{align*} y_1 & \sim N\left(0, \frac{\sigma^2}{1 - \rho^2}\right) \\ y_t - \rho y_{t-1} &\sim N(0, \sigma^2) \end{align*} \]

When observation error variances are known (as in here, provided in the dataset), this also works.

Non-Identifiability of Parameters

But when observation errors are also uncertain, run into non-identifiability:

\[N(0, \sigma^2) + N(0, \sigma_\text{obs}^2) = N(0, {\color{red}\sigma^2 + \sigma_\text{obs}^2})\]

Joint Likelihood for AR(1)

Could also use joint likelihood for residuals \(y_t - \text{EBM}(\theta; F_t)\):

\[ \begin{align*} \mathbf{y} &\sim \mathcal{N}(\mathbf{0}, \Sigma) \\ \Sigma &= \frac{\sigma^2}{1 - \rho^2} \begin{pmatrix}1 + \sigma_{\text{obs}, 1}^2 & \rho & \ldots & \rho^{T-1} \\ \rho & 1 + \sigma_{\text{obs}, 2}^2 & \ldots & \rho^{T-2} \\ \vdots & \vdots & \ddots & \vdots \\ \rho^{T-1} & \rho^{T-2} & \ldots & 1+ \sigma_{\text{obs}, T}^2\end{pmatrix} \end{align*} \]

Whitened Likelihood for AR(1) (in code)

# temp_obs: temperature data
# temp_err: standard deviations for temperatures
# m: model function
function ar_loglik(params, temp_obs, temp_err, m)
    S, γ, α, d, D, T₀, ρ, σ = params 
    ebm_sim = m((S, γ, α, d, D, T₀))
    residuals = temp_obs - ebm_sim
    # whiten residuals
    ll = 0
    # notice addition of observation errors
    for t = 1:length(temp_sd)
        if t == 1
            ll += logpdf(Normal(0, sqrt(σ^2 / (1 - ρ^2) + temp_err[1]^2)), residuals[1])
        else
            resid_wn = residuals[t] - ρ * residuals[t-1]
            ll += logpdf(Normal(0, sqrt(σ^2 + temp_err[t]^2)), resid_wn)
        end
    end
    return ll
end

AR(1) MLE (Code)

lower = [1.0, 0.5, 0.0, 50.0, 200.0, temp_lo[1], -1.0, 0.0]
upper = [5.0, 1.5, 2.0, 200.0, 1000.0, temp_hi[1], 1.0, 10.0]
p0 = [3.0, 1.0, 1.0, 100.0, 800.0,temp_obs[1], 0.0, 5.0]

result = Optim.optimize(params -> -ar_loglik(params, temp_obs, temp_sd, ebm_wrap), lower, upper, p0)
θ_ar = result.minimizer

Comparison of MLE

Parameters String	IID Any	AR Float64
S	3.4	3.3
γ	1.5	1.5
α	1.0	0.9
d	77.8	88.6

Parameters String	IID Any	AR Float64
D	623.7	649.1
T₀	-0.1	-0.1
ρ	-	0.5
σ	0.1	0.1

AR(1) Hindcast

n = 10_000

# get model hindcasts
temp_iid = ebm_wrap(θ_iid[1:end-1])
temp_ar = ebm_wrap(θ_ar[1:end-2])

# get iid and AR residuals from relevant processes
residuals_iid = stack(rand.(Normal.(0, sqrt.(temp_sd.^2 .+ θ_iid[end].^2)), n), dims=1)
residuals_ar = zeros(length(hind_idx), n)
for t = 1:length(temp_sd)
    if t == 1
        residuals_ar[t, :] = rand(Normal(0, sqrt(θ_ar[end]^2 / (1 - θ_ar[end-1]^2) + temp_sd[1]^2)), n)
    else
        residuals_ar[t, :] = θ_ar[end-1] * residuals_ar[t-1, :] + rand(Normal(0, sqrt(θ_ar[end]^2 + temp_sd[t]^2)), n)
    end
end

# add residuals back to model simulations
model_sim_iid = (residuals_iid .+ temp_iid)'
model_sim_ar = (residuals_ar .+ temp_ar)'

# get quantiles
q90_iid = mapslices(col -> quantile(col, [0.05, 0.5, 0.95]), model_sim_iid; dims=1) # compute 90% prediction interval
q90_ar = mapslices(col -> quantile(col, [0.05, 0.5, 0.95]), model_sim_ar; dims=1) # compute 90% prediction interval

p = scatter(time_obs, temp_obs, yerr=(temp_obs - temp_lo, temp_hi - temp_obs), color=:black, label="Observations", ylabel="(°C)", xlabel="Year", title="Temperature Anomaly", markersize=5)
plot!(p, hind_years, q90_iid[2, :], ribbon=(q90_iid[2, :] - q90_iid[1, :], q90_iid[3, :] - q90_iid[2, :]), fillalpha=0.2, linewidth=3, label="IID")
plot!(p, hind_years, q90_ar[2, :], ribbon=(q90_ar[2, :] - q90_ar[1, :], q90_ar[3, :] - q90_ar[2, :]), fillalpha=0.2, linewidth=3, label="AR")
plot!(size=(1200, 500))

Figure 3: Hindcast of the EBM with IID and AR(1) residuals.

Coverage Rates

Coverage Rates of 90% CIs:

• IID: 90.2%
• AR(1): 93.1%

Have We Captured Residual Autocorrelation?

Use simulations to look at distributions of residual probability assumptions.

resids_ar_sim = model_sim_ar .- temp_obs'
boxplot(mapslices(col -> pacf(col, 1:5), resids_ar_sim; dims=1)', label=:false, xlabel="Lag", ylabel="Partial Autocorrelation")
plot!(size=(500, 500))

Figure 4: Distribution of residual partial autocorrelations

Hindcasts May Not Differ Much…

Figure 5: Hindcast of the EBM with IID and AR(1) residuals.

…But Projections Can

ebm_sim_85(params) = ebm(forcing_non_aerosol_85[sim_idx], forcing_aerosol_85[sim_idx], p = params)
ebm_sim_26(params) = ebm(forcing_non_aerosol_26[sim_idx], forcing_aerosol_26[sim_idx], p = params)

# iid residuals
y_err = zeros(length(sim_idx))
y_err[1:length(hind_idx)] = temp_sd

residuals_iid = stack(rand.(Normal.(0, sqrt.(y_err.^2 .+ θ_iid[end]^2)), n), dims=1)
model_iid_85 = ebm_sim_85(θ_iid[1:end-1])
model_iid_26 = ebm_sim_26(θ_iid[1:end-1])
model_sim_iid_85 = (residuals_iid .+ model_iid_85)'
model_sim_iid_26 = (residuals_iid .+ model_iid_26)'
q90_iid_85 = mapslices(col -> quantile(col, [0.05, 0.5, 0.95]), model_sim_iid_85; dims=1) # compute 90% prediction interval```
q90_iid_26 = mapslices(col -> quantile(col, [0.05, 0.5, 0.95]), model_sim_iid_26; dims=1) # compute 90% prediction interval```

# AR residuals
residuals_ar = zeros(length(sim_idx), n)
for t = 1:length(sim_idx)
    if t == 1
        residuals_ar[t, :] = rand(Normal(0, sqrt(θ_ar[end]^2 / (1 - θ_ar[end-1]^2 + temp_sd[1]^2))), n)
    elseif t <= length(hind_idx)
        residuals_ar[t, :] = θ_ar[end-1] * residuals_ar[t-1, :] + rand(Normal(0, sqrt(θ_ar[end]^2 + temp_sd[t]^2)), n)
    else
        residuals_ar[t, :] = θ_ar[end-1] * residuals_ar[t-1, :] + rand(Normal(0, θ_ar[end]), n)
    end
end
model_ar_85 = ebm_sim_85(θ_ar[1:end-2])
model_ar_26 = ebm_sim_26(θ_ar[1:end-2])
model_sim_ar_26 = (residuals_ar .+ model_ar_26)'
model_sim_ar_85 = (residuals_ar .+ model_ar_85)'
q90_ar_26 = mapslices(col -> quantile(col, [0.05, 0.5, 0.95]), model_sim_ar_26; dims=1) # compute 90% prediction interval
q90_ar_85 = mapslices(col -> quantile(col, [0.05, 0.5, 0.95]), model_sim_ar_85; dims=1) # compute 90% prediction interval

p_sim= scatter(time_obs, temp_obs, color=:black, label="Data", ylabel="Temperature Anomaly (°C)", xlabel="Year", right_margin=5mm)
plot!(p_sim, sim_years, q90_iid_26[2, :], ribbon=(q90_iid_26[2, :] - q90_iid_26[1, :], q90_iid_26[3, :] - q90_iid_26[2, :]), fillalpha=0.2, linewidth=3, color=:royalblue, label="IID/SSP1-2.6")
plot!(p_sim, sim_years, q90_ar_26[2, :], ribbon=(q90_ar_26[2, :] - q90_ar_26[1, :], q90_ar_26[3, :] - q90_ar_26[2, :]), fillalpha=0.2, linewidth=3, color=:firebrick1, label="AR/SSP1-2.6")
plot!(p_sim, sim_years, q90_iid_85[2, :], ribbon=(q90_iid_85[2, :] - q90_iid_85[1, :], q90_iid_85[3, :] - q90_iid_85[2, :]), fillalpha=0.2, linewidth=3, color=:blue3, label="IID/SSP5-8.5")
plot!(p_sim, sim_years, q90_ar_85[2, :], ribbon=(q90_ar_85[2, :] - q90_ar_85[1, :], q90_ar_85[3, :] - q90_ar_85[2, :]), fillalpha=0.2, linewidth=3, color=:firebrick, label="AR/SSP5-8.5")
plot!(p_sim, size=(1200, 500))
xlims!((2000, 2100))
ylims!((0.75, 5.5))

Figure 6: Differences in the 90% confidence intervals

These Differences Could Be Decision-Relevant

p1 = histogram(model_sim_iid_26[:, end], color=:blue, xlabel="°C", ylabel="Count", label="IID", size=(600, 450), alpha=0.4, title="SSP1-2.6")
histogram!(p1, model_sim_ar_26[:, end], color=:red, label="AR", alpha=0.4)
p2 = histogram(model_sim_iid_85[:, end], color=:blue, xlabel="°C", ylabel="Count", label="IID", size=(600, 450), alpha=0.4, title="SSP5-8.5")
histogram!(p2, model_sim_ar_85[:, end], color=:red, label="AR", alpha=0.4)

display(p1)
display(p2)

(a) Projections of global mean temperature in 2100

(b)

Figure 7

Bootstrapping Residuals

Residual Resampling

• Hybrid between non-parametric and parametric bootstrapping.
• Use trend model but bootstrap only residuals, then refit.

resids_obs = temp_obs .- temp_ar
p = scatter(time_obs, resids_obs, color=:black, label="Residuals", ylabel="Temperature Anomaly (°C)", xlabel="Year", right_margin=5mm)
plot!(p, size=(600, 450))

Figure 8: Residuals of EBM

Residual Resampling

nboot = 1_000 # number of bootstrap replicates
resids_boot = sample(resids_obs, (length(resids_obs), nboot), replace=true) # bootstrap samples
ebm_boot = temp_ar .+ resids_boot
plot(time_obs, temp_ar, lw=3, color=:black, label="Modeled Trend")
plot!(time_obs, ebm_boot[:, 1], lw=2, label="Bootstrapped Realization")
plot!(time_obs, ebm_boot[:, 2], lw=2, label=false)
plot!(time_obs, ebm_boot[:, 3], lw=2, label=false)
plot!(size=(1200, 500))

Figure 9: Bootstrapped realizations from EBM

Refitting Model

temp_init = (temp_obs[1], temp_lo[1], temp_hi[1])

function fit_ebm_boot(temp_boot, temp_sd, ebm_wrap, temp_init)
    lower = [1.0, 0.5, 0.0, 50.0, 200.0, temp_init[2], -1.0, 0.0]
    upper = [5.0, 1.5, 2.0, 200.0, 1000.0, temp_init[3], 1.0, 10.0]
    p0 = [3.0, 1.0, 1.0, 100.0, 800.0,temp_init[1], 0.0, 5.0]

    result = Optim.optimize(params -> -ar_loglik(params, temp_boot, temp_sd, ebm_wrap), lower, upper, p0)
    θ = result.minimizer
    return θ
end

S = zeros(nboot)
d = zeros(nboot) 
α = zeros(nboot)

for i = 1:nboot
    S[i], α[i], d[i] = fit_ebm_boot(ebm_boot[:, i], temp_sd, ebm_wrap, temp_init)[[1,3,4]]
end

p1 = histogram(S, legend=false, ylabel="Count", xlabel=L"S ($^\circ$C)")
vline!(p1, [θ_ar[1]], color=:red, lw=3, linestyle=:dash)
q90 = quantile(S, [0.95, 0.05])
vspan!(p1, 2 * θ_ar[1] .- q90, linecolor=:gray, fillcolor=:gray, lw=3, alpha=0.3, fillalpha=0.3)
p2 = histogram(α, legend=false, ylabel="Count", xlabel=L"$\alpha$")
vline!(p2, [θ_ar[3]], color=:red, lw=3, linestyle=:dash)
q90 = quantile(α, [0.95, 0.05])
vspan!(p2, 2 * θ_ar[3] .- q90, linecolor=:gray, fillcolor=:gray, lw=3, alpha=0.3, fillalpha=0.3)
p3 = histogram(d, legend=false, ylabel="Count", xlabel="d (m)")
vline!(p3, [θ_ar[4]], color=:red, lw=3, linestyle=:dash)
q90 = quantile(d, [0.95, 0.05])
vspan!(p3, 2 * θ_ar[4] .- q90, linecolor=:gray, fillcolor=:gray, lw=3, alpha=0.3, fillalpha=0.3)

plot(p1, p2, p3, layout=(1, 3), size=(1200, 500))

Figure 10: EBM Bootstrap Parameter Distributions

Upcoming Schedule

Next Classes

Friday: Missing Data

Next Week: Multiple Imputation and Bayesian Statistics

Assessments

HW5 Due Friday (4/17)

References

References (Scroll for Full List)

Brynjarsdóttir, J., & O’Hagan, A. (2014). Learning about physical parameters: the importance of model discrepancy. Inverse Problems, 30, 114007. https://doi.org/10.1088/0266-5611/30/11/114007

Palmer, M. D., Harris, G. R., & Gregory, J. M. (2018). Extending CMIP5 projections of global mean temperature change and sea level rise due to thermal expansion using a physically-based emulator. Environ. Res. Lett., 13, 084003. https://doi.org/10.1088/1748-9326/aad2e4