Breaching 1.5°C: Give me the odds

Abstract

This is a companion notebook to the article “Breaching 1.5°C: Give me the odds” by J. Eduardo Vera-Valdés and Olivia Kvist. It contains the code for the additional exercise considering the Berkeley Earth dataset. The code is written in Julia and it is organized into sections follow those of the replication notebook.

0. Load Packages

In [1]:

Code

using Pkg
Pkg.activate(pwd())
using Plots, Dates, CSV, DataFrames, LongMemory, Statistics, MarSwitching, Random
include("TrendEstimators.jl")
using .TrendEstimators

Random.seed!(123456)
theme(:ggplot2)

  Activating project at `~/Library/CloudStorage/OneDrive-AalborgUniversitet/Research/CLIMATE/Paris Goal/Odds-of-breaching-1.5C`

1. Load Data

1.1 Temperature

In [2]:

Code

rawtemp = CSV.read("data/BerkeleyEarth.csv", DataFrame)
T = size(rawtemp, 1);
dates = collect(Date(1850, 1, 1):Month(1):(Date(1850, 1, 1)+Dates.Month(T - 1)));
rawtemp = DataFrame("Dates" => dates, "Temperature" => rawtemp[!, 3]);
last(rawtemp, 5)

5×2 DataFrame

Row	Dates	Temperature
	Date	Float64
1	2024-05-01	1.209
2	2024-06-01	1.227
3	2024-07-01	1.226
4	2024-08-01	1.379
5	2024-09-01	1.269

1.2 El Niño

Loading the data and removing the missing values. They appear only at the beginning of the time series.

The data has been neatly collected by https://bmcnoldy.earth.miami.edu/tropics/oni/

In [3]:

Code

rawnino = CSV.read("data/Nino_1854_2024.csv", DataFrame)
delete!(rawnino, rawnino.ONI .== -99.99)
first(rawnino, 5)

5×7 DataFrame

Row	YEAR	MON/MMM	NINO34_MEAN	NINO34_CLIM	NINO34_ANOM	ONI	PHASE
	Int64	Int64	Float64	Float64	Float64	Float64	String1
1	1871	2	25.59	26.03	-0.45	-0.26	N
2	1871	3	26.66	26.55	0.11	0.08	N
3	1871	4	27.57	26.99	0.58	0.31	N
4	1871	5	27.4	27.16	0.24	0.36	N
5	1871	6	27.17	26.93	0.24	0.26	N

In [4]:

Code

date_nino = Date.(rawnino[!, 1], rawnino[!, 2])
rawnino.Dates = date_nino
first(date_nino, 3)

3-element Vector{Date}:
 1871-02-01
 1871-03-01
 1871-04-01

1.3 Merging Data

Matching to first non-missing value of El Niño

In [5]:

Code

tempnino = rawnino[rawnino.Dates.>=rawtemp.Dates[1], [:Dates, :ONI]]
tempnino[!, :Temp] = rawtemp[rawtemp.Dates.>=rawnino.Dates[1], :Temperature]
T = length(tempnino.Dates)
first(tempnino, 5)

5×3 DataFrame

Row	Dates	ONI	Temp
	Date	Float64	Float64
1	1871-02-01	-0.26	-0.593
2	1871-03-01	0.08	-0.253
3	1871-04-01	0.31	-0.061
4	1871-05-01	0.36	-0.284
5	1871-06-01	0.26	-0.226

1.4 Rebaseline to pre-industrial levels (1850-1900)

In [6]:

Code

newbaseline = mean(tempnino[(tempnino.Dates.>Date(1850, 1, 1)).&(tempnino.Dates.<Date(1900, 1, 1)), :Temp]);
tempnino.Temp = tempnino.Temp .- newbaseline;

1.5 Temperature plots, pre-industrial baseline against 1951-1980 baseline

In [7]:

Code

theme(:ggplot2)
plot(tempnino.Dates, [rawtemp.Temperature[254:end] tempnino.Temp], label=["Baseline 1951-1980" "Baseline 1870-1900 (Pre-industrial)"], xlabel="", ylabel="", title = "", linewidth=[1 1.1], linestyle=[:dash :dot], xticks=(tempnino.Dates[1:240:end], Dates.format.(tempnino.Dates[1:240:end], "Y")), xlims=(Date(1880, 1, 1), Date(2021, 1, 1)))
plot!(size=(700, 400), fontfamily="Computer Modern", legendfontsize=12, tickfontsize=12, titlefontfamily="Computer Modern", legendfontfamily="Computer Modern", tickfontfamily="Computer Modern", ylabelfontsize=12, xlabelfontsize=14, titlefontsize=16, xlabel="", ylabel="")

