データ点検

直感 #

新しい時系列データを受け取ったら、いきなり ARIMA や Prophet を当てるのではなく、まず「健康診断」をします。人間の健康診断が体重・血圧・血液検査の順に進むように、時系列では列型と欠損（info()・isnull()）、全体の形（折れ線）、値の分布（ヒストグラム）、周期（自己相関）、定常性（ADF 検定）の順に見ていきます。この 5 つを 10 分で済ませておくと、後で「実は Date 列が文字列だった」「夏だけ欠損していた」と気づいて最初からやり直す事故を防げます。

1. 使用ライブラリ #

実験で使う標準ライブラリとサードパーティをまとめて import しておきます。stattools は ADF 検定、stats は後で分布変換に使います。

1
2
3
4
5
6
7
8

import datetime
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

from scipy import stats
from statsmodels.tsa import stattools

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

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

見るポイント：埋め方で結果が変わる。前方補完は「変わらない」仮定、線形補間は「なめらかに変わる」仮定。

2. データを読み込む #

最初の 10 行を表示し、列名や桁数を把握しておきます。Date と Temp の 2 列だけという単純な構成で、1 行が 1 日の気温です。

1
2

data = pd.read_csv("sample.csv")
data.head(10)

3. 日付列を datetime 型に変換する #

object のままでは並べ替えや差分計算、.dt.year の取り出しができません。datetime.datetime.strptime で 1 行ずつ日付に変換します。

1
2
3
4
5

data["Date"] = data["Date"].apply(
    lambda x: datetime.datetime.strptime(str(x), "%Y-%m-%d")
)

print(f"Date列の dtype: {data['Date'].dtype}")

4. 基本統計・欠損・重複を確認する #

まず data.info() で総行数と各列の型を見て、Temp が float64、Date が datetime64[ns] になっているかを確かめます。次に data.isnull().sum() で列ごとの欠損個数を数え、0 でなければ上のステップ図で見た前方補完や線形補間のどちらで埋めるかを決めます。最後に data.duplicated().sum() で同じ日付が 2 回入っていないかを調べ、1 以上なら drop_duplicates(subset="Date") で片方を落とします（コードは省略、Notebook では print して確認）。

5. 折れ線でざっと俯瞰する #

点検の最初の 1 枚は折れ線です。10 年分を 1 枚に描くと、季節の山谷と長期トレンドの有無が同時に見えます。

1
2
3
4
5
6
7

plt.figure(figsize=(12, 6))
plt.plot(data["Date"], data["Temp"])
plt.grid(axis="x")
plt.xlabel("Date")
plt.ylabel("Temp")
plt.title("Daily temperature")
plt.show()

図では、2〜3 日続けて外れた点がないか、観測間隔が飛んで線が途切れる区間がないかを確かめます。気温なら年ごとの山谷が 10 回繰り返され、全体がゆるく右上がりになっていれば温暖化のトレンドの疑いが出ます。

6. トレンド線を重ねる #

3 次多項式でトレンドを推定するヘルパー関数を作ります。numpy.polyfit で係数を求め、poly1d で曲線の値に戻す 2 段構えです。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

def get_trend(timeseries, deg=3):
    """時系列のトレンド（多項式近似）を返す"""
    x = list(range(len(timeseries)))
    y = timeseries.values
    coef = np.polyfit(x, y, deg)
    trend = np.poly1d(coef)(x)
    return pd.Series(data=trend, index=timeseries.index)

data["Trend"] = get_trend(data["Temp"])

plt.figure(figsize=(12, 6))
sns.lineplot(x=data["Date"], y=data["Temp"], alpha=0.5, label="Temp")
sns.lineplot(x=data["Date"], y=data["Trend"], label="Trend")
plt.grid(axis="x")
plt.legend()
plt.show()

メモ: 序盤でトレンドの有無を知っておくと、差分やスケーリングをどこで導入すべきか判断しやすくなります。トレンド線がほぼ水平なら差分は不要で、右上がりなら 1 次差分か Temp - Trend の残差を使う候補になります。

