Back to Article
Supplementary notebook (probabilities at the start of the Paris Agreement)
Download Notebook

Breaching 1.5°C: Give me the odds

Probability at the Paris Agreement notebook

Authors
Affiliation

J. Eduardo Vera-Valdés

Aalborg University

Olivia Kvist

Aalborg University

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 [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/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 [3]:
Code 
menstemp = reduce(+, eachcol(rawtemp[:, 4:203])) ./ ncol(rawtemp[:, 4:203]);
oldbase = mean(menstemp[(rawtemp.Time.>Date(1850, 1, 1)).&(rawtemp.Time.<Date(1900, 1, 1))]);
menstempind = menstemp .- oldbase;

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 [4]:
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 [5]:
Code 
date_nino = Date.(rawnino[!, 1], rawnino[!, 2])
tempnino = rawtemp[rawtemp.Time.>=date_nino[1], :]
tempnino[!, :ONI] = rawnino.ONI
rename!(tempnino, :Time => :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 [6]:
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 Data at the Paris Agreement

In [7]:
Code 
delete!(tempnino, tempnino.Date .> Date(2016, 12, 1));
T = length(tempnino.Date)
last(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 2016-08-01 0.991773 0.0172583 1.33266 1.25319 1.31662 1.30374 1.27611 1.32196 1.24759 1.36038 1.28391 1.40771 1.33004 1.22629 1.3325 1.30346 1.28883 1.28635 1.22357 1.36391 1.3142 1.26113 1.30701 1.39831 1.26307 1.35769 1.31247 1.2968 1.28484 1.33995 1.35125 1.31621 1.2541 1.31973 1.18768 1.24703 1.22309 1.37071 1.37493 1.34399 1.31083 1.25889 1.30028 1.35589 1.29191 1.28582 1.32764 1.22311 1.31315 1.33609 1.34624 1.33891 1.29256 1.29399 1.36801 1.32013 1.28255 1.21898 1.31823 1.34048 1.32233 1.28525 1.28646 1.34888 1.28476 1.23207 1.23326 1.35448 1.27182 1.2794 1.31055 1.25967 1.30247 1.2766 1.22596 1.36308 1.38134 1.29603 1.31084 1.34038 1.25577 1.30301 1.37435 1.24523 1.32459 1.28246 1.3187 1.33859 1.28796 1.23809 1.32422 1.32759 1.33448 1.31919 1.31166 1.37471 1.19676 1.25841 1.33656 ⋯
2 2016-09-01 0.992085 0.0125059 1.21306 1.13789 1.26487 1.24345 1.15731 1.24097 1.11536 1.29381 1.18664 1.28091 1.21797 1.11733 1.24952 1.18648 1.17698 1.14906 1.09883 1.25227 1.23092 1.15971 1.23811 1.29972 1.18045 1.2491 1.20763 1.17198 1.18514 1.25289 1.24516 1.21454 1.14067 1.22136 1.12567 1.151 1.14925 1.25519 1.25474 1.22822 1.20195 1.17245 1.20604 1.22517 1.17172 1.17761 1.23029 1.08533 1.20587 1.23285 1.24811 1.28831 1.16787 1.20784 1.26837 1.19041 1.19012 1.10914 1.23096 1.24595 1.20521 1.16566 1.21079 1.25704 1.19812 1.14436 1.15407 1.26781 1.17455 1.19429 1.19571 1.16783 1.22572 1.15929 1.14474 1.26758 1.26863 1.17672 1.23181 1.2409 1.12121 1.25339 1.20652 1.14071 1.21226 1.15486 1.2433 1.25263 1.16881 1.13893 1.23384 1.21181 1.24893 1.18184 1.19629 1.23976 1.10918 1.1791 1.24134 ⋯
3 2016-10-01 0.99364 0.0107343 1.19529 1.12868 1.20483 1.18544 1.11991 1.20642 1.14693 1.22894 1.11646 1.2732 1.20371 1.09856 1.22016 1.12659 1.14699 1.11759 1.06321 1.24578 1.20677 1.11234 1.20498 1.25947 1.1766 1.24593 1.21215 1.13225 1.17366 1.19871 1.25335 1.19569 1.1371 1.21732 1.08797 1.1259 1.13324 1.19519 1.22635 1.2272 1.1604 1.11121 1.1308 1.20722 1.13884 1.15458 1.19521 1.09225 1.18102 1.17486 1.22524 1.2218 1.15362 1.20275 1.26361 1.18742 1.14036 1.10379 1.22984 1.20985 1.21088 1.15484 1.16388 1.19329 1.15979 1.13084 1.11535 1.24272 1.12883 1.15281 1.15187 1.08469 1.19006 1.14408 1.14172 1.24383 1.24355 1.16757 1.19536 1.16816 1.12524 1.21414 1.17901 1.13037 1.21301 1.14989 1.20304 1.27112 1.16408 1.11688 1.2004 1.21578 1.17412 1.18425 1.17971 1.19463 1.11517 1.14761 1.19779 ⋯
4 2016-11-01 0.995774 0.00491544 1.23774 1.12936 1.25456 1.21385 1.17515 1.23793 1.15418 1.2455 1.14366 1.29669 1.2319 1.11695 1.221 1.17502 1.19923 1.1537 1.10613 1.25378 1.19249 1.16188 1.2483 1.30677 1.19418 1.27915 1.21797 1.16546 1.19637 1.23704 1.26138 1.23693 1.15549 1.23564 1.1495 1.15155 1.13994 1.22143 1.24159 1.21302 1.20162 1.15428 1.13786 1.25148 1.2175 1.20229 1.24933 1.11434 1.23311 1.2397 1.23588 1.26193 1.18479 1.2244 1.26599 1.21858 1.20018 1.12312 1.21589 1.20885 1.22464 1.16406 1.20789 1.23477 1.17979 1.15699 1.16252 1.27527 1.20841 1.15886 1.17992 1.15686 1.21292 1.1543 1.14247 1.26053 1.27224 1.19572 1.21854 1.21689 1.12568 1.2288 1.21574 1.1264 1.22056 1.17049 1.2359 1.24885 1.15435 1.1331 1.20594 1.20456 1.20625 1.22382 1.20115 1.20104 1.13389 1.15076 1.23062 ⋯
5 2016-12-01 0.998609 0.00178567 1.1571 1.08115 1.16442 1.17574 1.09624 1.16194 1.07736 1.20042 1.11312 1.22271 1.1443 1.04529 1.1617 1.12217 1.13677 1.09661 1.03635 1.19311 1.12665 1.11091 1.16238 1.19203 1.10492 1.20448 1.17843 1.09184 1.11814 1.10348 1.22308 1.15174 1.08113 1.14337 1.0527 1.07576 1.08675 1.19955 1.19554 1.14453 1.15584 1.09745 1.1124 1.16406 1.13651 1.10736 1.15856 1.10486 1.13542 1.16655 1.17029 1.20719 1.13051 1.14507 1.21617 1.1296 1.12877 1.07732 1.15468 1.17052 1.16429 1.1012 1.16071 1.20717 1.13191 1.08014 1.1017 1.19977 1.10354 1.11125 1.09564 1.09134 1.14638 1.08465 1.06844 1.17754 1.18863 1.13539 1.1546 1.14744 1.0759 1.17532 1.14751 1.07481 1.16183 1.11073 1.1581 1.19305 1.09472 1.08359 1.177 1.17465 1.13701 1.13232 1.13201 1.15723 1.07483 1.12847 1.1568 ⋯

2. 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 [8]:
Code 
nino_model3 = MSModel(tempnino[!, :ONI], 3);
summary_msm(nino_model3);
Markov Switching Model with 3 regimes
=================================================================
# of observations:         1751 AIC:                      2342.501
# of estimated parameters:   12 BIC:                      2408.116
Error distribution:    Gaussian Instant. adj. R^2:           0.709
Loglikelihood:          -1159.3 Step-ahead adj. R^2:        0.5987
-----------------------------------------------------------------
------------------------------
Summary of regime 1: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |      1.21  |       0.043  |   28.003  |    < 1e-3  
σ            |     0.525  |       0.014  |   37.692  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 9.98 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 2: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |    -0.745  |       0.033  |  -22.697  |    < 1e-3  
σ            |     0.417  |        0.01  |   41.138  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 15.94 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 3: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |     0.181  |       0.039  |    4.696  |    < 1e-3  
σ            |     0.271  |       0.014  |   19.665  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 7.68 periods
-------------------------------------------------------------------
left-stochastic transition matrix: 
          | regime 1   | regime 2   | regime 3
----------------------------------------------------
 regime 1 |   89.985%  |      0.0%  |     5.87%  |
 regime 2 |      0.0%  |   93.726%  |    7.146%  |
 regime 3 |   10.015%  |    6.274%  |   86.985%  |
In [9]:
Code 
nino_model5 = MSModel(tempnino[!, :ONI], 5);
summary_msm(nino_model5);
Markov Switching Model with 5 regimes
=================================================================
# of observations:         1751 AIC:                      2232.031
# of estimated parameters:   30 BIC:                       2396.07
Error distribution:    Gaussian Instant. adj. R^2:           0.627
Loglikelihood:          -1086.0 Step-ahead adj. R^2:        0.5488
-----------------------------------------------------------------
------------------------------
Summary of regime 1: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |    -0.092  |       0.019  |   -4.982  |    < 1e-3  
σ            |     0.117  |       0.017  |    6.668  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 1.50 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 2: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |     0.317  |       0.012  |   25.609  |    < 1e-3  
σ            |     0.191  |       0.012  |   15.939  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 4.87 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 3: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |    -0.616  |        0.04  |  -15.395  |    < 1e-3  
σ            |     0.319  |       0.007  |   47.125  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 11.30 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 4: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |    -0.478  |       0.064  |   -7.466  |    < 1e-3  
σ            |     1.547  |       0.021  |   73.803  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 5.35 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 5: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |     1.139  |        0.03  |   37.999  |    < 1e-3  
σ            |     0.598  |       0.012  |   47.915  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 12.24 periods
-------------------------------------------------------------------
left-stochastic transition matrix: 
          | regime 1   | regime 2   | regime 3   | regime 4   | regime 5
------------------------------------------------------------------------------
 regime 1 |    33.49%  |    5.387%  |    1.046%  |      0.0%  |      0.0%  |
 regime 2 |   25.878%  |   79.466%  |    0.446%  |      0.0%  |    8.171%  |
 regime 3 |   33.518%  |    1.536%  |   91.149%  |   18.706%  |      0.0%  |
 regime 4 |      0.0%  |      0.0%  |    2.966%  |   81.294%  |      0.0%  |
 regime 5 |    7.114%  |   13.611%  |    4.393%  |      0.0%  |   91.829%  |
In [10]:
Code 
nino_model7 = MSModel(tempnino[!, :ONI], 7);
summary_msm(nino_model7);
Markov Switching Model with 7 regimes
=================================================================
# of observations:         1751 AIC:                      1119.398
# of estimated parameters:   56 BIC:                      1425.603
Error distribution:    Gaussian Instant. adj. R^2:          0.9053
Loglikelihood:           -503.7 Step-ahead adj. R^2:        0.7855
-----------------------------------------------------------------
------------------------------
Summary of regime 1: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |     0.543  |       0.011  |   48.849  |    < 1e-3  
σ            |     0.132  |       0.006  |   22.206  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 3.51 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 2: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |    -1.183  |       0.028  |  -42.172  |    < 1e-3  
σ            |     0.323  |       0.011  |    29.35  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 8.91 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 3: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |    -0.466  |       0.015  |  -30.539  |    < 1e-3  
σ            |      0.22  |       0.005  |   46.246  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 5.93 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 4: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |     0.185  |       0.017  |   11.225  |    < 1e-3  
σ            |     0.129  |       0.017  |    7.379  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 5.74 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 5: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |     1.875  |       0.045  |   41.325  |    < 1e-3  
σ            |     0.448  |       0.019  |   23.053  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 5.63 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 6: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |     1.057  |       0.024  |    44.26  |    < 1e-3  
σ            |     0.226  |       0.009  |   24.771  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 4.11 periods
-------------------------------------------------------------------
------------------------------
Summary of regime 7: 
------------------------------
Coefficient  |  Estimate  |  Std. Error  |  z value  |  Pr(>|z|) 
-------------------------------------------------------------------
β_0          |    -0.012  |        0.04  |   -0.298  |     0.765  
σ            |       0.1  |       0.021  |    4.874  |    < 1e-3  
-------------------------------------------------------------------
Expected regime duration: 2.07 periods
-------------------------------------------------------------------
left-stochastic transition matrix: 
          | regime 1   | regime 2   | regime 3   | regime 4   | regime 5   | regime 6   | regime 7
--------------------------------------------------------------------------------------------------------
 regime 1 |   71.534%  |      0.0%  |      0.0%  |    6.923%  |      0.0%  |   10.078%  |      0.0%  |
 regime 2 |      0.0%  |   88.772%  |    5.664%  |      0.0%  |      0.0%  |      0.0%  |      0.0%  |
 regime 3 |      0.0%  |   11.228%  |   83.131%  |      0.0%  |      0.0%  |      0.0%  |   44.851%  |
 regime 4 |   11.859%  |      0.0%  |    9.635%  |   82.574%  |      0.0%  |    0.432%  |    3.519%  |
 regime 5 |      0.0%  |      0.0%  |      0.0%  |      0.0%  |   82.245%  |    7.776%  |      0.0%  |
 regime 6 |   13.577%  |      0.0%  |      0.0%  |      0.0%  |   12.086%  |   75.675%  |      0.0%  |
 regime 7 |    3.029%  |      0.0%  |     1.57%  |   10.503%  |    5.669%  |    6.039%  |    51.63%  |

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

3. 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 [11]:
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 [12]:
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 [13]:
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 [14]:
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 [15]:
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)

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.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 [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)
In [18]:
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 of temperature anomalies for realization 10")
plot!(tempnino.Date, bmodel.Yfit, linestyle=:dot, linewidth=2, label="Trend", color=2)
plot!(tempnino.Date, bmodel_exo.Yfit, linestyle=:dash, linewidth=2, label="Trend + El Niño", color=3)
plot!(date_for, bmodel_forecast.Yforecastmean, linestyle=:dot, linewidth=2, label="Trend", color=4)
plot!(date_for, bmodel_forecast.Yforecasterr, linestyle=:dash, linewidth=2, label="Trend + long-range dependence", color=5)
plot!(date_for, bmodel_exo_forecast.Yforecasterr, linestyle=:dashdot, linewidth=2, label="Trend + long-range dependence + 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 1: 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.

Selecting the best trend model

AIC

In [19]:
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.10016852003856833, 0.00031332907360425246, 0.0012464691349785165, 0.07936505677116559], σ² = 0.027850629049125584, break_point = 1267, rss = 48.6550489488224, aic = -6266.159995701619, bic = -6244.288226372475, betavar = [8.201209142226718e-5 -8.986437084944063e-8 1.884892591682906e-7 -1.0085878242786883e-6; -8.986437084944061e-8 1.3039803301218045e-10 -3.637020597852166e-10 7.49661044729996e-10; 1.884892591682905e-7 -3.637020597852166e-10 1.9414801862852096e-9 -1.2303692379564436e-9; -1.008587824278688e-6 7.496610447299961e-10 -1.2303692379564436e-9 2.2089998057643436e-5], res = [-0.14604616423706884, 0.3600929573871306, 0.26516412525615823, 0.0449171133439957, 0.15637112994750801, 0.3849515761967872, 0.1728424687429314, 0.0435968079661922, 0.13065666570512785, 0.03478152628262174  …  0.32824372447647576, 0.21050956260143905, 0.10083484561699718, 0.022383273523150837, 0.061242481778206015, 0.24412672378843325, 0.11748075068925523, 0.15225477588694247, 0.1563632165429364, 0.0700434126572369], Yfit = [-0.12049010572546713, -0.09319265734966657, -0.07462536521869423, -0.07034378330653171, -0.07796695991004401, -0.0697171251593232, -0.059086338705467416, -0.0563920579287282, -0.06560253566766382, -0.07878126624515772  …  1.1639140055609882, 1.114680167436025, 1.0725891844204667, 1.037641056514313, 1.016184988259258, 1.0034590762490307, 0.9978760193482088, 0.9946739141505215, 0.9978210134945277, 1.007317317380227]))

