In this post we are going to play with pyextremes a python library that does everything (and more) I have been talking about the two previous posts. After a bit of EDA, we are going to fit use BM and POT to compute the return values – could be interesting to know what magnitude to expect at a given horizon, no?

\(m\)-years return values/levels are the “high quantile for which the probability that the annual maximum exceeds this quantile is \(1/m\)” – see 100-year flood (wiki).

Setup and data - Earthquakes in CA

Setup

A bit of setup first…

!pip3.8 install -q -U emcee folium pyextremes statsmodels

import pandas as pd
import numpy as np

%matplotlib inline
import matplotlib as mpl
from matplotlib import pylab as plt
import matplotlib.dates as mdates
import seaborn as sns
sns.set(rc={'figure.figsize':(18,6)})
sns.set_style("ticks")

cmap = plt.get_cmap('plasma')

import statsmodels.api as sm

from pyextremes import EVA, __version__
import pyextremes

from pyextremes.plotting import (
    plot_extremes,
    plot_trace,
    plot_corner,
)

vmin=3.0
vmax=9.0

print("pyextremes", __version__)

pyextremes 2.2.1

Data

Data are from the Northern California Earthquake Data Center (NCEDC).

• Waldhauser, F. and D.P. Schaff, Large-scale relocation of two decades of Northern California seismicity using cross-correlation and double-difference methods, J. Geophys. Res.,113, B08311, doi:10.1029/2007JB005479, 2008.
• Waldhauser, F., Near-real-time double-difference event location using long-term seismic archives, with application to Northern California, Bull. Seism. Soc. Am., 99, 2736-2848, doi:10.1785/0120080294, 2009.

Go here if you want to download them.

df_ = pd.read_csv('ncedc.csv')
df_

magTypeColors = sns.color_palette("hls", 8)
palette = {
    'ML': magTypeColors[0],
    'Md': magTypeColors[1],
    'Mh': magTypeColors[2],
    'Mw': magTypeColors[3],
    'Mx': magTypeColors[4],
    'Unk': magTypeColors[5],
}
sns.histplot(data=df_.sort_values(by='magType'),
             x='magnitude',
             hue='magType',
             bins=20,
             multiple='stack',
             palette=palette)
sns.despine(left=True)
plt.show()

And more specifically above 5.

sns.histplot(data=df_[df_['magnitude'] > 5].sort_values(by='magType'),
             x='magnitude',
             hue='magType',
             bins=20,
             multiple='stack',
             palette=palette)
sns.despine(left=True)
plt.show()

The magniture are based on different scales.

df_.groupby(by='magType').size().to_frame('count')

We will stick to only one scale as mapping on to each other seems to be difficult. They basically measure different things… – See USGS - what-are-different-magnitude-scales.

We pick moment magnitude (Mw) scale as “[it is] based on the concept of seismic moment, is uniformly applicable to all sizes of earthquakes but is more difficult to compute than the other types” – See magnitude. Other scales are not covering the entire range of eathquakes creating truncated data. Note that this could introduce a selection bias as we are not sure why a particular events is using one scale or another.

df__ = df_[(df_['magType'] == 'Mw')]
df__.groupby(by='magType').size()

Mw    1456

df_all = df__[['time', 'latitude', 'longitude', 'magnitude']].copy()
df_all.rename({'magnitude': 'value',
               'time': 'obs'}, inplace=True, axis=1)
df_all = df_all.dropna(axis=0)
df_all = df_all.set_index(pd.DatetimeIndex(df_all['obs'])).drop('obs', axis=1)
df_all.head()

sns.histplot(data=df_all, x='value', bins=20, color=palette['Mw'])
sns.despine(left=True)
plt.show()

The data we are looking at cover the following time interval:

pd.DataFrame.from_dict(data={'min': [min(df_all.index)],
                             'max': [max(df_all.index)]},
                       orient='index',
                       columns=['obs'])

and range as follows:

pd.DataFrame(df_all['value'].describe())

The biggest earthquake in CA (in this dataset and using this scale…) are:

df_all.sort_values('value', ascending=False) \
      .nlargest(columns=['value'], n=10) \
      .style.background_gradient(subset=['value'],
                                 cmap=cmap,
                                 vmin=vmin,
                                 vmax=vmax)

import folium

def transparent(r, g, b, a, alpha=None):
    if alpha is None:
        return (r, g, b, a)
    return (r, g, b, alpha)
def rgb2hex(r, g, b, a=None):
    if a is None:
        return "#{:02x}{:02x}{:02x}".format(int(r*255),
                                            int(g*255),
                                            int(b*255))
    else:
        return "#{:02x}{:02x}{:02x}{:02x}".format(int(r*255),
                                                  int(g*255),
                                                  int(b*255),
                                                  int(a*255))

m = folium.Map(location=[df_all['latitude'].mean(),
                         df_all['longitude'].mean()],
                         tiles="Stamen Terrain", zoom_start=6)

