変化点の検出

直感 #

小売店の在庫率は、景気が良ければ商品がよく売れて 1.4 前後まで下がり、2008 年や 2020 年のように景気が悪化すると売れ残って跳ね上がります。折れ線を眺めれば「ここで水準が変わった」と感じる場所はありますが、人によって指す場所が違います。変化点検出は系列をいくつかの区間に切ります。各区間の中のばらつきの合計を最小にする区切り方を、Pelt などのアルゴリズムで探すわけです。区間を細かく切れば切るほどばらつきは減るので、区切り 1 本ごとに罰則 pen を払わせ、「ばらつきの減少が罰則を上回る場所だけ」を変化点として残します。

何が起きているか（ステップで見る） #

コードを読む前に、この手法が図の上で何をしているかを追ってみましょう。「再生」を押すか、スライダーを動かすと 1 手順ずつ進みます。

見るポイント：どこで区切ると左右のばらつきの合計が最小になるかを総当たりで探す。Pelt はこの探索を枝刈りで高速化し、pen で区切りの数を抑えます。

ruptures とは #

deepcharles/ruptures はオフライン（系列が全部そろった後の）変化点検出を集めた Python ライブラリです。アルゴリズムの詳細は次の論文で説明されています。

C. Truong, L. Oudre, N. Vayatis. Selective review of offline change point detection methods. Signal Processing, 167:107299, 2020.

ライセンスは BSD 2-Clause です。

BSD 2-Clause License

Copyright (c) 2017-2021, ENS Paris-Saclay, CNRS
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright notice, this
  list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

FRED からデータを取り出す準備 #

関数 get_fred_data() は FRED データベース のページと同じもので、系列 ID を渡すと日付索引の DataFrame を返します。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

import matplotlib.pyplot as plt
import os
import pandas as pd
import ruptures as rpt
import seaborn as sns
from full_fred.fred import Fred

# FRED_API_KEY = os.getenv('FRED_API_KEY')
fred = Fred()
print(f"FRED APIキーが環境変数に設定されている：{fred.env_api_key_found()}")

def get_fred_data(name, start="2013-01-01", end=""):
    df = fred.get_series_df(name)[["date", "value"]].copy()
    df["date"] = pd.to_datetime(df["date"])
    df["value"] = pd.to_numeric(df["value"], errors="coerce")
    df = df.set_index("date")

    if end == "":
        df = df.loc[f"{start}":]
    else:
        df = df.loc[f"{start}":f"{end}"]

    return df

FRED APIキーが環境変数に設定されている：True

小売業の在庫率（RETAILIRSA）を取得する #

Retailers: Inventories to Sales Ratio は小売業者の棚卸資産を月間売上高で割った比率で、『1 か月の売上高に対して手元にある在庫が何か月分あるか』の目安です。1.5 なら在庫が 1.5 か月分あることになります。2000 年以降の月次を取得します。

1

df_RETAILIRSA = get_fred_data("RETAILIRSA", start="2000-01-01")

Pelt で変化点を検出して描く #

Pelt に系列を渡して区切り位置を求め、元の折れ線に赤い縦線で重ねます。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

N = len(df_RETAILIRSA["value"])
signal = df_RETAILIRSA["value"].values.reshape((N, 1))

algo = rpt.Pelt(model="rbf").fit(signal)
result = algo.predict(pen=10)

sns.set(rc={"figure.figsize": (15, 5)})
sns.lineplot(data=df_RETAILIRSA, x="date", y="value")
for idx in result[:-1]:
    plt.axvline(x=df_RETAILIRSA.index[idx - 1], color="r")
plt.show()

実行結果の読み方 #

赤線は 2001 年、2008〜2009 年、2020 年前半の周辺に立ち、いずれも景気後退で在庫率が急上昇した時期と一致します。線と線の間では在庫率が 1.4〜1.6 の狭い範囲で安定しています。Pelt は「水準が切り替わる前後」を区間の境目として選んでいるわけです。2020 年の区切りの後、在庫率は 1.1 近くまで下がってから回復しており、1 本の線でもその前後で平均が大きく違うことが確認できます。

使いどころと落とし穴 #

pen に万能の既定値はなく、系列のスケールと長さで適正値が変わります。model="rbf" はカーネルの計算で値のスケールの影響を受けにくい一方、"l2" を使うときは値を標準化してから pen を決めないと、単位が変わるだけで変化点の数が変わります。Pelt はオフライン手法で、系列全体を見てから区切るので「今月が変化点か」をリアルタイムに判定することはできず、新しい月が追加されると過去の区切りが動くこともあります。区切りの数を自分で決めたい場合は Pelt ではなく rpt.Binseg(model="rbf").fit(signal).predict(n_bkps=3) のように n_bkps で本数を指定する方法に替えます。

まとめ #

変化点検出は、rpt.Pelt(model="rbf").fit(signal) と predict(pen=10) の 2 行で「性質が切り替わった位置」のリストを返します。signal は (N, 1) に整形した数値配列、pen は区切り 1 本あたりの罰則で、大きくすれば大きな変化だけが残ります。小売在庫率 RETAILIRSA では 2001 年・2008 年・2020 年の景気後退に赤線が立ち、折れ線を眺めて感じる「水準の切り替わり」が数値的に裏付けられました。pen を 3 と 30 で試して線の本数の変化を見るのが、この手法の感覚をつかむ一番の近道です。