BIC

In [20]:
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.10016852003856833, 0.00031332907360425246, 0.0012464691349785165, 0.07936505677116559], σ² = 0.027850629049125584, break_point = 1267, rss = 48.6550489488224, aic = -6266.159995701619, bic = -6244.288226372475, betavar = [8.201209142226718e-5 -8.986437084944063e-8 1.884892591682906e-7 -1.0085878242786883e-6; -8.986437084944061e-8 1.3039803301218045e-10 -3.637020597852166e-10 7.49661044729996e-10; 1.884892591682905e-7 -3.637020597852166e-10 1.9414801862852096e-9 -1.2303692379564436e-9; -1.008587824278688e-6 7.496610447299961e-10 -1.2303692379564436e-9 2.2089998057643436e-5], res = [-0.14604616423706884, 0.3600929573871306, 0.26516412525615823, 0.0449171133439957, 0.15637112994750801, 0.3849515761967872, 0.1728424687429314, 0.0435968079661922, 0.13065666570512785, 0.03478152628262174  …  0.32824372447647576, 0.21050956260143905, 0.10083484561699718, 0.022383273523150837, 0.061242481778206015, 0.24412672378843325, 0.11748075068925523, 0.15225477588694247, 0.1563632165429364, 0.0700434126572369], Yfit = [-0.12049010572546713, -0.09319265734966657, -0.07462536521869423, -0.07034378330653171, -0.07796695991004401, -0.0697171251593232, -0.059086338705467416, -0.0563920579287282, -0.06560253566766382, -0.07878126624515772  …  1.1639140055609882, 1.114680167436025, 1.0725891844204667, 1.037641056514313, 1.016184988259258, 1.0034590762490307, 0.9978760193482088, 0.9946739141505215, 0.9978210134945277, 1.007317317380227]))

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

