change_points.py

#!pip install ruptures
import ruptures as rpt
from ruptures.exceptions import BadSegmentationParameters

def change_points(data, model='l1', pen=3):
    signal = data.copy()
    # If there are no valid data points, returns one regime
    # that starts at 0, ends at the full length, and is nan
    if np.all(np.isnan(signal)):
        return np.array([0, len(signal)]), np.array([np.nan])

    # Gets the index of the first non-missing data
    ini = np.argmax(~np.isnan(signal))

    # Uses the Pelt method of the ruptures package
    # to find regime changes/breaks
    algo = rpt.Pelt(model=model).fit(signal[ini:])
    try:
        breaks = algo.predict(pen=pen)
    except BadSegmentationParameters:
        # If no changes/breaks were found, assumes
        # the regime goes up to the end of the series
        breaks = [len(signal)-ini]

    # Gets the indices of the breaks and 
    breaks = np.array([0, *breaks]) + ini
    # Computes the average of each regime
    averages = np.array([signal[breaks[i-1]:breaks[i]].mean() 
                         for i in range(1, len(breaks))])

    # First regime must be at least 10 years long
    i = np.argmax(np.diff(breaks) >= 10)
    return breaks[i:], averages[i:]