Figure 1: Temperature Anomalies (°C) [Berkeley Earth; Rohde and Hausfather (2020)]

In [8]:

Code

maximum(rawtemp.Temperature)

1.543

2. First Look at the Data

In [9]:

Code

plot(tempnino.Dates, [tempnino.Temp tempnino.ONI], label=["Temperature Anomalies (°C)" "El Niño"], xlabel="Date (Monthly)", ylabel="", legend=:topleft)

Figure 2: Raw Temperature Anomalies [Berkeley Earth; Rohde and Hausfather (2020)]

3. Markov Switching Model

Two specifications are considered: one with 3 regimes (El Niño, La Niña, and Neutral) and one with 7 regimes (Very Strong El Niño, Strong El Niño, Moderate El Niño, Neutral, Moderate La Niña, Strong La Niña, and Very Strong La Niña).

In [10]:

Code

nino_model3 = MSModel(tempnino[!, :ONI], 3);
summary_msm(nino_model3);

Markov Switching Model with 3 regimes
=================================================================
# of observations:         1844 AIC:                      2437.933
# of estimated parameters:   12 BIC:                       2504.17
Error distribution:    Gaussian Instant. adj. R^2:          0.7132
Loglikelihood:          -1207.0 Step-ahead adj. R^2:        0.6018
-----------------------------------------------------------------
------------------------------
Summary of regime 1: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |     1.207  |        0.04  |   30.086  |    < 1e-3  
σ            |     0.523  |       0.013  |   38.808  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 9.92 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 2: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |    -0.747  |       0.029  |  -26.072  |    < 1e-3  
σ            |      0.41  |        0.01  |   42.339  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 15.80 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 3: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |     0.183  |       0.033  |    5.493  |    < 1e-3  
σ            |      0.27  |       0.013  |   21.068  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 7.66 periods
-------------------------------------------------------------------
left-stochastic transition matrix: 
          | regime 1   | regime 2   | regime 3
----------------------------------------------------
 regime 1 |   89.918%  |      0.0%  |    5.904%  |
 regime 2 |      0.0%  |   93.671%  |    7.143%  |
 regime 3 |   10.082%  |    6.329%  |   86.953%  |

In [11]:

Code

nino_model7 = MSModel(tempnino[!, :ONI], 7);
summary_msm(nino_model7);

Markov Switching Model with 7 regimes
=================================================================
# of observations:         1844 AIC:                      1546.974
# of estimated parameters:   56 BIC:                      1856.076
Error distribution:    Gaussian Instant. adj. R^2:          0.9174
Loglikelihood:           -717.5 Step-ahead adj. R^2:        0.8096
-----------------------------------------------------------------
------------------------------
Summary of regime 1: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |    -0.035  |       0.015  |   -2.306  |     0.021  
σ            |     0.083  |       0.012  |    6.995  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 1.45 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 2: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |     -0.97  |       0.026  |  -37.467  |    < 1e-3  
σ            |     0.487  |       0.019  |    25.93  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 10.60 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 3: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |     1.615  |        0.04  |   40.215  |    < 1e-3  
σ            |     0.456  |       0.016  |   28.629  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 9.78 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 4: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |      0.32  |       0.021  |   15.508  |    < 1e-3  
σ            |     0.178  |       0.019  |    9.466  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 3.26 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 5: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |    -0.343  |       0.015  |  -23.276  |    < 1e-3  
σ            |     0.127  |        0.01  |   12.706  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 3.21 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 6: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |     0.284  |  221368.518  |      0.0  |       1.0  
σ            |    31.563  |     336.897  |    0.094  |     0.925  
-------------------------------------------------------------------
Expected regime duration: 1.00 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 7: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |     0.574  |       0.029  |   19.852  |    < 1e-3  
σ            |     0.226  |       0.019  |    11.95  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 4.51 periods
-------------------------------------------------------------------
left-stochastic transition matrix: 
          | regime 1   | regime 2   | regime 3   | regime 4   | regime 5   | regime 6   | regime 7