def add_to_map(row):
    folium.Circle(
        radius=100*row['value'],
        location=[row['latitude'], row['longitude']],
        popup="{obs} {value}".format(obs=row.name,
                                     value=row['value']),
        color=rgb2hex(*transparent(*cmap((row['value']-vmin)/(vmax - vmin)),
                                   alpha=0.5)),
        fill=True,
    ).add_to(m)

df_all.apply(add_to_map, axis=1)

m

from matplotlib.dates import YearLocator, DateFormatter

sns.scatterplot(data=df_all, x='obs', y='value',
                marker='.', palette=cmap, hue='value',
                hue_norm=(vmin, vmax), s=150, alpha=0.9)
ax = plt.gca()
ax.xaxis.set_major_locator(YearLocator(2))
ax.xaxis.set_major_formatter(DateFormatter('%Y'))
sns.despine(left=True)
plt.show()

The scale (Mw) seems to have been used throughout the period we are looking at.

data_df = df_all.groupby(by='obs')['value'].max()
data_df.index = pd.to_datetime(data_df.index)

How independent are these events?

sm.graphics.tsa.plot_acf(data_df, lags=20)
sns.despine(left=True)
plt.show()

There is a bit of autocorrelation. Not completely surprising as “Aftershocks can continue over a period of weeks, months, or years.” USGS glossary. To remove that we are going look at the max in 1-week windows. It is not perfect as aftershocks can span longer than that and there is no guarantee that all the aftershocks are going to fall in the same buckets. Let see how it does.

data_df_ = data_df.resample('1W').max().dropna()
sm.graphics.tsa.plot_acf(data_df_.squeeze(), lags=20)
sns.despine(left=True)
plt.show()

Looking better.

sns.scatterplot(data=data_df_.to_frame('value'), x='obs', y='value', marker='.', palette=cmap, hue='value', hue_norm=(vmin, vmax), s=150, alpha=0.9)
ax = plt.gca()
ax.xaxis.set_major_locator(YearLocator(2))
ax.xaxis.set_major_formatter(DateFormatter('%Y'))
sns.despine(left=True)
plt.show()

Let’s put the data in the form pyextremes expect them.

data=data_df_.to_frame('value').squeeze()
data

    obs
    1992-04-26    6.57
    1992-09-20    4.47
    1993-01-17    4.91
    1993-01-24    3.92
    1993-02-21    3.88
                  ...
    2020-12-06    4.96
    2020-12-13    4.77
    2020-12-20    4.81
    2020-12-27    3.86
    2021-01-03    3.63
    Name: value, Length: 777, dtype: float64

BM model

model = EVA(data=data)
model

                           Univariate Extreme Value Analysis
===========================================================================
                                      Source Data
---------------------------------------------------------------------------
Data label:                         value     Size:                     777
Start:                         April 1992     End:             January 2021
===========================================================================
                                     Extreme Values
---------------------------------------------------------------------------
Extreme values have not been extracted
===========================================================================
                                         Model
---------------------------------------------------------------------------
Model has not been fit to the extremes
===========================================================================

model.get_extremes(
    method="BM",
    extremes_type="high",
    block_size="365.2425D",
    errors="ignore",
)
model

                           Univariate Extreme Value Analysis
===========================================================================
                                      Source Data
---------------------------------------------------------------------------
Data label:                         value     Size:                     777
Start:                         April 1992     End:             January 2021
===========================================================================
                                     Extreme Values
---------------------------------------------------------------------------
Count:                                 29     Extraction method:         BM
Type:                                high     Block size: 365 days 05:49:12
===========================================================================
                                         Model
---------------------------------------------------------------------------
Model has not been fit to the extremes
===========================================================================

model.plot_extremes()
sns.despine(left=True)
plt.show()

MLE

model.fit_model()
model

                           Univariate Extreme Value Analysis
===========================================================================
                         Source Data
---------------------------------------------------------------------------
Data label:            value      Size:                                 777
Start:            April 1992      End:                         January 2021
===========================================================================
                        Extreme Values
---------------------------------------------------------------------------
Count:                    29      Extraction method:                     BM
Type:                   high      Block size:             365 days 05:49:12
===========================================================================
                            Model
---------------------------------------------------------------------------
Model:                   MLE      Distribution:                    gumbel_r
Log-likelihood:      -28.048      AIC:                               60.558
---------------------------------------------------------------------------
Free parameters:   loc=5.340      Fixed parameters: All parameters are free
                 scale=0.550
===========================================================================

model.plot_diagnostic(alpha=0.95)
sns.despine(left=True)
plt.show()

model.plot_return_values(
    return_period=np.logspace(0.01, 2, 100),
    return_period_size="365.2425D",
    alpha=0.95,
)
sns.despine(left=True)
plt.show()

summary_bm_mle = model.get_summary(
    return_period=[2, 5, 10, 25, 50, 100],
    alpha=0.95,
    n_samples=1000,
)
summary_bm_mle.style.background_gradient(cmap=cmap, vmin=vmin, vmax=vmax)

MCMC