7. 分布と外れ値を把握する #

値がどの範囲にどれだけ集まっているかをヒストグラムで見ます。気温なら 5 ℃ 付近と 25 ℃ 付近に 2 つの山ができる（冬と夏）ことが多く、これは正規分布ではありません。

1
2
3
4
5
6

plt.figure(figsize=(12, 6))
plt.hist(x=data["Temp"], rwidth=0.8)
plt.xlabel("Temp")
plt.ylabel("頻度")
plt.grid(axis="y")
plt.show()

ヒストグラムで偏りや多峰性を確認し、裾が長く偏っていれば Box-Cox やログ変換を検討します。気温の二峰性は夏と冬の混在が原因なので、変換よりも STL などの季節分解で成分を分けるほうが筋がよいです。

8. 自己相関で季節性を確認する #

pandas.plotting.autocorrelation_plot を使うと、ラグ 1 から 3,650 までの自己相関が 1 本の曲線で見えます。365 付近に垂線を引いておくと、年周期があるかどうかを一目で判断できます。

1
2
3
4
5
6
7

plt.figure(figsize=(12, 6))
pd.plotting.autocorrelation_plot(data["Temp"])
plt.grid()
plt.axvline(x=365)
plt.xlabel("Lag")
plt.ylabel("Autocorrelation")
plt.show()

9. 定常性テスト（ADF） #

Augmented Dickey-Fuller 検定で単位根（ランダムウォーク的な性質）の有無を判定します。帰無仮説は「単位根がある＝非定常」なので、p 値が 0.05 未満なら帰無仮説を棄却して定常とみなします。

1

stattools.adfuller(data["Temp"], autolag="AIC")

出力例:

(-4.44, 0.00024, 20, 3629, {'1%': -3.43, '5%': -2.86, '10%': -2.57}, 16642.82)

統計量 -4.44 は 1% 水準の臨界値 -3.43 より小さい（より負）ので、帰無仮説「単位根あり」は 1% 水準で棄却されます。p 値 0.00024 も 0.05 を大きく下回っており、この気温系列は定常とみなして AR 系のモデルに渡せます。

10. トレンド抽出の注意点 #

多項式トレンドには落とし穴が 3 つあります。最初の課題はスケールへの敏感さで、長期的に指数的な増加がある系列に 3 次多項式を当てると、終盤の大きな値に引きずられて序盤のトレンドが歪みます。この場合は対数変換や年代ダミーを事前に挟みます。二つ目は次数の上げ過ぎで、高次項は外れ値に引っぱられやすく端で跳ねるため、deg は 3〜4 で止めます。三つ目は備忘で、トレンド除去後の系列を Temp_detrended、元データを Temp_raw のように別名で持っておくと、後で「どちらを見ているか」を迷わずに比較できます。

次のコードは、元の気温に log(年 - 1980) * 10 という人工的な上昇トレンドを足し、同じ get_trend が右上がりの曲線を捉えるかを確かめる実験です。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

data_sub = data.copy()
data_sub["Temp"] = (
    data_sub["Temp"] + np.log(data_sub["Date"].dt.year - 1980) * 10
)
data_sub["Trend"] = get_trend(data_sub["Temp"])

plt.figure(figsize=(12, 6))
sns.lineplot(x=data_sub["Date"], y=data_sub["Temp"], alpha=0.5, label="Temp")
sns.lineplot(x=data_sub["Date"], y=data_sub["Trend"], label="Trend")
plt.grid(axis="x")
plt.legend()
plt.show()

11. まとめ #

データ点検は、info()・isnull().sum()・duplicated().sum() で型と欠損と重複を数える健康診断から始めます。次に折れ線・ヒストグラム・autocorrelation_plot の 3 枚で、トレンドの有無、分布の二峰性、365 日周期を目で確かめ、adfuller の p 値（この例では 0.00024）で定常性を数値で裏づけます。np.polyfit(deg=3) を包んだ get_trend は 7 行の小さな関数ですが、次のページのトレンド除去や季節分解でもそのまま使い回せるので、ここで一度書いておく価値があります。