--------------------------------------------------------------------------------------------------------
 regime 1 |   30.832%  |      0.0%  |      0.0%  |   13.485%  |   13.464%  |      0.0%  |      0.0%  |
 regime 2 |      0.0%  |   90.566%  |      0.0%  |      0.0%  |   10.758%  |    100.0%  |      0.0%  |
 regime 3 |      0.0%  |      0.0%  |   89.772%  |      0.0%  |      0.0%  |      0.0%  |    9.318%  |
 regime 4 |   37.278%  |      0.0%  |      0.0%  |   69.317%  |    1.588%  |      0.0%  |    7.994%  |
 regime 5 |     28.7%  |    7.804%  |      0.0%  |    1.712%  |   68.888%  |      0.0%  |    4.883%  |
 regime 6 |      0.0%  |      0.0%  |      0.0%  |      0.0%  |      0.0%  |      0.0%  |      0.0%  |
 regime 7 |     3.19%  |     1.63%  |   10.228%  |   15.485%  |    5.301%  |      0.0%  |   77.805%  |

Looking at both the AIC and BIC, the 7-regime model is preferred. Hence, we continue with the 7-regime model.

4. Trend Model

Fitting the trend model

One example, quadratic trend

Fitting the quadratic trend model to the temperature anomalies and to the temperature anomalies with the El Niño index as an exogenous variable.

In [12]:

Code

qmodel_exo = TrendEstimators.quad_trend_est(tempnino.Temp, tempnino.ONI)
qmodel = TrendEstimators.quad_trend_est(tempnino.Temp)
plot(tempnino.Dates, tempnino.Temp, linewidth=1, label="Temperature Anomalies", xlabel="", ylabel="", legend=:topleft)
plot!(tempnino.Dates, qmodel_exo.Yfit, linestyle=:dash, linewidth=2, label="Quadratic Trend + Niño", color=2)
plot!(tempnino.Dates, qmodel.Yfit, linestyle=:dot, linewidth=2, label="Quadratic Trend", color=3)

Forecasting the quadratic trend model in two ways: only the quadratic trend and the quadratic trend with the long memory component.

In [13]:

Code

h = 800
date_for = collect((tempnino.Dates[1]+Dates.Month(T)):Month(1):(tempnino.Dates[1]+Dates.Month(T - 1)+Dates.Month(h)));
qmodel_forecast = TrendEstimators.quad_trend_forecast(qmodel, h);
plot!(date_for, qmodel_forecast.Yforecastmean, linestyle=:dot, linewidth=2, label="Forecast Quadratic Trend", color=4)
plot!(date_for, qmodel_forecast.Yforecasterr, linestyle=:dash, linewidth=2, label="Forecast Quadratic Trend + LM", color=5)
#plot!(xlims=(Date(2016,1,1),Date(2050,1,1)))

Final forecasting adding the exogenous variable El Niño.

In [14]:

Code

simul_nino = generate_msm(nino_model7, h)[1]
qmodel_exo_forecast = TrendEstimators.quad_trend_forecast(qmodel_exo, h, simul_nino)
plot!(date_for, qmodel_exo_forecast.Yforecasterr, linestyle=:dash, linewidth=2, label="Forecast Quadratic Trend + LM + El Niño", color=6)

A second example, broken linear trend

Fitting the broken linear trend model to the temperature anomalies and to the temperature anomalies with the El Niño index as an exogenous variable.

In [15]:

Code

bmodel = TrendEstimators.broken_trend_est(tempnino.Temp)
bmodel_exo = TrendEstimators.broken_trend_est(tempnino.Temp, tempnino.ONI)
plot(tempnino.Dates, tempnino.Temp, linewidth=1, label="Temperature Anomalies", xlabel="", ylabel="", legend=:topleft)
plot!(tempnino.Dates, bmodel_exo.Yfit, linestyle=:dash, linewidth=2, label="Broken Trend + Niño", color=2)
plot!(tempnino.Dates, bmodel.Yfit, linestyle=:dot, linewidth=2, label="Broken Trend", color=3)

Forecasting the broken linear trend model in two ways: only the broken linear trend and the broken linear trend with the long memory component.

In [16]:

Code

date_for = collect((tempnino.Dates[1]+Dates.Month(T)):Month(1):(tempnino.Dates[1]+Dates.Month(T - 1)+Dates.Month(h)));
bmodel_forecast = TrendEstimators.broken_trend_forecast(bmodel, h);
plot!(date_for, bmodel_forecast.Yforecastmean, linestyle=:dot, linewidth=2, label="Broken Trend", color=4)
plot!(date_for, bmodel_forecast.Yforecasterr, linestyle=:dash, linewidth=2, label="Broken Trend + LM", color=5)

Final forecasting adding the exogenous variable El Niño.

In [17]:

Code

simul_nino = generate_msm(nino_model7, h)[1]
bmodel_exo_forecast = TrendEstimators.broken_trend_forecast(bmodel_exo, h, simul_nino)
plot!(date_for, bmodel_exo_forecast.Yforecasterr, linestyle=:dash, linewidth=2, label="Forecast Broken Trend + LM + El Niño", color=6)
plot!(fontfamily="Computer Modern", legendfontsize=12, tickfontsize=12, titlefontfamily="Computer Modern", legendfontfamily="Computer Modern", tickfontfamily="Computer Modern", ylabelfontsize=12, xlabelfontsize=12, titlefontsize=16, xlabel="", ylabel="", ylims=(-0.9, 2.3))

In [18]:

Code

plot!(fontfamily="Computer Modern", legendfontsize=12, tickfontsize=12, titlefontfamily="Computer Modern", legendfontfamily="Computer Modern", tickfontfamily="Computer Modern", ylabelfontsize=12, xlabelfontsize=12, titlefontsize=16, xlabel="", ylabel="", ylims=(-0.6, 2.8), xlims=(Date(2010, 1, 1), Date(2040, 1, 1)))

Selecting the best model

In [19]:

Code

thisseries = tempnino.Temp;
TrendEstimators.select_trend_model(TrendEstimators.trend_est(thisseries, tempnino.ONI), TrendEstimators.quad_trend_est(thisseries, tempnino.ONI), TrendEstimators.broken_trend_est(thisseries, tempnino.ONI))

("break", (β = [-0.10353755755164067, 0.00033556254446457545, 0.001390561401476513, 0.08465931023982343], σ² = 0.024066882169802795, break_point = 1244, rss = 44.28306319243714, aic = -6868.442182904796, bic = -6846.36341328833, betavar = [7.123035009773671e-5 -7.826445084279888e-8 1.435226571508739e-7 -8.348482209377505e-7; -7.826445084279889e-8 1.1353386006150297e-10 -2.7082182723749503e-10 6.054216120128408e-10; 1.4352265715087397e-7 -2.708218272374951e-10 1.0872951327886004e-9 -9.647761292462743e-11; -8.348482209377503e-7 6.054216120128405e-10 -9.647761292462656e-11 1.822193578978558e-5], res = [-0.17829667078580122, 0.13258360118819423, 0.30477639728857026, 0.07720786923211453, 0.14333823771163226, 0.19953674414318534, 0.15319547126754374, -0.04667987058411556, 0.11514368410019869, -0.16079979570349592  …  0.21751616566113974, 0.09318212445596896, 0.23485379358197522, 0.18175842821997223, 0.23143580767715477, 0.12195406990555924, 0.15939277351957415, 0.16597917370001336, 0.32995194629004576, 0.22753834647048565], Yfit = [-0.12521341566953018, -0.09609368764352565, -0.07628648374390166, -0.07171795568744593, -0.07984832416696369, -0.07104683059851678, -0.059705557722875155, -0.05683021587121588, -0.06665377055553011, -0.08071029075183553  …  1.4999737478835287, 1.4873077890886994, 1.4636361199626935, 1.4357314853246965, 1.4010541058675137, 1.3765358436391093, 1.3570971400250946, 1.349510739844655, 1.3385379672546227, 1.330951567074183]))

5. Forecasting paths

There is only one realization of the Berkeley Earth dataset, so we add the data uncertainty by sampling from the distribution of the estimated parameters.

In [20]:

Code

nsim = 1000
matforecasts = zeros(h, nsim)
uncertainty = 1.96 

thisseries = tempnino.Temp