In [21]:
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}:
 -5607.83  -6231.08  -6266.16
 -5591.42  -6209.21  -6244.29

Computing the breach point

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

anim = @animate for ii = 30:24:360
    plot(fulldate[1:T], thisreal[1:T], linewidth=1, label="Temperature anomalies", xlabel="Date (monthly)", ylabel="°C", legend=:bottomright, title="Breaching of the 1.5°C threshold for realization 10")
    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/breach_realization10_prePA.gif", fps = 5)
[ Info: Saved animation to /Users/jeddy/Library/CloudStorage/OneDrive-AalborgUniversitet/Research/CLIMATE/Paris Goal/Odds-of-breaching-1.5C/figures/breach_realization10_prePA.gif
Figure 2: 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 [23]:
Code 
breachreal10 = rolling_sample_mean(thisreal)
fulldate[breachreal10]
2040-05-01
In [24]:
Code 
ps = plot(fulldate[1:T], thisreal[1:T], linewidth=1, label="Temperature anomalies", xlabel="Date (monthly)", ylabel="°C", legend=:bottomright, title="Breaching of the 1.5°C threshold for realization 10")
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 3: 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.

4. 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 [25]:
Code 
nsim = 5
matforecasts = zeros(h, 200 * nsim)
uncertainty = 1.96

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)

    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