The library also supports fitting the distribution using a Bayesian approach. Wont go over this here. Check the tutorials that come with it.

POT

POT requires to select a threshold. A method to do that is to apply the fit the GPD for a range of thresholds and look at the stability of its parameters. This is what the following does:

pyextremes.plot_parameter_stability(ts=data, thresholds=np.linspace(4, 6, 50))
sns.despine(left=True)
plt.show()

We then pick a threshold in the range where the scale and shape parameters of the GPD are stable.

threshold = 4.25

model = EVA(data=data)
model

                           Univariate Extreme Value Analysis
===========================================================================
                                      Source Data
---------------------------------------------------------------------------
Data label:                         value      Size:                    777
Start:                         April 1992      End:            January 2021
===========================================================================
                                     Extreme Values
---------------------------------------------------------------------------
Extreme values have not been extracted
===========================================================================
                                         Model
---------------------------------------------------------------------------
Model has not been fit to the extremes
===========================================================================

sns.histplot(data=data[data>threshold], bins=20)
sns.despine(left=True)
plt.show()

model.get_extremes(
    method="POT",
    extremes_type="high",
    threshold=threshold)
model

                           Univariate Extreme Value Analysis
===========================================================================
                                  Source Data
---------------------------------------------------------------------------
Data label:                     value      Size:                       777
Start:                     April 1992      End:               January 2021
===========================================================================
                                 Extreme Values
---------------------------------------------------------------------------
Count:                            237      Extraction method:          POT
Type:                            high      Threshold:                 4.25
===========================================================================
                                     Model
---------------------------------------------------------------------------
Model has not been fit to the extremes
===========================================================================

model.plot_extremes()
sns.despine(left=True)
plt.show()

model.fit_model()
model

                           Univariate Extreme Value Analysis
===========================================================================
                                  Source Data
---------------------------------------------------------------------------
Data label:                     value      Size:                        777
Start:                     April 1992      End:                January 2021
===========================================================================
                                 Extreme Values
---------------------------------------------------------------------------
Count:                            237      Extraction method:           POT
Type:                            high      Threshold:                  4.25
===========================================================================
                                     Model
---------------------------------------------------------------------------
Model:                            MLE      Distribution:              expon
Log-likelihood:               -77.084      AIC:                     156.185
---------------------------------------------------------------------------
Free parameters:          scale=0.509      Fixed parameters:     floc=4.250
===========================================================================

model.plot_diagnostic(alpha=0.95)
sns.despine(left=True)
plt.show()

model.plot_return_values(
    return_period=np.logspace(0.01, 2, 100),
    return_period_size="365.2425D",
    alpha=0.95,
)
sns.despine(left=True)
plt.show()

summary_pot_mle = model.get_summary(
    return_period=[2, 5, 10, 25, 50, 100],
    alpha=0.95,
    n_samples=1000,
)
summary_pot_mle.style.background_gradient(cmap=cmap, vmin=vmin, vmax=vmax)

Comparison

Let’s compare what both methods.

df_bm = summary_bm_mle
df_pot = summary_pot_mle
# bm
plt.errorbar(df_bm.index,
             df_bm['return value'],
             yerr=[df_bm['return value'] - df_bm['lower ci'],
                   df_bm['upper ci'] - df_bm['return value']],
             fmt='-ob', label='return value (bm)');
plt.fill_between(df_bm.index,
                 df_bm['lower ci'],
                 df_bm['upper ci'],
                 color='b',
                 alpha=0.2)
# pot
plt.errorbar(df_pot.index,
             df_pot['return value'],
             yerr=[df_pot['return value'] - df_pot['lower ci'],
                   df_pot['upper ci'] - df_pot['return value']],
             fmt='-or',
             label='return value (pot)');
plt.fill_between(df_pot.index,
                 df_pot['lower ci'],
                 df_pot['upper ci'],
                 color='r',
                 alpha=0.2)
plt.grid(True, which="both", linestyle='--', linewidth=1)
plt.gca().set_xscale('log')
plt.title('global earthquakes - return values - bm and pot 95% CIs')
plt.legend()
sns.despine(left=True)
plt.show()

pd.concat([summary_bm_mle, summary_pot_mle],
          axis=1,
          keys=['bm', 'pot']).style.background_gradient(cmap=cmap,
                                                        vmin=vmin,
                                                        vmax=vmax)

Lets see what we had between 1992 and 2021:

df_all.sort_values('value', ascending=False) \
      .nlargest(columns=['value'], n=10) \
      .style.background_gradient(subset=['value'],
                                 cmap=cmap,
                                 vmin=vmin,
                                 vmax=vmax)

Quakes greater than 6.5 and 7.2 (Mw) have definitively been observed between 1992 and 2021 - looking at CI for the return value at 20 years - remember it doesn’t have to by the definition of the return value.

Conclusion

Cool piece of theory too look at! Nice lib! Nice tool!

References

• An Introduction to Statistical Modeling of Extremes Values, Stuart Coles
• pyextremes - python library