Breaching 1.5°C: Give me the odds

Abstract

This is the replication 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 used to generate the figures and tables in the article. The code is written in Julia and it is organized into sections that correspond to the sections of the paper.

0. Load Packages

In [2]:

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`
WARNING: replacing module TrendEstimators.

1. Load Data

1.1 Temperature

In [3]:

Code

rawtemp = CSV.read("data/HadCRUT.5.0.2.0.analysis.ensemble_series.global.monthly.csv", DataFrame)
first(rawtemp, 5)

5×203 DataFrame

103 columns omitted

Row	Time	Fraction of area represented	Coverage uncertainty (1 sigma)	Realization 1	Realization 2	Realization 3	Realization 4	Realization 5	Realization 6	Realization 7	Realization 8	Realization 9	Realization 10	Realization 11	Realization 12	Realization 13	Realization 14	Realization 15	Realization 16	Realization 17	Realization 18	Realization 19	Realization 20	Realization 21	Realization 22	Realization 23	Realization 24	Realization 25	Realization 26	Realization 27	Realization 28	Realization 29	Realization 30	Realization 31	Realization 32	Realization 33	Realization 34	Realization 35	Realization 36	Realization 37	Realization 38	Realization 39	Realization 40	Realization 41	Realization 42	Realization 43	Realization 44	Realization 45	Realization 46	Realization 47	Realization 48	Realization 49	Realization 50	Realization 51	Realization 52	Realization 53	Realization 54	Realization 55	Realization 56	Realization 57	Realization 58	Realization 59	Realization 60	Realization 61	Realization 62	Realization 63	Realization 64	Realization 65	Realization 66	Realization 67	Realization 68	Realization 69	Realization 70	Realization 71	Realization 72	Realization 73	Realization 74	Realization 75	Realization 76	Realization 77	Realization 78	Realization 79	Realization 80	Realization 81	Realization 82	Realization 83	Realization 84	Realization 85	Realization 86	Realization 87	Realization 88	Realization 89	Realization 90	Realization 91	Realization 92	Realization 93	Realization 94	Realization 95	Realization 96	Realization 97	⋯
	Date	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	⋯
1	1850-01-01	0.53036	0.12529	-0.648141	-0.687369	-0.659397	-0.795922	-0.575723	-0.872706	-0.697845	-0.938079	-0.714489	-0.693846	-0.801877	-0.601654	-0.745235	-0.688707	-0.51959	-0.580684	-0.733635	-0.775802	-0.685629	-0.750843	-0.593316	-0.796245	-0.681935	-0.709096	-0.673113	-0.627142	-0.87027	-0.734045	-0.739397	-0.733723	-0.722726	-0.733075	-0.704837	-0.646199	-0.63273	-0.687739	-0.686664	-0.766938	-0.660519	-0.500427	-0.761633	-0.593218	-0.673419	-0.633932	-0.53979	-0.534376	-0.730824	-0.672585	-0.709929	-0.782612	-0.587371	-0.775013	-0.669417	-0.738468	-0.65991	-0.641778	-0.705	-0.728567	-0.854561	-0.622728	-0.666254	-0.762724	-0.654979	-0.698282	-0.566776	-0.641214	-0.608169	-0.721251	-0.643558	-0.737449	-0.638663	-0.610622	-0.526234	-0.738695	-0.819199	-0.678843	-0.683284	-0.784564	-0.613478	-0.755483	-0.777482	-0.624742	-0.627155	-0.62747	-0.688723	-0.854917	-0.592061	-0.583304	-0.820669	-0.859069	-0.661617	-0.548561	-0.743047	-0.796279	-0.746007	-0.975871	-0.626403	⋯
2	1850-02-01	0.471513	0.15873	-0.198692	-0.379979	-0.353424	-0.499135	-0.307929	-0.437846	-0.226076	-0.577648	-0.405545	-0.392731	-0.597314	-0.396609	-0.332673	-0.431029	-0.186967	-0.323219	-0.320514	-0.585462	-0.386137	-0.31327	-0.374146	-0.393756	-0.377606	-0.270782	-0.281058	-0.297886	-0.440194	-0.272709	-0.375954	-0.316262	-0.395188	-0.45457	-0.310405	-0.380035	-0.414555	-0.284795	-0.391644	-0.302887	-0.429657	-0.208564	-0.397717	-0.427922	-0.29149	-0.208944	-0.232016	-0.34283	-0.303425	-0.159931	-0.409084	-0.379061	-0.441772	-0.384035	-0.168103	-0.374769	-0.331487	-0.281074	-0.473583	-0.308901	-0.519065	-0.37427	-0.373581	-0.392881	-0.259453	-0.358303	-0.189542	-0.268094	-0.19915	-0.513153	-0.427336	-0.380227	-0.389459	-0.32174	-0.139402	-0.278627	-0.342224	-0.304565	-0.209553	-0.431765	-0.262323	-0.418606	-0.363407	-0.292048	-0.373431	-0.360892	-0.476022	-0.406109	-0.389504	-0.239712	-0.329308	-0.478804	-0.303512	-0.305987	-0.337012	-0.644136	-0.404883	-0.496024	-0.366001	⋯
3	1850-03-01	0.443069	0.146062	-0.559715	-0.58127	-0.613536	-0.770914	-0.649228	-0.674574	-0.422177	-0.66901	-0.633439	-0.52613	-0.742085	-0.673435	-0.500491	-0.750384	-0.441357	-0.66936	-0.684599	-0.871666	-0.528859	-0.450323	-0.676369	-0.747881	-0.515235	-0.564603	-0.517356	-0.65774	-0.799563	-0.562602	-0.7321	-0.708126	-0.601581	-0.750345	-0.491204	-0.521375	-0.63429	-0.650293	-0.695271	-0.567457	-0.661922	-0.402203	-0.729825	-0.630483	-0.591241	-0.595347	-0.520872	-0.534415	-0.675771	-0.521642	-0.579516	-0.71082	-0.538983	-0.637005	-0.613244	-0.527521	-0.518286	-0.549176	-0.646482	-0.474831	-0.749775	-0.593508	-0.65672	-0.616531	-0.641366	-0.579132	-0.353847	-0.593789	-0.379657	-0.577604	-0.665159	-0.552994	-0.700469	-0.621582	-0.519006	-0.692925	-0.651994	-0.63646	-0.483905	-0.629932	-0.463239	-0.664408	-0.699684	-0.507791	-0.731221	-0.577734	-0.688031	-0.647881	-0.681086	-0.607305	-0.700216	-0.774448	-0.562562	-0.489596	-0.625427	-0.771725	-0.605319	-0.642851	-0.665689	⋯
4	1850-04-01	0.477059	0.116737	-0.582562	-0.630568	-0.667448	-0.662525	-0.535898	-0.639695	-0.530702	-0.777102	-0.586981	-0.560966	-0.728786	-0.682735	-0.672664	-0.766407	-0.453648	-0.560511	-0.602442	-0.735504	-0.63805	-0.540513	-0.739323	-0.585814	-0.718682	-0.603517	-0.559599	-0.594016	-0.883967	-0.482028	-0.72607	-0.738021	-0.665697	-0.593461	-0.588228	-0.420146	-0.570728	-0.525597	-0.762746	-0.651996	-0.67501	-0.645	-0.770939	-0.635276	-0.575136	-0.474133	-0.439831	-0.45318	-0.473134	-0.518875	-0.570078	-0.59986	-0.661941	-0.67698	-0.512728	-0.550911	-0.497807	-0.615813	-0.453946	-0.432336	-0.790955	-0.860717	-0.531102	-0.672874	-0.730638	-0.601917	-0.429681	-0.665624	-0.535866	-0.594661	-0.656524	-0.685286	-0.49377	-0.605827	-0.551834	-0.556864	-0.688364	-0.615699	-0.659199	-0.586911	-0.499775	-0.698917	-0.642701	-0.524081	-0.586376	-0.729377	-0.701497	-0.737507	-0.544223	-0.426124	-0.786389	-0.706019	-0.53647	-0.61126	-0.617072	-0.755983	-0.614728	-0.72169	-0.544256	⋯
5	1850-05-01	0.487233	0.0959004	-0.405015	-0.493155	-0.550864	-0.544453	-0.47662	-0.672841	-0.49744	-0.602508	-0.587425	-0.501821	-0.602731	-0.661768	-0.543273	-0.610651	-0.416928	-0.468399	-0.599917	-0.601003	-0.517457	-0.532942	-0.559357	-0.613025	-0.540389	-0.424855	-0.623499	-0.473184	-0.582298	-0.488883	-0.597144	-0.580549	-0.596275	-0.515859	-0.479709	-0.451944	-0.450419	-0.513041	-0.654336	-0.641418	-0.526376	-0.440389	-0.501845	-0.485711	-0.625926	-0.433709	-0.496531	-0.421105	-0.541942	-0.478641	-0.559323	-0.551123	-0.677664	-0.617393	-0.537903	-0.498957	-0.40096	-0.499489	-0.400485	-0.450071	-0.640857	-0.553566	-0.649817	-0.483871	-0.56162	-0.587049	-0.408533	-0.486738	-0.515691	-0.4793	-0.607593	-0.468211	-0.523822	-0.490501	-0.496283	-0.508472	-0.526525	-0.486144	-0.449892	-0.609757	-0.281833	-0.624157	-0.528959	-0.380567	-0.620688	-0.520108	-0.608274	-0.551075	-0.464028	-0.444602	-0.423735	-0.538749	-0.574072	-0.397351	-0.445786	-0.760367	-0.51173	-0.664402	-0.439849	⋯

Saving the mean of the ensemble to be used later.

In [4]:

Code

menstemp = reduce(+, eachcol(rawtemp[:, 4:203])) ./ ncol(rawtemp[:, 4:203]);

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 [5]:

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

1.3 Merging Data

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

In [6]:

Code

date_nino = Date.(rawnino[!, 1], rawnino[!, 2])
tempnino = rawtemp[rawtemp.Time.>=date_nino[1], :]
tempnino[!, :ONI] = rawnino.ONI
rename!(tempnino, :Time => :Date)
T = length(tempnino.Date)
first(tempnino, 5)

5×204 DataFrame

104 columns omitted

Row	Date	Fraction of area represented	Coverage uncertainty (1 sigma)	Realization 1	Realization 2	Realization 3	Realization 4	Realization 5	Realization 6	Realization 7	Realization 8	Realization 9	Realization 10	Realization 11	Realization 12	Realization 13	Realization 14	Realization 15	Realization 16	Realization 17	Realization 18	Realization 19	Realization 20	Realization 21	Realization 22	Realization 23	Realization 24	Realization 25	Realization 26	Realization 27	Realization 28	Realization 29	Realization 30	Realization 31	Realization 32	Realization 33	Realization 34	Realization 35	Realization 36	Realization 37	Realization 38	Realization 39	Realization 40	Realization 41	Realization 42	Realization 43	Realization 44	Realization 45	Realization 46	Realization 47	Realization 48	Realization 49	Realization 50	Realization 51	Realization 52	Realization 53	Realization 54	Realization 55	Realization 56	Realization 57	Realization 58	Realization 59	Realization 60	Realization 61	Realization 62	Realization 63	Realization 64	Realization 65	Realization 66	Realization 67	Realization 68	Realization 69	Realization 70	Realization 71	Realization 72	Realization 73	Realization 74	Realization 75	Realization 76	Realization 77	Realization 78	Realization 79	Realization 80	Realization 81	Realization 82	Realization 83	Realization 84	Realization 85	Realization 86	Realization 87	Realization 88	Realization 89	Realization 90	Realization 91	Realization 92	Realization 93	Realization 94	Realization 95	Realization 96	Realization 97	⋯
	Date	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	⋯
1	1871-02-01	0.515577	0.123416	-0.771081	-0.521396	-0.723983	-0.735996	-0.597696	-0.685311	-0.591748	-0.684326	-0.648538	-0.784344	-0.617085	-0.664948	-0.721234	-0.58498	-0.557159	-0.72782	-0.635699	-0.581992	-0.597385	-0.547334	-0.682821	-0.877503	-0.7258	-0.667874	-0.721747	-0.618982	-0.674885	-0.623378	-0.782285	-0.683221	-0.573295	-0.53616	-0.502875	-0.679611	-0.629361	-0.72961	-0.785473	-0.742797	-0.857847	-0.620454	-0.66871	-0.80094	-0.649195	-0.637652	-0.747482	-0.633486	-0.712828	-0.56038	-0.792881	-0.693514	-0.677521	-0.723436	-0.773365	-0.53446	-0.635264	-0.564067	-0.600464	-0.802482	-0.750094	-0.5885	-0.60942	-0.833858	-0.717688	-0.535839	-0.671908	-0.774701	-0.779586	-0.673174	-0.61587	-0.781161	-0.725901	-0.637947	-0.656255	-0.610496	-0.682024	-0.700044	-0.766065	-0.895671	-0.741039	-0.762777	-0.675761	-0.599172	-0.71617	-0.690623	-0.73101	-0.837565	-0.691019	-0.552972	-0.57577	-0.702674	-0.731318	-0.695985	-0.67621	-0.745337	-0.422871	-0.597319	-0.692075	⋯
2	1871-03-01	0.561066	0.134799	-0.168048	-0.123813	-0.206045	-0.187313	-0.181054	-0.242656	-0.0583111	-0.146779	-0.147942	-0.307023	-0.110226	-0.170375	-0.265095	-0.185197	-0.135571	-0.181729	-0.10528	-0.0818921	-0.0610045	-0.17454	-0.305941	-0.150049	-0.270539	-0.152655	-0.293326	-0.140365	-0.1669	-0.197977	-0.356856	-0.241114	0.0143284	-0.0311026	0.0111244	-0.194278	-0.166907	-0.254994	-0.271976	-0.123126	-0.3393	-0.148702	-0.16758	-0.221218	-0.170746	-0.152673	-0.159621	-0.245974	-0.133716	-0.12461	-0.252393	-0.130006	-0.264026	-0.152203	-0.30441	-0.0116849	-0.170787	-0.097697	-0.125526	-0.173275	-0.361964	-0.128557	-0.176711	-0.370949	-0.171573	-0.076722	-0.139446	-0.282856	-0.226627	-0.325808	-0.0656022	-0.180601	-0.167828	-0.171125	-0.0849002	-0.194156	-0.222889	-0.234873	-0.339556	-0.304459	-0.163763	-0.268489	-0.222932	-0.182512	-0.203238	-0.137038	-0.157819	-0.298329	-0.188355	-0.158444	-0.214296	-0.170551	-0.163948	-0.13393	-0.062116	-0.247188	0.00207629	-0.141833	-0.206706	⋯
3	1871-04-01	0.649393	0.0914514	-0.252919	-0.243722	-0.274539	-0.197268	-0.233354	-0.317776	-0.134673	-0.204517	-0.231595	-0.370374	-0.198386	-0.223869	-0.328581	-0.276489	-0.179403	-0.274998	-0.0689078	-0.22583	-0.201089	-0.173408	-0.244664	-0.231315	-0.244879	-0.336012	-0.288421	-0.223314	-0.188126	-0.299007	-0.359414	-0.18775	-0.106765	-0.205048	-0.0566086	-0.151868	-0.234503	-0.282121	-0.388991	-0.260582	-0.336696	-0.182403	-0.15996	-0.146656	-0.263237	-0.183453	-0.218428	-0.24586	-0.151798	-0.129426	-0.381452	-0.219743	-0.283651	-0.327327	-0.361677	-0.135703	-0.154314	-0.131061	-0.159177	-0.183503	-0.26722	-0.112593	-0.19456	-0.303903	-0.174837	-0.146504	-0.168074	-0.279533	-0.24531	-0.254839	-0.173626	-0.307909	-0.289818	-0.26642	-0.189888	-0.203676	-0.237189	-0.26805	-0.322924	-0.367398	-0.212629	-0.297684	-0.214526	-0.306628	-0.334873	-0.272869	-0.23807	-0.383339	-0.189262	-0.202215	-0.186874	-0.329263	-0.295474	-0.18122	-0.167671	-0.231628	-0.0259454	-0.197327	-0.263817	⋯
4	1871-05-01	0.569534	0.100992	-0.390965	-0.358866	-0.3872	-0.412647	-0.289133	-0.444212	-0.350638	-0.411561	-0.314236	-0.409788	-0.293004	-0.423972	-0.352542	-0.352773	-0.324253	-0.361603	-0.212133	-0.319255	-0.243624	-0.37716	-0.417869	-0.422049	-0.414478	-0.469378	-0.358084	-0.424911	-0.38943	-0.368253	-0.622229	-0.384898	-0.217172	-0.320815	-0.179714	-0.381533	-0.336759	-0.324799	-0.467845	-0.43728	-0.485023	-0.291502	-0.362939	-0.404393	-0.334562	-0.339195	-0.397332	-0.224567	-0.353203	-0.318229	-0.39959	-0.323261	-0.399928	-0.247869	-0.400006	-0.279719	-0.352996	-0.326335	-0.404642	-0.374809	-0.419938	-0.301217	-0.392191	-0.457006	-0.280861	-0.312897	-0.298531	-0.445634	-0.504316	-0.39271	-0.366774	-0.312235	-0.350363	-0.365224	-0.25892	-0.422766	-0.326767	-0.456135	-0.382438	-0.484093	-0.378428	-0.53167	-0.401953	-0.386346	-0.543575	-0.329703	-0.406563	-0.612432	-0.487918	-0.345462	-0.330318	-0.512348	-0.489328	-0.403575	-0.346914	-0.384715	-0.147214	-0.288468	-0.412974	⋯
5	1871-06-01	0.581247	0.0946694	-0.28967	-0.135594	-0.332532	-0.360571	-0.204	-0.314435	-0.246807	-0.332408	-0.265587	-0.310019	-0.31175	-0.356402	-0.295284	-0.275481	-0.11862	-0.359438	-0.192357	-0.22494	-0.184043	-0.218272	-0.347667	-0.389377	-0.26043	-0.392255	-0.245163	-0.351099	-0.18067	-0.355544	-0.456089	-0.300251	-0.19407	-0.221448	-0.124136	-0.22837	-0.304018	-0.319256	-0.389229	-0.287898	-0.447419	-0.274095	-0.361008	-0.279793	-0.292989	-0.221987	-0.326887	-0.143718	-0.255443	-0.194447	-0.353247	-0.284191	-0.309094	-0.343588	-0.331752	-0.170288	-0.425651	-0.222224	-0.225238	-0.31519	-0.325241	-0.248408	-0.208976	-0.371835	-0.279356	-0.214969	-0.334997	-0.328621	-0.394733	-0.30453	-0.299655	-0.217936	-0.271893	-0.184974	-0.188288	-0.341147	-0.361841	-0.250004	-0.311449	-0.479353	-0.295131	-0.451731	-0.300321	-0.275753	-0.429856	-0.265486	-0.358396	-0.463469	-0.306513	-0.254791	-0.28826	-0.28799	-0.382331	-0.295289	-0.304173	-0.455725	-0.124766	-0.190635	-0.327341	⋯

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

In [7]:

Code

for ii = 1:200
    newbaseline = mean(tempnino[(tempnino.Date.>=Date(1850, 1, 1)).&(tempnino.Date.<Date(1900, 1, 1)), 3+ii])
    tempnino[!, 3+ii] = tempnino[!, 3+ii] .- newbaseline
end

1.5 Plotting baselines on mean ensemble temperature

In [8]:

Code

oldbase = mean(menstemp[(rawtemp.Time.>=Date(1850, 1, 1)).&(rawtemp.Time.<Date(1900, 1, 1))]);
menstempind = menstemp .- oldbase;

In [9]:

Code

plot(rawtemp.Time, [menstemp menstempind], label=["Baseline 1961-1990" "Baseline 1870-1900 (pre-industrial)"], xlabel="Date (monthly)", ylabel="°C", linewidth=[1 1], linestyle=[:dash :solid], xticks=(rawtemp.Time[372:240:end], Dates.format.(rawtemp.Time[372:240:end], "Y")), xlims=(Date(1880, 1, 1), Date(2021, 1, 1)), title="Temperature anomalies", color=[2 1])
plot!(fontfamily="Computer Modern", legendfontsize=12, tickfontsize=12, titlefontfamily="Computer Modern", legendfontfamily="Computer Modern", tickfontfamily="Computer Modern", ylabelfontsize=12, xlabelfontsize=12, titlefontsize=12)

Figure 1: Temperature anomalies (°C) in the HadCRUT5 dataset (Morice et al. 2021). The dashed line presents the data according to the baseline period in HadCRUT5 (1961-1990). The solid line represents the temperature anomalies with the pre-industrial baseline period (1870-1900).

In [10]:

Code

savefig("figures/BaselinesHadCRUT5.png")

"/Users/jeddy/Library/CloudStorage/OneDrive-AalborgUniversitet/Research/CLIMATE/Paris Goal/Odds-of-breaching-1.5C/figures/BaselinesHadCRUT5.png"

2. First Look at the Data

In [11]:

Code

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

Figure 2: Raw Temperature Anomalies [HadCRUT5; Morice et al. (2021)]

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 [12]:

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 [13]:

Code

nino_model5 = MSModel(tempnino[!, :ONI], 5);
summary_msm(nino_model5);

Markov Switching Model with 5 regimes
=================================================================
# of observations:         1844 AIC:                      2342.062
# of estimated parameters:   30 BIC:                      2507.653
Error distribution:    Gaussian Instant. adj. R^2:          0.6366
Loglikelihood:          -1141.0 Step-ahead adj. R^2:         0.538
-----------------------------------------------------------------
------------------------------
Summary of regime 1: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |    -0.613  |         NaN  |      NaN  |       NaN  
σ            |     0.327  |         NaN  |      NaN  |       NaN  
-------------------------------------------------------------------
Expected regime duration: 11.65 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 2: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |     1.169  |         NaN  |      NaN  |       NaN  
σ            |     0.453  |         NaN  |      NaN  |       NaN  
-------------------------------------------------------------------
Expected regime duration: 11.07 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 3: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |     0.254  |         NaN  |      NaN  |       NaN  
σ            |      0.35  |         NaN  |      NaN  |       NaN  
-------------------------------------------------------------------
Expected regime duration: 11.97 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 4: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |    -1.039  |         NaN  |      NaN  |       NaN  
σ            |     0.383  |         NaN  |      NaN  |       NaN  
-------------------------------------------------------------------
Expected regime duration: 9.40 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 5: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |    -0.041  |         NaN  |      NaN  |       NaN  
σ            |     0.737  |         NaN  |      NaN  |       NaN  
-------------------------------------------------------------------
Expected regime duration: 1.00 periods
-------------------------------------------------------------------
left-stochastic transition matrix: 
          | regime 1   | regime 2   | regime 3   | regime 4   | regime 5
------------------------------------------------------------------------------
 regime 1 |   91.413%  |      0.0%  |      0.0%  |     2.28%  |   65.755%  |
 regime 2 |      0.0%  |   90.969%  |      0.0%  |      0.0%  |   34.245%  |
 regime 3 |    4.178%  |    7.626%  |   91.648%  |      0.0%  |      0.0%  |
 regime 4 |    1.213%  |      0.0%  |      0.0%  |   89.364%  |      0.0%  |
 regime 5 |    3.196%  |    1.405%  |    8.352%  |    8.356%  |      0.0%  |

In [14]:

Code

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

Markov Switching Model with 7 regimes
=================================================================
# of observations:         1844 AIC:                      1393.647
# of estimated parameters:   56 BIC:                       1702.75
Error distribution:    Gaussian Instant. adj. R^2:          0.8172
Loglikelihood:           -640.8 Step-ahead adj. R^2:         0.725
-----------------------------------------------------------------
------------------------------
Summary of regime 1: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |     0.723  |       0.004  |  170.438  |    < 1e-3  
σ            |     0.133  |       0.019  |    6.918  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 1.93 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 2: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |    -0.973  |       0.023  |  -43.151  |    < 1e-3  
σ            |     0.343  |       0.008  |   42.215  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 10.30 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 3: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |    -0.383  |       0.023  |   -16.52  |    < 1e-3  
σ            |     0.195  |       0.011  |   17.499  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 5.75 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 4: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |      0.24  |       0.019  |    12.89  |    < 1e-3  
σ            |     0.242  |       0.008  |   30.771  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 1.00 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 5: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |     0.223  |       0.022  |   10.314  |    < 1e-3  
σ            |     0.238  |       0.012  |   20.184  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 4.91 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 6: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |     1.222  |       0.017  |   70.628  |    < 1e-3  
σ            |     0.181  |       0.015  |   12.413  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 2.72 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 7: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |     1.841  |       0.043  |   42.585  |    < 1e-3  
σ            |     0.429  |       0.011  |   39.005  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 3.16 periods
-------------------------------------------------------------------
left-stochastic transition matrix: 
          | regime 1   | regime 2   | regime 3   | regime 4   | regime 5   | regime 6   | regime 7
--------------------------------------------------------------------------------------------------------
 regime 1 |   48.261%  |      0.0%  |      0.0%  |      0.0%  |   11.016%  |   22.165%  |      0.0%  |
 regime 2 |      0.0%  |   90.291%  |    5.144%  |      0.0%  |      0.0%  |      0.0%  |      0.0%  |
 regime 3 |      0.0%  |    9.709%  |   82.616%  |   12.296%  |    9.367%  |      0.0%  |      0.0%  |
 regime 4 |    26.83%  |      0.0%  |   11.157%  |      0.0%  |      0.0%  |      0.0%  |      0.0%  |
 regime 5 |      0.0%  |      0.0%  |      0.0%  |    83.78%  |   79.617%  |      0.0%  |      0.0%  |
 regime 6 |   15.444%  |      0.0%  |      0.0%  |      0.0%  |      0.0%  |   63.231%  |   31.689%  |
 regime 7 |    9.465%  |      0.0%  |    1.083%  |    3.925%  |      0.0%  |   14.604%  |   68.311%  |

Looking at the 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 [15]:

Code

qmodel_exo = TrendEstimators.quad_trend_est(tempnino[:, 10], tempnino.ONI)
qmodel = TrendEstimators.quad_trend_est(tempnino[:, 10])
plot(tempnino.Date, tempnino[:, 10], linewidth=1, label="Temperature Anomalies", xlabel="", ylabel="", legend=:topleft)
plot!(tempnino.Date, qmodel_exo.Yfit, linestyle=:dash, linewidth=2, label="Quadratic Trend + Niño", color=2)
plot!(tempnino.Date, 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 [16]:

Code

h = 800
date_for = collect((tempnino.Date[1]+Dates.Month(T)):Month(1):(tempnino.Date[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 [17]:

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)

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)), legend=:topleft)

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 [19]:

Code

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

Date of the break

In [20]:

Code

tempnino.Date[bmodel.break_point]

1976-08-01

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 [21]:

Code

date_for = collect((tempnino.Date[1]+Dates.Month(T)):Month(1):(tempnino.Date[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 [22]:

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)

In [23]:

Code

fulldate = collect((tempnino.Date[1]):Month(1):date_for[end])

plot(tempnino.Date, tempnino[:, 10], linewidth=1, label="Temperature anomalies", xlabel="Date (monthly)", ylabel="°C", legend=:bottomright, title="Forecast for realization 10 of HadCRUT5")
plot!(tempnino.Date, bmodel.Yfit, linestyle=:dot, linewidth=2, label="Fitted trend", color=2)
plot!(tempnino.Date, bmodel_exo.Yfit, linestyle=:dash, linewidth=2, label="   ' ' + El Niño", color=3)
plot!(date_for, bmodel_forecast.Yforecastmean, linestyle=:dot, linewidth=2, label="Forecasted trend", color=4)
plot!(date_for, bmodel_forecast.Yforecasterr, linestyle=:dash, linewidth=2, label="   ' ' + long-range dependence", color=5)
plot!(date_for, bmodel_exo_forecast.Yforecasterr, linestyle=:dashdot, linewidth=2, label="   ' ' +  ' ' + El Niño", color=6)
vline!([tempnino.Date[end]], label="Last observation", color=:gray, linewidth=2, linestyle=:dash)
plot!(fontfamily="Computer Modern", legendfontsize=10, tickfontsize=12, titlefontfamily="Computer Modern", legendfontfamily="Computer Modern", tickfontfamily="Computer Modern", ylabelfontsize=12, xlabelfontsize=12, titlefontsize=14, ylims=(0.66, 1.8), xlims=(Date(2010, 1, 1), Date(2045, 1, 1)), yticks=0.75:0.25:1.75,  xticks=(fulldate[1660:60:end], Dates.format.(fulldate[1660:60:end], "Y")) )

Figure 3: Forecast of temperature anomalies for realization 10 of the HadCRUT5 dataset. The forecast is based on the broken trend model with long-range dependence and El Niño as an exogenous variable. A simulated El Niño series using a Markov-switching model with 7 states was used to generate the forecast.

In [24]:

Code

savefig("figures/Forecast-HadCRUT5-Realization10.png")

"/Users/jeddy/Library/CloudStorage/OneDrive-AalborgUniversitet/Research/CLIMATE/Paris Goal/Odds-of-breaching-1.5C/figures/Forecast-HadCRUT5-Realization10.png"

Selecting the best trend model

AIC

In [25]:

Code

thisseries = tempnino[:, 10];
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.10064851870194579, 0.0003142358906008469, 0.0013711294769147653, 0.0821471222095873], σ² = 0.02742877070247553, break_point = 1285, rss = 50.46893809255498, aic = -6627.3285309908915, bic = -6605.249761374425, betavar = [7.905587713826891e-5 -8.469404729710748e-8 1.644954055553638e-7 -8.60674321622298e-7; -8.469404729710748e-8 1.19976350771688e-10 -3.062215896694223e-10 5.206050039230095e-10; 1.6449540555536383e-7 -3.062215896694224e-10 1.3899969031667437e-9 6.445658024938706e-10; -8.606743216222978e-7 5.206050039230093e-10 6.445658024938709e-10 2.0767636360366003e-5], res = [-0.14484373537669837, 0.3603485771814411, 0.26477896318263516, 0.04439194118155497, 0.15612325751191286, 0.38442459040035326, 0.17195290762250606, 0.042622878065617605, 0.13001567684016727, 0.034612581725196254  …  0.13299221765294367, 0.0666538630610578, 0.21174433435641848, 0.21834136176225827, 0.19303945894486518, 0.10606840146232188, 0.177100816647203, 0.17435913472274223, 0.29678643768666446, 0.22408215576220347], Yfit = [-0.12169253458583763, -0.0934482771439771, -0.07424020314517119, -0.06981861114409098, -0.07771908747444885, -0.06919013936288929, -0.05819677758504208, -0.05541812802815361, -0.06496154680270325, -0.07861232168773223  …  1.3911039123845204, 1.3788242669764061, 1.3558654956810454, 1.3287993682752057, 1.2951614710925987, 1.271381228575142, 1.252529813390261, 1.2451789953147219, 1.2345422923507996, 1.2271914742752605]))

BIC

In [26]:

Code

TrendEstimators.select_trend_model_bic(TrendEstimators.trend_est(thisseries, tempnino.ONI), TrendEstimators.quad_trend_est(thisseries, tempnino.ONI), TrendEstimators.broken_trend_est(thisseries, tempnino.ONI))

("break", (β = [-0.10064851870194579, 0.0003142358906008469, 0.0013711294769147653, 0.0821471222095873], σ² = 0.02742877070247553, break_point = 1285, rss = 50.46893809255498, aic = -6627.3285309908915, bic = -6605.249761374425, betavar = [7.905587713826891e-5 -8.469404729710748e-8 1.644954055553638e-7 -8.60674321622298e-7; -8.469404729710748e-8 1.19976350771688e-10 -3.062215896694223e-10 5.206050039230095e-10; 1.6449540555536383e-7 -3.062215896694224e-10 1.3899969031667437e-9 6.445658024938706e-10; -8.606743216222978e-7 5.206050039230093e-10 6.445658024938709e-10 2.0767636360366003e-5], res = [-0.14484373537669837, 0.3603485771814411, 0.26477896318263516, 0.04439194118155497, 0.15612325751191286, 0.38442459040035326, 0.17195290762250606, 0.042622878065617605, 0.13001567684016727, 0.034612581725196254  …  0.13299221765294367, 0.0666538630610578, 0.21174433435641848, 0.21834136176225827, 0.19303945894486518, 0.10606840146232188, 0.177100816647203, 0.17435913472274223, 0.29678643768666446, 0.22408215576220347], Yfit = [-0.12169253458583763, -0.0934482771439771, -0.07424020314517119, -0.06981861114409098, -0.07771908747444885, -0.06919013936288929, -0.05819677758504208, -0.05541812802815361, -0.06496154680270325, -0.07861232168773223  …  1.3911039123845204, 1.3788242669764061, 1.3558654956810454, 1.3287993682752057, 1.2951614710925987, 1.271381228575142, 1.252529813390261, 1.2451789953147219, 1.2345422923507996, 1.2271914742752605]))

Table with the AIC and BIC values for the different trend models.

In [27]:

Code

[TrendEstimators.trend_est(thisseries, tempnino.ONI).aic TrendEstimators.quad_trend_est(thisseries, tempnino.ONI).aic TrendEstimators.broken_trend_est(thisseries, tempnino.ONI).aic;
TrendEstimators.trend_est(thisseries, tempnino.ONI).bic TrendEstimators.quad_trend_est(thisseries, tempnino.ONI).bic TrendEstimators.broken_trend_est(thisseries, tempnino.ONI).bic]

2×3 Matrix{Float64}:
 -5613.2   -6551.17  -6627.33
 -5596.64  -6529.09  -6605.25

Computing the breach point

In [28]:

Code

thisreal = [tempnino[:, 10]; bmodel_exo_forecast.Yforecasterr]

anim = @animate for ii = 30:6:160
    plot(fulldate[1:T], thisreal[1:T], linewidth=1, label="Temperature anomalies", xlabel="Date (monthly)", ylabel="°C", legend=:bottomright, title="")
    plot!(fulldate[T+1:end], thisreal[T+1:end], linewidth=1, label="Predicted path", color=2, linestyle=:dash)
    vspan!(fulldate[[T-119+ii,T+120+ii]], linecolor = :darkgray, fillcolor = :darkgray, label = "20-year period", alpha = 0.5)
    plot!(fulldate[[T-119+ii,T+120+ii]], [ mean(thisreal[T-119+ii:T+120+ii]), mean(thisreal[T-119+ii:T+120+ii])], label="20-year average", linewidth=4, linestyle=:dot, color = :black)
    plot!(fontfamily="Computer Modern", legendfontsize=10, tickfontsize=12, titlefontfamily="Computer Modern", legendfontfamily="Computer Modern", tickfontfamily="Computer Modern", ylabelfontsize=12, xlabelfontsize=12, titlefontsize=14, ylims=(0.66, 1.8), xlims=(Date(2010, 1, 1), Date(2045, 1, 1),), yticks=0.75:0.25:1.75,  xticks=(fulldate[1660:60:end], Dates.format.(fulldate[1660:60:end], "Y")) )
end

gif(anim, "figures/Breaching-Realization10.gif", fps = 5)

[ Info: Saved animation to /Users/jeddy/Library/CloudStorage/OneDrive-AalborgUniversitet/Research/CLIMATE/Paris Goal/Odds-of-breaching-1.5C/figures/Breaching-Realization10.gif

Figure 4: Breaching of the 1.5°C threshold for realization 10 of the HadCRUT5 dataset. The figure shows the temperature anomalies and the forecasted path for the next several months. The forecasted path is the same as in the previous figure. The 20-year period is highlighted in gray, and the 20-year average is shown as a black dashed line.

Finding the breach point of the broken linear trend model.

In [29]:

Code

breachreal10 = TrendEstimators.rolling_sample_mean(thisreal)
fulldate[breachreal10]

2033-03-01

In [30]:

Code

ps = plot(fulldate[1:T], thisreal[1:T], linewidth=1, label="Temperature anomalies", xlabel="Date (monthly)", ylabel="°C", legend=:bottomright, title="")
plot!(fulldate[T+1:end], thisreal[T+1:end], linewidth=1, label="Predicted path", color=2, linestyle=:dash)
vspan!(fulldate[[breachreal10-119,breachreal10+120]], linecolor = :darkgray, fillcolor = :darkgray, label = "20-year period", alpha = 0.5)
plot!(fulldate[[breachreal10-119,breachreal10+120]], [ mean(thisreal[breachreal10-119:breachreal10+120]), mean(thisreal[breachreal10-119:breachreal10+120])], label="20-year average", linewidth=4, linestyle=:dot, color = :black)
plot!(fontfamily="Computer Modern", legendfontsize=10, tickfontsize=12, titlefontfamily="Computer Modern", legendfontfamily="Computer Modern", tickfontfamily="Computer Modern", ylabelfontsize=12, xlabelfontsize=12, titlefontsize=14, ylims=(0.66, 1.8), xlims=(Date(2010, 1, 1), Date(2045, 1, 1),), yticks=0.75:0.25:1.75,  xticks=(fulldate[1660:60:end], Dates.format.(fulldate[1660:60:end], "Y")) )

Figure 5: Breaching of the 1.5°C threshold for realization 10 of the HadCRUT5 dataset. The figure shows the temperature anomalies and the forecasted path for the next several months. The 20-year period is highlighted in gray, and the 20-year average is shown as a black dashed line.

In [31]:

Code

savefig("figures/Breaching-Realization10.png")

"/Users/jeddy/Library/CloudStorage/OneDrive-AalborgUniversitet/Research/CLIMATE/Paris Goal/Odds-of-breaching-1.5C/figures/Breaching-Realization10.png"

In [32]:

Code

breakdates = Vector{Date}()

Date[]

5. Forecasting paths

Columns 4 to 203 are time series of the temperature anomalies for each realization. We fit the trend model to each realization. Column 204 is the El Niño index.

In [33]:

Code

# label: forecast-simulations

nsim = 5
matforecasts = zeros(h, 200 * nsim)
uncertainty = 1.96 #1 #2.326

models = Vector{String}()
breakdates = Vector{Date}()

for ii = 1:200

    thisseries = tempnino[!, 3+ii]

    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)
    push!(models, model_selection[1])

    if model_selection[1] == "break"
        push!(breakdates, fulldate[model_selection[2].break_point])
    end

    for jj = 1:nsim

        simul_nino = generate_msm(nino_model7, h)[1]

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

    end
end

Models selected

In [34]:

Code

number_break_models = sum(models.== "break")
number_quad_models = sum(models.== "quad")
number_trend_models = sum(models.== "trend")

display([number_break_models number_quad_models number_trend_models])

1×3 Matrix{Int64}:
 190  10  0

Break dates for the broken linear trend model

In [35]:

Code

plot(breakdates, seriestype=:scatter, ylabel="Date", xlabel="Number of realization", title="", legend=:topleft)

Plotting the forecasts

In [36]:

Code

plot(rawtemp.Time, menstempind, linewidth=1, label="Temperature anomalies", legend=:topleft, title="", xlabel="Date (monthly)", ylabel="°C")
plot!(date_for, matforecasts[:,rand(1:200*nsim,100)], linestyle=:dot, linewidth=0.8, label="", xticks=(fulldate[108:240:end], Dates.format.(fulldate[108:240:end], "Y")) )
plot!(date_for, zeros(size(date_for,1)), color=nothing, label="Simulated forecast paths")
plot!(fontfamily="Computer Modern", legendfontsize=12, tickfontsize=12, titlefontfamily="Computer Modern", legendfontfamily="Computer Modern", tickfontfamily="Computer Modern", ylabelfontsize=12, xlabelfontsize=12, titlefontsize=14, ylims=(-0.65, 2.8), xlims=(Date(1870, 1, 1), Date(2060, 1, 1)))

Figure 6: Simulated forecast paths for HadCRUT5 temperature anomalies. One hundred paths of a total of 1000 paths are shown to ease visualization. The forecasts are based on the best-fitting model for each realization, with El Niño as an exogenous variable. For ease of visualization, the mean of all temperature anomaly realizations is shown as a solid line.

In [37]:

Code

savefig("figures/Paths-HadCRUT5-All-Realizations.png")

"/Users/jeddy/Library/CloudStorage/OneDrive-AalborgUniversitet/Research/CLIMATE/Paris Goal/Odds-of-breaching-1.5C/figures/Paths-HadCRUT5-All-Realizations.png"

Coverage of the predictions

In [38]:

Code

quantilesforecasts = zeros(h, 7)

for ii = 1:h
    quantilesforecasts[ii, :] = quantile(matforecasts[ii, :], [0.005, 0.025, 0.05, 0.5, 0.95, 0.975, 0.995 ])
end

In [39]:

Code

plot(rawtemp.Time, menstempind, linewidth=1, label="Temperature anomalies", legend=:topleft, title="Simulated forecast paths for temperature anomalies", xlabel="Date (monthly)", ylabel="°C")
plot!(date_for, matforecasts[:,rand(1:200*nsim,100)], linestyle=:dot, linewidth=0.5, label="", xticks=(fulldate[108:240:end], Dates.format.(fulldate[108:240:end], "Y")) )
plot!(date_for, quantilesforecasts[:, 2], fillrange=quantilesforecasts[:, 6], fillalpha = 0.35,label="Prediction interval 95%", color=5, linestyle=:auto)
plot!(date_for, quantilesforecasts[:, 6], label="", color=5, linestyle=:auto)
plot!(date_for, quantilesforecasts[:, 1], fillrange=quantilesforecasts[:, 7], fillalpha = 0.35,label="Prediction interval 99%", color=4, linestyle=:dot)
plot!(date_for, quantilesforecasts[:, 7], label="", color=4, linestyle=:dot)
plot!(fontfamily="Computer Modern", legendfontsize=12, tickfontsize=12, titlefontfamily="Computer Modern", legendfontfamily="Computer Modern", tickfontfamily="Computer Modern", ylabelfontsize=12, xlabelfontsize=12, titlefontsize=14, ylims=(-0.65, 2.55), xlims=(Date(1870, 1, 1), Date(2060, 1, 1)))

Figure 7: Simulated forecast paths for HadCRUT5 temperature anomalies. One hundred paths of a total of 1000 paths are shown to ease visualization. The forecasts are based on the best-fitting model for each realization, with El Niño as an exogenous variable. For ease of visualization, the mean of all temperature anomaly realizations is shown as a solid line.

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 [40]:

Code

PA_start = Date("2016-11-01") # PA enters into force
inicio = findfirst(isequal(PA_start),tempnino.Date)-119
fin = T + h - 120 # 10 years

datejoin = collect(PA_start-Month(119):Month(1):date_for[h-120])
meantemp = mean(Matrix(tempnino[inicio:T, 4:203]), dims=2)
tsize = fin-inicio+1
dummies = zeros(tsize, 200 * nsim)

for jj = 1:(200*nsim)
    completo = [meantemp; matforecasts[:, jj]]
    for ii = 120:tsize
        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 [41]:

Code

xls = (Date(2020, 11, 1), Date(2080, 1, 1))

plot(datejoin, pa15, label="1.5C", color=:darkorange, legend=:bottomright, xticks=(collect(xls[1]:Dates.Month(120):xls[2]), Dates.format.(collect(xls[1]:Dates.Month(120):xls[2]), "Y")), linewidth=4, linestyle=:dash, title="Probability of breaching the 1.5°C and 2°C thresholds", xlabel="Date (monthly)", ylabel="Probability")
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=14, ylims=(0, 1), xlims=xls)

Figure 8: Proportion of scenarios that breach the 1.5°C and 2°C thresholds for the HadCRUT5 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 [42]:

Code

savefig("figures/Coverage-HadCRUT5-All-Realizations.png")

"/Users/jeddy/Library/CloudStorage/OneDrive-AalborgUniversitet/Research/CLIMATE/Paris Goal/Odds-of-breaching-1.5C/figures/Coverage-HadCRUT5-All-Realizations.png"

Data frame to save probability paths

In [43]:

Code

results_probability_paths = DataFrame("Date (month)" => datejoin, "1.5°C Threshold" => vec(pa15), "2°C Threshold" => vec(pa20))
CSV.write("tables/ProbabilityPathsHadCRUT5.csv", results_probability_paths)

"tables/ProbabilityPathsHadCRUT5.csv"

Data frame to store the probabilities

In [44]:

Code

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

push!(results_probabilities, ["Above 1%, 20-years avg.", datejoin[findfirst(pa15 .> 0.01)], datejoin[findfirst(pa20 .> 0.01)]])
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 1%, 20-years avg.	2025-04-01	2044-04-01
2	Above 50%, 20-years avg.	2030-05-01	2055-09-01
3	Above 99%, 20-years avg.	2040-10-01	2067-04-01

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

In [45]:

Code

PA_start = Date("2016-11-01") # PA enters into force
inicio30 = findfirst(isequal(PA_start),tempnino.Date)-179
fin30 = T + h - 180 # 15 years

datejoin30 = collect(PA_start-Month(179):Month(1):date_for[h-180])
meantemp30 = mean(Matrix(tempnino[inicio30:T, 4:203]), dims=2)
tsize30 = fin30-inicio30+1
dummies30 = zeros(tsize30, 200 * nsim)

for jj = 1:(200*nsim)
    completo30 = [meantemp30; matforecasts[:, jj]]
    for ii = 180:h
        dummies30[ii, jj] = mean(completo30[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 [46]:

Code

push!(results_probabilities, ["Above 1%, 30-years avg.", datejoin[findfirst(pa15_30 .> 0.01)], datejoin[findfirst(pa20_30 .> 0.01)]])
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 1%, 20-years avg.	2025-04-01	2044-04-01
2	Above 50%, 20-years avg.	2030-05-01	2055-09-01
3	Above 99%, 20-years avg.	2040-10-01	2067-04-01
4	Above 1%, 30-years avg.	2031-02-01	2049-04-01
5	Above 50%, 30-years avg.	2035-07-01	2060-09-01
6	Above 99%, 30-years avg.	2045-08-01	2072-08-01

Saving the results to a csv file.

In [47]:

Code

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

"tables/ResultsHadCRUT5.csv"

Full table with the results.

In [48]:

Code

results_probabilities

In [48]:

6×3 DataFrame

Row	Probability level and period	1.5°C Threshold	2°C Threshold
	String	Date	Date
1	Above 1%, 20-years avg.	2025-04-01	2044-04-01
2	Above 50%, 20-years avg.	2030-05-01	2055-09-01
3	Above 99%, 20-years avg.	2040-10-01	2067-04-01
4	Above 1%, 30-years avg.	2031-02-01	2049-04-01
5	Above 50%, 30-years avg.	2035-07-01	2060-09-01
6	Above 99%, 30-years avg.	2045-08-01	2072-08-01

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

Compare with status at start of the Paris Agreement

Loading the data from the start of the Paris Agreement.

In [49]:

Code

beforePApaths = CSV.read("tables/ProbabilityPathsHadCRUT5BeforePA-Extensive.csv", DataFrame);
afterPApaths = CSV.read("tables/ProbabilityPathsHadCRUT5-Extensive.csv", DataFrame);

In [50]:

Code

xls = (Date(2015, 11, 1), Date(2090, 1, 1))

ppa = plot(afterPApaths[!, 1], afterPApaths[!, 2], label="1.5C", color=:darkorange, legend=:bottomright, xticks=(collect(xls[1]:Dates.Month(120):xls[2]), Dates.format.(collect(xls[1]:Dates.Month(120):xls[2]), "Y")), linewidth=4, linestyle=:dashdot, title="", xlabel="Date (monthly)", ylabel="Probability")
plot!(beforePApaths[!, 1], beforePApaths[!, 2], label="1.5C (At PA start)", color=:orange, linewidth=4, linestyle=:dot)
plot!(afterPApaths[!, 1], afterPApaths[!, 3], label="2.0C", color=:red3, linewidth=4, linestyle=:dash)
plot!(beforePApaths[!, 1], beforePApaths[!, 3], label="2.0C (At PA start)", color=:red2, linewidth=4, linestyle=:solid)
plot!(fontfamily="Computer Modern", legendfontsize=12, tickfontsize=12, titlefontfamily="Computer Modern", legendfontfamily="Computer Modern", tickfontfamily="Computer Modern", ylabelfontsize=12, xlabelfontsize=12, titlefontsize=14, ylims=(0, 1), xlims=xls)

Figure 9: Proportion of scenarios that breach the 1.5°C and 2°C thresholds for the HadCRUT5 temperature anomalies for each month at the start of the Paris Agreement and at the last observation. 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 [51]:

Code

savefig(ppa,"figures/Coverage-HadCRUT5-PrePostPA-Extensive.png")

"/Users/jeddy/Library/CloudStorage/OneDrive-AalborgUniversitet/Research/CLIMATE/Paris Goal/Odds-of-breaching-1.5C/figures/Coverage-HadCRUT5-PrePostPA-Extensive.png"

Morice, Colin P, John J Kennedy, Nick A Rayner, JP Winn, Emma Hogan, RE Killick, RJH Dunn, TJ Osborn, PD Jones, and IR Simpson. 2021. “An Updated Assessment of Near-Surface Temperature Change from 1850: The HadCRUT5 Data Set.” Journal of Geophysical Research: Atmospheres 126 (3): e2019JD032361.

Replication notebook

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 Plotting baselines on mean ensemble temperature

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 trend model

Computing the breach point

5. Forecasting paths

Plotting the forecasts

Coverage of the predictions

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

Data frame to save probability paths

Data frame to store the probabilities

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

Compare with status at start of the Paris Agreement