Plotting the forecasts

In [26]:
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.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 4: 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.

5. Coverage of the prediction intervals

We compute the coverage of the prediction intervals for the predictions of the temperature anomalies at the Paris Agreement.

In [27]:
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

IPCC estimates

In [28]:
Code 
ipcc = CSV.read("data/ipcc-projections-sr15.csv", DataFrame)
first(ipcc, 5)
5×3 DataFrame
Row year min_t max_t
Int64 Float64 Float64
1 1850 -0.023 -0.035
2 1851 -0.022 -0.033
3 1852 -0.021 -0.032
4 1853 -0.02 -0.03
5 1854 -0.019 -0.029 
In [29]:
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, 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)
vline!([tempnino.Date[end]], label="Paris Agreement in effect", color=:gray, linewidth=2, linestyle=:dash)
plot!(Date.(ipcc.year), [ipcc.min_t ipcc.max_t], linewidth=2, linestyle=[:dashdot :dash], label=["IPCC minimum" "IPCC maximum"], xlabel="Date", ylabel="°C", legend=:topleft, color = [6 7])
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 5: Simulated forecast paths for HadCRUT5 temperature anomalies. The 95% and 99% prediction intervals are shown as shaded areas. The IPCC projections for the minimum and maximum temperature anomalies are shown as dashed lines (Allen et al. 2018). 
In [30]:
Code 
savefig("figures/Paths-HadCRUT5-PrePA-Coverage.png")
"/Users/jeddy/Library/CloudStorage/OneDrive-AalborgUniversitet/Research/CLIMATE/Paris Goal/Odds-of-breaching-1.5C/figures/Paths-HadCRUT5-PrePA-Coverage.png"

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 [31]:
Code 
inicio = T - 119 # 10 years starting 2004
fin = T + h - 120 # 10 years

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