test = TrendEstimators.trend_est(thisseries, tempnino.ONI)
qtest = TrendEstimators.quad_trend_est(thisseries, tempnino.ONI)
btest = TrendEstimators.broken_trend_est(thisseries, tempnino.ONI)

model_selection = TrendEstimators.select_trend_model(test, qtest, btest)

for jj = 1:nsim

    simul_nino = generate_msm(nino_model7, h)[1]

    if model_selection[1] == "trend"
        matforecasts[:, jj] = TrendEstimators.trend_forecast(test, h, simul_nino, uncertainty).Yforecasterr
    elseif model_selection[1] == "quad"
        matforecasts[:, jj] = TrendEstimators.quad_trend_forecast(qtest, h, simul_nino, uncertainty).Yforecasterr
    elseif model_selection[1] == "break"
        matforecasts[:, jj] = TrendEstimators.broken_trend_forecast(btest, h, simul_nino, uncertainty).Yforecasterr
    else
        @warn "No model selected"
    end

end

Plotting the forecasts

In [21]:

Code

 p1 = plot(tempnino.Dates, tempnino.Temp, linewidth=1, label="Temperature Anomalies (BEST)", xlabel="", ylabel="", legend=:topleft)
 plot!(date_for, matforecasts[:,rand(1:nsim,100)], linestyle=:dot, linewidth=2, label="")
 plot!(fontfamily="Computer Modern", legendfontsize=12, tickfontsize=12, titlefontfamily="Computer Modern", legendfontfamily="Computer Modern", tickfontfamily="Computer Modern", ylabelfontsize=12, xlabelfontsize=12, titlefontsize=16, xlabel="", ylabel="", ylims=(-0.6, 2.8), xlims=(Date(1890, 1, 1), Date(2090, 1, 1)))

Figure 3: Simulated forecast paths for Berkeley Earth temperature anomalies.

6. Estimating the probability of exceeding 1.5°C

Using the middle point of the first 20 years period where the mean temperature exceeds 1.5°C

In [22]:

Code

inicio = T - 119 # 10 years starting 2004
fin = T + h - 120 # 10 years

datejoin = collect((tempnino.Dates[inicio]):Month(1):date_for[h-120])
meantemp = tempnino[inicio:T, :Temp]
dummies = zeros(h, nsim)

for jj = 1:nsim
    completo = [meantemp; matforecasts[:, jj]]
    for ii = 120:h
        dummies[ii, jj] = mean(completo[ii-119:ii+120])
    end
end

pa15 = mean(dummies[:, :] .>= 1.5, dims=2);
pa20 = mean(dummies[:, :] .>= 2, dims=2);

In [23]:

Code

plot(datejoin, pa15, label="1.5C", color=:darkorange, legend=:bottomright, xlims=(datejoin[120], datejoin[end-48]), xticks=(datejoin[1:120:end], Dates.format.(datejoin[1:120:end], "Y")), linewidth=4, linestyle=:dash)
plot!(datejoin, pa20, label="2.0C", color=:red3, linewidth=4, linestyle=:dashdot)
plot!(fontfamily="Computer Modern", legendfontsize=12, tickfontsize=12, titlefontfamily="Computer Modern", legendfontfamily="Computer Modern", tickfontfamily="Computer Modern", ylabelfontsize=12, xlabelfontsize=12, titlefontsize=16, xlabel="", ylabel="", ylims=(0, 1))

Figure 4: Proportion of scenarios that breach the 1.5°C and 2°C thresholds for the Berkeley Earth temperature anomalies for each month. The figure considers 1000 scenarios, each based on the best-fitting model for each realization, with five simulations for El Niño as an exogenous variable each.

In [24]:

Code

results_probabilities = DataFrame("Probability level and period" => String[], "1.5°C Threshold" => Date[], "2°C Threshold" => Date[])

push!(results_probabilities, ["Above 0%, 20-years avg.", datejoin[findfirst(pa15 .> 0)], datejoin[findfirst(pa20 .> 0)]])
push!(results_probabilities, ["Above 50%, 20-years avg.", datejoin[findfirst(pa15 .>= 0.5)], datejoin[findfirst(pa20 .>= 0.5)]])
push!(results_probabilities, ["Above 99%, 20-years avg.", datejoin[findfirst(pa15 .>= 0.99)], datejoin[findfirst(pa20 .>= 0.99)]])