for jj = 1:(200*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 [32]:
Code 
plot(datejoin, pa15, label="1.5C", color=:darkorange, legend=:bottomright, xlims=(datejoin[120], datejoin[end-48]), xticks=(datejoin[120:120:end], Dates.format.(datejoin[120:120:end], "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=(Date(2024, 1, 1), Date(2070, 1, 1)))
Figure 6: 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. 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 [33]:
Code 
savefig("figures/Coverage-HadCRUT5-PrePA.png")
"/Users/jeddy/Library/CloudStorage/OneDrive-AalborgUniversitet/Research/CLIMATE/Paris Goal/Odds-of-breaching-1.5C/figures/Coverage-HadCRUT5-PrePA.png"

Determning the first month where the probability of exceeding 1.5°C and 2°C is greater than 0%.

In [34]:
Code 
display(datejoin[findfirst(pa15 .> 0)])
display(datejoin[findfirst(pa20 .> 0)])
2022-07-01
2042-09-01

Determining the first month that the probability of exceeding 1.5°C and 2°C is greater than 99%.

In [35]:
Code 
display(datejoin[findfirst(pa15 .>= 0.99)])
2049-09-01
In [36]:
Code 
date_for[end]
2083-08-01
Allen, Myles, Opha Pauline Dube, William Solecki, Fernando Aragón-Durand, Wolfgang Cramer, Stephen Humphreys, Mikiko Kainuma, et al. 2018. “Special Report: Global Warming of 1.5 c.” Intergovernmental Panel on Climate Change (IPCC) 27: 677.