3×3 DataFrame

Row	Probability level and period	1.5°C Threshold	2°C Threshold
	String	Date	Date
1	Above 0%, 20-years avg.	2024-09-01	2041-10-01
2	Above 50%, 20-years avg.	2028-06-01	2053-08-01
3	Above 99%, 20-years avg.	2036-01-01	2063-08-01

Using the middle point of the first 30 years period where the mean temperature exceeds 1.5°C

In [25]:

Code

inicio30 = T - 179 # 15 years starting 2004
fin30 = T + h - 180 # 15 years
datejoin30 = collect((tempnino.Dates[inicio30]):Month(1):date_for[h-180])
meantemp30 = tempnino[inicio30:T, :Temp]
dummies30 = zeros(h, nsim)

for jj = 1:nsim
    completo = [meantemp30; matforecasts[:, jj]]
    for ii = 180:h
        dummies30[ii, jj] = mean(completo[ii-179:ii+180])
    end
end

pa15_30 = mean(dummies30[:, :] .>= 1.5, dims=2);
pa20_30 = mean(dummies30[:, :] .>= 2, dims=2);

Updating table with the results for the 30 years period.

In [26]:

Code

push!(results_probabilities, ["Above 0%, 30-years avg.", datejoin[findfirst(pa15_30 .> 0)], datejoin[findfirst(pa20_30 .> 0)]])
push!(results_probabilities, ["Above 50%, 30-years avg.", datejoin[findfirst(pa15_30 .>= 0.5)], datejoin[findfirst(pa20_30 .>= 0.5)]])
push!(results_probabilities, ["Above 99%, 30-years avg.", datejoin[findfirst(pa15_30 .>= 0.99)], datejoin[findfirst(pa20_30 .>= 0.99)]])

6×3 DataFrame

Row	Probability level and period	1.5°C Threshold	2°C Threshold
	String	Date	Date
1	Above 0%, 20-years avg.	2024-09-01	2041-10-01
2	Above 50%, 20-years avg.	2028-06-01	2053-08-01
3	Above 99%, 20-years avg.	2036-01-01	2063-08-01
4	Above 0%, 30-years avg.	2029-09-01	2046-05-01
5	Above 50%, 30-years avg.	2033-12-01	2058-10-01
6	Above 99%, 30-years avg.	2040-01-01	2068-11-01

Saving the results to a csv file.

In [27]:

Code

CSV.write("tables/ResultsBerkeley.csv", results_probabilities)

"tables/ResultsBerkeley.csv"

Full table with the results.

In [28]:

Code

results_probabilities

In [28]:

6×3 DataFrame

Row	Probability level and period	1.5°C Threshold	2°C Threshold
	String	Date	Date
1	Above 0%, 20-years avg.	2024-09-01	2041-10-01
2	Above 50%, 20-years avg.	2028-06-01	2053-08-01
3	Above 99%, 20-years avg.	2036-01-01	2063-08-01
4	Above 0%, 30-years avg.	2029-09-01	2046-05-01
5	Above 50%, 30-years avg.	2033-12-01	2058-10-01
6	Above 99%, 30-years avg.	2040-01-01	2068-11-01

Table 1: Months to breach the 1.5°C and 2°C thresholds for the Berkeley Earth temperature anomalies at a given probability level.

Rohde, Robert A, and Zeke Hausfather. 2020. “The Berkeley Earth Land/Ocean Temperature Record.” Earth System Science Data 12 (4): 3469–79. https://doi.org/10.5194/essd-12-3469-2020.

Supplementary notebook (alternative data source Berkeley Earth)

0. Load Packages

1. Load Data

1.1 Temperature

1.2 El Niño

1.3 Merging Data

1.4 Rebaseline to pre-industrial levels (1850-1900)

1.5 Temperature plots, pre-industrial baseline against 1951-1980 baseline

2. First Look at the Data

3. Markov Switching Model

4. Trend Model

Fitting the trend model

One example, quadratic trend

A second example, broken linear trend

Selecting the best model

5. Forecasting paths

Plotting the forecasts

6. Estimating the probability of exceeding 1.5°C

Using the middle point of the first 20 years period where the mean temperature exceeds 1.5°C

Using the middle point of the first 30 years period where the mean temperature exceeds 1.5°C