ビジネスデータの予測に効く時系列分析の実践ガイド

序章 時間の順序を持つデータを分析する前に

順序そのものが情報を持つ

通常のデータ分析では、データの行を並べ替えても結論は変わりません。顧客1000人のアンケートは、どの順番で集計しても同じ結果になります。しかし時系列データは違います。順序そのものが情報であり、昨日と今日はつながっています。この「つながり」は予測を可能にする一方で、独立性を前提とする多くの統計手法の前提を崩す原因にもなります。

本ガイドの前半は、こうした前提の崩れを避けるための知識に充てています。自己相関、定常性、単位根、見せかけの回帰といった概念は、いずれも派手さはありませんが、これらを踏まえずに時系列データへ回帰分析や機械学習を適用すると、もっともらしく見えて実は誤った結論に至りやすくなります。月次データの分析でつまずく原因の多くは、モデルの高度さの不足ではなく、この基礎の見落としにあると考えています。たとえば前年同月比が改善したという報告も、季節性を取り除いて実力の推移を見なければ、改善なのか季節の綾なのかを判断できません。

予測の価値は不確実性を定量化することにある

本ガイドの後半は予測に充てています。ビジネスにおける予測の価値は、未来をぴたりと当てることよりも、不確実性を定量化して備えられるようにすることにあると考えています。来月の需要を「1万個」と1点で当てにいくのではなく、「9千から1万1千個の間に9割の確率で収まる」と幅で語れれば、在庫も人員もその幅に合わせて設計できます。だからこそ本ガイドでは、点予測だけでなく予測区間、そして予測を業務に組み込む設計（第12章）までを扱います。

古典的手法と機械学習の使い分け

時系列予測では、「ARIMAのような古典的手法と機械学習のどちらが優れているか」という議論が繰り返されてきました。本ガイドの立場は、どちらが勝つかはデータと課題の条件次第であり、必ずベンチマークとバックテストで検証して選ぶ、というものです。系列が少なく構造の解釈が要る場面では古典的手法が扱いやすく、大量の関連系列と豊富な特徴量がある場面では機械学習が有利になりやすいという傾向があります。両方を学ぶ利点は、この使い分けの判断軸を持てることにあると考えています。

本ガイドが一貫して重視する点は2つです。時間の順序がもたらす前提の崩れ（自己相関・非定常性）を避けること。そして予測を点ではなく幅で捉え、検証と運用まで設計すること。全12章はこの2つを軸に積み上げていきます。

Anagraftでは、AIプロジェクトの構想・課題設計から、データ分析・機械学習モデルの開発、AI人材の育成まで一貫したご支援を行っています。ご相談は、以下よりお問い合わせください。

お問い合わせ

全体の地図

本ガイドの構成です。第1部から順に積み上がりますが、各章は参照用に独立して読める形にもしています。

第1章 時系列データはなぜ特別なのか

こんにちは、Anagraftの伊藤です。本ガイドは「時間の順序を持つデータ」、すなわち時系列データを対象にした実践ガイドです。本シリーズではこれまで、統計学基礎編にあたる『データから正しく結論を導く統計学の実践ガイド』で記述統計や仮説検定の基礎を、機械学習編にあたる『機械学習の仕組みと使いどころがわかる実践ガイド』で予測モデリングの考え方を、回帰編にあたる『ビジネスの要因分析に効く回帰分析と統計モデリングの実践ガイド』で要因分析のための回帰分析を扱ってきました。回帰編の最終章では、時系列データに通常の回帰分析をそのまま当てはめると見せかけの回帰と呼ばれる誤った結論に導かれかねないことに触れ、時系列分析という専門領域への入口を示しました。本ガイドはその続きとして、時系列データそのものの構造や振る舞いを体系的に扱います。

月次の売上、日次のWebサイトのアクセス数、店舗ごとの在庫水準、工場の設備に取り付けたセンサーの計測値。経営層やDX推進担当者が日常的に向き合う数字の多くは、こうした時間軸に沿って並んだ時系列データです。ところが、これまでの統計学基礎編や機械学習編で前提としてきた「データの並び順には意味がない」という考え方は、時系列データには通用しません。同じ手元のデータであっても、時間の順序を無視して扱うか、順序を主役に据えて扱うかで、得られる結論も予測の精度も大きく変わってきます。本章では、時系列データが通常のデータ分析の常識からどのようにはみ出しているのかを整理し、本ガイド全体で扱う内容の見取り図を示します。

本ガイドが出発点に置く問題意識は単純です。「先月に比べて売上が10%増えた」という報告を受けたとき、それが事業の構造的な改善なのか、単に毎年この時期に売上が伸びる季節要因なのか、あるいはたまたま数か月だけ良い数字が続いただけなのかを、多くの現場では感覚と経験則で切り分けています。この切り分けを勘に頼らずデータで裏付けるためには、時系列データ特有の性質を踏まえた分析の枠組みが必要になります。同様に、需要予測の担当者が来月の発注量を決める際にも、過去の推移をどう読み解けば精度の高い見積もりになるのかは、通常の回帰分析や機械学習の知識だけでは十分に答えられません。時系列データを正しく扱うための最初の一歩は、まさにこの「なぜ通常のやり方が通用しないのか」を理解することにあります。

時系列データの4つの構成要素：トレンド・季節性・循環・不規則変動

時系列データを眺めるとき、最初に押さえておきたいのが、1本の折れ線グラフの動きは通常、性質の異なる複数の成分が重なり合ってできているという見方です。伝統的に、時系列データはトレンド、季節性、循環、不規則変動という4つの成分に分解して理解されます。

小売業の月次売上を例に考えると、この4成分の重なりが具体的にイメージしやすくなります。まず事業全体が成長していれば、売上は数年単位で緩やかに右肩上がりのトレンドを描きます。そこに、12月の年末商戦や夏のセールで売上が跳ね上がるという、1年周期で毎年繰り返される季節性が重なります。さらに、景気の良し悪しによって数年単位で売れ行きの波が生じる循環成分が加わり、最後に、特定の週にだけ発生した突発的なキャンペーンの反響のような、パターン化できない不規則変動が残ります。日次データであれば、平日と週末で来店客数が大きく変わる曜日効果も、季節性の一種として扱われる代表的な例です。

循環と季節性は混同されやすいため、区別のポイントを補足しておきます。季節性は「1年」「1週間」のように周期の長さがあらかじめ決まっている変動であるのに対し、循環は周期の長さが一定でない、より緩やかで不規則なうねりを指します。景気循環にともなう企業の設備投資の増減や、業界特有の数年単位の受注サイクルなどが典型例です。実務データでは循環と長期トレンドの境目があいまいになりがちで、両者を厳密に切り分けることよりも、「1年や1週間といった決まった周期で規則的に繰り返す部分（季節性）」と「それより長く、周期の長さも一定しない緩やかな動き（トレンドと循環をまとめたもの）」という大掴みの区別で扱うことが多くなります。不規則変動についても補足すると、これは理想的にはどの時点とも相関を持たない、いわゆるホワイトノイズに近い性質を持つことが期待される成分です。裏を返せば、トレンドと季節性を取り除いた後の残差になお強い自己相関やパターンが残っているとすれば、それはモデル化しきれていない構造がまだ残っているという診断のサインになります。

実務でこの4成分を意識する意味は、変化の「原因」を大まかに切り分けられる点にあります。前月比で売上が落ちたとき、それが事業構造そのものの悪化（トレンドの転換）なのか、単に季節的に売上が落ちる月に入っただけなのか、それとも一時的な特殊要因（不規則変動）なのかを区別できなければ、対策の方向性を誤ります。この4成分への分解を体系的に行う手法は第2章で詳しく扱い、季節性を取り除いた季節調整済みの系列を作る方法まで踏み込みます。

自己相関：「昨日と今日は似ている」をどう測るか

時系列データの最大の特徴は、隣り合う時点の値どうしが独立ではなく、互いに似通う傾向を持つという点にあります。今日の売上が良ければ明日の売上もある程度良い可能性が高く、今月のアクセス数が多ければ来月のアクセス数もある程度多い可能性が高いというのは、多くのビジネスデータで直感的に成り立つ経験則です。この「昨日と今日は似ている」という関係を数値として測る道具が自己相関です。

自己相関を測るには、まずラグという考え方を理解する必要があります。ラグとは、ある時点から何時点分さかのぼるかを表す整数で、ラグ1は1時点前、ラグ12は12時点前（月次データであれば1年前）を指します。自己相関係数（ACF、Autocorrelation Function）は、ある系列とそれ自身をラグ分だけずらした系列との間の相関係数であり、ラグごとに1つの値が定まります。ラグ1の自己相関係数が高ければ「隣り合う時点どうしがよく似ている」ことを、ラグ12の自己相関係数が高ければ「1年前の同じ時期とよく似ている（季節性が強い）」ことを意味します。

複数のラグにわたる自己相関係数を棒グラフとして並べたものがコレログラム（自己相関プロット）です。コレログラムには通常、統計的に有意な自己相関かどうかを判断するための信頼区間の帯が併せて描かれ、棒がこの帯の外に突き出ているラグは、偶然とは考えにくい自己相関があると判断できます。月次売上のコレログラムを描いたとき、ラグ1の棒が高く、そこから緩やかに減衰していく形はトレンドの影響を、ラグ12の付近で棒が突出する形は季節性の影響を示唆します。この読み方は、後の章で扱うARモデルやMAモデルにおいて、モデルの次数（過去何時点分をモデルに取り込むか）を決める手がかりとしても使われます。

自己相関の強さは、業種や指標によっても大きく異なります。在庫水準のように前日の値がそのまま翌日の出発点になる指標は、ラグ1の自己相関が非常に高くなる傾向があります。一方、日々のキャンペーン当選者数のように毎回の抽選で値が決まる指標は、自己相関がほとんど見られないことも珍しくありません。自己相関の強さそのものが、その指標がどれだけ「慣性を持って動くか」を映す1つの手がかりになる、という見方も実務では有効です。

独立性が崩れると何が起きるか：見かけの情報量の水増し

統計学基礎編で扱った仮説検定や区間推定の多くは、手元のデータが互いに独立に得られたという前提の上に成り立っています。標準誤差の計算式も、p値の計算式も、100個のデータ点があれば100個分の独立な情報が手に入っているという想定で組み立てられています。ところが時系列データでは、前節で見た通り隣り合う時点の値どうしが強く似通っているため、この前提が崩れます。

この崩れが引き起こす実害が、有効サンプルサイズの減少です。自己相関の強いデータでは、見かけ上のデータ数（観測点の数）ほどには、独立な情報が含まれていません。たとえば、隣り合う月の売上が非常によく似た値を取り続けるデータが72か月分（6年分）あったとしても、実質的に独立とみなせる情報の量は72よりもずっと少なくなります。にもかかわらず、通常の検定や区間推定の計算式は観測点の数をそのまま独立な情報の数として扱うため、実際の不確実性よりも標準誤差を小さく見積もり、本来は誤差の範囲内に収まる程度の違いを「統計的に有意な差」と誤って判定してしまう危険が生まれます。

自己相関があるデータでは、観測点の数がそのまま独立な情報の数にはなりません。100点のデータを集めても、自己相関が強ければ実質的な情報量は数十点分程度にとどまることがあり、この水増しされた情報量をもとに検定やモデルの当てはまりを評価すると、実態よりも自信過剰な結論に至りやすくなります。時系列データを扱う際は、常にこの「見かけの情報量」と「実質的な情報量」のずれを意識する姿勢が欠かせません。

この問題は、回帰編で扱った見せかけの回帰とも根が同じです。トレンドを持つ2つの無関係な系列を回帰にかけると決定係数やp値が実態以上に強い関係を示してしまう現象も、突き詰めれば独立性の前提が崩れていることに由来します。この独立性の崩れをどう診断し、どう対処するかは、第3章で扱う定常性という概念と、第4章で扱う見せかけの回帰・共和分の議論において、本ガイドを通じて繰り返し立ち返る土台になります。

なぜランダム分割の交差検証は使えないのか

機械学習編では、モデルの汎化性能を評価する手段として交差検証を扱い、データをランダムにシャッフルして訓練用とテスト用に分割するk分割交差検証（k-Fold CV）を紹介しました。この手法は、データの各行が互いに独立であるという前提のもとでは合理的な検証方法です。しかし時系列データにそのまま適用すると、深刻な問題を引き起こします。

ランダム分割を時系列データに適用すると、時間的に未来にあたる観測点が訓練データに、過去にあたる観測点がテストデータに入り込むという組み合わせが起こり得ます。これは、実運用では絶対に手に入らない「未来の情報」を使って「過去」を予測しているのと同じ状況であり、リーク（データ漏洩）の一種です。自己相関が強いデータでは、訓練データに含まれる未来の点とテストデータに含まれる過去の点が時間的に近接しているケースも多く、モデルは実質的に「ほぼ答えが見えている」状態で評価されてしまいます。その結果、ランダム分割による交差検証のスコアは実態よりも大幅に楽観的な数値を示し、実運用に投入した途端に精度が大きく低下するという事態を招きます。

この問題への対処が、時間の順序を守った分割です。scikit-learnのTimeSeriesSplitのように、常に時間的に前のデータを訓練用、それより後のデータをテスト用に割り当てる分割方法を使うことで、モデルは実運用と同じ「過去から未来を予測する」条件のもとで評価されます。時系列データを対象にモデルの精度を検証する際は、通常のk-Fold CVではなく時間順を守った分割を用いているかを必ず確認する必要があります。この論点は、第10章で機械学習による時系列予測を扱う際にあらためて具体的な実装とともに整理します。

時系列分析の2つの目的：理解と予測

時系列分析に取り組む目的は、大きく理解と予測の2つに整理できます。理解を目的とする分析は、系列そのものの構造（トレンドの強さ、季節性の大きさ、他の系列との関係性など）を把握し、意思決定の材料とすることに主眼を置きます。たとえば、売上のトレンドが本当に上向いているのか、それとも季節性による一時的な押し上げに過ぎないのかを見極めることや、広告費の系列と売上の系列の間に、単なる見せかけではない実質的な関係があるのかを検証することが、この目的にあたります。

一方、予測を目的とする分析は、これまでの値の推移から、まだ観測されていない将来の値を数値として当てにいくことに主眼を置きます。来月の売上がいくらになるか、来週の在庫がどの水準まで減るかを見積もる需要予測は、この目的の典型例です。理解と予測は互いに独立した作業ではなく、多くの場合、系列の構造を正しく理解できていることが、精度の高い予測を組み立てるための土台になります。

本ガイドの構成は、この2つの目的を意識して組み立てられています。第2章から第4章にかけては、季節調整、定常性、見せかけの回帰と共和分という、系列の構造を正しく理解するための基礎理論を扱います。第5章から第9章にかけては、ARモデル・MAモデル・ARIMA・SARIMA・VARモデル・状態空間モデルという、理解と予測の両方に使える代表的な統計モデル群を扱います。第10章では機械学習を用いた時系列予測、第11章では異常検知と変化点検知という、予測とは異なる角度から時系列データを活用する手法を扱い、最終章となる第12章では、これらすべてを踏まえた需要予測の実務設計を扱います。

Pythonによる時系列データの基本操作

ここからは、本章の締めくくりとして、以降の章で繰り返し使う時系列データの基本的な扱い方をpandasとstatsmodelsで確認します。時系列データを扱う実務の出発点は、日付や時刻の情報をDatetimeIndex（日時型のインデックス）としてデータフレームに持たせることです。これにより、日付をキーにした抽出や、後述するresample・shift・rollingといった時系列専用の操作が可能になります。

まず、トレンドと季節性を持つ月次売上の合成データを作成し、そのコレログラムをstatsmodelsのplot_acf関数で描画します。合成データを使うのは、真の構造（トレンドの傾き、季節性の強さ）があらかじめ分かっている状態でACFの読み方を確認するためです。

なお、本ガイドのコード例では図に日本語のラベルを付けています。matplotlibは既定のフォントで日本語を表示できないため、そのまま実行するとラベルが文字化けします。japanize-matplotlibを導入するか、plt.rcParams["font.family"] = "Meiryo"のようにフォントを指定してから実行してください。

import numpy as np
import pandas as pd
from statsmodels.graphics.tsaplots import plot_acf

np.random.seed(0)

# 2018年1月から2024年12月までの月次売上を合成する
dates = pd.date_range("2018-01-01", "2024-12-01", freq="MS")
n = len(dates)

trend = np.linspace(100, 220, n)  # 緩やかな右肩上がりのトレンド
month = dates.month
seasonality = 15 * np.sin(2 * np.pi * (month - 3) / 12) + np.where(month == 12, 25, 0)
noise = np.random.normal(0, 5, n)

sales = trend + seasonality + noise
df = pd.DataFrame({"sales": sales}, index=dates)
df.index.name = "date"

plot_acf(df["sales"], lags=24)

freq=”MS”は月初（Month Start）を意味し、date_rangeで生成される各日付が各月の1日になるよう指定しています。このデータのACFプロットを確認すると、ラグ1から緩やかに減衰していく形が見られ、これはトレンドの影響によるものです。トレンドを持つ系列は、直前の値との相関が常に高くなる性質があるため、こうした緩やかな減衰はトレンドの存在を疑うサインになります。加えて、ラグ12の付近で減衰の傾きが緩む、あるいは棒が再び高くなる様子が見られれば、それは1年周期の季節性が残っていることを示唆します。

次に、このトレンドを大まかに取り出すための移動平均（rolling）と、月次データを四半期単位に集約するresampleを確認します。

# 12か月移動平均でトレンド成分だけを大まかに取り出す
df["rolling_12"] = df["sales"].rolling(window=12, center=True).mean()

# 月次データを四半期ごとの合計値に集約する
quarterly = df["sales"].resample("QE").sum()
print(quarterly.head())

rolling(window=12, center=True)は、各時点のおおむね中央に値を置く12か月の移動平均になります（偶数窓のため厳密な中心化には2段階の移動平均を使います。詳しくは第2章で扱います）。centerをTrueにすることで、平均の対象期間の中央にその値が対応するように揃えられ、トレンドの転換点をより実態に近いタイミングで捉えられます。resample("QE")は、月次データを四半期末（Quarter End）を基準に集約する処理で、sum()と組み合わせることで四半期ごとの売上合計を得られます。resampleはこのほか、日次データを週次や月次に集約する場合や、逆に月次データを日次に補間する場合など、時間の粒度を変換するあらゆる場面で使う基本操作です。

最後に、前月比・前年同月比といった変化量を作るdiffと、過去の値をそのまま参照するshift、そして自己相関係数を数値として直接計算する方法を確認します。

# 前月比・前年同月比を計算する
df["mom_diff"] = df["sales"].diff(1)     # 前月との差(前月比の変化量)
df["yoy_diff"] = df["sales"].diff(12)    # 12か月前との差(前年同月比の変化量)
df["prev_month"] = df["sales"].shift(1)  # 前月の値そのもの

# 自己相関係数を数値でも確認する(ラグ1とラグ12)
from statsmodels.tsa.stattools import acf

acf_values = acf(df["sales"], nlags=13)
print(f"ラグ1の自己相関係数: {acf_values[1]:.3f}")
for lag in (11, 12, 13):
    print(f"ラグ{lag}の自己相関係数: {acf_values[lag]:.3f}")

diff(1)は各時点の値から1時点前の値を引いた差分を、diff(12)は12時点前（月次データであれば1年前の同じ月）との差分を計算します。shift(1)はdiffとは異なり、値そのものを1時点分後ろにずらして持ってくる操作で、「前月の値を今月の行に並べて比較する」といった用途に使われます。diffとshiftはいずれも、次章以降で扱う季節調整や定常性の確認、ARIMAモデルの構築において土台となる基本操作です。acf関数は、コレログラムの元になっている自己相関係数を数値として直接取得するための関数です。このデータではラグ1の自己相関係数が0.90と高く、ラグ12は0.60で、隣接するラグ11（0.57）やラグ13（0.52）より高い局所的な盛り上がりになります。季節性のサインとして読み取るべきなのは、絶対値の高さではなく、この1年周期での盛り上がりのほうです。グラフで大まかな傾向をつかんだうえで、重要なラグについては数値でも裏付けを取るという進め方が実務では有効です。

本章で扱った内容を、新しい時系列データに向き合う際の確認事項として整理すると、次のようなチェック項目にまとめられます。

• データにDatetimeIndexを持たせ、日付・時刻の抜け漏れや重複がないかを確認したか
• 折れ線グラフを描き、トレンド・季節性・循環・不規則変動のどれが目立って見えるかを大まかに把握したか
• ACF（コレログラム）を確認し、どのラグで自己相関が強く出ているかを確認したか
• 自己相関が強いデータに対して、独立性を前提とする通常の検定や信頼区間をそのまま当てはめていないか
• モデルの検証で、ランダム分割ではなく時間順を守った分割（TimeSeriesSplitなど）を使っているか

本章では、時系列データがトレンド・季節性・循環・不規則変動という4成分の重なりとして理解できること、そして自己相関という性質が独立性を前提とする通常の統計的手法の土台を崩し、有効サンプルサイズの減少や交差検証の設計ミスといった具体的な問題を引き起こすことを確認しました。次章では、この4成分をデータから実際に取り出す分解の手法と、季節性を取り除いた系列を作る季節調整の具体的な方法を扱います。

この章を深めたい方への参考書籍

『現場ですぐ使える時系列データ分析』（横内大介・青木義充、技術評論社）：Pythonの実装コードを豊富に交えながら、自己相関やトレンド・季節性の見方から状態空間モデルまでを実務目線で解説した一冊です。本章で扱ったACFやresample・rollingの読み方を、より多くの実データの例で確認したい読者に向いています。

第2章 トレンドと季節性を分解して時系列を理解する

前章では、時系列データが自己相関を持つこと、そしてトレンドや季節性といった性質を備えていることを確認しました。この章では、その性質のうち「トレンド」と「季節性」を実際に分けて取り出す方法、すなわち時系列分解を扱います。

ビジネスの現場で売上や来店客数を分析する際、最も広く使われている比較方法は前年同月比でしょう。しかしこの方法には見落とされがちな弱点があります。比較の基準にしている「前年の同じ月」が、たまたま一時的な要因で異常に高い、あるいは低い数値だった場合、その異常を引きずったまま今年の数字を評価することになるためです。今月の前年同月比がマイナスだったとして、それが事業の実力の低下を意味するのか、単に前年の基準が高すぎただけなのかを、前年同月比という1つの数字だけから判別することはできません。この区別をつけるための道具が、次節以降で扱う時系列分解です。

この章で扱う内容は次の通りです。

• 「分解して見る」ことの価値：前年同月比だけに頼る分析の限界と、成分分解によって実力と季節要因を切り分ける考え方
• 移動平均によるトレンド抽出：中心化移動平均の仕組み、窓幅の選び方（12か月・7日）、端点で生じる欠損の問題
• 古典的分解：加法モデルと乗法モデルの使い分け、対数変換によって乗法モデルを加法モデルに直す方法
• STL分解：Loessに基づく頑健な分解の実務標準、statsmodelsでの実装と外れ値への強さ
• 季節調整の考え方：季節調整済み系列とは何か、GDPなど公的統計での位置づけ
• 分解結果の読み方の実務：トレンドの転換点、季節パターンの変化、残差に残るイベント効果
• 曜日効果・祝日効果・営業日数：月次データの営業日調整と、日次データの曜日ダミーという発想
• ビジネス適用：「売上が下がったのは季節のせいか実力か」を分解で答える方法と、KPIモニタリングでの活用

「分解して見る」ことの価値

前年同月比は、比較の基準を「前年の同じ月」というただ1点の観測値に置いています。この1点が、大型のプロモーションや悪天候による臨時休業といった一時的な要因で平年並みの水準からずれていた場合、前年同月比の数字はそのずれをそのまま引き継ぎます。今月の前年同月比がマイナスだったとしても、それが事業の実力の低下によるものか、前年の基準が特殊に高かっただけなのかを、前年同月比という指標そのものから読み取ることはできません。

時系列分解は、観測された系列を複数の成分に分けて捉えることで、この区別を可能にします。売上の系列であれば、長期的な水準の動きを表す「トレンド成分」、1年周期などで繰り返すパターンを表す「季節成分」、そしてトレンドと季節では説明しきれない残りの変動を表す「残差成分」に分けて考えます。季節成分を取り除いた後のトレンド成分の動きを見れば、事業の「実力」に近い部分がどう変化しているかを、前年の1点に左右されずに確認できます。分解は観測期間全体のパターンから季節成分を推定する手法であるため、1年前のたった1つの月が特殊な事情を抱えていたとしても、その影響を受けにくいという利点があります。

「今月の売上は前年同月比でマイナスだった、事業の実力が落ちたのだろうか」という問いに、前年同月比という1つの数字だけで答えを出すことはできません。トレンド成分が横ばいか上向きであれば、マイナスの主因は季節要因や前年の基準の特殊性である可能性が高く、トレンド成分自体が下向きに転じていれば、実力の低下を疑う根拠になります。分解によって初めて、この2つの可能性を切り分けられます。

移動平均によるトレンド抽出

トレンド成分を取り出すための最も基本的な道具が移動平均です。ある時点の前後一定期間の値を平均することで、季節性や短期的なノイズをならし、緩やかな水準の動きだけを残す考え方です。実務でよく使われるのは、月次データに対する12か月移動平均と、日次データに対する7日移動平均です。窓幅を12か月にとると、1年周期の季節パターンをちょうど1周期分含んだ平均になるため季節変動が打ち消され、窓幅を7日にとると、曜日による1週間周期の変動が打ち消されます。窓幅を短くしすぎると季節性が残ったままトレンドに混入し、長くしすぎるとトレンドの転換に反応するまでに時間がかかるという、感度と滑らかさのトレードオフがあります。

ここで実務上つまずきやすいのが「中心化」の扱いです。ある時点\(t\)のトレンドを求めるには、その時点を中心にして前後対称に値を集めた平均、すなわち中心化移動平均を使うのが基本です。窓幅が7のように奇数であれば、前後3つずつと自分自身を合わせた7点の平均が、そのまま\(t\)を中心とした対称な平均になります。ところが窓幅が12のように偶数の場合、前後にちょうど半分ずつという対称な区切り方ができません。このため古典的な分解では、まず12か月移動平均を2回、位置を1つずらして計算し、その2つを平均するという「2×12移動平均」の手続きを踏んで、\(t\)を中心とする値に合わせ直します。

もう1つの実務上の制約が端点の問題です。移動平均は前後の値を使って計算するため、系列の最初と最後の数点では、窓に必要な数のデータがそろわず、トレンドの値を計算できません。窓幅12か月の中心化移動平均であれば、最初と最後のそれぞれ6か月分がこの理由で欠損します。直近の数か月のトレンドを移動平均だけで評価しようとすると、まさに知りたい「直近」の部分がデータの欠損によって埋まらないという、実務上悩ましい制約に直面することになります。この端点問題は、後述するSTL分解でも形を変えて残る論点です。

import numpy as np
import pandas as pd

rng = np.random.default_rng(0)
n = 60  # 5年分の月次データ

trend_true = 100 + 1.2 * np.arange(n)
seasonal_true = 15 * np.sin(2 * np.pi * np.arange(n) / 12)
noise = rng.normal(0, 5, n)
sales = trend_true + seasonal_true + noise

idx = pd.date_range("2021-01-01", periods=n, freq="MS")
ts = pd.Series(sales, index=idx, name="sales")

# 単純な12か月移動平均(center=Trueでも偶数窓は厳密には中心化されない)
ma12 = ts.rolling(window=12, center=True).mean()

# 2×12移動平均による中心化(古典的分解と同じ手続き)
ma_2x12 = ma12.rolling(window=2).mean().shift(-1)

print(ma12.head(7))       # 先頭6か月はNaN(端点の問題)
print(ma_2x12.head(7))
print(ma_2x12.tail(7))    # 末尾6か月もNaN。直近側にも同じ制約が出る

このコードでは、まず12か月移動平均（ma12）を求め、それをさらに2期分平均して1つ分だけ位置をずらす（ma_2x12）ことで、偶数窓を中心化しています。出力を確認すると、どちらの系列も先頭6か月分がNaNになっており、末尾（直近側）にも同じくNaNが並びます。これが端点の問題であり、まさに知りたい直近数か月分のトレンド値が、移動平均だけでは得られないことの具体的な現れです。なお2×12移動平均（ma_2x12）は先頭・末尾ともに6か月分が欠損しますが、単純な12か月移動平均（ma12）はpandasの偶数窓の位置合わせの都合で末尾の欠損が5か月分になります。

古典的分解：加法モデルと乗法モデル

移動平均でトレンドを取り出した後、観測値からトレンドを差し引く（あるいは割る）ことで、季節性と残差の部分を取り出せます。この考え方を体系化したものが古典的分解（classical decomposition）です。観測値\(y_t\)を、トレンド\(T_t\)・季節\(S_t\)・残差\(R_t\)の3成分にどう組み合わせるかによって、加法モデルと乗法モデルの2通りの定式化があります。

加法モデルは \( y_t = T_t + S_t + R_t \) という形で、季節変動の大きさ（振幅）が、トレンドの水準にかかわらずほぼ一定である場合に適しています。例えば毎年12月に売上が「一定額（例えば200万円）だけ」上乗せされるような系列です。一方、乗法モデルは \( y_t = T_t \times S_t \times R_t \) という形で、季節変動の振幅がトレンドの水準に比例して大きくなる場合に適しています。事業が成長してトレンドの水準そのものが年々高くなっているのに合わせて、12月の季節的な上乗せも「一定額」ではなく「一定の割合（例えば通常月の1.3倍）」で大きくなっていくような系列がこれにあたります。売上や来店客数のように、水準が大きく成長する系列では、乗法モデルの方が実態に合うことが多くなります。

どちらのモデルを使うべきか迷う場合の簡便な判断材料は、系列をプロットして季節変動の振幅がトレンドの水準と一緒に大きくなっているかどうかを目で確認することです。振幅が水準に比例して広がっているように見えれば乗法モデル、振幅がほぼ一定に見えれば加法モデルを選びます。もう1つの実務的な解決策は、観測値を対数変換してしまうことです。乗法モデルの両辺の対数をとると、次の関係が成り立ちます。

\( \ln(y_t) = \ln(T_t) + \ln(S_t) + \ln(R_t) \)

対数を取った系列に対しては、加法モデルとして分解を実行できます。つまり乗法モデルは、対数変換という一手間を挟むことで、加法モデルの枠組みに帰着させることができます。statsmodelsの実装では、model引数に”additive”か”multiplicative”を指定するだけで乗法モデルを直接扱うこともできますが、対数変換を経由する方法は、後の章で扱う回帰やARIMAといった他の手法とも組み合わせやすいという利点があります。

STL分解：頑健な分解の実務標準

古典的分解には実務上の弱点があります。移動平均に基づく手続きのため、端点でトレンドが欠損すること、季節成分が観測期間全体を通じて完全に一定であることを前提としているため季節パターンがゆっくり変化するような系列にうまく対応できないこと、そして外れ値の影響を強く受けやすいことです。これらの弱点を補う手法として実務で広く使われているのがSTL（Seasonal and Trend decomposition using Loess）分解です。

STLは、Loess（局所重み付き回帰）と呼ばれる平滑化手法を繰り返し適用することで、トレンド成分と季節成分を推定します。Loessは、注目する時点の近傍のデータだけを使って局所的な回帰を行う手法で、遠くのデータの影響を受けにくいという性質があります。この性質のおかげで、STLは季節成分が年ごとに緩やかに変化することを許容でき、また外れ値の影響を抑える頑健（robust）なオプションも備えています。robustオプションを有効にすると、大きく外れた残差を持つ観測値の重みを繰り返し計算の中で小さくしていくため、一時的なイベントによる異常値がトレンドや季節成分の推定を歪める度合いが抑えられます。

statsmodelsでは、statsmodels.tsa.seasonal.STLクラスでSTL分解を実行します。周期を表すperiod引数（月次データであれば12）と、季節成分の滑らかさを制御するseasonal引数（奇数を指定する必要があり、値が大きいほど季節成分は年をまたいでより一定になります）、そして頑健性を有効にするrobust引数が主な設定項目です。従来型のseasonal_decompose関数と比較すると、STLの方が外れ値に対して頑健であることが多くの実務例で確認されています。

import numpy as np
import pandas as pd
from statsmodels.tsa.seasonal import seasonal_decompose, STL

rng = np.random.default_rng(1)
n = 84  # 7年分の月次データ
idx = pd.date_range("2018-01-01", periods=n, freq="MS")

trend_true = 200 + 1.0 * np.arange(n)
seasonal_true = 20 * np.sin(2 * np.pi * np.arange(n) / 12)
noise = rng.normal(0, 6, n)
sales = trend_true + seasonal_true + noise

# 50番目の月に一時的な外れ値(入力ミスや突発的な特需を想定)を混入
sales[50] += 150
ts = pd.Series(sales, index=idx, name="sales")

result_classic = seasonal_decompose(ts, model="additive", period=12)
result_stl = STL(ts, period=12, seasonal=13, robust=True).fit()

print("外れ値月の残差(古典的分解):", round(result_classic.resid.iloc[50], 1))
print("外れ値月の残差(STL, robust=True):", round(result_stl.resid.iloc[50], 1))
print("外れ値月周辺のトレンド(古典的分解):", result_classic.trend.iloc[48:53].round(1).tolist())
print("外れ値月周辺のトレンド(STL):", result_stl.trend.iloc[48:53].round(1).tolist())

このコードでは、50番目の月に一時的な外れ値を混入させた系列に対し、古典的分解（seasonal_decompose）とSTL(robust=True)を並べて実行しています。古典的分解では、外れ値の影響が移動平均を通じてトレンド成分にまで及びます。混入させた+150のうち12か月移動平均の重み分（150÷12＝12.5）がそのままトレンドに乗るため、外れ値の前後およそ11か月にわたって、トレンドが真の水準（248〜252）より12から13ポイント高い260から265へ押し上げられます。一方STLでは、robustオプションによって外れ値の重みが繰り返し計算の中で小さくなるため、トレンドは247.8から252.1と真の水準をほぼそのまま保ち、外れ値の影響のほぼ全量が残差成分（152.3）に吸収されます。実務データには、キャンペーンの一時的な効果や、システム上の計上ミスなど、この種の外れ値がしばしば混ざるため、頑健性を備えたSTLを標準の分解手法として使うことが実務上の定石になっています。

季節調整の考え方

分解によって得られた季節成分を観測値から取り除いた系列を、季節調整済み系列（seasonally adjusted series）と呼びます。加法モデルであれば「観測値マイナス季節成分」、乗法モデルであれば「観測値割る季節成分」によって作成します。対数変換してから加法モデルで分解した場合は、対数の世界で「観測値マイナス季節成分」を計算し、元の単位に戻すときに指数を取ることになるため、結果としては乗法モデルと同じ「割り算」に相当します。季節調整済み系列は、トレンド成分と残差成分だけを残した系列であり、季節による周期的な変動を取り除いた上で、期間同士を単純に比較できるようにする目的で使われます。

季節調整は、企業のKPI分析だけで使われる特殊な操作ではありません。GDP（国内総生産）の四半期速報値をはじめとする公的統計の多くは、季節調整済みの数値を前期比という形で公表しています。これは、四半期ごとに存在する構造的な季節パターン（年度末にあたる四半期に公共投資が集中しやすい、など）を取り除かないまま前期と比較すると、実力の変化ではなく季節的な変動の影響を受けた比較になってしまうためです。統計を作成する機関では、X-12-ARIMAやX-13ARIMA-SEATSと呼ばれる季節調整専用の手法が標準的に用いられており、季節調整という操作自体は、公的統計の作成実務においても広く確立された標準的な手続きです。企業のKPI分析における季節調整済み系列の利用は、この公的統計と同じ発想を、自社のデータに適用したものと捉えることができます。

import matplotlib.pyplot as plt

# STLの季節成分を使って季節調整済み系列を作成(加法モデルなので引き算)
seasonally_adjusted = ts - result_stl.seasonal

fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(ts.index, ts, label="原系列", alpha=0.5)
ax.plot(ts.index, seasonally_adjusted, label="季節調整済み系列", linewidth=2)
ax.plot(ts.index, result_stl.trend, label="STLトレンド成分", linestyle="--")
ax.legend()
ax.set_title("原系列・季節調整済み系列・トレンド成分の比較")
plt.tight_layout()
plt.savefig("seasonal_adjustment.png")

このコードでは、先ほど求めたSTLの季節成分（result_stl.seasonal）を原系列から差し引くことで、季節調整済み系列を作成しています。原系列は季節ごとの上下動を伴ってジグザグに動きますが、季節調整済み系列はその上下動がならされ、トレンド成分に沿った動きに近づきます。ただしこの系列には、前節で頑健性を確かめるために混入させた外れ値がそのまま残っており、その月だけがトレンドから150ポイント以上跳ね上がったスパイクとして残ります。季節調整は季節性を取り除く操作であって、外れ値を取り除く操作ではないという点も、この図から読み取れます。KPIダッシュボードで前月比を語る際には、原系列の前月比ではなく、この季節調整済み系列の前月比を使うことで、季節的な要因による見かけの増減に振り回されずに済みます。

分解結果の読み方の実務

分解を実行して満足するのではなく、得られたトレンド・季節・残差の3成分を、それぞれ意味のある形で読み解くことが実務では重要です。

トレンド成分を読む際にまず注目すべきは転換点です。トレンドが右肩上がりから横ばいに転じた、あるいは横ばいから右肩下がりに転じた時点がどこにあるかを確認することで、事業の実力がいつから変化し始めたのかを、季節要因に惑わされずに特定できます。ある月の実績が悪化していても、トレンド成分がまだ上向きを保っていれば、それは一時的な季節要因やイベントの影響である可能性が高いと判断できますし、逆に良好な実績が続いていてもトレンド成分がすでに横ばいに転じていれば、季節要因による下支えがなくなった途端に数字が崩れる可能性を警戒できます。

季節成分についても、一定不変のものとして扱わない視点が必要です。STLのように季節成分が緩やかに変化することを許容する手法を使うと、過去数年間で季節パターンそのものが変わってきたかどうかを確認できます。例えば、新型コロナウイルス流行の前後で、外出やイベント関連の消費データにおいて、それまでの季節パターン（繁忙期・閑散期の位置づけ）が変化した企業は少なくありません。生活様式や消費行動の変化によって季節パターンが変わるのであれば、古いデータで固定的に推定した季節成分を使い続けることは、将来の季節調整を誤らせる原因になります。季節成分の変化そのものも、定点観測しておくべき対象です。

残差成分は「トレンドでも季節性でも説明できない部分」であり、その中には偶発的なノイズだけでなく、特定できるイベントの効果が残っていることがよくあります。大型キャンペーンの実施月、想定外の悪天候、システム障害による営業停止といった出来事は、トレンドにも季節性にも当てはまらないため、残差成分に大きなスパイクとして現れます。残差成分の大きな値を放置せず、その月に何が起きていたかを事業側の記録と突き合わせる作業は、分解結果を単なる統計的な副産物で終わらせず、意思決定に活かすための重要な一手間です。

曜日効果・祝日効果・営業日数

ここまでの分解は、主に1年周期の季節性を念頭に置いてきましたが、実務のデータにはもう1つ、営業日数や曜日構成に由来する変動が存在します。月次データでは、月によって営業日数（平日の日数から祝日を除いた実質的な営業可能日数）が異なることが、売上や来店客数の月次比較を歪める要因になります。2月は他の月より日数が少なく、ゴールデンウィークを含む月は祝日で営業日が減るといった具合に、月ごとの営業日数の違いは無視できません。これに対する実務上の対処が営業日調整で、月間の合計値を営業日数で割って「1営業日あたりの平均」に直してから月同士を比較する、あるいは営業日数そのものを説明変数として回帰モデルに組み込む、という2通りの方法があります。

日次データでは、営業日数という単位では粗すぎるため、曜日そのものを説明変数として扱う発想が使われます。平日と土日祝日とで来店客数の水準が大きく異なる業態であれば、月曜から日曜までの各曜日を表すダミー変数（その曜日であれば1、それ以外は0）を作り、回帰モデルに組み込むことで、曜日による変動を明示的に推定できます。これは分解によって季節成分を推定する考え方と、回帰モデルによって説明変数の効果を推定する考え方が、根っこでは同じ発想を共有していることを示しています。曜日や祝日を、単に「ノイズとして均してしまう対象」ではなく、「効果を推定して活用できる説明変数」として扱えることは、需要予測やシフト計画の精度を上げる上で実務的な価値があります。

ビジネス適用：季節のせいか実力か、KPIモニタリングでの活用

この章の内容を実務で最も端的に活かせる場面が、「売上が下がったのは季節のせいか、それとも実力が落ちたのか」という、経営会議で頻繁に問われる質問への回答です。分解を行っていれば、この問いに対して次のような手順で答えられます。まず観測値の減少幅のうち、季節成分の変化で説明できる部分がどれだけあるかを確認します。次に、季節成分を除いた季節調整済み系列やトレンド成分が、実際にどちらの方向を向いているかを確認します。トレンド成分が横ばいまたは上向きを保っているのであれば、今回の落ち込みは季節要因による一時的なものである可能性が高いと説明できますし、トレンド成分自体が下向きに転じているのであれば、季節要因だけでは説明のつかない実力の低下が起きていると答えられます。前年同月比という1つの数字だけで答えていた頃と比べて、根拠を持って説明できる幅が格段に広がります。

KPIモニタリングの運用面では、原系列と季節調整済み系列を並べてダッシュボードに表示する運用が有効です。原系列は現場の実感（今月は繁忙期だから数字が良い、閑散期だから悪い）と一致しやすく、季節調整済み系列は経営判断に必要な「実力の趨勢」を映します。両方を並べて見ることで、現場の実感と経営判断の材料の両方を1つのダッシュボードでカバーできます。特に、月次の予実管理やアラート設計においては、原系列の対前月比だけで異常を検知しようとすると季節性による誤検知が多発するため、季節調整済み系列やトレンド成分の変化率を監視対象に加えることで、誤検知を減らしながら本当に注意すべき変化を捉えやすくなります。

時系列分解は、統計的な手続きである以前に、「今起きている変化のうち、毎年繰り返す季節要因で説明できる部分」と「それ以外の、注意を払うべき変化」を切り分けるための実務上の道具です。分解結果をレポートに添えることで、経営会議での議論を、感覚的な印象論から、成分ごとの根拠に基づいた議論へと引き上げることができます。

この章を深めたい方への参考書籍

『基礎からわかる時系列分析』（萩原淳一郎・瓜生真也・牧山幸史、技術評論社）：Rを用いた時系列分析の入門書で、古典的分解や状態空間モデルによる季節調整の考え方を、実データに近い例題を通じて丁寧に解説しています。本章で扱ったトレンドと季節性の分解を、異なる実装（Rのdecompose関数やstlに相当する処理）の視点から確認し、理解を補強したい場合に参考になります。

第3章 定常性と単位根、時系列分析の土台を固める

前章では、時系列データをトレンド・季節性・残差という3つの成分に分解し、それぞれの動きを可視化する方法を見てきました。分解によって系列の見た目の構造は整理できましたが、そこから先、ARIMAのような統計モデルを当てはめたり、複数の系列を回帰で結びつけたりするためには、もう1つ越えなければならない関門があります。それが本章のテーマである「定常性（stationarity）」です。定常性は地味な概念に見えますが、時系列分析の手法の大半は、この前提が成り立っていることを暗黙のうちに要求しています。前提を確認せずにモデルを当てはめると、統計的には有意に見える結果が実は無意味な偶然の産物だった、という事態を招きかねません。本章では定常性の定義から、代表的な非定常パターンであるランダムウォーク、単位根という考え方、そしてADF検定・KPSS検定による定常性の確認方法までを、実務で使える形で整理します。

定常性とは何か

時系列 \( y_t \) が弱定常（weakly stationary、または共分散定常）であるとは、次の3つの性質が時間によらず一定であることを指します。1つ目は平均 \( E[y_t] \) が時点 \( t \) によらず一定であること、2つ目は分散 \( \mathrm{Var}(y_t) \) が時点によらず一定であること、3つ目は自己共分散、つまり \( y_t \) と \( y_{t-k} \) の共分散が、時点 \( t \) そのものには依存せず、時点の間隔 \( k \) だけに依存することです。実務で単に「定常性」と言う場合、多くはこの弱定常性を指しています。

この性質が重視される理由は、統計的なモデル化という営みそのものの前提に関わっています。ARIMAにせよ回帰分析にせよ、モデルは「過去のデータから読み取ったパターンが、将来のデータにも同じように当てはまる」ことを暗黙の前提としています。しかし平均や分散が時点によって変わり続ける系列では、去年のデータから推定した平均や分散のパラメータが、来年のデータにそのまま当てはまる保証がありません。定常性は、いわば「過去の統計的な性質を将来に外挿してよい」という許可証のようなものであり、これが崩れている系列に通常の時系列モデルを適用すると、推定されたパラメータの意味自体が不安定になります。

非定常の典型パターン

非定常な系列にはいくつかの典型的な現れ方があります。実務のデータで頻繁に出会うパターンを整理します。

• トレンドのある系列：売上や会員数のように、時間とともに平均水準そのものが右肩上がり（または右肩下がり）に変化していく系列。平均が一定という条件を満たしません。
• 分散が拡大する系列：株価やアクセス数のように、水準が大きくなるほど変動の幅そのものも大きくなる系列。分散が一定という条件を満たしません。対数変換によって緩和できることが多いパターンです。
• ランダムウォーク：次節で詳しく扱う、前期の値にランダムなショックが加わり続けることで生じる非定常系列。見た目にはトレンドがあるように見えることもありますが、トレンドのある系列とは性質が本質的に異なります。

実務でこれらのパターンを見分ける第一歩は、まず系列をそのままプロットして目で確認することです。統計的な検定はあくまで判断を補強する材料であり、グラフを見ずに検定結果の数字だけで判断すると、季節性や構造変化といった検定では拾いにくい特徴を見落とすことがあります。

ランダムウォークを理解する

ランダムウォーク（random walk）は、非定常な時系列の中でも特に重要な性質を持つモデルです。定義は単純で、次の式で表されます。

\( y_t = y_{t-1} + \varepsilon_t \)

ここで \( \varepsilon_t \) は平均0、分散 \( \sigma^2 \) の互いに独立なノイズ（ホワイトノイズ）です。つまりランダムウォークとは、「今日の値は、昨日の値にランダムなショックを1つ加えただけのもの」という、極めて単純な仕組みで動く系列です。株価の日次終値がしばしばこのランダムウォークに近い動きを見せることは古くから知られており、市場価格には入手可能な情報が瞬時に織り込まれるという効率的市場仮説の議論とも関連の深いモデルです。ある日の株価が前日の株価にランダムな値動きを加えただけで決まるという考え方は、直感的にも受け入れやすいものでしょう。

ランダムウォークには「ドリフト付き（with drift）」と呼ばれる派生形もあります。式で表すと次のようになります。

\( y_t = c + y_{t-1} + \varepsilon_t \)

定数項 \( c \) が加わっただけの違いですが、この \( c \) が系列の見た目を大きく変えます。\( c \) がゼロでなければ、系列は毎期平均的に \( c \) だけ増加（または \( c \) が負なら減少）し続けるため、グラフに描くと右肩上がり（または右肩下がり）のトレンドがあるように見えます。ここで注意すべき点があります。この「トレンドのように見える動き」は、後の章で扱うような決定的トレンド（deterministic trend、時間 \( t \) の関数として書けるトレンド）とは本質的に異なるものだという点です。

ドリフト付きランダムウォークが持つトレンドは、確率的トレンド（stochastic trend）と呼ばれます。決定的トレンドの場合、トレンドの形は時間の関数としてあらかじめ決まっており、そこからのズレ（残差）だけが確率的に変動します。これに対しドリフト付きランダムウォークでは、過去に加わったすべてのショック \( \varepsilon_1, \varepsilon_2, \ldots, \varepsilon_t \) の累積そのものが系列の水準を決めています。つまりトレンドの経路自体がランダムに変化し続けるのです。グラフを1枚見ただけでは両者はよく似た右肩上がりの曲線に見えることがありますが、この違いは後述する定常化の方法（トレンド除去か、差分か）を左右する、実務上とても重要な区別になります。

ドリフト付きランダムウォークが見せる右肩上がりの動きは、決まった軌道からのランダムなブレではなく、過去のショックがすべて積み上がった結果としての「確率的トレンド」です。同じ右肩上がりのグラフでも、背後にある仕組みが決定的トレンドか確率的トレンドかによって、適切な定常化の方法（トレンド除去か差分か）が変わってきます。

次のコードは、ドリフトなしのランダムウォークと、平均0・分散一定の定常な系列（ホワイトノイズ）を同じ本数だけ生成し、見た目を比較する例です。乱数のシードを固定しているため、実行するたびに同じ系列が得られます。

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)
n = 300

# 定常な系列(ホワイトノイズ):平均・分散が時間によらず一定
stationary_series = rng.normal(loc=0, scale=1, size=n)

# ランダムウォーク:前期の値にノイズを加えた累積和
noise = rng.normal(loc=0, scale=1, size=n)
random_walk = np.cumsum(noise)

fig, axes = plt.subplots(2, 1, figsize=(9, 6), sharex=True)
axes[0].plot(stationary_series)
axes[0].set_title('定常な系列(ホワイトノイズ)')
axes[1].plot(random_walk)
axes[1].set_title('ランダムウォーク')
plt.tight_layout()
plt.show()

print(f'定常系列の前半・後半の平均: {stationary_series[:150].mean():.3f} / {stationary_series[150:].mean():.3f}')
print(f'ランダムウォークの前半・後半の平均: {random_walk[:150].mean():.3f} / {random_walk[150:].mean():.3f}')

ホワイトノイズは前半・後半とも平均が-0.051と-0.032で、0の近傍に留まります。これに対しランダムウォークは、前半7.583・後半4.853と、どちらの区間でも出発点の0から大きく離れた場所にいて、しかも区間によってその水準自体が変わります。定常な系列のように「平均はここ」と1つの数字で言い切れないことが、非定常な系列の特徴です。累積和という仕組み上、加わったショックが消えずに水準へ積み上がり続け、平均に引き戻される力がどこにも働いていないためです。この「引き戻す力の有無」こそが、次に説明する単位根の考え方の核心になります。

単位根という概念

単位根（unit root）は、ランダムウォークを一般化した言い方だと理解しておくと扱いやすくなります。時系列を1期前の自分自身に回帰する次のようなモデルを考えます。

\( y_t = \phi y_{t-1} + \varepsilon_t \)

この係数 \( \phi \) がちょうど1のとき、この式はランダムウォークそのものになります。この \( \phi = 1 \) という状態を指して「単位根を持つ」と表現します。一方 \( \phi \) の絶対値が1より小さければ、系列は平均へと回帰する力を持つ定常過程になります。単位根を持つかどうかを見分けることが、実務上なぜそれほど重要なのかを、直感的な言葉で言い換えると次のようになります。単位根を持つ系列では、ある時点で加わったショックの影響が、時間が経ってもまったく減衰せずに将来までずっと残り続けます（ショックの永続性）。これに対し単位根を持たない定常な系列では、ショックの影響は時間とともに徐々に薄れていき、系列は元の平均水準へと回帰していきます。

このショックの永続性という性質は、ビジネスの意思決定に直結する含意を持ちます。たとえば、あるキャンペーンの失敗や外部環境の悪化によって売上が一時的に落ち込んだとします。この売上の系列が単位根を持つ（非定常な）系列であれば、落ち込んだ水準そのものが新しい基準点になってしまい、何か新しい施策を打たない限り、自然に元の水準へ戻ることは期待できません。逆に売上の系列が定常であれば、ショックの影響は時間とともに減衰し、平均的な水準へと自然に回帰していくと期待できます。「一度落ち込んだ数字は放っておいても戻るのか、それとも戻らないのか」という、経営判断に直結する問いに対して、単位根の有無は1つの手がかりを与えてくれます。

単位根検定：ADF検定とKPSS検定

単位根を持つかどうかをデータから判断するための代表的な統計的検定が、拡張ディッキー・フラー検定（Augmented Dickey-Fuller test、ADF検定）です。ADF検定を使ううえで、実務担当者が必ず押さえておくべき注意点が1つあります。それは、ADF検定の帰無仮説が「単位根が存在する（系列は非定常である）」という側に置かれている点です。つまりADF検定は、「非定常である」ことを出発点として、それを否定できるだけの十分な証拠がデータにあるかどうかを調べる検定です。検定の結果、p値が有意水準（たとえば0.05）を下回れば帰無仮説を棄却でき、「単位根はない、系列は定常である」と判断します。逆にp値が大きければ、帰無仮説を棄却できないというだけであり、これは「定常であることが証明された」わけではなく、「非定常でないとは言い切れない」という消極的な結論に留まる点に注意が必要です。

この帰無仮説の向きを補う目的でよく併用されるのが、KPSS検定（Kwiatkowski-Phillips-Schmidt-Shin test）です。KPSS検定の帰無仮説は、ADF検定とはちょうど逆向きに設定されています。KPSS検定の帰無仮説は「系列は定常である」という側に置かれており、p値が有意水準を下回れば帰無仮説を棄却し、「定常ではない（非定常である）」と判断します。この帰無仮説の向きの違いを取り違えると、検定結果の解釈がまるきり逆になってしまうため、実務で扱う際には特に注意が必要です。

2つの検定の帰無仮説が逆向きであることを逆手にとって、実務では両方の検定を組み合わせて結論の確からしさを高める割り切った使い方がよく採られます。ADF検定で単位根が棄却され（定常と判断）、かつKPSS検定でも定常が棄却されない（定常と判断）場合は、両検定の結論が一致しているため、定常であるという判断にある程度の自信を持てます。逆にADF検定で単位根が棄却されず、KPSS検定でも定常が棄却される場合は、両検定とも非定常を支持していることになります。一方で2つの検定の結論が食い違う場合は、系列がトレンド定常（後述）に近い性質を持っている可能性や、サンプルサイズが不十分である可能性が考えられ、どちらか一方の検定結果だけを鵜呑みにせず、実際のグラフや業務知識と照らし合わせて総合的に判断する姿勢が求められます。

次のコードは、statsmodelsのadfullerとkpssを使って、前節で生成したランダムウォークと定常な系列それぞれに対して単位根検定を実行し、結果を解釈する例です。

from statsmodels.tsa.stattools import adfuller, kpss

def report_unit_root_tests(series, name):
    adf_stat, adf_p, _, _, _, _ = adfuller(series, autolag='AIC')
    kpss_stat, kpss_p, _, _ = kpss(series, regression='c', nlags='auto')
    print(f'--- {name} ---')
    print(f'ADF検定: 統計量={adf_stat:.3f}, p値={adf_p:.3f}')
    print(f'KPSS検定: 統計量={kpss_stat:.3f}, p値={kpss_p:.3f}')

report_unit_root_tests(stationary_series, '定常な系列(ホワイトノイズ)')
report_unit_root_tests(random_walk, 'ランダムウォーク')

定常な系列に対しては、ADF検定のp値が有意水準を大きく下回り単位根が棄却される一方、KPSS検定のp値は有意水準を上回り定常性が棄却されない、という両検定が一致した結果になるのが典型的です。これに対しランダムウォークに対しては、ADF検定のp値が有意水準を上回り単位根を棄却できず、KPSS検定のp値は有意水準を下回り定常性が棄却される、という非定常を支持する結果になるのが典型的です。なお、kpss関数を実行するとサンプルサイズによっては警告（InterpolationWarning）が表示されることがありますが、これはp値が用意されている参照表の範囲を超えていることを知らせる警告であり、検定結果の解釈自体を妨げるものではありません。

定常化の技法

系列が非定常だと判断された場合、そのままの形でARMAモデルや通常の回帰に投入することはできません。定常性を満たす形に変換してからモデル化するのが基本的な流れになります。なお第6章で扱うARIMAは、この変換のうち差分をモデルの中に組み込んだもので、差分の回数を次数として指定すれば非定常な系列をそのまま渡せます。代表的な定常化の技法を整理します。

• 差分（1階差分）：\( \Delta y_t = y_t – y_{t-1} \) を取ることで、単位根を持つ系列（確率的トレンドを持つ系列）を定常化する、最も基本的な方法です。ランダムウォークの1階差分は、定義からそのままホワイトノイズになります。
• 季節差分：月次データで前年同月との差分 \( y_t – y_{t-12} \) を取るなど、周期を1つ前と比較する差分です。季節性に由来する非定常性を取り除く際に使います。
• 対数変換：分散が水準とともに拡大していく系列に対して、対数を取ることで変動の幅を安定させる方法です。差分と組み合わせて対数差分（近似的に前期比の変化率）として使われることも多くあります。
• トレンド除去（デトレンド）：時間 \( t \) を説明変数とする回帰を当てはめ、その予測値を元の系列から差し引く方法です。決定的トレンドを持つ系列に対して用います。

差分を取る際に注意したいのが、必要以上に差分を繰り返してしまう「過差分（over-differencing）」の弊害です。すでに1階差分で定常になっている系列にさらに2階差分を取ってしまうと、系列は定常のままですが、本来存在しなかった不必要な自己相関構造を人為的に作り出してしまいます。技術的には、過差分を行った系列ではMA多項式の根が単位円上に乗ってしまう（反転可能性の条件を満たさなくなる）ことが知られていますが、実務上は「差分を取りすぎると、モデルが表現すべきでない見かけの構造が生まれ、推定が不安定になったり予測の幅がかえって広がったりする」という結論を押さえておけば十分です。差分は「定常性を得るために必要な最小限の階数」に留めるのが原則であり、何回差分を取ったら定常になったかを、単位根検定とグラフの両方で確認しながら1段階ずつ進めることが推奨されます。

次のコードは、前節のドリフトなしランダムウォークに1階差分を適用し、差分を取る前と後でADF検定の結果がどう変わるかを確認する例です。

import pandas as pd
from statsmodels.tsa.stattools import adfuller

rw_series = pd.Series(random_walk)
rw_diff = rw_series.diff().dropna()

adf_before = adfuller(rw_series, autolag='AIC')
adf_after = adfuller(rw_diff, autolag='AIC')

print(f'差分前のADF検定 p値: {adf_before[1]:.3f}')
print(f'1階差分後のADF検定 p値: {adf_after[1]:.3f}')

差分を取る前はp値が有意水準を上回り単位根を棄却できませんが、1階差分を取った後はp値が大きく下がり、単位根が棄却され定常であると判断できる水準になります。これはランダムウォークの1階差分が定義上ホワイトノイズになるという性質と整合しており、単位根検定が想定どおりに機能していることの確認にもなります。実務のデータでも、多くの場合1階差分、季節性が強い場合はこれに季節差分を組み合わせる程度で定常性が得られることが多く、2階差分以上が必要になるケースはそれほど多くありません。

トレンド定常と差分定常の違い

非定常な系列は、大きく2種類に分けて考えることができます。1つはトレンド定常（trend stationary）な系列で、決定的トレンドを取り除けば定常になる系列です。この場合の適切な定常化の方法は、時間を説明変数とした回帰でトレンド成分を除去するデトレンドです。もう1つは差分定常（difference stationary）な系列で、単位根を持ち、確率的トレンドによって非定常になっている系列です。この場合の適切な定常化の方法は、デトレンドではなく差分を取ることになります。

この2つを取り違えると、定常化がうまくいきません。差分定常な系列（ランダムウォークなど）に対してデトレンドだけを行っても、単位根そのものは残ったままのため非定常性は解消されません。逆にトレンド定常な系列に対して差分を取ると、過差分と同様の弊害が生じ、本来不要な自己相関構造を作り出してしまいます。したがって非定常な系列に出会ったときは、それがトレンド定常なのか差分定常なのかを、単位根検定とグラフの両方から見極めたうえで、デトレンドか差分かという対処法を選ぶことが重要になります。

次章への接続

本章では、1本の時系列がそれ自体として定常か非定常かを見極める方法を扱いました。しかし実務では、複数の時系列の間の関係を調べたい場面が数多くあります。売上と広告費、為替レートと輸出額といった具合です。ここで注意しなければならないのが、非定常な系列同士を、定常性を確認しないまま単純に回帰分析にかけてしまうと何が起こるか、という問題です。本来まったく無関係な2つの非定常系列同士を回帰しても、決定係数が非常に高く、回帰係数の統計的有意性を示すt値も大きくなるという、見かけ上は極めて良好な回帰結果が得られてしまうことがあります。これは「見せかけの回帰（spurious regression）」と呼ばれる現象であり、次章ではこの現象がなぜ起こるのか、そしてこれを避けるための共和分という考え方を扱います。

この章を深めたい方への参考書籍

『計量経済学』（西山慶彦・新谷元嗣・川口大司・奥井亮、有斐閣）：定常性・単位根・ADF検定やKPSS検定の理論的な背景を、経済データへの応用例とともに体系的に解説しており、本章で直感的に説明した内容を数式レベルで裏付けたい読者に適しています。単位根検定の検出力や過差分の問題についても踏み込んで扱われています。

第4章 見せかけの回帰を見破り、共和分を活かす

前章では、時系列データの定常性と単位根という概念を扱いました。ある系列の平均や分散、自己相関の構造が時間によらず一定であれば定常、そうでなければ非定常であり、ランダムウォークのように「1つ前の値に乱数の増分を加えていく」だけの系列は単位根を持つ非定常過程の代表例だという整理をしています。本章では、この非定常性が回帰分析にどのような罠を仕掛けるかを本格的に掘り下げます。

回帰編でも触れたとおり、無関係な2つの非定常系列を回帰にかけると、決定係数が高く、係数も統計的に有意という結果が出てしまう現象があります。これを見せかけの回帰（spurious regression）と呼びます。回帰編ではこの現象への入口として簡単に触れましたが、本章ではなぜこれが起きるのかというメカニズムを掘り下げたうえで、非定常同士でも本物の長期的な関係が存在する場合がある共和分という概念、そしてその関係を実務でどう検証し、どう活用するかまでを扱います。

無関係な系列同士に高い決定係数が出る割合を数える

見せかけの回帰がどの程度の規模で起きる現象なのかを、シミュレーションで確認します。互いに完全に独立な乱数から作った2本のランダムウォークを用意し、一方をもう一方に回帰にかけるという操作を、乱数シードを変えながら1,000回繰り返し、p値が0.05を下回る（統計的に有意と判定される）割合を数えます。本来、両者の間には何の関係もないため、有意と判定される割合は理論上の水準である5%程度に収まるはずです。

import numpy as np
import statsmodels.api as sm

def random_walk(n, rng):
    # 乱数の累積和でランダムウォークを生成する
    return np.cumsum(rng.normal(0, 1, n))

def spurious_rate(n, n_trials=1000, seed_offset=0):
    count_significant = 0
    r2_list = []
    for trial in range(n_trials):
        x = random_walk(n, np.random.default_rng(trial + seed_offset))
        y = random_walk(n, np.random.default_rng(trial + seed_offset + 50000))
        model = sm.OLS(y, sm.add_constant(x)).fit()
        r2_list.append(model.rsquared)
        if model.pvalues[1] < 0.05:
            count_significant += 1
    return count_significant / n_trials, np.mean(r2_list)

for n in [50, 200, 1000]:
    rate, r2_mean = spurious_rate(n)
    print(f"サンプルサイズn={n}: 有意判定の割合={rate:.2f}, 平均R2={r2_mean:.3f}")

このコードを実行すると、有意と判定される割合はn=50で68%、n=200で83%、n=1,000で92%となり、いずれも本来あるべき5%を大幅に超えます。互いに何の関係もない2本の乱数列を回帰にかけているだけであるにもかかわらず、です。平均R2も0.23から0.25と、無関係な系列同士とは思えない水準になります（個別の試行では0.9近くに達するものもあります）。さらに注目すべきは、サンプルサイズnを50から1,000まで増やすと、この割合が改善するどころか68%から92%へ単調に悪化する点です。通常、統計的な検定は標本サイズを増やすほど誤判定が減っていくものですが、見せかけの回帰ではこの直感が通用しません。データを増やせば増やすほど、誤った有意判定に自信を持ってしまう危険すらあるということです。

見せかけの回帰はなぜ起きるのか：平均回帰しない単位根系列の性質

この現象の根っこにあるのは、単位根を持つ系列が平均回帰（mean reversion）しないという性質です。平均回帰する系列であれば、値が平均から離れるほど元の水準へ引き戻す力が働くため、長期間にわたって一方向に動き続けることはありません。ところがランダムウォークのような単位根過程には、この引き戻す力がありません。ある時点で受けたショック（ランダムな増分）は消えることなく、その後もずっと系列の水準に足し込まれ続けます。この性質のため、単位根系列は「たまたま」数十期にわたって上昇し続けたり下降し続けたりする局面を生みやすく、無関係な2本の系列であっても、たまたま似た期間を切り取ると、片方が上がるときにもう片方も上がって見えるという偶然の一致が起こりやすくなります。

この偶然の一致を、通常の回帰分析は本物の関係と区別する仕組みを持っていません。むしろ問題はもう一段深いところにあります。無関係な非定常系列同士を回帰にかけると、その残差（誤差項）自体も非定常になってしまうのです。通常の最小二乗法におけるt検定やp値の計算は、残差が定常であり、かつ観測値どうしが実質的に独立な情報を持つという前提のもとで組み立てられています。残差が非定常になるということは、この前提が根底から崩れているということであり、t値は本来従うはずの分布には収束せず、標本サイズを増やしても発散していく方向に振れることが知られています。先ほどのシミュレーションでサンプルサイズを増やしても有意判定の割合が改善しなかったのは、この技術的な事情を反映した結果です。

見せかけの回帰は、統計ソフトの計算間違いでも、たまたま運が悪かった一度きりの偶然でもありません。単位根を持つ系列同士を回帰にかけると、残差が非定常になり、t検定やp値が前提としている条件そのものが崩れてしまうという、構造的に繰り返し起こる現象です。データの点数を増やしても解決しない（むしろ悪化しうる）という点が、通常の統計的な誤りとの大きな違いです。

右肩上がりの指標同士は相関の解釈に注意する

実務で扱う指標の多くは、トレンドを持つ非定常な系列です。事業が成長を続けていれば自社の月次売上は右肩上がりになりますし、経済全体が拡大している局面ではGDPのような集計指標も右肩上がりになります。この2つを並べて回帰にかければ、高い決定係数と有意な係数が出ることが多くなりますが、それが「GDPが伸びたから自社の売上が伸びた」という因果関係を意味するとは限りません。両者がたまたま同じ時期に右肩上がりだったという、トレンドどうしの偶然の一致である可能性を常に疑う必要があります。

同様の注意は、SNSのフォロワー数と株価、Webサイトの累積訪問者数と売上高、従業員数と経費総額など、時間とともに一貫して増え続ける指標どうしの組み合わせ全般に当てはまります。こうした指標はいずれも非定常であることが多く、単純に並べて相関係数や回帰係数を計算すると、実態以上に強い関係があるかのような数字が出やすくなります。経営会議やレポートでこの種の相関を目にしたときは、両方の系列がトレンドを持つ非定常な系列ではないか、その関係は差分を取っても残るのかを、まず確認する姿勢が実務上の防衛線になります。

BIツールのダッシュボードが自動で算出する相関係数にも、同じ注意が当てはまります。ダッシュボード上の数字は「統計的に正しい手続きを経て出てきた値」であるとは限らず、単に2つの右肩上がりの指標を並べただけで高い相関係数が表示されているだけ、という場合が少なくありません。相関係数や回帰係数が経営判断の根拠として使われる場面ほど、その背後にある系列が非定常かどうかを事前に確認しておく価値があります。

見せかけの回帰を見破る4つの確認項目

非定常な系列同士の回帰結果を見せかけの回帰と疑うべきかどうかは、いくつかの診断項目を順番に確認することで、ある程度機械的に見破ることができます。以下では回帰分析の基本的な出力を使います。決定係数（R2）はモデルがデータのばらつきをどれだけ説明できているかを0から1で表す指標、t値は推定した係数が0からどれだけ離れているかを標準誤差で割った値、p値はその係数が偶然0から離れて見える確率です。いずれも回帰編で扱った内容ですが、ここでは「値が大きいほど関係が強く見える」という程度の理解があれば読み進められます。

• 回帰の残差に対して単位根検定（ADF検定など、前章で扱った手法）を行い、残差に単位根が残っていないか
• ダービン・ワトソン統計量（DW統計量）が0に極端に近い、非常に小さい値になっていないか。DW統計量は回帰の残差に自己相関が残っていないかを測る指標で、値はおおむね0から4の範囲を取り、2の付近であれば自己相関はほとんどないと判断します。0に近いほど残差が前の時点を強く引きずっている状態を意味します
• 決定係数（R2）がDW統計量の値そのものより大きくなっていないか（グレンジャー・ニューボルドの目安として知られる経験則で、見せかけの回帰の可能性を示唆する）
• 2つの系列を差分系列（前期との差）に変換してから回帰し直したとき、有意だった関係が消えてしまわないか

この4項目のうち、最初の2つと3つ目は、回帰を1回実行した後の残差やDW統計量を確認するだけで済む簡便なチェックです。とりわけR2がDW統計量を上回っているという状態は、見せかけの回帰の古典的な兆候として知られています。最後の1項目である「差分を取ると関係が消えるか」は、次に説明する対処の基本そのものでもあり、最も実務的な判断材料になります。

差分系列で回帰するという対処とその代償

見せかけの回帰への最も基本的な対処は、水準（値そのもの）ではなく、差分系列（前期との差）に変換してから回帰することです。前章で扱ったとおり、ランダムウォークの増分は本来独立な乱数そのものであるため、差分を取った時点で「無関係な系列」という実態が正しく数値に反映されるようになります。無関係な2本のランダムウォークについて、水準の回帰と差分の回帰を比較すると、水準の回帰では試行によってR2が0.9近くまで達し有意なp値も出やすいのに対し、差分の回帰では平均R2が0.02以下、有意と判定される割合も5%程度と、名目の有意水準どおりに戻ります。

ただし、この対処には見過ごせない代償があります。差分を取るという操作は、系列が持っていた水準の情報、つまり「長期的にどのくらいの関係で動くか」という情報をそぎ落としてしまうのです。仮に2つの系列の間に本物の長期的な均衡関係が存在していたとしても、差分系列だけを見ていると、その長期的な関係はデータから消え去り、短期的なノイズどうしの関係しか見えなくなります。差分を取れば見せかけの回帰は避けられますが、それと引き換えに、本物の長期的な関係を分析する手段まで手放してしまうことになりかねません。この代償を払わずに、非定常な系列どうしの本物の関係を扱う方法が、次に説明する共和分という考え方です。

共和分：非定常同士でも差が定常なら長期的な関係がある

共和分（cointegration）とは、それぞれ単独では非定常な系列であっても、それらを適切な比率で組み合わせた線形結合（差や比率）が定常になる場合、その関係を指す概念です。2つの系列が共和分の関係にあるとき、両者は短期的にはそれぞれ好き勝手にふらふらと動くように見えても、長期的には一定の関係を保ち続けるという性質を持ちます。

この関係をイメージするうえでよく使われるのが、酔っ払いとその飼い犬の散歩というたとえ話です。酔っ払いは足取りがおぼつかず、次にどちらへ歩き出すか予測がつきません。犬もまた、あちこちの匂いにつられて自由気ままに動きます。それぞれの歩みだけを見れば、どちらも「次にどこへ向かうか分からない」非定常な動きです。しかし2人はリードでつながれているため、2人の距離が際限なく開いていくことはなく、ある範囲内に収まり続けます。それぞれの位置（水準）は非定常でも、2人の距離（線形結合）は定常だというのが、共和分の本質です。

ビジネスの文脈でも、こうした関係は珍しくありません。ガソリン価格と灯油価格は、どちらも原油価格という共通の要因に連動して動くため、それぞれの水準は時間とともに大きく変動する非定常な系列でありながら、両者の価格差はある範囲に収まり続ける傾向があります。近隣に立地する2店舗の売上高も、地域全体の経済状況という共通の要因を受けて動くため、それぞれの売上高は非定常でも、両者の差や比率は比較的安定した範囲で推移することがあります。こうした「共通の要因に連動して動く」関係にある系列のペアは、共和分の候補として検討する価値があります。

共和分検定：Engle-Granger法の手順とJohansen検定

2つの系列が本当に共和分の関係にあるかどうかを確かめる代表的な手法が、Engle-Granger法です。手順は次の2段階です。まず、一方の系列をもう一方に回帰し、通常の最小二乗法で係数を推定します。次に、その回帰から得られた残差に対して単位根検定を行い、残差が定常である（単位根を持たない）と判定されれば、両者は共和分の関係にあると結論づけます。残差が非定常なままであれば、両者の間に長期的な均衡関係はなく、見せかけの回帰だったということになります。

ここで注意が必要なのは、回帰残差に対する単位根検定は、通常のADF検定の臨界値をそのまま使えないという点です。残差はすでに推定されたパラメータから計算された量であるため、通常のADF検定より厳しい（より大きなマイナス側の値が必要な）臨界値を使う必要があります。この調整済みの臨界値はEngle-GrangerあるいはMacKinnonの臨界値と呼ばれ、statsmodelsのcoint関数を使えば、この調整を自動的に行ったうえで共和分の有無を検定できます。

import numpy as np
import statsmodels.api as sm
from statsmodels.tsa.stattools import coint, adfuller

rng = np.random.default_rng(42)
n = 300

# 共通の確率的トレンドを持つ、本物の共和分関係にある2系列を作る
common_trend = np.cumsum(rng.normal(0, 1, n))
y = common_trend + rng.normal(0, 0.5, n)        # 共通トレンド+定常な誤差
x = 0.7 * common_trend + rng.normal(0, 0.5, n)  # 同じ共通トレンドの影響を受ける

# 互いに無関係な2系列(見せかけの回帰の例、共和分なし)
y_indep = np.cumsum(rng.normal(0, 1, n))
x_indep = np.cumsum(rng.normal(0, 1, n))

score1, pvalue1, _ = coint(y, x)
score2, pvalue2, _ = coint(y_indep, x_indep)

print(f"共通トレンドを持つペア: coint検定p値={pvalue1:.4f}")
print(f"無関係な2系列のペア: coint検定p値={pvalue2:.4f}")

# Engle-Granger法の手順を手動でなぞる場合
model = sm.OLS(y, sm.add_constant(x)).fit()
resid = model.resid
adf_stat, adf_p, _, _, crit, _ = adfuller(resid)
print(f"回帰残差のADF統計量={adf_stat:.3f}, p値(参考値)={adf_p:.4f}")

実行すると、共通のトレンドを持つように作った1組目のペアではcoint検定のp値が小さく（共和分ありと判定されやすく）、互いに完全に無関係な2組目のペアではp値が大きく（共和分なしと判定されやすい）なる傾向が確認できます。なお、手動でadfullerを残差にかけた場合のp値は、通常のADF検定用の臨界値をもとに計算されているため参考値にとどまり、正式な判定にはcoint関数が返す（MacKinnonの調整済み臨界値に基づく）p値を使うのが適切です。

2系列の組み合わせであればEngle-Granger法で十分ですが、3つ以上の系列の間で共和分関係を検討したい場合や、複数の共和分関係が同時に存在しうる場合には、Johansen検定と呼ばれる、より一般化された手法が使われます。Engle-Granger法は「どちらの系列を説明変数に、どちらを目的変数にするか」で結果が変わりうるという弱点を持ちますが、Johansen検定はこの非対称性を解消し、複数の系列をまとめてベクトル自己回帰（VAR）モデルとして扱ったうえで、その枠組みの中で共和分関係の本数まで推定できる点が大きな違いです。VARモデルの枠組みそのものは第7章で扱いますが、Johansen検定の具体的な手順は専門書に譲り、本ガイドでは考え方の紹介にとどめます。実務では、まず2系列でEngle-Granger法による簡易な当たりをつけ、3系列以上の分析が必要になった段階でJohansen検定に進むという順序で十分な場合がほとんどです。

誤差修正モデル：長期均衡からの乖離が短期の調整を生む

2つの系列が共和分の関係にあると分かった場合、その関係を活かして予測や分析に使う枠組みが誤差修正モデル（Error Correction Model、ECM）です。ECMの背後にある発想はシンプルです。長期的には2つの系列が一定の均衡関係を保つとしても、短期的にはその均衡から外れることがあります。ECMは、この「前期に均衡からどれだけ乖離していたか」という情報を、当期の変化量を説明する要因の1つとして組み込みます。

具体的には、yの当期の変化量（差分）を、xの当期の変化量と、前期における長期均衡からの乖離（共和分の回帰から得られる残差の1期前の値）の2つで説明します。数式で表すと、おおよそ\(\Delta y_t = \beta \Delta x_t + \gamma (y_{t-1} – \theta x_{t-1}) + \varepsilon_t\)という形になり、\(\gamma\)が乖離解消の速さを表す係数です。長期均衡からの乖離を表す項の係数\(\gamma\)は、通常マイナスの値を取ります。これは「前期に均衡より上振れしていれば、当期は下落方向に調整され、逆に下振れしていれば上昇方向に調整される」という、均衡へ戻ろうとする力を表しています。この係数の絶対値が大きいほど、乖離が速く解消される（調整速度が速い）ことを意味し、絶対値が小さければ、均衡から外れた状態が長く続きやすいことを意味します。

import numpy as np
import statsmodels.api as sm

# 前セルのcommon_trendから作ったy, xが共和分の関係にあるとして、
# 長期均衡関係(共和分回帰)とECMを推定する
long_run = sm.OLS(y, sm.add_constant(x)).fit()
ecm_term = long_run.resid[:-1]  # 前期の均衡からの乖離(1期ラグ)

dy = np.diff(y)
dx = np.diff(x)

X_ecm = sm.add_constant(np.column_stack([dx, ecm_term]))
model_ecm = sm.OLS(dy, X_ecm).fit()

print(f"短期の関係(dxの係数): {model_ecm.params[1]:.3f}")
print(f"調整速度(乖離解消の係数): {model_ecm.params[2]:.3f}")

このコードを実行すると、乖離解消の係数はマイナスの値になり、前期に生じた均衡からのずれが、当期にある程度解消される方向に働いていることが確認できます。ECMの強みは、短期の変化（差分）と長期の均衡関係（水準）の両方の情報を、1つのモデルの中に同時に取り込める点にあります。差分だけを見る回帰では失われてしまう長期的な関係の情報を、乖離解消の項というかたちで取り戻せるというのが、共和分を確認したうえでECMを使う実務上の利点です。

非定常×非定常の回帰を見たら、まず「差分を取るべきか、共和分の関係を疑うべきか」を問うことが実務の第一歩です。差分だけを見て終わらせると、本物の長期的な関係を見落とすかもしれません。逆に共和分を確認せずに水準のまま回帰を続けると、見せかけの回帰を本物の関係と誤認するリスクを抱えたままになります。

本章では、見せかけの回帰が起きるメカニズムを、単位根系列の平均回帰しないという性質と、残差の非定常性という技術的な事情の両面から掘り下げ、それを見破るための診断項目、そして非定常同士でも本物の長期的な関係を扱う共和分とECMという枠組みまでを整理しました。次章では、こうした非定常性の議論から離れ、時系列データ自身の過去の値をどう使ってモデル化するか、AR（自己回帰）モデルとMA（移動平均）モデルという、時系列予測の基本的な構成要素を扱います。

この章を深めたい方への参考書籍

『経済・ファイナンスデータの計量時系列分析』（沖本竜義、朝倉書店）は、単位根、見せかけの回帰、共和分、誤差修正モデルという本章の中核テーマを、経済・金融データを題材に数理的な背景とともに体系立てて解説した一冊です。本章で扱ったEngle-Granger法やECMの導出をより厳密に追いたい読者、また続く第7章のVARモデルの理解を深めておきたい読者に向いています。

第5章 過去から未来を語るARモデルとMAモデル

前章までで、時系列データが持つ自己相関という性質（第1章）、トレンドと季節性を取り除く分解の考え方（第2章）、そして定常性という土台となる概念と、その土台が崩れたときに起きる見せかけの回帰や単位根の問題（第3章、第4章）を確認してきました。ここからは、いよいよ時系列データそのものを使って未来の値を語るための、最も基本的なモデル群に入ります。本章で扱うのは、ARモデル（自己回帰モデル）とMAモデル（移動平均モデル）、そしてその組み合わせであるARMAモデルです。

これらのモデルは、ニューラルネットワークを使った複雑な予測手法が広まった現在でも、時系列分析の基礎体力として欠かせない存在です。理由は単純で、複雑なモデルの多くも、根っこの部分ではARやMAが捉えている「過去の値や過去の誤差が現在にどう影響するか」という発想を土台にしているためです。この土台を理解しないまま高度な手法に進むと、モデルが出す予測の意味を説明できなくなってしまいます。

過去の値で未来の値を説明するという発想

通常の回帰分析では、売上を説明するために広告費や気温といった、目的変数とは別の変数を用意します。ところが時系列予測の多くの場面では、そうした説明変数を用意しなくても、対象となる系列そのものの過去の値だけを手がかりに、ある程度筋の通った予測を組み立てることができます。これが時系列モデリングの最も基本的な発想であり、「過去の自分」を使って「未来の自分」を説明するという考え方です。

この発想が成り立つのは、時系列データの多くが自己相関、すなわち近い時点同士が似た値を取りやすいという性質を持っているためです。今月の売上は先月の売上と無関係な数字ではなく、何らかの慣性や余韻を引き継いでいます。ARモデルとMAモデルは、この慣性や余韻の伝わり方を、それぞれ異なる角度から定式化したものだと理解しておくと、以降の内容が整理しやすくなります。

ARモデルは「値そのものの慣性」に着目し、MAモデルは「誤差（ショック）の余韻」に着目します。この違いを念頭に置きながら、それぞれのモデルを見ていきます。

ARモデル：前期の値をどれだけ引き継ぐか

自己回帰モデル（AutoRegressive Model、ARモデル）は、現在の値を、過去の自分自身の値の一次結合として説明するモデルです。最も基本的な形であるAR（1）モデルは、次の式で表されます。

\( y_t = c + \phi y_{t-1} + \varepsilon_t \)

ここで\( y_t \)は時点\( t \)における観測値、\( c \)は定数項、\( \phi \)（ファイと読みます）は自己回帰係数、\( \varepsilon_t \)は平均0の白色雑音（過去の情報とは無関係な誤差項）です。この式が意味しているのは、「今日の値は、昨日の値の\( \phi \)倍に、いくらかのランダムな変動が加わったものである」ということです。

\( \phi \)は、前期の値のうち何割を今期に引き継ぐかを表す係数だとイメージすると直感的です。たとえば\( \phi = 0.8 \)であれば、前期からの水準の8割程度が今期に持ち越され、\( \phi = 0.2 \)であれば2割程度しか持ち越されません。\( \phi \)が大きいほど、系列は前期の値に強く引きずられ、変化がゆっくりとしたなだらかな動きになります。逆に\( \phi \)が小さいほど、系列は前期の値への依存が弱くなり、平均的な水準へ素早く戻ろうとする動きになります。

この「平均的な水準へ戻ろうとする速さ」は平均回帰（mean reversion）と呼ばれ、\( \phi \)の大きさによって決まります。AR（1）モデルの理論上の平均は\( c / (1 – \phi) \)で与えられ、系列は何らかのショックで一時的にこの平均から離れても、時間が経つにつれてこの水準に戻っていきます。\( \phi \)が0に近いほど平均への回帰は速く、1に近いほど回帰は遅く、ショックの影響が長く尾を引きます。

ここで重要になるのが、\( \phi \)の値とモデルの定常性の関係です。AR（1）モデルが定常であるための条件は\( |\phi| < 1 \)、すなわち\( \phi \)の絶対値が1未満であることです。この条件が満たされる限り、ショックの影響は時間とともに指数関数的に減衰し、系列は一定の平均と分散を持つ状態に落ち着きます。

一方で\( \phi = 1 \)になると、この式はそのまま第3章で扱ったランダムウォークと一致します。\( y_t = y_{t-1} + \varepsilon_t \)という形になり、前期の値をそのまま100%引き継ぎ、平均に回帰する力が働かなくなります。第3章で確認したとおり、この状態は単位根を持つ非定常過程であり、ショックの影響は減衰せずに系列に半永久的に残り続けます。AR（1）モデルは、\( \phi \)を1未満に保つことで定常性を確保しているモデルであり、\( \phi = 1 \)という境界線の向こう側が、まさに単位根の世界だという理解が、第3章と本章をつなぐ橋渡しになります。

実務上は、AR（1）よりも高次のAR（p）モデル、すなわち\( p \)期前までの値を使って現在の値を説明するモデルもよく使われます。式は次のように一般化されます。

\( y_t = c + \phi_1 y_{t-1} + \phi_2 y_{t-2} + \cdots + \phi_p y_{t-p} + \varepsilon_t \)

次数\( p \)が大きいほど、より遠い過去の情報まで現在の値の説明に使うことになりますが、パラメータ数が増える分、推定の不安定さや過学習のリスクも高まります。次数をどう決めるかは、後述するACF・PACFやAICを使った判断が関わってきます。

MAモデル：過去のショックが有限期間だけ残る

移動平均モデル（Moving Average Model、MAモデル）は、ARモデルとは異なる角度から時系列を説明します。MAモデルは、現在の値を、現在および過去の誤差項（ショック）の一次結合として表現します。最も基本的なMA（1）モデルは、次の式で表されます。

\( y_t = \mu + \varepsilon_t + \theta \varepsilon_{t-1} \)

ここで\( \mu \)は系列の平均、\( \varepsilon_t \)は現時点の白色雑音、\( \varepsilon_{t-1} \)は1期前の白色雑音、\( \theta \)（シータと読みます）はその余韻の大きさを表す移動平均係数です。この式が意味しているのは、「今日の値は、平均を中心に、今日起きたランダムなショックと、昨日起きたショックの余韻を\( \theta \)倍だけ引きずったものである」ということです。

MA（1）モデルの直感で押さえておきたいのは、ショックの影響がちょうど1期分だけ残り、2期先には完全に消えてなくなるという点です。より一般化したMA（q）モデルでは、ショックの影響は\( q \)期分だけ残り、\( q + 1 \)期以降にはきっぱりと消えます。これは、時間とともに影響が指数関数的に減衰していくものの、理論上は永遠にゼロにならないARモデルの慣性とは、性質が明確に異なります。

この違いを整理すると、ARモデルの慣性は「じわじわと薄まりながらも消えきらない余韻」であるのに対し、MAモデルの慣性は「決まった期間だけ完全に残り、その後はきっぱり消える余韻」だといえます。両者はともに「過去が現在に影響する」という点では似ていますが、影響の残り方の形が根本的に異なるモデルです。

ARMAモデル：ARとMAを組み合わせて倹約的に表現する

ARモデルとMAモデルは、それぞれ単独でも使われますが、両者を組み合わせたARMAモデル（AutoRegressive Moving Average Model）もよく使われます。ARMA（p, q）モデルは、次数\( p \)の自己回帰項と、次数\( q \)の移動平均項をあわせ持つモデルで、式は次のようになります。

\( y_t = c + \phi_1 y_{t-1} + \cdots + \phi_p y_{t-p} + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \cdots + \theta_q \varepsilon_{t-q} \)

ARMAモデルが使われる理由は、少ないパラメータ数で複雑な自己相関の構造を表現できる、いわゆる倹約性（パーシモニー）にあります。ある系列の挙動を、次数の高いARモデルだけで無理に近似しようとすると、多くのパラメータが必要になることがありますが、AR項とMA項をあわせ持つARMAモデルであれば、より低い次数の組み合わせで同等の当てはまりを実現できる場合が少なくありません。パラメータ数を抑えることは、推定の安定性を高め、モデルの過学習を防ぐことにもつながります。

次数\( p \)と\( q \)をどう選ぶかは、実務上の重要な論点です。次節では、この次数選びの古典的な手がかりであるACF（自己相関関数）とPACF（偏自己相関関数）の使い方を整理します。

ACFとPACFによる次数の見極め

ARモデルとMAモデル、それぞれの次数を見極める古典的な手がかりとして使われてきたのが、ACF（自己相関関数）とPACF（偏自己相関関数）です。ACFは、系列とそのラグ（何期前の値か）との単純な相関を表す指標で、PACFは、間のラグの影響を除いたうえでの、そのラグ固有の相関を表す指標です。

この2つの指標には、モデルの種類によって特徴的な現れ方があることが知られており、次数の見当をつける古典的な型として整理されています。

この表を整理すると、「PACFが特定のラグでスパッと切れていればARモデルらしく、その切れた場所が次数\( p \)の候補になる」「ACFが特定のラグでスパッと切れていればMAモデルらしく、その切れた場所が次数\( q \)の候補になる」という古典的な読み方が成り立ちます。ARMAのようにAR項とMA項の両方を含む場合は、どちらの指標もはっきりとは打ち切られず、緩やかに減衰する形になるため、この読み方だけでは正確な次数を一意に特定するのは難しくなります。

実務では、ACFとPACFのグラフを目視で確認することは、系列の大まかな性質（強い自己回帰的な慣性を持つのか、限られた期間だけ余韻が残るショック型の動きなのか）をつかむうえで今でも有用な最初のステップです。ただし、現実のデータではACFやPACFの減衰パターンが教科書どおりにきれいに現れることは少なく、目視だけで次数を一意に決め切るのは難しい場合がほとんどです。そこで実務では、ACF・PACFで大まかな見当をつけたうえで、AIC（赤池情報量規準）などの情報量規準を使い、複数の次数の組み合わせを自動的に試して最も当てはまりの良いものを選ぶという、機械的な自動選択を併用するのが一般的です。ACF・PACFによる目視の型は「見立てのための道具」、AICによる自動選択は「絞り込みのための道具」と役割を分けて理解しておくと、実務での使い分けがしやすくなります。

なお、モデルのパラメータ（\( \phi \)や\( \theta \)）そのものの推定方法については、本章では立ち入りません。かつてはユール・ウォーカー方程式のような解析的な手法が使われていましたが、現在はstatsmodelsをはじめとする多くの統計ソフトウェアで、最尤推定（尤度を最大にするパラメータを数値的に探す手法）が標準的に用いられています。実務でモデルを当てはめる際には、この推定の内部処理を意識する必要はほとんどなく、ライブラリに任せてしまって差し支えありません。

残差が白色雑音になっているかで妥当性を確認する

ARモデルやARMAモデルを当てはめた後、次数の選び方が適切だったかどうかを確認する作業が欠かせません。その基本的な考え方は、「モデルが系列の持つ自己相関の構造をきちんと捉えられていれば、当てはめた後に残る残差には、もはや自己相関が残っていないはずだ」というものです。残差に自己相関が残っているということは、モデルがまだ捉えきれていない構造がデータの中に残っている、というサインになります。

この確認によく使われるのが、残差のACFプロットと、Ljung-Box検定（リュング・ボックス検定）です。残差のACFプロットは、残差を系列とみなしてACFを計算し、各ラグの自己相関がほぼ0の周辺に収まっているか（統計的に有意な自己相関が残っていないか）を視覚的に確認するものです。

Ljung-Box検定は、この確認を統計的な仮説検定の形で行うもので、「残差にはあるラグまでの範囲で自己相関が存在しない」という帰無仮説を検定します。p値が小さい（一般的な有意水準である0.05を下回るなど）場合は、残差にまだ有意な自己相関が残っている、すなわちモデルの次数や構造がデータに対して不十分である可能性を示唆します。逆にp値が大きい場合は、残差が白色雑音に近い状態にあると判断し、モデルの当てはまりに大きな問題はないと考えることができます。

ここで注意しておきたいのは、Ljung-Box検定で残差に有意な自己相関が見当たらないからといって、そのモデルが唯一の正解だと確定するわけではないという点です。あくまで「まだ捉えられていない自己相関という明白な不備は見当たらない」ことを確認しているに過ぎず、他の次数の組み合わせでも同様に問題のない結果が得られることは珍しくありません。残差診断は、モデルを積極的に選び取るための道具というよりも、明らかに不十分なモデルを除外するための道具だと位置づけておくと、過信を防ぎやすくなります。

慣性型の現象とショック残存型の現象

ARモデルとMAモデルの違いは、抽象的な数式としてだけでなく、ビジネス上の現象を見立てる視点としても役立ちます。ある現象がどちらの性質に近いかを考えることは、モデルを選ぶ前の仮説作りとして有用です。

ARモデルの発想に近いのは、値そのものに慣性が働く現象です。たとえば顧客数の推移は、多くの場合AR的な性質を持ちます。今月の顧客数は、来店や利用の習慣を持つ既存顧客の存在によって、先月の顧客数からある程度の割合を引き継ぎます。急激な離反やキャンペーンがない限り、顧客数は前月の水準に強く引きずられながら、緩やかに増減していく傾向があります。この「前月の何割を引き継ぐか」という発想は、そのままAR（1）モデルの\( \phi \)に対応します。定数項がある場合、厳密に引き継がれるのは値そのものではなく、系列の平均からのずれです。また\( \phi \)が負の値のときは、平均より高かった翌期は低くなるという形で、符号を反転させながら平均へ戻っていきます。

一方、MAモデルの発想に近いのは、特定のイベントによる影響が一定期間だけ残り、その後きっぱりと消える現象です。たとえば「効果が翌月まで残るプロモーション」は、MA的な性質を持つ典型例です。ある月に大型のキャンペーンを実施すると、その月だけでなく翌月まで購買が押し上げられる効果が観察されることがありますが、キャンペーンの記憶や在庫の押し上げ効果は、多くの場合いつまでも続くわけではなく、数か月を過ぎればほぼ解消されます。この「効果が何か月分だけ残るか」という発想が、MAモデルの次数\( q \)に対応する考え方です。

ある現象をモデル化する前に、「これは前期の水準そのものが慣性的に引き継がれる現象か（AR的）、それとも特定のショックの影響が一定期間だけ残って消える現象か（MA的）」という問いを立ててみることは、次数選びの機械的な手続きに入る前の、実務的に有用な仮説作りになります。現実の多くの現象は両方の性質を併せ持つため、最終的にはARMAモデルとしてデータに語らせる場面が多くなりますが、この問いを立てておくこと自体が、モデルの結果を解釈するときの土台になります。

Pythonで確認する

ここまでの内容を、Pythonで確認します。まず、AR（1）モデルの\( \phi \)を変えたときに、系列の挙動がどう変わるかをシミュレーションで比較します。

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(0)
n = 200
phis = [0.2, 0.6, 0.95]  # 引き継ぐ割合を変えて比較

fig, axes = plt.subplots(len(phis), 1, figsize=(8, 7), sharex=True)
for ax, phi in zip(axes, phis):
    eps = rng.normal(0, 1, n)
    y = np.zeros(n)
    for t in range(1, n):
        y[t] = phi * y[t - 1] + eps[t]
    ax.plot(y)
    ax.set_title(f"AR(1)シミュレーション: phi={phi}")
    ax.axhline(0, color="gray", linewidth=0.8)
plt.tight_layout()
plt.show()

このコードを実行すると、\( \phi = 0.2 \)の系列は0の周辺で細かく上下しながらすぐに水準を戻すのに対し、\( \phi = 0.95 \)の系列はいったん水準がずれると、なかなか0に戻らず、なだらかな山や谷を描きながら長く漂う様子が確認できます。\( \phi \)を1に近づけるほど、系列の見た目がランダムウォークに近づいていく様子も合わせて観察できます。

次に、ACFとPACFのプロットを使って、次数の見当をつける手順を確認します。

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# y: 分析対象の時系列(1次元配列またはpandas.Series)。手元のデータに置き換えて使うテンプレートです
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
plot_acf(y, lags=20, ax=axes[0])
axes[0].set_title("ACF")
plot_pacf(y, lags=20, method="ywm", ax=axes[1])
axes[1].set_title("PACF")
plt.tight_layout()
plt.show()

PACFが特定のラグを境にほぼ0の範囲（信頼区間の帯）に収まっていれば、そのラグまでがAR項の次数\( p \)の候補になります。逆にACFが特定のラグで打ち切られていれば、そのラグまでがMA項の次数\( q \)の候補です。両方が緩やかに減衰している場合は、ARMAモデルとしてAR項とMA項の両方を含めることを検討します。

最後に、statsmodelsのARIMAクラスを使ってARMAモデルを当てはめ、残差診断まで一通り確認します。ARMAモデルはARIMAモデルの特殊ケース（差分の次数を0としたもの）にあたるため、statsmodelsでは共通のARIMAクラスを使って推定します。

from statsmodels.tsa.arima.model import ARIMA
from statsmodels.stats.diagnostic import acorr_ljungbox

# series: 分析対象の時系列(pandas.Series、日次や月次のインデックスを想定)。
# 手元のデータに置き換えて使うテンプレートです
# order=(p, d, q): pはAR次数、dは差分の次数(定常な系列なら0)、qはMA次数
model = ARIMA(series, order=(1, 0, 1))
result = model.fit()
print(result.summary())

resid = result.resid

# 残差のACFを目視で確認
plot_acf(resid, lags=20)
plt.title("残差のACF")
plt.show()

# Ljung-Box検定で残差に自己相関が残っていないか確認
lb_test = acorr_ljungbox(resid, lags=[10, 20], return_df=True)
print(lb_test)

result.summary()の出力では、AR項とMA項それぞれの係数、標準誤差、p値に加えて、AIC・BIC（ベイズ情報量規準。AICと同じくモデルの当てはまりと複雑さのバランスを測る指標で、AICよりもパラメータ数の多さに厳しい）といった情報量規準の値も確認できます。複数の次数の組み合わせでモデルを当てはめ、AICが最も小さい組み合わせを比較検討することで、ACF・PACFの目視による見立てを、数値的な基準で補強できます。Ljung-Box検定の結果は、lb_stat列とlb_pvalue列を確認し、lb_pvalueが0.05を下回るラグがなければ、残差に明白な自己相関は残っていないと判断できます。もし小さいp値のラグが見つかった場合は、次数を見直すか、季節性など他の構造が残っていないかを確認する必要があります。季節性を伴う系列への拡張は、次章のSARIMAモデルで扱います。

この章を深めたい方への参考書籍

『時系列解析 自己回帰型モデル・状態空間モデル・異常検知』（島田直希、共立出版）：自己回帰型モデルから状態空間モデル、異常検知までを一貫した枠組みで扱っており、本章で最小限に留めたAR・MAモデルの推定理論（最尤推定の数理的な背景）を、より厳密に掘り下げたい読者に向いています。次章以降で扱うARIMAや状態空間モデルへの橋渡しとしても参照しやすい一冊です。

第6章 ARIMAとSARIMAで実務の予測を組み立てる

前章では、時系列データを自己回帰（AR）モデル、移動平均（MA）モデル、そして両者を組み合わせたARMAモデルという枠組みで捉える方法を確認しました。ARMAモデルが前提としているのは、対象のデータがすでに定常であるという条件です。しかし第3章で見たとおり、実務で扱う売上や在庫、アクセス数といったデータの多くはトレンドや季節性を含んでおり、そのままでは定常性の条件を満たしません。

この章で扱うARIMAモデルは、第3章で学んだ「差分を取って定常化する」という発想と、前章で学んだARMAモデルを1つの手続きに統合したものです。さらに、季節性を持つデータに対応させたSARIMA、外生変数を組み込んだSARIMAXへと拡張し、予測区間の考え方、バックテストによる検証、そしてARIMA系と並ぶもう1つの定番であるExponential Smoothing（指数平滑法）との使い分けまでを、月次売上の合成データを使った一連の実務例を通して確認します。

ARIMAへの統合：差分（I）がARとMAをつなぐ

ARIMAは、AutoRegressive Integrated Moving Averageの略で、名前が示すとおり3つの要素から構成されています。1つ目はAR（自己回帰）、2つ目はI（和分、Integrated）、3つ目はMA（移動平均）です。このうちAR成分とMA成分は前章で扱った内容と同じで、過去の値そのものを使って現在を説明するのがAR、過去の予測誤差（ショック）を使って現在を説明するのがMAでした。

ARIMAが加えるI成分は、モデルを推定する前にデータを何回差分するかという操作を表します。第3章で確認したとおり、単位根を持つデータ（非定常なデータの代表的なケース）は、1階差分を取ることで定常なデータに変換できることが多くあります。ARIMAは、この「差分を取って定常化してからARMAを当てはめる」という2段階の作業を、1つのモデルの中に組み込んだものだと理解すると見通しがよくなります。差分を何回取るか、ARの次数をいくつにするか、MAの次数をいくつにするかという3つの数字を並べて、ARIMA（p, d, q）と表記します。pがAR次数、dが差分の回数、qがMA次数です。d=0のときのARIMA（p, 0, q）は、前章のARMA（p, q）そのものになります。ARIMAはARMAを包含する、より広い枠組みだと捉えることができます。

（p, d, q）の決め方：dは単位根検定、p・qはAICで

ARIMAを実務で使ううえで最初につまずきやすいのが、（p, d, q）という3つの数字をどう決めるかという点です。経験的には、この3つを別々の手順で決めるのが分かりやすい進め方です。

dの決定には、第3章で扱った単位根検定（ADF検定など）を使います。元データに対してADF検定を行い、帰無仮説（単位根がある、すなわち非定常）が棄却されなければ1階差分を取り、差分後のデータに対して再度ADF検定を行うという手順を、帰無仮説が棄却されるまで繰り返します。多くの実務データでは1階差分、季節性が強い場合や成長率が変化し続けるデータではまれに2階差分まで必要になりますが、3階以上の差分が必要になるケースは稀で、差分を重ねすぎるとかえってノイズを増幅させてしまう点には注意が必要です。差分後のデータの自己相関（ACF）を確認し、急激な減衰が見られるかどうかも、差分回数が適切かどうかの補助的な判断材料になります。

dが決まったら、p（AR次数）とq（MA次数）を決めます。前章ではACF・PACFの形状からARかMAかを見分ける考え方を扱いましたが、実データではAR成分とMA成分が混在していることが多く、目視だけで正確な次数を絞り込むのは簡単ではありません。実務でよく使われるのは、複数の（p, q）の組み合わせでモデルを推定し、赤池情報量規準（AIC）が最小になる組み合わせを選ぶという総当たりの方法です。AICは、モデルの当てはまりの良さと複雑さ（パラメータ数）のバランスを評価する指標で、値が小さいほど良いモデルとされます。次のコードは、pとqの候補範囲を総当たりし、AICが小さい順に並べる例です。

# y_train は学習用の時系列(手元のデータに置き換えてください)
import itertools
import warnings
import pandas as pd
from statsmodels.tsa.arima.model import ARIMA

warnings.filterwarnings("ignore")

def grid_search_arima(y, p_range, d_range, q_range):
    records = []
    for p, d, q in itertools.product(p_range, d_range, q_range):
        try:
            fitted = ARIMA(y, order=(p, d, q)).fit()
            records.append({
                "order": (p, d, q),
                "aic": fitted.aic,
                "bic": fitted.bic,
            })
        except Exception:
            # 収束しない組み合わせは無視して次に進む
            continue
    return pd.DataFrame(records).sort_values("aic").reset_index(drop=True)

# dは単位根検定で先に決めておき、ここでは固定してp・qだけを探索する
d_fixed = 1
order_table = grid_search_arima(
    y_train, p_range=range(0, 4), d_range=[d_fixed], q_range=range(0, 4)
)
print(order_table.head(10))

ここでdをAICの探索対象に含めていない点には理由があります。差分回数dが異なるモデルどうしのAICは、それぞれ別の変換を施したデータに対して尤度を計算しているため、直接比較してはいけません。dは単位根検定で先に決めて固定し、AICで探索するのはpとqだけにする、という順序を守る必要があります。

この総当たりはAICのみを機械的に比較する方法であり、次数を決める作業のすべてを代替するものではありません。AICが最小の組み合わせが必ずしも実務上の最良解とは限らず、上位数件の候補について残差の自己相関（Ljung-Box検定）やパラメータの有意性もあわせて確認したうえで、最終的な次数を選ぶという姿勢が欠かせません。

SARIMA：季節差分と（P, D, Q, s）

ARIMAは、トレンドを持つデータには対応できますが、季節性、たとえば毎年12月に売上が伸びる、毎週末にアクセス数が増えるといった周期的な変動には、そのままでは十分に対応できません。この課題に対応するのがSARIMA（Seasonal ARIMA）です。SARIMAは、通常の（p, d, q）に加えて、季節成分の（P, D, Q, s）という4つのパラメータを持ちます。sは季節周期の長さで、月次データで年周期を扱う場合はs=12、四半期データならs=4です。P、D、Qはそれぞれ季節AR次数、季節差分の回数、季節MA次数で、通常の（p, d, q）が1期前・2期前といった直近のラグとの関係を捉えるのに対し、（P, D, Q）はsか月前、2sか月前といった、周期を1つ分・2つ分さかのぼった関係を捉えます。表記はSARIMA（p, d, q）（P, D, Q, s）となります。

季節差分は、時点tの値から、s期前の値を差し引く操作です。通常の1階差分が「前日・前月との差」を取るのに対し、季節差分は「1年前の同じ月との差」を取るイメージで、これによって毎年繰り返す季節パターンを除去できます。実務では、通常の差分（d）と季節差分（D）を両方使うケースが多く、たとえばSARIMA（1, 1, 1）（1, 1, 1, 12）であれば、トレンド成分を通常の1階差分で、年次の季節成分を季節差分で、それぞれ定常化したうえでARMA構造を当てはめるモデルになります。

データの周期をどう扱うかは、集計単位によって変わります。月次データであればs=12が定番ですが、週次データで年周期を扱う場合はs=52、日次データで週周期を扱う場合はs=7です。なお1年は正確には約52.18週であり、sは整数でしか指定できないため、週次データの年周期はs=52による近似になります。この端数を厳密に扱いたい場合は、フーリエ項を外生変数として与える方法や、第8章で扱う状態空間モデルを検討します。日次データに年周期（s=365）を持たせるSARIMAは、季節ラグの数が非常に大きくなり計算負荷も推定の安定性も厳しくなるため、実務では日次データの年周期は第8章・第9章で扱う状態空間モデルや、第10章の機械学習ベースの手法に任せることが多く、SARIMAは月次・週次・四半期といった、季節周期がそれほど大きくないデータに向いていると考えておくとよいでしょう。

SARIMAX：プロモーションや祝日を説明変数に組み込む

実務の売上データには、価格改定やプロモーション、祝日といった、時系列モデルの自己回帰構造だけでは説明しきれない要因が数多く影響します。回帰編第12章では、時系列データに通常の回帰分析をそのまま当てはめると、残差に自己相関が残ってしまい、標準誤差や有意性の判断が歪むという問題を扱いました。SARIMAXは、この問題への1つの回答です。SARIMAの枠組みに外生変数（Exogenous variables）を追加し、時系列としての自己回帰構造（SARIMA部分）と、価格やプロモーションといった説明変数の効果（回帰部分）を同じモデルの中で同時に推定します。

SARIMAXに組み込む代表的な外生変数には、プロモーション実施の有無を示す0/1のフラグ、価格や割引率といった連続値、祝日や大型連休を示すダミー変数などがあります。これらの変数は、モデルの自己回帰構造では説明できない「その時点固有の事情」を捉える役割を果たします。将来の予測を行う際には、予測対象の期間についても外生変数の値が必要になる点に注意してください。プロモーション予定や祝日カレンダーのように事前に判明している変数であれば問題ありませんが、価格のように将来値が確定していない変数を使う場合は、想定シナリオごとに予測を出す、あるいは価格自体を別途予測するといった工夫が必要になります。

点予測だけでなく予測区間で語る

統計学基礎編で扱った区間推定の考え方は、時系列の予測にもそのまま引き継がれます。SARIMA・SARIMAXが出力するのは、単一の予測値（点予測）だけではなく、その予測値がどの範囲に収まりそうかを示す予測区間です。経営会議で「来月の売上は120としか言えません」と報告するよりも、「来月の売上は120で、95%予測区間は105から135です」と報告するほうが、意思決定者にとって不確実性の大きさを踏まえた判断がしやすくなります。

ARIMA系モデルの予測区間には、先の期間になるほど区間が広がるという構造的な特徴があります。1か月先の予測は直近の実績に近い情報をもとに計算できるため誤差が小さく収まりますが、12か月先の予測になると、途中の各時点での誤差が積み重なっていくため、区間は大きく広がります。この「遠い将来ほど自信を持って言い切れなくなる」という性質は、予測モデルの欠陥ではなく、時系列予測が本質的に持つ性質だと理解しておく必要があります。逆に、遠い先の予測がやけに狭い区間で示されているモデルを見かけたときは、モデルの前提（誤差の分散が一定であるという仮定など）が実データに合っていない可能性を疑う材料になります。

時系列予測の報告では、点予測の数字だけを1つ示すのではなく、予測区間もあわせて示すことで、どこまでが確度の高い予測でどこからが不確実性の大きい領域なのかを、意思決定者と共有できます。特に先の期間の予測ほど区間が広がるという構造そのものを説明に含めることで、数字を過信するリスクを減らせます。

バックテスト：ホールドアウトとローリング予測

ARIMA・SARIMAX・後述するETSのいずれを使う場合も、モデルを選んだ根拠を、学習データへの当てはまりの良さ（AICなど）だけで済ませず、実際に予測させてみて検証するバックテストの手続きが欠かせません。もっとも基本的な方法はホールドアウト検証で、手元のデータの末尾の一定期間（たとえば直近12か月）を検証用に取り分け、それより前のデータだけでモデルを学習させ、取り分けた期間を予測させて実績値と突き合わせます。

ホールドアウト検証は手軽ですが、検証に使える期間が1回分しかないため、たまたまその期間の実績が特殊だった場合に評価が歪みやすいという弱点があります。これを補うのがローリング予測（逐次再学習）で、学習期間を少しずつ後ろにずらしながら、複数回にわたって「学習してn期先を予測する」という手続きを繰り返し、そのつど生じる予測誤差を積み上げて評価します。ローリング予測は、単発のホールドアウト検証よりも季節や景気循環の異なる複数のタイミングでの予測精度を確認できるため、より頑健な評価が可能になりますが、モデルを繰り返し学習させる計算コストが増える点はトレードオフになります。ローリング予測の具体的な設計（検証窓の幅や再学習の頻度をどう決めるか）は、需要予測の実務設計を扱う第12章で詳しく取り上げます。

Exponential Smoothing（指数平滑法）とARIMAの使い分け

ARIMA系モデルと並ぶ、時系列予測のもう1つの定番がExponential Smoothing（指数平滑法、ETS）です。ETSの基本的な発想は単純で、直近の実績ほど大きな重みを、過去にさかのぼるほど小さな（指数的に減衰する）重みを与えて平均を取り、その加重平均を将来の予測値とするというものです。ARIMAが「差分による定常化とAR・MA構造」という視点からデータを捉えるのに対し、ETSは「レベル（水準）・トレンド・季節性という3つの成分に分解し、それぞれを指数平滑で更新する」という視点からデータを捉えます。

ETSファミリーの中でも実務でよく使われるのがHolt-Winters法です。Holt法がレベルとトレンドの2つの成分を指数平滑で追跡するのに対し、Holt-Winters法はさらに季節成分を加え、レベル・トレンド・季節性の3つを同時に更新しながら予測します。季節成分の入れ方には加法（季節変動の大きさが水準によらずほぼ一定）と乗法（水準が大きいほど季節変動の振れ幅も大きくなる）の2種類があり、月次売上のように、売上規模が大きくなるほど季節ごとの変動幅も大きくなるデータでは乗法モデルが向いていることが多いとされます。

ARIMA（SARIMAX含む）とETSのどちらを使うべきかについては、一方が常に優れているという単純な優劣関係はありません。ETSはパラメータの数が少なく、モデルの解釈（レベル・トレンド・季節性の分解）が直感的で、データ量が少なくても比較的安定して動く傾向があります。一方でETSは、外生変数を組み込む標準的な枠組みを持たず（状態空間モデルとして拡張する方法はありますが）、プロモーションや価格の効果を明示的に扱いたい場面ではSARIMAXに分があります。ARIMA系はACF・PACFに基づく次数選択や外生変数の追加など、モデルの構造をより細かく作り込める柔軟性を持ちますが、そのぶんパラメータ数が増え、データ量が少ないと過学習気味になることもあります。実務での結論は、どちらが理論的に優れているかを先に決めるのではなく、両方を実際にバックテストにかけて予測精度を比較し、対象のデータで実際に精度が高い方を採用するという進め方になります。

auto-arimaによる自動化とその注意点

（p, d, q）や季節次数の探索を自動化するツールとして、pmdarimaライブラリのauto_arima関数がよく使われます。auto_arimaは、単位根検定によるdの自動判定、AIC（またはBIC）を基準にしたp・qの段階的な探索、季節次数の探索までを1つの関数呼び出しでまとめて実行してくれる、実務上非常に便利なツールです。特に、対象とする時系列の数が数十から数百に及び、1本ずつ手作業で次数を検討する時間が取れないような場面では、auto_arimaによる一次選定は有効な出発点になります。

ただし、auto_arimaが選んだ次数をそのまま最終結果として採用することには注意が必要です。auto_arimaはAIC・BICという統計的な基準だけをもとに探索するため、実務上の制約（たとえば季節性の周期は業務カレンダー上12か月で固定したい、プロモーションの効果は必ず正の係数になるはずといった知識）までは反映されません。また探索範囲を広く取りすぎると、見かけ上AICが最も良い次数が、検証データに対しては必ずしも良い予測を出さないという事態も起こり得ます。auto_arimaはあくまで候補を効率よく絞り込むための一次選定のツールと位置づけ、選ばれた次数についても、前節までで扱った残差診断とバックテストによる検証を必ず行うという姿勢が欠かせません。

月次売上の合成データでSARIMAXを当てはめる

ここまでの内容を、月次売上を模した合成データを使って一連の流れとして確認します。トレンド、年次の季節性、プロモーションによる押し上げ効果を持つ60か月分のデータを作り、直近12か月をホールドアウト検証用に取り分けます。

ここで、プロモーションの実施月を「毎年6月・11月・12月」のように固定していない点には理由があります。外生変数が季節周期と完全に一致していると、季節差分（D=1、s=12）を取った時点でその変数は恒等的に0になり、係数を推定する手がかりが消えてしまいます。実務でも、毎年必ず同じ月に同じ施策を打っている場合、その効果は季節成分と区別できません。外生変数として意味を持つのは、季節パターンから外れた動きを含む変数だという点は、SARIMAXを設計するうえでの重要な前提です。

import numpy as np
import pandas as pd
from statsmodels.tsa.statespace.sarimax import SARIMAX

rng = np.random.default_rng(0)
n_months = 60
date_index = pd.date_range("2020-01-01", periods=n_months, freq="MS")

# トレンド + 年次季節性 + プロモーション効果 + ノイズで合成売上データを作成
trend = np.linspace(100, 160, n_months)
seasonal = 15 * np.sin(2 * np.pi * date_index.month / 12)
# プロモーションは不定期に実施される(毎年同じ月に固定すると季節成分と区別できない)
promo_flag = rng.binomial(1, 0.3, n_months)
noise = rng.normal(0, 5, n_months)
sales = trend + seasonal + 20 * promo_flag + noise

df = pd.DataFrame({"sales": sales, "promo": promo_flag}, index=date_index)
train, test = df.iloc[:-12], df.iloc[-12:]

model = SARIMAX(
    train["sales"],
    exog=train[["promo"]],
    order=(1, 1, 1),
    seasonal_order=(1, 1, 1, 12),
    enforce_stationarity=False,
    enforce_invertibility=False,
)
result = model.fit(disp=False)
print(result.summary())

# 残差の自己相関をLjung-Box検定で確認(モデルが情報を取り切れているかの診断)
# lagsを省略すると観測数に応じた少ないラグまでしか見ない。季節周期12まで指定する
lb = result.test_serial_correlation(method="ljungbox", lags=12)
print(f"Ljung-Box(ラグ12) 統計量={lb[0, 0, -1]:.3f}, p値={lb[0, 1, -1]:.3f}")

# 12か月先の予測と95%予測区間
forecast = result.get_forecast(steps=12, exog=test[["promo"]])
pred_mean = forecast.predicted_mean
conf_int = forecast.conf_int(alpha=0.05)
print(pd.concat([pred_mean.rename("予測値"), conf_int], axis=1))

result.summaryでは各次数の係数とp値、AIC・BICが確認できます。この例ではpromo変数の係数が19.38（p値0.000）と推定され、データ生成時に設定した押し上げ幅20をほぼ回収できています。プロモーションの効果がモデルに適切に取り込めているかどうかは、この係数の符号と大きさで判断します。test_serial_correlationによるLjung-Box検定は、ラグ12まで見て統計量9.719、p値0.641となり、残差に自己相関は残っていないと判断できます。なおlags引数を省略した場合に何ラグまで検定するかは、有効な観測数から自動的に決まります。この例では4ラグまでしか見てくれないため、季節周期を持つモデルでは、周期までのラグを自分で明示する必要があります。get_forecastによる予測では、predicted_meanが点予測、conf_intが予測区間に相当します。予測区間の幅は1か月先の24.0から12か月先の25.0へと少しずつ広がっていきます。この例で広がり方が緩やかなのは、推定されたMA項が1階差分の効果をほぼ打ち消しているためで、系列やモデルの次数によっては、区間はもっと急速に広がります。

最後に、同じホールドアウト期間に対してETS（Holt-Winters法）を当てはめ、SARIMAXとの予測精度を比較します。

from statsmodels.tsa.holtwinters import ExponentialSmoothing
from sklearn.metrics import mean_absolute_percentage_error

ets_model = ExponentialSmoothing(
    train["sales"],
    trend="add",
    seasonal="add",
    seasonal_periods=12,
).fit()
ets_forecast = ets_model.forecast(12)

sarimax_mape = mean_absolute_percentage_error(test["sales"], pred_mean)
ets_mape = mean_absolute_percentage_error(test["sales"], ets_forecast)

print(f"SARIMAX(外生変数あり)のMAPE: {sarimax_mape:.3f}")
print(f"ETS(Holt-Winters)のMAPE: {ets_mape:.3f}")

実行すると、SARIMAXのMAPE（平均絶対パーセント誤差）が0.028、ETSが0.058となり、SARIMAXが優位でした。このデータではプロモーションが売上を不定期に押し上げる構造を組み込んで生成しているため、その情報を外生変数として明示的に使えるSARIMAXが有利になります。裏を返せば、ETSの誤差の多くはプロモーション実施月を予測しきれないことから生じており、外生変数が本当に効く場面ではこれだけの差がつくということでもあります。一方でETSは、外生変数を持たずトレンドと季節性だけで平均5.8%の誤差に収めており、外生変数の情報が入手しづらい、あるいはモデルをシンプルに保ちたい場面での有力な選択肢であることも読み取れます。実務では、この章で見た項目、次数選択の根拠、残差診断、予測区間の妥当性、バックテストでの精度比較を一通りそろえたうえで、SARIMAXとETSのどちらを採用するかを判断する流れが基本になります。

この章を深めたい方への参考書籍

『Forecasting: Principles and Practice（第3版）』（Rob J. Hyndman、George Athanasopoulos、Web公開・英語）は、ETS（指数平滑法）とARIMA系モデルの双方を体系的にカバーし、本章で扱った次数選択やバックテストの考え方をより詳しく確認できる定番書です。原著者自身がforecastパッケージ・fableパッケージの開発者でもあり、理論と実装の橋渡しが丁寧な点が特徴で、無料でWeb公開されているため、章末に出てきた用語をさらに深掘りしたい読者には最初にあたる1冊としておすすめできます。

第7章 複数の時系列の関係を読むVARモデルとGranger因果

1本の系列だけでは答えられない、系列間の関係

第5章・第6章で扱ってきたAR・MA・ARIMA・SARIMAは、いずれも1本の時系列を、その系列自身の過去の値や過去の誤差から説明しようとするモデルでした。売上という1本の系列があれば、過去の売上の動きから将来の売上を予測する、という枠組みです。この考え方は、季節性やトレンドを含む単一の指標の先行きを見通す場面では強力ですが、実務で直面する問いの多くは、実は「1本の系列の中で完結しない」形をしています。

広告費を増やした翌週、翌月の売上はどう動くのか。店頭価格を変更したとき、販売数量への影響は何週間続くのか。Webサイトへの流入数の増減は、実店舗への来客数にどれだけ先行して現れるのか。これらはいずれも、2つ以上の系列の間の時間的な関係を問うものであり、片方の系列だけを見ていては答えが出ません。広告費だけをARIMAで分析しても売上への波及は見えず、売上だけを分析しても、その背後に広告費の動きがあったのかは分かりません。

本章で扱うVAR（Vector Autoregression、ベクトル自己回帰）モデルは、こうした「複数の系列が互いにどう影響し合っているか」を捉えるための枠組みです。あわせて、VARの推定結果を使って「ある系列の過去の情報を加えると、別の系列の予測が統計的に改善するか」を検定するGranger因果という考え方、そして「ある系列に一時的なショックが起きたとき、他の系列がその後どう反応するか」を可視化するインパルス応答関数を扱います。これらはいずれも、広告と売上、価格と数量、Web流入と来店数といった、複数の指標を横断してビジネスの因果構造を探るための実務的な道具です。

VARモデル: 全変数の過去で全変数を説明する連立の自己回帰

VARモデルの発想は、AR（自己回帰）モデルの考え方をそのまま複数の系列に拡張したものだと捉えると理解しやすくなります。ARモデルでは、ある系列の現在の値を、その系列自身の過去の値だけで説明しました。VARモデルでは、対象とする複数の系列それぞれの現在の値を、自分自身の過去の値に加えて、他のすべての系列の過去の値でも説明します。つまり、変数の数だけ回帰式を用意し、それぞれの式の右辺に全変数の過去の値を並べた、連立方程式の形を取ります。

具体的なイメージをつかむために、広告費\( x_t \)と売上\( y_t \)という2つの系列を対象にした、最も単純なVAR（1）（ラグ1期のVAR）を考えます。この場合、モデルは次の2本の式で構成されます。

\( y_t = c_1 + a_{11} y_{t-1} + a_{12} x_{t-1} + e_{1,t} \)
\( x_t = c_2 + a_{21} y_{t-1} + a_{22} x_{t-1} + e_{2,t} \)

1本目の式は、今期の売上\( y_t \)を、前期の売上\( y_{t-1} \)と前期の広告費\( x_{t-1} \)の両方で説明しています。2本目の式は逆に、今期の広告費\( x_t \)を、前期の広告費\( x_{t-1} \)と前期の売上\( y_{t-1} \)の両方で説明しています。それぞれの式は見た目こそ単なる重回帰式ですが、両方の式が同時に成り立つ連立の体系として推定される点がVARの特徴です。係数\( a_{12} \)が統計的に意味のある大きさを持てば、「前期の広告費が今期の売上に影響を与えている」という関係を数値として捉えられますし、\( a_{21} \)に意味があれば、逆に「前期の売上の水準が今期の広告費の決め方に影響している」（たとえば好調な月の翌月に広告予算を積み増す、といった意思決定のパターン）を捉えられます。

変数の数が3つ、4つと増えても考え方は同じで、それぞれの変数を目的変数とする回帰式を1本ずつ用意し、右辺にはすべての変数の過去の値（ラグ）を並べます。ラグの次数を\( p \)期まで広げたものをVAR（\( p \)）と呼びます。個々の式だけを見ると通常の重回帰と変わりませんが、変数どうしが互いの過去を説明変数として参照し合う「相互作用のネットワーク」を1つのモデルとして推定する点が、単変量のARIMAにはない発想です。

VARの前提と手順: 定常性の確認とラグ次数の選択

VARモデルを実際に当てはめる前に、いくつかの手順を踏む必要があります。最初の関門は、対象とする各系列が定常であることの確認です。第3章で扱った通り、VARのような自己回帰の枠組みは、系列が定常である（平均・分散・自己相関構造が時間を通じて安定している）ことを前提にしており、非定常な系列をそのままVARに投入すると、係数の推定量や有意性検定の結果が信頼できなくなります。実務では、拡張ディッキー・フラー検定（ADF検定）などを用いて各系列の単位根の有無を確認し、非定常であれば階差を取って定常化してからVARを組むのが基本の対応です。

ただし、ここで第4章の内容が関わってきます。複数の非定常な系列どうしが、それぞれの階差では捉えきれない長期的な均衡関係（共和分関係）で結ばれている場合、単純にすべての系列を階差にしてVARを組んでしまうと、系列間の長期的なつながりに関する情報を捨ててしまうことになります。この場合には、階差を取ったVARではなく、共和分関係を明示的に組み込んだVECM（Vector Error Correction Model、ベクトル誤差修正モデル）を使うのが適切な選択です。つまり、複数系列を分析する際の定常性の確認は、単に「差分を取ればよい」という単純な話ではなく、「非定常だが共和分しているかどうか」まで踏み込んで見極める必要があり、この判断こそが第3章・第4章とVARの分析の橋渡しになります。本章では、各系列が定常である（あるいは階差を取ることで定常化できる）ケースを前提に、標準的なVARの推定と解釈を扱います。

2つ目の手順は、ラグ次数\( p \)の選択です。ラグを増やすほどモデルは過去の情報をより多く取り込めますが、その分だけ推定すべきパラメータの数が増え、過学習気味になったり、推定が不安定になったりするリスクも高まります。実務では、赤池情報量規準（AIC）をはじめとする情報量規準を、候補となる複数のラグ次数それぞれについて計算し、値が最小になるラグ次数を選ぶのが標準的な手順です。AICに加えて、BIC（ベイズ情報量規準）やHQIC（ハナン・クイン情報量規準）、FPE（最終予測誤差）といった規準が同時に算出されることも多く、いずれも当てはまりと複雑さのバランスを別々の重みで評価するものですが、規準ごとに推奨するラグ次数が食い違う場合には、実務上の解釈のしやすさやサンプルサイズとのバランスも踏まえて最終的なラグ次数を決めます。

Granger因果: 「予測に役立つか」という限定された意味での因果

VARを推定できると、「ある変数の過去の情報を加えることで、別の変数の予測が統計的に改善するかどうか」を検定できるようになります。これがGranger因果（グレンジャー因果）と呼ばれる考え方です。計量経済学者クライブ・グレンジャーが定式化したもので、「Xの過去の値を、Y自身の過去の値に加えて予測モデルに含めたとき、Yの予測誤差が統計的に有意に小さくなるならば、XはYをGranger因果する」という、純粋に予測の改善度合いに基づく定義です。

ここで強調しておきたいのは、Granger因果はあくまで「予測に役立つかどうか」を測る統計的な概念であり、私たちが日常的に使う意味での「真の因果関係」を証明するものではないという点です。この違いを丁寧に押さえておかないと、分析結果を読み違えてしまいます。Granger因果が確認された場合でも、次のような可能性を排除できていません。

• 第3の変数による見せかけの関係: XとYの両方が、観測されていない共通の要因（季節性、景気動向、競合の動きなど）によって動かされている場合、XとYの間に直接の因果関係がなくても、Xの過去がYの予測に役立つという統計的な関係だけは検出されることがあります。
• 同時性・双方向のフィードバック: 広告費が売上を動かすと同時に、売上の見込みが広告予算の意思決定に影響する、という双方向の関係がある場合、どちらが「原因」でどちらが「結果」かを、Granger因果検定だけで一方向に切り分けることはできません。
• 因果の向きの取り違え: 統計的にはXがYをGranger因果するという結果が出ても、それは「Xの過去の情報がYの予測精度を上げる」という事実を示しているに過ぎず、Xを操作すればYが変化するという介入の効果を保証するものではありません。

真に「Xを変化させればYがどう変わるか」という介入の効果を明らかにしたい場合には、ランダム化比較試験（A/Bテストなど）や、それが難しい場面での差分の差分法、傾向スコア、操作変数法といった、因果推論の専門的な手法が必要になります。Granger因果は、そうした因果推論の枠組みに進む前段階として、「どの変数とどの変数の間に、時間的な先行関係を示す統計的な手がかりがありそうか」を探索するスクリーニングの道具として位置づけるのが実務的な使い方だと言えます。

Granger因果検定の実行と解釈: 双方向・一方向・無関係のパターン

Granger因果検定を実際に2系列に適用すると、結果は大きく3つのパターンに分かれます。それぞれのパターンが、ビジネス上どのような状況に対応するかを整理すると次のようになります。

実務でこの検定結果を報告する際には、単に「有意だった/なかった」という二値の結論だけでなく、有意になったラグの範囲（何期前の情報が効いているか）や係数の符号（プラスかマイナスか）まで確認することが重要です。同じ「XがYをGranger因果する」という結果でも、1期前だけが効くのか3期前まで効くのかで、施策の効果が波及するまでの時間差についての解釈が変わってくるためです。

インパルス応答関数: ショックの波及とその持続期間を可視化する

Granger因果検定が「関係の有無」を判定する検定であるのに対し、インパルス応答関数（Impulse Response Function、IRF）は、「関係の大きさと時間的な広がり」を可視化するための道具です。具体的には、ある変数に1単位分の一時的なショック（インパルス）が加わったとき、他の変数（および自分自身）がその後何期にもわたってどのように反応するかを、VARの推定結果から計算して折れ線グラフとして描き出します。

ビジネスの文脈で言えば、「広告費に一時的な追加投下（ショック）があったとき、売上はその翌週、翌々週、それ以降と、どのように反応するか」をグラフとして確認できます。ショックの直後に売上が跳ね上がり、数週間かけてゆっくり元の水準に戻っていくのか、それとも効果はごく短期間で消えてしまうのかといった、施策の効果の「持続期間」を定量的に把握できる点が、インパルス応答関数の実務的な価値です。広告施策の予算配分やキャンペーンの実施間隔を検討する際に、この持続期間の情報は重要な判断材料になります。

ただし、インパルス応答関数の解釈には注意が必要な点が2つあります。1つ目は、識別のための仮定です。VARの推定誤差（残差）は、通常、複数の変数の間で同時点の相関を持っています。ある期に広告費と売上の両方に影響する予期しない出来事（たとえば競合の値下げ）が起きれば、両方の残差が同時に動くためです。この同時点の相関をそのままにしては「Xだけに1単位のショックが起きた」という状況を作れないため、多くの場合、残差をコレスキー分解などの手法で直交化（お互いに無相関な成分に分解）してからインパルス応答を計算します。2つ目は、この直交化の結果が、変数を並べる順序に依存するという点です。コレスキー分解では、並びの先頭に置いた変数が、同時点でも他の変数から影響を受けない（他の変数への影響のみを及ぼす）という仮定が暗黙のうちに置かれるため、変数の順序を変えると、計算されるインパルス応答の形も変わり得ます。順序をどう決めるかは、どの変数がより「外生的」に動くかという業務知識に基づく判断であり、機械的に決められるものではありません。この点を踏まえずにインパルス応答の図だけを見て結論を出すと、順序の置き方次第で結果が変わるという落とし穴にはまりやすいため、複数の妥当な順序で結果を確認する、あるいは順序に依存しない識別手法（構造VARなど）を検討するといった慎重さが求められます。

分散分解: 変動の何割が他の変数のショックで説明されるか

VARに関連するもう1つの分析手法として、分散分解（Forecast Error Variance Decomposition、FEVD）を簡単に紹介します。分散分解は、ある変数の将来の予測誤差の分散のうち、何パーセントが自分自身のショックによるもので、何パーセントが他の変数のショックによるものかを、予測のホライズン（何期先か）ごとに分解して示すものです。

たとえば、売上の予測誤差の分散のうち、1期先ではほとんどが売上自身のショックで説明されるものの、8期先になると3割程度が広告費のショックで説明されるようになる、といった結果が得られれば、広告費の変動が売上に影響を及ぼすまでにはある程度の時間を要し、その影響は短期よりも中長期の予測誤差に強く表れる、という解釈につながります。インパルス応答関数が「時間を追った反応の形」を見るのに対し、分散分解は「どの変数がどれだけ変動の説明力を持つか」を割合で見る、相互に補完的な整理の仕方だと捉えておくとよいと思います。分散分解もインパルス応答関数と同様、コレスキー分解による直交化を前提にしており、変数の順序への依存という同じ注意点を抱えている点は押さえておく必要があります。

実務の使いどころと限界: 変数を絞り込むことの重要性

VARは複数系列の関係を体系的に扱える強力な枠組みですが、実務での適用にはいくつかの現実的な制約が伴います。最も大きな制約は、変数の数とラグ次数を増やすほど、推定すべきパラメータの数が急激に膨らむという点です。変数の数を\( k \)、ラグ次数を\( p \)とすると、各方程式には定数項を含めて\( kp+1 \)個の係数があり、これが\( k \)本の方程式それぞれに存在するため、推定すべき係数の総数はおよそ\( k^2 p + k \)個になります。変数が2つでラグが2期であれば係数はおよそ10個程度で済みますが、変数が6つに増えるとラグ2期だけで70個を超える係数を推定する必要があり、月次データのような限られたサンプルサイズでは、推定量が不安定になったり、過学習によって将来の予測性能がかえって悪化したりするリスクが高まります。

もう1つの制約は、変数が増えるほど結果の解釈が難しくなるという点です。組み合わせごとにGranger因果検定を行うと検定回数が変数の数の2乗のオーダーで増え、多重比較の問題（偶然に有意な結果が紛れ込む確率の上昇）への配慮も必要になります。インパルス応答関数や分散分解も、変数が増えるほど図表が膨れ上がり、重要な関係を見失いやすくなります。

こうした事情を踏まえると、VARを実務で活用する際は、あらかじめ答えたい問いを明確にしたうえで、その問いに関わる2から4変数程度に絞り込んで分析するのが現実的なアプローチです。「広告費と売上の関係を見たい」のであればこの2変数に価格や競合動向といった重要な調整変数を1つか2つ加える程度にとどめ、手元の指標を全部投入するような進め方は避けるべきです。焦点を絞ったVARのほうが、パラメータ数を抑えて推定を安定させられるだけでなく、結果も意思決定者に説明しやすい形で提示できます。

• 分析に含める変数は、明確なビジネス上の問いに基づいて2から4本程度に絞り込む
• 各変数の定常性を事前に確認し、非定常な場合は差分か共和分（VECM）かを見極める
• ラグ次数はAICなど複数の情報量規準を比較しながら選び、規準どうしが食い違う場合は解釈のしやすさも加味する
• Granger因果の結果は「予測に役立つか」の指標であり、真の因果関係の証明ではないことを常に明記して報告する
• インパルス応答関数・分散分解は、変数の順序に依存する識別の仮定のもとで計算されていることを踏まえて解釈する

Pythonコード例: VARの推定とラグ次数の選択

statsmodelsライブラリは、VARモデルの推定、ラグ次数の選択、Granger因果検定、インパルス応答関数の計算までを一貫して扱えるインターフェイスを提供しています。まず、広告費と売上を想定した2変数のデータに対して、VARモデルを推定し、複数の情報量規準に基づいてラグ次数を選択する例を確認します。

import numpy as np
import pandas as pd
from statsmodels.tsa.api import VAR

# 広告費と売上を想定したダミーデータ(いずれも定常であることを事前に確認済みとする)
rng = np.random.default_rng(0)
n = 200
ad_spend = np.zeros(n)
sales = np.zeros(n)
for t in range(2, n):
    ad_spend[t] = 0.5 * ad_spend[t-1] + 0.1 * sales[t-1] + rng.normal(0, 1.0)
    sales[t] = 0.3 * sales[t-1] + 0.4 * ad_spend[t-1] + rng.normal(0, 1.0)

df = pd.DataFrame({"ad_spend": ad_spend, "sales": sales}).iloc[10:].reset_index(drop=True)

model = VAR(df)

# 複数のラグ次数についてAIC/BIC/HQIC/FPEを比較
lag_order_result = model.select_order(maxlags=8)
print(lag_order_result.summary())

# AICが示す次数でVARを推定
selected_lag = lag_order_result.aic
results = model.fit(selected_lag)
print(results.summary())

model.select_order()は、指定した最大ラグまでの各次数についてAIC・BIC・HQIC・FPEを一括で計算し、それぞれの規準が推奨するラグ次数をまとめて表示します。lag_order_result.aicのようにすると、AIC基準で選ばれたラグ次数を直接取得でき、そのままモデルの推定に使えます。results.summary()の出力では、各方程式（この例では広告費の式と売上の式）ごとに、すべての変数の各ラグに対応する係数と有意性が並び、どのラグのどの変数が、どちらの式に効いているかを確認できます。

Pythonコード例: Granger因果検定とインパルス応答関数

次に、推定したVARを使ってGranger因果検定を行い、あわせてインパルス応答関数を計算する例を確認します。Granger因果検定には、statsmodels.tsa.stattoolsモジュールのgrangercausalitytests関数を使う方法と、推定済みのVARの結果オブジェクトが持つtest_causalityメソッドを使う方法があります。ここでは後者を使い、双方向それぞれについて検定します。

# 「ad_spendの過去がsalesの予測に有意に寄与するか」を検定
gc_ad_to_sales = results.test_causality(caused="sales", causing=["ad_spend"], kind="f")
print(gc_ad_to_sales.summary())

# 逆方向: 「salesの過去がad_spendの予測に有意に寄与するか」を検定
gc_sales_to_ad = results.test_causality(caused="ad_spend", causing=["sales"], kind="f")
print(gc_sales_to_ad.summary())

# インパルス応答関数の計算(10期先までのショックの波及を確認)
irf = results.irf(10)
irf.plot(orth=True)  # コレスキー分解による直交化インパルス応答をプロット

# 分散分解(10期先までの予測誤差分散の内訳)
fevd = results.fevd(10)
fevd.summary()   # summary()自身が表を出力する(printで包むとNoneが余分に出る)

test_causalityメソッドのcaused引数には予測の対象となる変数を、causing引数には「その過去が予測に寄与するかどうかを検定したい変数」を指定します。kind="f"はF検定を、kind="wald"はワルド検定を指定するオプションで、いずれも帰無仮説は「causingで指定した変数の過去のラグ項の係数がすべてゼロである（=予測に寄与しない）」というものです。p値が慣例的な有意水準（0.05など）を下回れば、その方向のGranger因果が統計的に支持されたと判断します。irf.plot(orth=True)は、先述のコレスキー分解による直交化を行ったうえでのインパルス応答関数をプロットするもので、orth=Falseにすると直交化を行わない生の応答を確認することもできますが、実務では変数間の同時点の相関を踏まえたorth=Trueの結果を基本として解釈するのが一般的です。

この合成データでの実行結果は、広告費から売上への方向がF値36.97・p値0.000で有意、売上から広告費への方向がF値0.65・p値0.420で有意ではなく、一方向のGranger因果と判定されました。ここで注意したいのは、データを作る際には売上から広告費への影響（係数0.1）も入れてあるという点です。効果がノイズに対して小さいため、検定では検出されませんでした。Granger因果が検出されなかったことは、その方向に関係がないことの証明にはならないという実例でもあります。

インパルス応答関数の出力を見ると、広告費に1単位のショックを与えたとき、売上の反応は1期目に0.42でピークを迎えます。ここでの「1単位のショック」は直交化した後のショックであり、広告費そのものを1だけ増やした場合の効果とは必ずしも一致しない点には注意が必要です。ピークの後はそこから0.32、0.21、0.12と減衰し、6期目には0.04まで下がります。効果は即座に消えるのでも永続するのでもなく、数期にわたって尾を引くという読み方になります。分散分解の出力では、売上の予測誤差分散のうち広告費のショックで説明される割合が、1期先で13.4%、3期先で21.9%、そこから徐々に頭打ちになって10期先で23.1%に落ち着きます。つまり売上の変動のうち広告費で説明できるのは2割強で、残りは売上自身のショックだということが数値で読み取れます。

本章のまとめと次章予告

本章では、1本の系列の中では閉じない「系列間の関係」を捉える枠組みとして、VARモデル、Granger因果、インパルス応答関数、分散分解を扱いました。VARは各変数を全変数の過去で説明する連立の自己回帰であり、適用にあたっては第3章の定常性の確認と第4章の共和分の考え方が土台になります。Granger因果は「予測に役立つか」を測る限定的な因果概念であり真の因果関係の証明とは別物である点、インパルス応答関数と分散分解は変数の順序という識別の仮定に依存する点は、結果を正しく読み解くうえで欠かせません。実務では、変数を絞り込んだ焦点の明確なVARのほうが、安定性と説明のしやすさの両面で有利になります。

次章では、時系列モデルのもう1つの大きな体系である状態空間モデルに話を進めます。トレンドや季節性が時間とともに緩やかに変化していく状況や、観測できない潜在的な水準を推定したい状況を、より柔軟に表現できる枠組みとして、状態空間モデルの基本的な考え方とローカルレベルモデルを扱います。

この章を深めたい方への参考書籍

『実践 時系列解析』（Aileen Nielsen、オライリー・ジャパン）：単変量のARIMAから多変量のVARまでを、Pythonの実装例とともに幅広く扱っており、本章で扱ったVARの推定やインパルス応答関数の考え方を、手を動かしながら確認したい読者に適した1冊です。理論的な導出よりも実務での使い方に重心が置かれている点が、本章の実務志向の内容とよく噛み合うと思います。

第8章 見えない状態を推定する状態空間モデル入門

第5章から第7章では、AR・MAといった基本モデルから出発し、ARIMA・SARIMAによる実践的な予測、さらにVARモデルとGranger因果性による複数系列の分析までを見てきました。これらのモデルはいずれも、観測されたデータそのものの自己回帰的な構造や系列間の関係を、数式として直接組み立てるアプローチでした。

本章では、これとは異なる発想に立つ時系列モデルの枠組み、状態空間モデルを取り上げます。状態空間モデルは、目の前に見えているデータの背後に、直接は観測できない「真の状態」が存在し、その状態が確率的に変化しながらデータを生み出している、という2層構造でデータを捉え直す考え方です。ARIMAが「観測値の系列そのもの」を数式化するのに対し、状態空間モデルは「観測値を生み出している見えない裏側」を数式化します。この視点の転換によって、トレンドや季節性といった構成要素を直接推定できるようになり、欠測値の扱いや将来予測の不確実性の表現にも柔軟性が生まれます。

本章では、状態空間モデルの基本的な考え方と、その代表例であるローカルレベルモデル・ローカル線形トレンドモデル・季節成分を含むモデルを順に確認し、状態の推定に使われるカルマンフィルタの直感的な仕組みにも触れます。具体的な推定と実務での使い方については、次の第9章で改めて扱います。

発想の転換：見えている数字と、見えていない状態

日々の売上データを例に考えます。ある日の売上が急に落ち込んでも、それが「事業の実力が本当に落ちた」ことを意味するとは限りません。天候や曜日、需要のばらつきによって、観測される数字は日々上下に振れます。その裏側には、急には変化しない「その時点での本当の需要水準」があり、観測される売上はそこにノイズが乗って現れたものだと考えるのが自然です。

状態空間モデルは、この直感をそのまま数式にした枠組みです。データの背後に、直接は観測できない「状態」があると仮定し、その状態が時間とともにどのように変化していくかと、状態から観測値がどのように生成されるかを、それぞれ別の式として書き分けます。この状態には、水準（レベル）だけでなく、トレンドの傾きや季節性の大きさなど、複数の成分を同時に持たせることができます。

ARIMAとの世界観の違いを整理すると、次のようになります。

どちらが優れているという単純な優劣関係ではなく、目的に応じた使い分けが重要です。系列全体の自己相関構造を柔軟に捉えたい場合はARIMA・SARIMAが扱いやすく、水準やトレンド、季節性といった構造そのものに意味を持たせて解釈したい場合や、欠測を含むデータを扱う場合には状態空間モデルが向いています。

状態空間モデルの2つの式：状態方程式と観測方程式

状態空間モデルは、次の2本の式で構成されます。1つ目は、見えない状態が時間とともにどう変化するかを表す状態方程式です。

\( \alpha_{t+1} = T_t \alpha_t + R_t \eta_t \)

ここで\( \alpha_t \)が時点\( t \)における状態（水準・トレンド・季節成分などをまとめたもの）、\( \eta_t \)が状態に加わるノイズです。\( T_t \)は状態が次の時点へどう引き継がれるかを決める行列、\( R_t \)はどの状態にノイズが乗るかを選ぶ行列で、後で扱う具体的なモデルではいずれも単純な形になります。まずは「状態が前の時点の値を引き継ぎながら、ランダムウォークのように少しずつ確率的に変化していく」という関係を表した式だと捉えれば十分です。

2つ目は、見えない状態から実際に観測されるデータがどう生成されるかを表す観測方程式です。

\( y_t = Z_t \alpha_t + \varepsilon_t \)

ここで\( y_t \)が実際に観測される値、\( \varepsilon_t \)が観測時に乗る誤差（観測ノイズ）、\( Z_t \)は状態のうちどれを観測値として取り出すかを決める行列です。状態\( \alpha_t \)そのものは直接観測できず、そこにノイズが加わった\( y_t \)だけが手元のデータとして得られる、という構図です。

日々の売上を例にすると、状態方程式は「その日の真の需要水準は、前日の水準を引き継ぎながら、じわじわと変化していく」ことを表し、観測方程式は「実際に記録される売上高は、その日の真の需要水準に、天候や曜日といった細かなノイズが加わったものである」ことを表します。分析者が本当に知りたいのは日々の売上の記録そのものよりも、その裏にある真の需要水準の推移であることが多く、状態空間モデルはこの「知りたいもの」を直接推定の対象にできる点に強みがあります。

ローカルレベルモデル：最もシンプルな状態空間モデル

状態空間モデルの中で最も単純な形が、ローカルレベルモデルです。状態として水準\( \mu_t \)だけを持ち、次の2式で表されます。

状態方程式：\( \mu_{t+1} = \mu_t + \eta_t, \quad \eta_t \sim N(0, \sigma_\eta^2) \)
観測方程式：\( y_t = \mu_t + \varepsilon_t, \quad \varepsilon_t \sim N(0, \sigma_\varepsilon^2) \)

水準\( \mu_t \)は、前の時点の水準にノイズ\( \eta_t \)が加わって少しずつランダムウォークのように動いていくと仮定し、実際に観測される値\( y_t \)は、その水準に観測ノイズ\( \varepsilon_t \)が加わったものとして表現されます。パラメータとして推定が必要なのは、状態が変化する大きさを決める\( \sigma_\eta^2 \)と、観測時に乗るノイズの大きさを決める\( \sigma_\varepsilon^2 \)の2つだけです。

このモデルの実務上の価値は、「観測値をそのまま鵜呑みにせず、平滑化された水準を見る」ことにあります。日々の売上や日次のアクセス数、月次の受注件数といった指標は、たとえ事業の実力に変化がなくても日々細かく上下します。この細かな上下動をそのまま追いかけて一喜一憂するのではなく、ローカルレベルモデルによって推定された水準\( \mu_t \)の推移を見ることで、ノイズに惑わされない「実力の変化」を把握しやすくなります。\( \sigma_\eta^2 \)が\( \sigma_\varepsilon^2 \)に比べて小さいほど、推定される水準は滑らかになり、観測値の細かな変動の大部分がノイズとして吸収されます。

ローカルレベルモデルは「今日の数字が良かった・悪かった」という表面的な変動と、「事業の実力そのものが変化したのか」という本質的な変化を切り分けるための道具です。KPIのモニタリングにおいて、平滑化された水準の推移を追うことは、日々の数字に振り回されない意思決定につながります。

ローカル線形トレンドモデル：傾きの変化も追う

ローカルレベルモデルは水準の変化だけを捉えますが、実際のビジネスデータでは、水準が緩やかに上昇・下降していくトレンドそのものも時間とともに変化します。この傾きの変化まで状態として持たせたものが、ローカル線形トレンドモデルです。水準\( \mu_t \)に加えて、傾き（トレンド）\( \beta_t \)も状態として持ち、次のように表されます。

\( \mu_{t+1} = \mu_t + \beta_t + \eta_t, \quad \beta_{t+1} = \beta_t + \zeta_t \)

水準は「前の水準に、その時点の傾きを足したもの」として更新され、傾き自体もまた、前の傾きにノイズ\( \zeta_t \)が加わって緩やかに変化していきます。観測方程式は、ローカルレベルモデルと同じく\( y_t = \mu_t + \varepsilon_t \)です。

このモデルの実務上の利点は、トレンドの転換を追跡できる点にあります。単純な直線回帰でトレンドを当てはめた場合、期間全体を通じて傾きは一定と仮定されてしまい、途中で成長が鈍化した、あるいは減少から増加に転じたといった変化を表現できません。ローカル線形トレンドモデルでは、傾き\( \beta_t \)自体が時間とともに変化することを許容しているため、成長が鈍化し始めたタイミングや、下降トレンドが下げ止まったタイミングを、推定された\( \beta_t \)の推移から読み取ることができます。売上や会員数の中長期的な伸び率の変化を追いたい場面で特に有効です。

季節成分の追加：分解を「モデル」として行う

第2章では、時系列データを水準・トレンド・季節性・残差に分解する古典的な手法を扱いました。状態空間モデルでは、季節成分\( \gamma_t \)を状態の一部として組み込むことで、この分解を「決め打ちの計算」ではなく「確率モデルによる推定」として行うことができます。ローカル線形トレンドモデルに季節成分を加えると、観測方程式は次のようになります。

\( y_t = \mu_t + \gamma_t + \varepsilon_t \)

季節成分\( \gamma_t \)自身も、一定の周期性を保ちながら緩やかに変化していく状態として推定されます。第2章の古典分解やSTL分解と目的は似ていますが、状態空間モデルとして季節成分を扱うことには、次のような利点があります。

• 成分ごとの不確実性を評価できる：分解手法の多くは季節性の推定値だけを返しますが、状態空間モデルでは水準・トレンド・季節性それぞれについて、推定の不確実性（標準誤差）を評価できます
• 欠測処理を組み込める：分解手法の多くは連続したデータを前提としますが、状態空間モデルでは欠測期間があっても季節成分を含めた推定を継続できます
• 季節パターンの時間変化を許容できる：季節パターンが固定ではなく年々少しずつ変化していく場合、状態空間モデルはその変化の大きさ自体もパラメータとして推定します

つまり、状態空間モデルにおける季節成分の追加は、単にトレンドと季節を足し算しているのではなく、「季節パターンは毎年まったく同じではなく、緩やかに変化しうる」という前提を組み込んだ、より柔軟な分解だと捉えることができます。

カルマンフィルタの直感：予測と観測の重み付き平均

状態空間モデルにおいて、見えない状態\( \alpha_t \)を実際にどう推定するかを担うのが、カルマンフィルタと呼ばれるアルゴリズムです。数式の詳細には立ち入らず、直感的な仕組みを確認します。

カルマンフィルタは、各時点で次の2段階の処理を繰り返します。まず、前の時点までの情報から「この時点の状態はおそらくこのくらいだろう」という予測を立てます（予測ステップ）。次に、実際にその時点の観測値が得られたら、事前の予測と新しい観測値を、それぞれの不確実性の大きさに応じた重みで足し合わせ、状態の推定値を更新します（更新ステップ）。

この「重み付き平均」の重みは、カルマンゲインと呼ばれる値によって決まります。観測ノイズが大きく信頼できない状況では、新しい観測値への重みを小さくし、これまでの予測を重視します。逆に観測ノイズが小さく信頼できる状況では、新しい観測値への重みを大きくします。この重み自体も、データを処理するたびに自動的に更新されていきます。

状態の推定には、目的に応じて3つの異なる操作があります。

• フィルタ化（filtering）：時点\( t \)までの情報だけを使って、時点\( t \)の状態を推定する。リアルタイムでの現状把握に対応します
• 平滑化（smoothing）：系列全体（未来の情報も含む）を使って、過去の各時点の状態を改めて推定し直す。過去を振り返って「あの時点の真の水準はこうだった」と評価する際に使います
• 予測（forecasting）：直近までの情報から、まだ観測されていない将来の状態と観測値を推定する

平滑化はフィルタ化よりも多くの情報（未来のデータ）を使えるため、一般に推定の精度が高くなります。過去の需要水準やトレンドの推移を振り返って分析する場合には平滑化された値を、直近の状態をリアルタイムで把握したい場合にはフィルタ化された値を、それぞれ使い分けます。

欠測値への自然な対応

状態空間モデルの実務上の大きな強みが、欠測値への対応の自然さです。ARIMAの古典的な推定手法の多くは、連続した観測値の並びを前提としており、欠測があると補完（前後の値からの穴埋めなど）をしてからモデルに投入する必要がありました。ただし、この点は実装によって事情が異なります。後述するとおり、statsmodelsのARIMA・SARIMAXは内部を状態空間表現に変換して推定しているため、欠測を含む系列をそのまま渡しても推定が進みます。

状態空間モデルでは、ある時点の観測値\( y_t \)が欠けている場合、その時点の更新ステップ（観測値との重み付き平均）を単純に飛ばし、予測ステップだけを進めます。つまり「観測がない期間は、状態がどう変化していくかの予測だけを続ける」という扱いになり、欠測を挟んでも前後のデータを矛盾なくつなげて状態を推定できます。店舗の臨時休業でPOSデータが記録されなかった日や、システム障害でログが欠落した期間があっても、その前後のデータから水準やトレンドの推定を継続できるという点は、実データを扱ううえで実務的な価値が大きい特性です。

補足：ARIMAも状態空間表現で推定されている

ARIMA・SARIMAは、状態空間モデルとは異なる発想のモデルとして紹介しましたが、実は多くの統計ソフトウェアの内部では、ARIMA・SARIMAのパラメータ推定自体が状態空間モデルの枠組みを使って行われています。statsmodelsのSARIMAXクラスも例外ではなく、ARIMA・SARIMAの構造を状態空間表現に変換したうえで、カルマンフィルタを使って尤度を計算し、最尤推定によってパラメータを求めています。

つまり、状態空間モデルは「ARIMAとは別の、もう一つの手法」というだけでなく、「ARIMAを含む、より一般的で柔軟な推定の枠組み」だと捉えることもできます。第6章で扱ったARIMA・SARIMAは、状態空間モデルという大きな枠組みの中の、特定の構造を持つ一事例だと理解しておくと、両者の関係がすっきりと整理できます。

Pythonで状態空間モデルを実装する

statsmodelsのUnobservedComponentsクラスを使うと、ローカルレベルモデルやローカル線形トレンドモデル、季節成分を含むモデルを簡潔に指定できます。まずはローカルレベルモデルを当てはめ、平滑化された水準を確認します。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.statespace.structural import UnobservedComponents

# ダミーの日次売上データ（真の水準がゆっくり変化し、観測ノイズが乗っている想定）
np.random.seed(0)
n = 200
true_level = np.cumsum(np.random.normal(0, 0.5, n)) + 100
y = true_level + np.random.normal(0, 3.0, n)
dates = pd.date_range("2025-01-01", periods=n, freq="D")
sales = pd.Series(y, index=dates)

# ローカルレベルモデル：状態は水準のみ
model_level = UnobservedComponents(sales, level='local level')
result_level = model_level.fit(disp=False)

print(result_level.summary().tables[1])

# 平滑化された水準を取り出してプロット
smoothed_level = result_level.level.smoothed

plt.figure(figsize=(10, 4))
plt.plot(sales.index, sales, label="観測値（売上）", alpha=0.4)
plt.plot(sales.index, smoothed_level, label="平滑化された水準", linewidth=2)
plt.legend()
plt.title("ローカルレベルモデルによる水準の推定")
plt.show()

level='local level'と指定すると、状態が水準のみのローカルレベルモデルになります。result_level.level.smoothedには、系列全体の情報を使って推定し直した平滑化後の水準が格納されており、観測値のノイズを取り除いた滑らかな推移を確認できます。フィルタ化された値を見たい場合はresult_level.level.filteredを参照します。

続いて、トレンドと季節成分を持つより実務的なモデルを当てはめます。ここでは月次データを想定し、傾きの変化と12カ月周期の季節性を同時に推定します。

import numpy as np
import pandas as pd
from statsmodels.tsa.statespace.structural import UnobservedComponents

# ダミーの月次データ（トレンド＋季節性＋ノイズ）
np.random.seed(1)
n_months = 96
t = np.arange(n_months)
trend = 0.3 * t
season = 8 * np.sin(2 * np.pi * t / 12)
noise = np.random.normal(0, 2.5, n_months)
y_month = 50 + trend + season + noise
idx = pd.date_range("2018-01-01", periods=n_months, freq="MS")
demand = pd.Series(y_month, index=idx)

# ローカル線形トレンド＋季節成分（周期12）を持つ状態空間モデル
model_trend_season = UnobservedComponents(
    demand,
    level='local linear trend',
    seasonal=12
)
result_ts = model_trend_season.fit(disp=False)

# 平滑化された水準・トレンド・季節成分を取り出す
level_hat = result_ts.level.smoothed
trend_hat = result_ts.trend.smoothed
season_hat = result_ts.seasonal.smoothed

# statsmodels標準の可視化メソッドでも各成分をまとめて確認できる
fig = result_ts.plot_components(figsize=(10, 8), legend_loc='lower right')
plt.show()

level='local linear trend'によって水準と傾きの両方が状態として推定され、seasonal=12によって周期12（月次データにおける年周期）の季節成分が追加されます。result_ts.plot_components()を使うと、観測値・水準・トレンド・季節成分・残差を1つの図にまとめて可視化でき、第2章で扱った古典的な分解結果と似た見た目でありながら、各成分の推定の不確実性まで踏まえた出力を得られます。

最後に、欠測値を含むデータでの挙動を確認します。状態空間モデルでは、欠測期間について特別な前処理を行わなくても、そのままNaNを含む系列を渡すだけで推定が継続されます。

import numpy as np
import pandas as pd
from statsmodels.tsa.statespace.structural import UnobservedComponents

# 先ほどの月次データに、意図的に欠測期間を作る（システム障害等を想定）
demand_missing = demand.copy()
demand_missing.iloc[40:44] = np.nan  # 4カ月分のデータが欠落

model_missing = UnobservedComponents(
    demand_missing,
    level='local linear trend',
    seasonal=12
)
result_missing = model_missing.fit(disp=False)

# .level.smoothed はNumPy配列で返るため、日付インデックスを付け直して確認する
level_hat_missing = pd.Series(
    result_missing.level.smoothed, index=demand_missing.index
)

print("欠測期間の水準推定値:")
print(level_hat_missing.iloc[38:46])

# 欠測を含まない場合の推定と比較しても、前後のデータから
# 妥当な水準が補われていることが確認できる

出力を見ると、欠測させた2021年5月から8月までの水準推定値は62.09、62.40、62.71、63.02と並び、その前後の2021年4月（61.78）から2021年9月（63.33）へ滑らかにつながっています。欠測期間中は観測値との照合を行わず予測だけが進められるため、その部分の状態の推定はやや不確実性が大きくなりますが、前後のトレンドや季節パターンを踏まえた妥当な値が得られていることが数値で確認できます。欠測を含む区間だけを事前に削除したり、平均値で補完したりする必要がない点は、実データを扱う分析業務において作業負荷を軽減する実務的な利点です。

本章では、状態空間モデルの基本的な考え方として、観測データの背後にある見えない状態を状態方程式と観測方程式によって表現する枠組みを確認し、ローカルレベルモデル、ローカル線形トレンドモデル、季節成分を含むモデルへと発展させながら、カルマンフィルタによる状態推定の直感と、欠測値への自然な対応、そしてARIMAとの関係を見てきました。次の第9章では、これらのモデルを実際の需要予測やモニタリングにどう適用していくか、パラメータの解釈やモデル診断も含めて具体的に扱います。

この章を深めたい方への参考書籍

『時系列分析と状態空間モデルの基礎 RとStanで学ぶ理論と実装』（馬場真哉、プレアデス出版）：状態空間モデルをカルマンフィルタの導出から丁寧に解説しており、本章で最小限にとどめた数式の背景を体系的に補うのに適しています。RとStanによる実装例も豊富で、本章のローカルレベル・ローカル線形トレンドモデルをさらに深く理解したい場合の一冊としておすすめします。

第9章 状態空間モデルを実務に活かす：説明変数と介入効果

前章では、状態空間モデルという枠組みを使うと、水準・トレンド・季節性といった時系列の構成要素を、それぞれ独立に動く「状態」として明示的に組み立てられることを確認しました。ローカルレベルモデルやローカル線形トレンドモデル、構造時系列モデルは、いずれも系列の内部構造を分解して理解するための土台であり、カルマンフィルタという逐次的な推定の仕組みが、この土台を支えていました。

本章は、その応用編にあたります。状態空間モデルの本当の強みは、水準やトレンドといった内生的な構成要素だけでなく、価格やプロモーション、気温といった外生的な説明変数を同じ枠組みの中に自然に組み込める点にあります。さらに、説明変数の効きそのものを時間とともに変化させる時変係数モデルを使えば、「広告の効果が年々弱まっている」といった構造の変化そのものを分析の成果物として取り出せます。単一の効果の大きさを1つの数値として求めるだけでは見えてこない、効きの推移そのものを可視化できる点は、本章で扱う応用の中でも実務的な価値が大きい部分です。加えて、状態空間モデルは施策やイベントの効果を測る介入分析の入口としても使われます。本章では、これらの応用を、月次売上に価格改定という介入があったケースを想定した合成データを通しで扱いながら確認していきます。

説明変数付きの状態空間モデル：回帰成分を組み込む

状態空間モデルに説明変数を組み込む発想は、通常の回帰分析に季節性や水準の変動という「時間的に動く土台」を追加する、という見方をすると理解しやすくなります。通常の重回帰分析では、切片は一定の値として固定されますが、状態空間モデルでは、この切片に相当する部分（水準）がローカルレベルとして時間とともに緩やかに動くことを許容します。そこに、価格・プロモーションの実施有無・気温といった説明変数による回帰成分を加えることで、水準や季節性では説明できない変動のうち、既知の要因に由来する部分を切り分けられるようになります。

statsmodelsのUnobservedComponentsクラスは、levelやseasonalといった構造成分と並べて、exog引数に説明変数のデータフレームを渡すだけで、この回帰成分を追加できます。まずは、係数を時間によらず一定と仮定した、最も基本的な形から確認します。

import numpy as np
import pandas as pd
from statsmodels.tsa.statespace.structural import UnobservedComponents

np.random.seed(0)

# 2019年1月から2024年12月までの月次売上を合成する
dates = pd.date_range("2019-01-01", "2024-12-01", freq="MS")
n = len(dates)
month = dates.month

trend = np.linspace(300, 380, n)
season = 20 * np.sin(2 * np.pi * (month - 3) / 12)
intervention = (dates >= "2023-07-01").astype(int)  # 2023年7月の価格改定
price = 100 + 8 * intervention + np.random.normal(0, 3, n)  # 販売価格(単位:百円)。改定後は水準が切り替わる
promo = np.random.binomial(1, 0.25, n)            # プロモーション実施の有無
temperature = 15 + 10 * np.sin(2 * np.pi * (month - 4) / 12) + np.random.normal(0, 2, n)

sales = (
    trend + season
    - 1.2 * (price - 100)   # 価格が上がると売上が減る
    + 15 * promo            # プロモーション実施月は売上が上乗せされる
    + 0.8 * temperature
    + np.random.normal(0, 6, n)
)

df = pd.DataFrame(
    {"sales": sales, "price": price, "promo": promo, "temperature": temperature},
    index=dates,
)

model_fixed = UnobservedComponents(
    endog=df["sales"],
    level="local level",
    seasonal=12,
    exog=df[["price", "promo", "temperature"]],
)
result_fixed = model_fixed.fit(disp=False)
print(result_fixed.summary())

levelに”local level”を指定することで水準が確率的に変動する成分として、seasonal=12で12か月周期の季節成分が組み込まれます。exogに渡した価格・プロモーション・気温は、通常の回帰分析と同じように係数（beta.price、beta.promo、beta.temperatureという名前でパラメータ推定値に現れます）を持ち、その符号と大きさから、価格が1単位上がると売上がどれだけ減るか、プロモーションを実施した月には売上がどれだけ上乗せされるかを読み取れます。この段階では係数は期間全体を通じて一定と仮定されており、通常の重回帰分析との違いは、残差の中に含まれていたはずの水準の緩やかな動きや季節性を、あらかじめ別の状態として切り分けている点にあります。この切り分けによって、価格やプロモーションの効果を、トレンドや季節性による見かけの変動と混同せずに推定しやすくなります。

係数を固定するか時間変化させるかという判断の前段階として、そもそもどのような水準成分を採用するかという設計も欠かせません。ローカルレベル（水準のみが確率的に変動する）で十分なのか、ローカル線形トレンド（水準に加えて傾きも確率的に変動する）まで必要なのかは、AIC（赤池情報量規準）やBIC（ベイズ情報量規準）といった指標でモデル同士を比較しながら判断します。あわせて、残差にコレログラムやLjung-Box検定で周期的な自己相関が残っていないかを点検する作業も実務では欠かせません。水準や季節性の設計を誤ったまま回帰成分を追加すると、価格やプロモーションの係数の推定自体が歪んでしまうためです。実際、この合成データでAICを比較すると、ローカルレベルが431.3、ローカル線形トレンドが406.9となり、上昇トレンドを持つこの系列では後者が支持されます。

もう1つ、この出力から読み取っておきたい点があります。3つの共変量のうち、気温の係数だけが0.57（p値0.193）と有意になりません。原因は、この合成データの気温が年周期の正弦波として作られており、モデルの季節成分と強く重なっている（両者の相関は約0.82）ことにあります。気温が持つ効果の大半を季節成分が吸収してしまうため、気温そのものの寄与を切り出せなくなります。年周期を持つ共変量を季節成分と同時に投入すると、両者は互いに識別しづらくなるという点は、状態空間モデルで説明変数を扱う際の実務的な注意点です。気温のように季節と連動する変数の効果を本当に取り出したい場合は、季節成分を落として気温に説明を委ねるか、気温から季節の平年値を引いた偏差（平年より暑いか寒いか）を変数として使う、といった設計の工夫が必要になります。

時変係数モデル：効果の推移そのものを分析の成果物にする

係数を固定するという仮定は、扱う期間が短ければ実務上大きな問題になりません。しかし、数年にわたる月次データを扱う場合、「プロモーションの効果が最初のうちは大きかったが、競合の追随や消費者の慣れによって年々弱まっている」といった、係数そのものが時間とともに変化する状況は珍しくありません。この「効きの変化」を捉えるための拡張が時変係数モデルです。

UnobservedComponentsでは、mle_regressionという引数をFalseに設定すると、回帰係数を固定パラメータとして最尤推定するのではなく、状態ベクトルの一部として扱い、カルマンフィルタで逐次的に推定させることができます。この設定では係数そのものは一定という仮定のままですが、フィルタ化推定値を見ると「各時点までの情報だけで見たときの係数の推定値」が時間とともに更新されていく様子を確認できます。係数自体がランダムウォークに従って確率的に変化する、文字どおりの時変係数モデルを推定したい場合は、SARIMAXクラスでtime_varying_regression=True（mle_regression=Falseと併用）を指定する方法が用意されています。

# プロモーション係数を状態として逐次推定するモデル
model_tv = UnobservedComponents(
    endog=df["sales"],
    level="local level",
    seasonal=12,
    exog=df[["promo"]],
    mle_regression=False,
)
result_tv = model_tv.fit(disp=False)

# フィルタ化された状態から、各時点までの情報で見た係数の逐次推定値を取り出す
beta_promo = result_tv.states.filtered["beta.promo"]
print(beta_promo.iloc[[0, 11, 23, -1]])

# 係数自体をランダムウォークとして変化させる時変回帰はSARIMAXで指定できる
# トレンドと季節性は差分で先に処理しておく(これを怠るとAR項がそれらを吸収し、
# 係数の軌跡が信用できなくなる)
from statsmodels.tsa.statespace.sarimax import SARIMAX
model_tvr = SARIMAX(
    df["sales"],
    exog=df[["promo"]],
    order=(0, 1, 1),
    seasonal_order=(0, 1, 1, 12),
    time_varying_regression=True,
    mle_regression=False,
)
result_tvr = model_tvr.fit(disp=False)

# 係数の軌跡は平滑化推定値から取り出す
beta_path = result_tvr.states.smoothed["beta.promo"]
print(beta_path.resample("YE").mean())

result_tv.states.filteredは、カルマンフィルタで推定された各時点の状態（水準・季節成分・回帰係数）をまとめたデータフレームで、beta.promoという列に、その月までの情報で見たプロモーション係数の逐次推定値が並びます。この例では初期値の0.00から始まり、12か月目で26.26、24か月目で16.30、最終的に15.99へと落ち着きます。

ここで、この逐次推定値の読み方には決定的な注意が要ります。この推移は「効果が変化した」ことを表しているのではなく、「データが増えるにつれて推定が定まっていく」学習の過程を表しています。mle_regression=Falseの設定では、係数は状態ベクトルに入るものの、時間を通じて一定という仮定は保たれたままです。実際、同じモデルの平滑化推定値（全期間の情報を使った推定値）を見ると、どの時点も15.99で一定になります。初期の大きな振れは、状態の初期値から推定が収束していく過程であり、これを効果の推移として折れ線に描くと「プロモーションの効きが年々弱まっている」という、まったく事実に反する読み方を招きます。

係数そのものが時間とともに動くことを許すには、SARIMAXのtime_varying_regressionを使います。ここで大切なのは、トレンドと季節性をあらかじめ差分で処理しておくことです。これを怠ってorder=(1, 0, 0)のようなモデルに当てはめると、AR項が単位根に張り付いてトレンドと季節性を丸ごと吸収してしまい、係数の軌跡は実態を反映しないものになります。差分を入れたうえでの推定結果では、係数の年平均は2019年の14.9から2024年の18.0までの範囲で緩やかに上下しました。

この軌跡をどう読むかが、実務では最も重要な判断になります。今回の合成データはプロモーション効果を一定（15）として作っているため、観測された14.2から18.1という振れ幅は、効果が実際に変化したことを示すものではなく、推定の不確実性そのものです。プロモーションを実施した月は72か月中16か月しかなく、係数の変化を語るにはそもそも情報が足りません。時変係数モデルを実務で使うときは、係数がどれだけ機敏に動くかを決める分散パラメータの推定値と、軌跡の不確実性の幅を必ず併せて確認し、その幅を明確に超える動きがあるときだけ「効きが変わった」と判断するのが安全です。あわせて、対象の施策が十分な回数実施されているか（施策を打った期間の数が足りているか）も、分析に着手する前に確認しておく必要があります。

通常の回帰分析であれば「期間全体を通じた平均的な効果」という1つの数値しか得られませんが、条件が整えば、時変係数モデルではこの係数の推移そのものが分析の成果物になります。広告費の投下効果を継続的にモニタリングしたい場合や、値引き施策の効果が消費者の慣れによって目減りしていないかを確認したい場合に有効な見方ですが、その判断には、施策の実施回数と観測期間の長さという前提条件がついて回ります。

時変係数モデルの価値は、単一の「効果の大きさ」を1つの数値として出すことではなく、その数値が時間とともにどう推移しているかという軌跡そのものを可視化できる点にあります。プロモーションや広告、価格施策の効きが年々変わっているかどうかを継続的に把握したい場面では、固定係数の回帰分析では見えてこない情報を引き出せます。

介入分析の入口：反実仮想と実績を比較する発想

状態空間モデルは、施策やイベントの効果を測る介入分析の入口としても使われます。基本的な発想は、介入が起きる前のデータだけを使ってモデルを作り、そのモデルが「もし介入が起きなかったとしたら、その後どう推移していたはずか」という反実仮想の予測を行い、実際に観測された値との差を効果とみなす、というものです。この考え方は、Googleが公開したCausalImpactというツールで広く知られるようになったアプローチで、状態空間モデルによる予測と信頼区間をそのまま因果効果の推定に転用する発想に基づいています。

この発想が成立するためには、いくつかの前提を満たす必要があります。第一に、モデルに使う共変量（価格以外の説明変数など）が、介入そのものの影響を受けていないことです。介入によって変化してしまう変数を予測に使うと、反実仮想の予測自体が歪んでしまいます。第二に、介入前に成立していた変数間の関係が、介入後も安定して成り立ち続けるという仮定です。介入前後で市場環境や消費者の行動様式が大きく変わってしまえば、介入前のデータで学習したモデルによる予測は当てになりません。第三に、介入と同じタイミングで、モデルが捉えていない別の要因（競合の値下げ、突発的な景気変動など）が重なっていないことです。こうした前提が崩れている疑いがある場合、算出された効果は介入そのものの効果ではなく、複数の要因が混ざった結果である可能性が残ります。

• 介入前後を通じて、モデルに使う共変量自体が介入の影響を受けていないか
• 介入前に見られた変数間の関係が、介入後も継続すると考えてよいか
• 介入と同時期に、他の無視できない要因（競合の動き、外部環境の急変など）が重なっていないか
• 反実仮想の予測区間（信頼区間）を踏まえても、実績との差が誤差の範囲を超えていると言えるか

ここで得られる効果の推定値は、あくまで単一の反実仮想モデルに基づく参考値である点にも注意が必要です。水準成分の選び方や、どの共変量を含めるかといったモデルの構造を変えると、反実仮想の予測自体も変わり、結果として推定される効果の大きさも動きます。実務では、水準成分や共変量の組み合わせを複数パターン試し、効果の推定値がモデルの選び方によってどの程度変動するかを確認したうえで、結論の頑健性を判断する姿勢が求められます。

本章で扱うのは、この介入分析の入口となる考え方と、状態空間モデルを使った実装の基本形までです。効果の統計的な有意性の検証や、複数の対照群を組み合わせた設計、傾向スコアなどを用いたより厳密な因果推論の方法論については、本シリーズの因果推論と効果検証編で、施策効果検証の実務プロセスとあわせて詳しく扱う予定です。

構造変化・レベルシフトをダミー変数で吸収する

介入の影響は、必ずしも「反実仮想との比較」という形でしか扱えないわけではありません。介入によって系列の水準が一段階だけ切り替わることがあらかじめ分かっている場合には、ダミー変数を説明変数として組み込むことで、その段差を直接モデルに吸収させる方法も実務ではよく使われます。

このうちステップ型のダミー変数は、価格改定のように、ある時点を境に水準そのものが切り替わり、その後も切り替わった水準が続くと想定される状況に対応します。exogに他の説明変数と並べてこのダミー列を加えるだけで、UnobservedComponentsは水準の恒久的な段差を、季節性やその他の回帰効果と切り離して推定します。前節で扱った反実仮想による効果測定と、このステップ型ダミーによる効果測定は、同じ介入の効果を異なる角度から捉える手法であり、両方の結果を突き合わせることで、推定の頑健性を確認できます。

ケーススタディ：価格改定の効果を状態空間モデルで読み解く

ここまでの内容を、月次売上に価格改定という介入があったケースの合成データで通して確認します。2019年から2024年までの月次売上データのうち、2023年7月に価格改定が実施され、それ以降、価格の水準そのものが切り替わったという設定です。介入前の54か月分（2019年1月から2023年6月まで）を使ってモデルを学習し、介入後の期間について「価格改定が起きなかったとしたら」という反実仮想の予測を行い、実際の売上の推移と比較します。

# 価格改定(2023年7月)の前後でデータを分割する
intervention_date = "2023-07-01"
pre = df.loc[:intervention_date].iloc[:-1]   # 介入前のデータ(学習用)
post = df.loc[intervention_date:]            # 介入後のデータ(検証用)

# 介入前のデータだけでモデルを学習する
# 価格自体は介入によって変化する変数なので、ここではプロモーションと気温のみを共変量に使う
# この系列は上昇トレンドを持つため、水準だけでなく傾きも状態に持つローカル線形トレンドを使う
model_pre = UnobservedComponents(
    endog=pre["sales"],
    level="local linear trend",
    seasonal=12,
    exog=pre[["promo", "temperature"]],
)
result_pre = model_pre.fit(disp=False)

# 介入後の期間について、介入が起きなかった場合の反実仮想を予測する
forecast = result_pre.get_forecast(steps=len(post), exog=post[["promo", "temperature"]])
counterfactual = forecast.predicted_mean
ci = forecast.conf_int(alpha=0.05)

# 実績との差を効果として算出する
effect = post["sales"] - counterfactual
cumulative_effect = effect.cumsum()
print(f"介入後1か月目の推定効果: {effect.iloc[0]:.1f}")
print(f"介入後12か月間の累積推定効果: {cumulative_effect.iloc[11]:.1f}")

get_forecastは、学習済みモデルから将来期間（ここでは介入後の期間）の予測分布を返すメソッドで、predicted_meanが点予測、conf_intが指定した信頼水準での予測区間にあたります。counterfactualは、価格改定が起きず、それまでの水準・季節性・プロモーションや気温との関係がそのまま続いた場合に予想される売上であり、これと実際の売上post[“sales”]との差effectが、価格改定による効果の推定値になります。実行すると、介入後1か月目の推定効果は-7.1、介入後12か月間の累積推定効果は-144.6となりました。データ生成時に設定した価格改定の影響は月あたり-9.6（単価8×係数-1.2）なので、12か月では-115.2にあたります。推定値がこれをやや上回っているのは、介入後の期間に含まれる通常の変動が上乗せされているためです。この差の推移を月ごとに追うことで、効果が改定直後に大きく出るのか、時間とともに拡大あるいは縮小するのかを確認できます。cumulative_effectは、この月ごとの効果を累積した系列で、介入後1年間でどれだけの売上への累積的な影響があったかを1つの数値にまとめる際に使われます。

この結果を解釈する際は、effectの推移が予測区間ciの幅と比べてどの程度大きいかを必ず確認する必要があります。実績が予測区間の中に収まっている月については、価格改定の効果とは言い切れず、通常の予測誤差の範囲内である可能性が残ります。反対に、実績が予測区間を明確に超えて推移している期間については、価格改定以外の要因が同時に重なっていないかを確認したうえで、効果として扱うのが妥当な進め方です。この例では、介入後18か月のうち5か月で実績が予測区間の外に出ており、効果が予測誤差では説明しきれない大きさであることを裏づけています。

なお、ここで水準成分にローカル線形トレンドを選んでいる点は、この分析の成否を分ける設計判断です。この系列は上昇トレンドを持つため、水準だけを状態に持つローカルレベルを使うと、反実仮想の予測が最後の水準のまま横ばいに固定されます。すると、時間が経つほど「介入がなかった場合の水準」を過小に見積もることになり、本来マイナスであるはずの効果が打ち消され、後半では符号が逆転してプラスに見えるという結果を招きます。前節で述べたAICによる水準成分の比較は、こうした事態を避けるための手順です。

この反実仮想による推定は、前節で扱ったステップ型ダミー変数による推定と突き合わせることで検証できます。介入前後を通じた全期間のデータに対して、promoとtemperatureに加えて2023年7月以降を1とするダミー変数を組み込んだモデルを別途学習し、そのダミーの係数が、cumulative_effectを介入後の月数で割った平均的な効果とおおむね近い水準になっているかを確認します。この例では、ステップダミーの係数が-12.57、反実仮想による累積効果を12か月で割った値が-12.05となり、2つの異なる角度からの推定がよく一致しました。両者が大きく食い違う場合は、水準成分の選び方や共変量の設定など、モデルの構造のどこかに見直すべき点が残っている可能性が高く、単一の手法だけで結論を急がない姿勢が実務では重要になります。

状態空間モデルが力を発揮するのは、系列を水準・季節性・回帰効果といった成分に分解して理解したい場面、欠測値を含むデータをそのまま扱いたい場面、係数が時間とともに変化する構造を捉えたい場面、そして施策やイベントの介入効果を測りたい場面です。反対に、成分の解釈や介入効果の測定を必要とせず、純粋に予測精度だけを追求したいのであれば、第6章で扱ったSARIMAや、第10章で扱う機械学習による時系列予測との比較検証を行い、目的に合った手法を選ぶ姿勢が実務では重要になります。

ベイズ的な拡張への案内

本章で扱ったUnobservedComponentsは、最尤推定によって単一のパラメータ値を求める枠組みですが、状態空間モデルはベイズ統計の枠組みとも相性が良く、StanやPyMCといった確率的プログラミング言語を使うことで、より柔軟な拡張が可能になります。たとえば、店舗や商品カテゴリごとに水準やトレンドを個別に推定しながら、全体で情報を共有する階層的な状態空間モデルや、売上のように非負で裾の重い分布に従うデータを正規分布ではなくポアソン分布や負の二項分布で扱う非ガウス型の状態空間モデルは、ベイズ的な推定の枠組みでこそ実装しやすくなります。こうしたベイズ統計と状態空間モデルを組み合わせた発展的な内容は、本シリーズのベイズ統計編で章を割いて扱う予定です。最尤推定による本章の枠組みでも、多くの実務課題には十分対応できますが、店舗数や商品カテゴリ数が多く、個別の系列ごとにデータ数が少ない状況では、ベイズ的な階層構造によって系列間で情報を共有する発想が精度の面で有利に働く場面が少なくありません。

本章では、状態空間モデルに説明変数を組み込む回帰成分の拡張、係数の推移そのものを分析の成果物とする時変係数モデル、そして介入分析の入口となる反実仮想の考え方とダミー変数による構造変化の扱いを確認しました。次章では、視点を変えて、勾配ブースティングやニューラルネットワークといった機械学習の手法を時系列予測にどう応用するか、そして状態空間モデルや古典的な統計モデルとどう使い分けるかを扱います。

この章を深めたい方への参考書籍

『カルマンフィルタ Rを使った時系列予測と状態空間モデル』（野村俊一、共立出版）：状態空間モデルの理論的な背景とカルマンフィルタの導出を丁寧に扱いながら、説明変数の組み込みや時変係数への拡張についても踏み込んで解説している一冊です。本章で扱った回帰成分や時変係数モデルの数理的な裏付けを、より深く追いたい読者に向いています。

第10章 機械学習による時系列予測と古典手法の使い分け

ここまでの章では、ARモデル・MAモデル・ARIMA・SARIMA・VARモデル・状態空間モデルという、統計学の伝統に連なる時系列モデル群を扱ってきました。これらの手法を学ぶと、「結局、時系列予測は機械学習で行えばよいのではないか」という疑問が自然に出てきます。GBDT（勾配ブースティング決定木）やニューラルネットワークは、機械学習編第8章で見た通りテーブルデータの予測において高い性能を発揮します。であれば時系列予測でも、機械学習が古典的な統計モデルを置き換えてよいはずだ、という発想も一理あります。

この問いへの答えは、「場合による」というのが本ガイドの立場です。機械学習は時系列予測において確かに強力な選択肢であり、特定の条件下では古典的な統計モデルを上回る精度を示すことが少なくありません。しかし、あらゆる状況で機械学習が古典手法に勝るわけではなく、系列の数や長さ、外部データの豊富さ、解釈可能性や予測区間の要否といった条件によって、有利な手法は入れ替わります。本章では機械学習を時系列予測に適用するための定石を確認したうえで、どのような場面でどちらの手法が有利になりやすいかを整理し、深層学習や時系列基盤モデルの現在地にも触れます。最終的には、手法への愛着ではなく検証によって選ぶという実務上の結論に至ります。

機械学習で時系列を扱う定石：教師あり学習への変換

機械学習の多くの手法は、本来「時間の順序」という概念を持ちません。ランダムフォレストやGBDT、線形回帰といった手法は、行と列からなる通常のテーブルデータを入力として受け取り、各行を互いに独立なサンプルとして扱います。したがって時系列データにこれらの手法を適用するには、時系列予測問題を教師あり学習の問題へと変換する作業が必要になります。

この変換の基本的な考え方は、ある時点の目的変数（たとえば来月の売上）を予測するために、その時点までに観測できたはずの情報を特徴量として横に並べ直すというものです。時間軸に沿って1列に並んでいたデータを、行ごとに「予測に使ってよい過去の情報」と「予測したい将来の値」の組み合わせへと組み替え、通常の回帰問題と同じ形式に落とし込みます。この変換で使われる代表的な特徴量には、次のような種類があります。

ラグ特徴量と移動平均特徴量は、ARモデルが「過去の値の線形結合で現在の値を説明する」という発想を特徴量エンジニアリングの形に置き換えたものと理解すると見通しがよくなります。カレンダー特徴量は、SARIMAが季節差分によって扱っていた季節性を、月や曜日といった特徴量として機械学習モデルに直接与える役割を果たします。イベントフラグは状態空間モデルの説明変数に近く、モデルが自力では学習できない「なぜその日だけ値が跳ねたのか」という背景知識を人間の側からモデルに教え込む手段です。

この変換作業を終えたあとの実務では、モデルのアルゴリズム選びよりも、どのようなラグを取るか、どの期間の移動統計量を作るか、どのカレンダー・イベント情報を用意するかという特徴量エンジニアリングそのものが主戦場になる傾向があります。同じGBDTのアルゴリズムを使っていても特徴量の設計次第で精度が大きく変わることは珍しくなく、時系列予測に機械学習を持ち込む作業の大半は、この特徴量表をどれだけ丁寧に作り込めるかにかかっています。

ここで、簡単な小売店舗の日次売上データを例に、ラグ特徴量とカレンダー特徴量を作成するコードを確認します。

import numpy as np
import pandas as pd

np.random.seed(0)

dates = pd.date_range("2022-01-01", "2024-12-31", freq="D")
n = len(dates)
trend = np.linspace(50, 90, n)
weekday_effect = np.where(dates.weekday >= 5, 15, 0)  # 週末に売上が増える効果
noise = np.random.normal(0, 5, n)
sales = trend + weekday_effect + noise

df = pd.DataFrame({"sales": sales}, index=dates)
df.index.name = "date"

# ラグ特徴量(1日前・7日前・14日前)
for lag in [1, 7, 14]:
    df[f"lag_{lag}"] = df["sales"].shift(lag)

# 移動平均・移動標準偏差特徴量(直近7日・28日、shiftで当日の値を含めない)
df["rolling_mean_7"] = df["sales"].shift(1).rolling(window=7).mean()
df["rolling_mean_28"] = df["sales"].shift(1).rolling(window=28).mean()
df["rolling_std_7"] = df["sales"].shift(1).rolling(window=7).std()

# カレンダー特徴量
df["weekday"] = df.index.weekday
df["month"] = df.index.month
df["is_weekend"] = (df["weekday"] >= 5).astype(int)
df["is_month_end"] = df.index.is_month_end.astype(int)

df = df.dropna()
print(df.head())

このコードで注意すべき点は、rolling_mean_7やrolling_mean_28を計算する前に必ずshift(1)を挟んでいることです。shiftを挟まずにdf["sales"].rolling(window=7).mean()とすると、当日の売上自体が移動平均の計算に含まれてしまい、予測時点ではまだ手に入らない「未来の情報」を特徴量に混入させることになります。この点は次節で扱う検証設計とも密接に関わる、時系列における特徴量エンジニアリングの最重要注意点です。dropna()は、ラグや移動平均の計算によって生じる欠損行（データの先頭部分）を取り除く処理です。

GBDTによる時系列予測とグローバルモデルという発想

特徴量への変換を終えたあとの時系列予測で広く使われるのが、LightGBMをはじめとするGBDT（勾配ブースティング決定木）です。機械学習編第8章で確認した通り、GBDTは欠損値の扱いやすさ、特徴量間の非線形な関係やカテゴリ変数との相性、チューニングの手間の少なさから、テーブルデータに対する予測モデルとして高い実績を持ちます。時系列予測をラグ・移動平均・カレンダー特徴量を持つテーブルデータへの回帰問題として組み替えれば、この強みをそのまま持ち込めます。

GBDTを時系列予測に用いる際に特に威力を発揮する発想が、グローバルモデルと呼ばれる考え方です。従来のARIMAや状態空間モデルの多くは、1つの系列につき1つのモデルを個別に当てはめる「ローカルモデル」の発想に立っています。店舗が100店舗、商品が50品目あれば、組み合わせごとに5,000本の売上系列が存在し、ローカルモデルの発想ではこの5,000系列それぞれに個別のARIMAモデルを構築することになります。これに対してグローバルモデルは、店舗や商品の識別子そのものを特徴量として持たせたうえで、5,000系列すべてを1つのモデルに学習させるという発想です。

グローバルモデルが威力を発揮しやすいのは、次のような場面です。系列の数が非常に多い場合、個別にモデルを構築・保守するコストが現実的でなくなるため、1つのモデルで全系列をまとめて扱える運用上の利点が大きくなります。また、系列どうしに関連する情報がある場合、たとえば同じ地域の店舗どうしや似た性質を持つ商品カテゴリどうしでは、ある系列で学習したパターンを他の系列の予測に活用できます。特に観測期間の短い新商品や新規出店の系列は、単独では十分なデータ量を持ちませんが、類似する既存系列の情報と合わせることで予測精度を底上げできる場合があります。

グローバルモデルの考え方は、系列が少数であればあまり恩恵がありません。系列が数本程度であれば、それぞれに丁寧なローカルモデルを当てはめるほうが、系列固有の癖を細かく捉えられることが多いためです。グローバルモデルが優位になりやすいのは、系列数が数十から数千、数万という規模に達し、かつ系列間で何らかの構造的な類似性が期待できる場面だと整理しておくとよいでしょう。

実装の面では、店舗と商品の組み合わせで複数系列を持つデータフレームに対し、店舗ID・商品IDそのものをカテゴリ特徴量として追加したうえで、次節で扱うLightGBMやHistGradientBoostingRegressorに1つのモデルとして学習させるだけで、グローバルモデルの構成は実現できます。scikit-learnのHistGradientBoostingRegressorであればcategorical_featuresパラメータに識別子の列名を指定することで、ダミー変数への変換なしに水準の違いをモデルへ直接取り込めます。ただしこの機能で扱えるカテゴリの数には上限（既定では列あたり255種類）があるため、商品SKUのように種類が桁違いに多い識別子をそのまま渡すことはできません。その場合は、カテゴリを粒度の粗い区分にまとめる、あるいは過去実績から作った統計量（商品ごとの平均売上など）に置き換えるといった前処理が必要になります。実務では、季節性の強さに応じてラグの種類を増やす、店舗の売場面積や商品の価格帯といった外部情報を加えるといった拡張がよく行われます。

検証設計：時間順を守った分割とリークへの警戒

機械学習モデルを時系列予測に適用する際、精度そのものと同じくらい重要なのが検証設計です。第1章では、通常のk-Fold交差検証がデータの各行の独立性を前提としており、時系列データに適用すると訓練データに未来の情報が入り込むリークを引き起こすことを確認しました。この問題への対処が、scikit-learnのTimeSeriesSplitに代表される、時間順を守った分割です。

TimeSeriesSplitは、常に時間的に前のデータを訓練用、それより後のデータを検証用に割り当てる分割方法で、分割を重ねるごとに訓練データの期間が伸びていく形を取ります。これにより、モデルは実運用時と同じ条件のもとで評価されます。

from sklearn.model_selection import TimeSeriesSplit
from sklearn.metrics import mean_absolute_error
from lightgbm import LGBMRegressor

feature_cols = [
    "lag_1", "lag_7", "lag_14",
    "rolling_mean_7", "rolling_mean_28", "rolling_std_7",
    "weekday", "month", "is_weekend", "is_month_end",
]
X = df[feature_cols]
y = df["sales"]

tscv = TimeSeriesSplit(n_splits=5, test_size=30)  # 各foldのテスト期間を30日に固定
mae_scores = []

for train_idx, test_idx in tscv.split(X):
    X_train, X_test = X.iloc[train_idx], X.iloc[test_idx]
    y_train, y_test = y.iloc[train_idx], y.iloc[test_idx]

    model = LGBMRegressor(max_depth=6, learning_rate=0.05, n_estimators=300, random_state=0)
    model.fit(X_train, y_train)
    pred = model.predict(X_test)
    mae_scores.append(mean_absolute_error(y_test, pred))

print(f"各foldのMAE: {[round(s, 2) for s in mae_scores]}")
print(f"平均MAE: {np.mean(mae_scores):.2f}")

test_size=30を指定することで、各foldのテスト期間を30日分に固定しています。ここで注意が必要なのは、テスト期間の長さと予測ホライズンは別物だという点です。この特徴量には前日の実測値（lag_1）が含まれているため、この評価は1日先予測を30日分繰り返した結果であり、30日先を一度に予測する条件とは異なります。数期先を予測する運用であれば、予測値を順に入力へ戻す再帰予測か、ホライズンごとに別のモデルを学習するdirect方式に切り替えたうえで、同じ形で検証する必要があります。分割数を増やすほど、異なる時期における予測精度のばらつきを確認でき、特定の期間だけ精度が悪化していないかを点検する材料にもなります。

検証設計において機械学習ならではの注意点となるのが、リークしやすい特徴量への警戒です。特に危険なのが、集計処理の対象期間に予測対象の時点自体が含まれてしまう特徴量です。前節のコード例でrolling_mean_7の計算前にshift(1)を挟んだのはこのためであり、この一手間を怠ると、モデルは学習時にだけ「本来知り得ない当日の情報」を利用でき、検証スコアは実態よりも大幅に楽観的な値を示します。同様の注意は、店舗全体の月間平均売上のような集計特徴量にも当てはまり、月がまだ終わっていない時点での予測に「未来を含む集計値」をそのまま使ってしまうと、同じ種類のリークが発生します。予測時点で確実に判明している情報だけを特徴量として使うという原則を、特徴量を追加するたびに確認する姿勢が欠かせません。

古典的手法が勝ちやすい場面、機械学習が勝ちやすい場面

ここまで機械学習による時系列予測の定石を見てきましたが、これはARIMAやETS、状態空間モデルといった古典的な統計モデルが不要になったことを意味しません。両者には得意な条件があり、実務ではこの条件を踏まえて手法を選び分ける必要があります。

まず、古典的な統計モデルが優位になりやすい場面を整理します。

• 予測対象の系列の数が少ない、あるいは1系列のみである。ラグや移動平均を作り込む学習データが十分に確保しづらく、機械学習モデルの複雑さを支えるだけのデータ量が不足しやすい
• 観測期間が短い。季節周期を十分な回数観測できていない場合、機械学習モデルは季節パターンをうまく学習できないが、SARIMAや状態空間モデルは季節成分を構造として明示的に組み込めるため、短い履歴でも季節性を捉えやすい
• トレンドや季節性、他の系列との関係といった構造そのものを解釈し、意思決定の根拠として説明する必要がある。ARIMAの係数や状態空間モデルの状態変数はその意味を言葉で説明しやすい
• 予測区間（不確実性の幅）が業務上必須である。ARIMAや状態空間モデルは理論的な裏付けを持つ予測区間を自然に算出できる
• データを生み出す背後の構造が比較的安定しており、急激な構造変化が起きにくい

一方、機械学習が優位になりやすい場面は、次のように整理できます。

• 予測対象の系列が大量にある。前節のグローバルモデルの発想が生きる場面であり、店舗×商品のような膨大な組み合わせを1つのモデルで扱える運用上の利点も大きい
• 気象データ、競合の価格、マーケティング施策の実施状況、経済指標といった外部特徴量が豊富に用意でき、多数の外部変数を柔軟に取り込める
• 目的変数と特徴量の関係が非線形である。ARIMAや状態空間モデルの多くは線形関係を前提とするのに対し、GBDTや決定木ベースの手法は非線形な交互作用を自動的に学習できる
• 長い履歴データが蓄積されており、ラグや移動平均を十分な期間にわたって作成できる

実務でよく起こるのは、この2種類の条件が混在するケースです。系列数は多いが一部の系列は観測期間が極端に短い、あるいは外部特徴量は豊富だが予測区間の提示も求められる、といった状況では、どちらか一方を機械的に選ぶのではなく、両者を組み合わせる、あるいは対象系列のグループごとに手法を使い分ける判断が必要になります。第12章では、こうした実務上の判断の設計をあらためて扱います。

深層学習と時系列基盤モデルの現在地

時系列予測における機械学習の話題として避けて通れないのが、深層学習の応用と、近年注目を集める時系列基盤モデルの動向です。深層学習による時系列予測は、大きくRNN（再帰型ニューラルネットワーク）やその発展形であるLSTM（Long Short-Term Memory）を用いる系統と、自然言語処理の分野で成果を上げたTransformerを時系列向けに応用する系統に分かれます。RNN・LSTM系は系列を時間方向に沿って逐次処理しながら内部状態を更新する構造を持ち、比較的長い依存関係を捉える設計として発展してきました。Transformer系は注意機構（Attention）によって系列内の離れた時点どうしの関係を直接学習でき、複数系列や長期の予測ホライズンを扱う研究で応用が進んでいます。

さらに近年の潮流として注目されているのが、時系列基盤モデルと呼ばれる考え方です。これは、大量かつ多様な時系列データであらかじめ大規模に事前学習を行ったモデルを用意しておき、個別のタスクに合わせた追加学習（ファインチューニング）をほとんど行わずに、初めて見る系列に対してもある程度の精度で予測を行う、ゼロショット予測を志向するアプローチです。自然言語処理における大規模言語モデルの発展に触発される形で、この種の基盤モデルの研究開発が活発化しています。

ただし、深層学習モデルや時系列基盤モデルが、あらゆる場面で古典的な統計モデルやGBDTを上回るという理解は正確ではありません。時系列予測の大規模な比較の場として知られるMコンペティションの結果は、この点を端的に示しています。10万系列を対象とした2018年のM4では、上位を占めたのは統計的手法の組み合わせと、統計モデルとニューラルネットワークを組み合わせたハイブリッド手法であり（優勝したのは後者です）、純粋な機械学習手法の成績は振るいませんでした。ところが小売の階層データを対象とした2020年のM5では、上位の手法がいずれも機械学習ベース（多くはLightGBMを含む構成）となり、統計的なベンチマークとその組み合わせを明確に上回りました。同じ「時系列予測」という課題でも、データの性質と評価条件が変われば順位は入れ替わるということです。深層学習モデルは一般に多くの学習データと計算資源を必要とし、系列数が少ない、あるいは履歴が短いデータでは、その複雑さに見合うだけの恩恵を得られないことも多いと報告されています。時系列基盤モデルについても、特定分野への適応やドメイン固有の季節性・イベントへの対応という点では、なお発展途上にあるという評価が一般的です。有力な選択肢の1つとして押さえつつ、常に古典手法を上回るという前提には立たず、個別のデータで実際に比較検証するという姿勢が求められます。

Prophetの位置づけ：手軽さと限界

時系列予測のライブラリとして名前が挙がることの多いProphetにも簡単に触れておきます。Prophetは、トレンドと季節性、休日効果を加法的に組み合わせるモデルをあらかじめ用意し、比較的少ないパラメータ設定で扱えるように設計されたライブラリです。日付とイベント情報さえ整えれば手早く予測モデルを構築できる手軽さが特長で、分析の初期段階でのたたき台として活用されてきました。

一方で、Prophetが前提とするトレンドと季節性の加法的な構造は系列によっては実際のデータ生成過程と合わず、複雑な非線形パターンや外部特徴量との入り組んだ相互作用を捉える力は、ARIMA・状態空間モデルやGBDTに比べて限定的だとされています。手軽に試せる選択肢の1つと位置づけつつ、精度を追求する局面では他の手法とのバックテスト比較を必ず行うべき点は変わりません。

実務の結論：ベンチマークと検証で決める

本章で見てきた通り、機械学習は時系列予測において強力な選択肢であり、系列数の多さや外部特徴量の豊富さといった条件が揃えば、古典的な統計モデルを上回る精度を発揮する場面が少なくありません。しかし、系列が少なく短い、構造の解釈や予測区間が求められるといった条件では、古典的な統計モデルが依然として有力です。深層学習や時系列基盤モデルも、常に他の手法を上回るとは限らないという慎重な立場が、現時点での適切な理解だと考えられます。

手法を選ぶ実務上の結論は明快です。まず、前年同月や前週同曜日のように、系列の周期に合わせて過去の実測値をそのまま予測値とする単純なナイーブ予測や、ETSのような素朴な古典的手法を必ずベンチマークとして用意します。そのうえで、ARIMA・状態空間モデル・GBDT・必要であれば深層学習モデルといった複数のアプローチを、時間順を守ったバックテストで横並びに比較し、対象のデータで実際に精度と運用性に優れている手法を選び取ります。特定の手法への愛着や流行で選ぶのではなく、検証によって決めるという姿勢が、時系列予測を実務に組み込むうえでの一貫した原則になります。

最後に、ナイーブ予測とGBDTによる予測を同じ検証データで比較するコードを確認します。

from lightgbm import LGBMRegressor

from sklearn.metrics import mean_absolute_error

# 検証期間(直近90日)を分離する
train_df = df.iloc[:-90]
test_df = df.iloc[-90:]

feature_cols = [
    "lag_1", "lag_7", "lag_14",
    "rolling_mean_7", "rolling_mean_28", "rolling_std_7",
    "weekday", "month", "is_weekend", "is_month_end",
]

model = LGBMRegressor(max_depth=6, learning_rate=0.05, n_estimators=300, random_state=0)
model.fit(train_df[feature_cols], train_df["sales"])
pred_gbdt = model.predict(test_df[feature_cols])

# ナイーブ予測: 7日前の実測値をそのまま予測値とする(週次周期を考慮した季節ナイーブ)
pred_naive = test_df["lag_7"].values

mae_gbdt = mean_absolute_error(test_df["sales"], pred_gbdt)
mae_naive = mean_absolute_error(test_df["sales"], pred_naive)

print(f"GBDTのMAE: {mae_gbdt:.2f}")
print(f"季節ナイーブ予測のMAE: {mae_naive:.2f}")
print(f"改善率: {(1 - mae_gbdt / mae_naive) * 100:.1f}%")

ここでのpred_naiveは、7日前の実測値（lag_7列）をそのまま今日の予測値とみなす季節ナイーブ予測で、週次の周期性を持つデータのベンチマークとしてよく使われます。実行結果は、GBDTのMAE（平均絶対誤差）が4.38、季節ナイーブ予測が4.86で、改善率は9.9%でした。この合成データはトレンドと曜日効果に標準偏差5の独立なノイズを加えたものなので、どれほど良いモデルでもMAEはおよそ4.0が下限であり、9.9%という改善幅は、縮められる余地のおよそ半分を回収した水準にあたります。改善率の絶対値だけを見るのではなく、データが持つノイズの大きさから決まる「これ以上は縮まらない誤差」と照らして読むことが大切です。GBDTがナイーブ予測を上回る改善を示せなければ、特徴量エンジニアリングの見直しや、そもそも機械学習モデルを採用する妥当性の再検討が必要になります。逆に明確な改善が確認できれば、その改善幅を根拠に運用への採用を判断できます。単純な手法との比較を出発点に据える習慣が、時系列予測を実務で運用するうえでの安全網になります。

本章では、機械学習を時系列予測に適用する際の基本である教師あり学習への変換、GBDTとグローバルモデルの発想、時間順を守った検証設計とリークへの警戒、古典的手法と機械学習それぞれが優位になりやすい条件、深層学習・時系列基盤モデルの現在地とProphetの位置づけを確認しました。次章では「値そのものを予測する」という視点から離れ、時系列データにおける異常検知と変化点検知という、予測とは異なる角度からの活用法を扱います。

この章を深めたい方への参考書籍

『Kaggleで勝つデータ分析の技術』（門脇大輔ほか、技術評論社）：特徴量エンジニアリングとモデル検証の実践知を、コンペティションの現場で磨かれた具体的な手法とともに体系立てて解説した一冊です。本章で扱ったラグ・移動平均特徴量の設計や、時間順を守った検証設計、リークへの警戒といった論点を、より広い応用例とともに深めたい読者に向いています。

第11章 時系列の異常検知と変化点検知

機械学習編第11章「実務の壁を越える、不均衡データと異常検知」では、Isolation ForestやOne-Class SVM、LOFといった手法を取り上げ、正常データの分布から外れたサンプルを検出する考え方を扱いました。あの章の手法は、データの並び順を前提としない、いわば「1件1件が独立したサンプルの集まり」を対象にした異常検知です。これに対して本章が扱うのは、時間の順序そのものが意味を持つ時系列データにおける異常検知です。売上やアクセス数、センサー値のように時間とともに水準や季節パターンが変化していく系列では、「今この値が正常か異常か」を判断するために、その時点の値そのものだけでなく、トレンドや季節性を踏まえた「その時点であるべき値」との比較が欠かせません。本章では、この時系列特有の異常検知と、系列の構造そのものが切り替わる変化点の検知を、実務での運用設計まで含めて整理します。

時系列の異常の3類型

時系列データに現れる「異常」は、性質の異なる3つのパターンに大別できます。この分類を意識せずに1つの手法だけで対応しようとすると、検出できる異常とできない異常の偏りが生じやすくなります。

外れ値型は第6章・第8章で扱った予測区間との比較で比較的検出しやすい一方、変化点型は「異常な1点」が存在するわけではないため、点ごとの外れ値検出では見逃されやすいという特徴があります。パターン異常はさらに厄介で、各時点の値だけを見ても平均や分散に大きな変化がないことが多く、周期構造そのものに着目した分析が必要になります。この3類型のどれを主に検出したいのかを最初に整理しておくことが、手法選定の出発点になります。

なぜ時系列の異常検知は難しいか

時系列の異常検知が一般のサンプルデータの異常検知より難しい最大の理由は、「正常な値の範囲」そのものが時間とともに動き続ける点にあります。第1章・第2章で見たとおり、多くのビジネス系列にはトレンドと季節性が含まれており、ある時点で「大きい値」であることが、その系列にとって異常を意味するとは限りません。

この難しさを最も分かりやすく示すのが、系列全体の平均と標準偏差から3シグマ（平均から標準偏差の3倍以上離れた値を異常とみなす）を計算し、それを閾値として使う単純な方法の失敗です。年末商戦で売上が跳ね上がる小売の月次データを例に考えます。12月の売上は他の月に比べて明らかに高い水準になりますが、これは毎年繰り返される季節性であり、事業にとってはむしろ「あるべき正常な状態」です。ところが系列全体の平均・標準偏差から3シグマの閾値を機械的に計算すると、この12月の値が閾値を超えてしまい、正常な季節ピークが異常として誤検知されることになります。逆に、本当に注意すべき小さな変化、たとえば通常月における前年同月比のわずかな落ち込みは、季節変動の大きさに埋もれてしまい、3シグマの閾値には引っかからないまま見過ごされることもあります。

さらに第3章で扱った非定常性の問題も重なります。トレンドを持つ系列では平均そのものが時間とともに変化し続けるため、「全期間を通じた平均」という基準自体が意味を持ちにくくなります。したがって時系列の異常検知では、系列全体に対して一律の閾値を当てはめるのではなく、その時点のトレンドと季節性を踏まえた「その時点であるべき値」を何らかの形で推定し、そこからのズレを異常度として評価するというアプローチが基本になります。次節以降で扱う予測ベースのアプローチと分解ベースのアプローチは、いずれもこの考え方に立脚しています。

時系列の異常検知における最大の落とし穴は、季節性のある系列に対して素朴な3シグマ法をそのまま適用してしまうことです。季節ピークが異常として誤検知され、逆に季節変動に埋もれた本当の異常が見逃されます。「異常」とは値そのものの大小ではなく、トレンドと季節性を踏まえた期待値からのズレとして定義する必要があります。

基本アプローチ1：予測ベースの異常検知

予測ベースのアプローチは、第6章で扱ったARIMA・SARIMAや第8章・第9章で扱った状態空間モデル、あるいは第10章で扱った機械学習による時系列予測モデルを使って「あるべき値」を予測し、実際の観測値が予測区間から外れているかどうかで異常を判定する方法です。これらのモデルはトレンドと季節性を明示的に組み込んで予測を行うため、季節ピークのような正常な変動を異常と誤認しにくいという利点があります。

ARIMA・SARIMAモデルやSARIMAXモデルは、点予測に加えて予測区間（たとえば95%予測区間）を出力できます。実際の観測値がこの予測区間の外側に出た場合に異常の候補としてフラグを立てる、というのが最も直接的な実装です。状態空間モデルの場合はさらに好都合な性質があります。第8章・第9章で扱ったとおり、状態空間モデルはカルマンフィルタによって1期先の予測値と、その予測がどれだけ外れる可能性があるかを示す予測誤差の分散を、逐次的に計算します。この1期先予測誤差（イノベーション）を分散で標準化した値は、系列全体が定常か非定常かによらず、モデルが正しく機能していれば理論上は平均0・分散1に近い分布に従うため、リアルタイムの異常判定に使いやすい指標になります。

次のコードは、SARIMAXモデルで学習期間のデータから予測区間を計算し、その後のテスト期間の観測値が予測区間を外れているかどうかで異常候補を判定する例です。テスト期間には、意図的にスパイク型の異常を1点だけ混入させてあります。

import numpy as np
import pandas as pd
from statsmodels.tsa.statespace.sarimax import SARIMAX

rng = np.random.default_rng(0)
n = 200

# 季節性(周期12)とトレンドを持つ系列を生成し、後半に外れ値を混入
t = np.arange(n)
seasonal = 10 * np.sin(2 * np.pi * t / 12)
trend = 0.05 * t
noise = rng.normal(0, 1, n)
values = 50 + trend + seasonal + noise
values[150] += 25  # スパイク型の異常を1点だけ混入

series = pd.Series(values)
train, test = series[:140], series[140:]

model = SARIMAX(train, order=(1, 1, 1), seasonal_order=(1, 1, 1, 12))
fitted = model.fit(disp=False)

forecast = fitted.get_forecast(steps=len(test))
pred_mean = forecast.predicted_mean
conf_int = forecast.conf_int(alpha=0.05)

is_anomaly = (test.values < conf_int.iloc[:, 0].values) | (test.values > conf_int.iloc[:, 1].values)
anomaly_points = test.index[is_anomaly]
print(f"予測区間から外れた時点: {list(anomaly_points)}")

実行すると、異常候補として3時点（146、150、159）が挙がります。実際に混入させたスパイクは150の1点だけで、146と159は予測区間の下限をそれぞれ0.1程度下回っただけの誤検知です。95%予測区間を60時点に適用すれば、まったく正常なデータであっても平均して3時点前後は区間の外に出る計算になります。この誤検知は手法の欠陥ではなく、有意水準の設定そのものが持つ性質です。異常度の大きさ（150では予測区間の上限を20以上超えています）まで併せて見て、閾値ぎりぎりの検出と明確な逸脱を区別する運用が、実務では欠かせません。

このアプローチの弱点は、モデル自体の当てはまりが悪いと、モデルの誤差なのか本物の異常なのかの切り分けが難しくなる点です。モデルの次数選定や残差診断（第6章）を怠ったまま予測区間だけを機械的に使うと、モデルが季節パターンをうまく捉えられていない箇所を大量に異常として誤検知することがあります。予測ベースの異常検知を導入する前に、まずモデル自体の予測精度と残差の妥当性を確認しておく手順が欠かせません。

基本アプローチ2：分解ベースの異常検知

分解ベースのアプローチは、第2章で扱ったSTL分解（Seasonal-Trend decomposition using Loess）を使い、系列をトレンド・季節性・残差の3成分に分解したうえで、残差成分に対して外れ値検出をかける方法です。トレンドと季節性をあらかじめ取り除いているため、残差はおおむね平均0で安定した振れ幅を持つ成分になっており、単純な閾値による異常検知がしやすくなります。

残差に対して閾値を設定する際には、単純な平均・標準偏差ではなく、中央値と絶対偏差の中央値（MAD、Median Absolute Deviation）を用いた頑健な指標を使うことが推奨されます。残差の中に外れ値がすでに含まれている場合、通常の平均・標準偏差はその外れ値自体に強く引っ張られてしまい、閾値そのものが歪んでしまうためです。中央値とMADは外れ値の影響を受けにくいという性質があり、頑健なZスコア（Modified Z-score）として次のように計算されます。

import numpy as np
import pandas as pd
from statsmodels.tsa.seasonal import STL

# STL分解(前節のvaluesを再利用し、周期12で分解)
stl = STL(pd.Series(values, index=pd.RangeIndex(n)), period=12, robust=True)
result = stl.fit()
residual = result.resid

# 頑健なZスコア(中央値とMADを用いる)による異常判定
median = np.median(residual)
mad = np.median(np.abs(residual - median))
modified_z = 0.6745 * (residual - median) / mad

threshold = 3.5
anomaly_idx = np.where(np.abs(modified_z) > threshold)[0]
print(f"残差ベースで異常と判定された時点: {anomaly_idx.tolist()}")

0.6745という係数は、残差が正規分布に従うと仮定した場合にMADを標準偏差と同等のスケールに変換するための定数です。閾値を3.5前後に置くのは統計の分野でよく使われる経験的な目安であり、実務では自社データの分布を確認しながら調整する必要があります。分解ベースのアプローチは、STL分解自体がモデルのパラメータ推定を必要とせず計算が軽いため、多数の系列を一括で監視するような場面で扱いやすいという利点があります。一方で、STL分解は季節周期をあらかじめ固定した値として与える必要があるため、周期が不安定な系列や季節性がそもそも明確でない系列には向きません。

変化点検知：系列の構造が切り替わる時点をとらえる

予測ベース・分解ベースのアプローチは、いずれも1時点ごとの値が異常かどうかを判定するのに適していますが、変化点型の異常、つまり系列の平均やトレンドそのものが持続的に切り替わる現象の検出には別の考え方が必要です。変化点検知（Change Point Detection）は、この「系列の統計的な性質がどこで切り替わったか」を特定することを目的とした手法群です。

変化点検知の直感を最もつかみやすい古典的な手法が、累積和（CUSUM、Cumulative Sum）です。CUSUMの考え方は単純で、各時点の観測値とある基準値（たとえば直近の平均）との差を、プラス方向とマイナス方向それぞれについて累積していきます。系列が基準値の周りで安定して推移している間は、この累積和はゼロ付近を行き来するだけですが、ある時点を境に平均水準が上（または下）にシフトすると、差が同じ符号で積み重なり続けるため、累積和が一方向に大きく伸び始めます。この累積和があらかじめ決めた閾値を超えた時点を、変化が起きたシグナルとして検出します。CUSUMは1時点ごとの小さなブレには反応せず、持続的な水準シフトが積み重なって初めて反応するという性質を持つため、変化点型の異常を検出するのに適しています。

CUSUMをより柔軟にした考え方として、ベイズ的な変化点検知（Bayesian Online Change Point Detection）があります。これは各時点において「直近に変化点が起きてから何時点が経過しているか」の確率分布を逐次的に更新していく手法で、変化点が起きた確信度合いを確率として扱えるため、閾値を1つの数値で固定するCUSUMよりも柔軟な運用がしやすくなります。ただし実装や解釈にはある程度の統計的な素養が必要で、実務で手軽に使うにはハードルが高い面もあります。

実務でより手軽に使える選択肢が、rupturesのような変化点検知に特化したPythonライブラリです。こうしたライブラリでは、系列全体をあらかじめコストの観点から最適に分割するアルゴリズム（PELT法やBinary Segmentation法など）が用意されており、変化点の候補を自動的に列挙できます。次のコードは、平均水準が途中で切り替わる合成データに対して、rupturesのようなライブラリでPELT法を使い変化点を検出する例です。なおrupturesは、最新のPythonへの対応が本体のリリースより遅れることがあります（執筆時点の対応は3.9以上3.14未満）。導入できない場合は、対応済みのバージョンの仮想環境を用意して実行してください。

import numpy as np
import ruptures as rpt

rng = np.random.default_rng(1)

# 前半150点は平均0、後半150点は平均4にシフトする系列
segment1 = rng.normal(loc=0, scale=1.0, size=150)
segment2 = rng.normal(loc=4, scale=1.0, size=150)
signal = np.concatenate([segment1, segment2])

# PELT法(コスト関数はrbf)で変化点を検出。penは変化点の数に対する罰則の強さ
algo = rpt.Pelt(model="rbf").fit(signal)
change_points = algo.predict(pen=10)

print(f"検出された変化点(系列末尾を含む): {change_points}")

penパラメータは、変化点を1つ増やすごとに課されるペナルティの大きさを表しており、値を小さくするほど変化点が検出されやすくなる（検出感度が上がる一方で誤検知も増える）、値を大きくするほど検出されにくくなるという関係にあります。この値の選び方に唯一の正解はなく、対象とする業務データの性質を踏まえて試行錯誤で調整する必要があります。またPELT法は基本的にオフライン（過去データをまとめて分析する）手法であり、リアルタイムに逐次データが届く監視業務にそのまま使う場合は、一定期間ごとにバッチで再計算する運用や、CUSUMのような逐次型の手法との組み合わせを検討する必要があります。

リアルタイム監視の設計：検知の速さと誤報の少なさのトレードオフ

異常検知や変化点検知の手法をどれだけ精緻に実装しても、実際の監視ダッシュボードに組み込む段階で必ず直面するのが、検知の速さと誤報の少なさのトレードオフです。閾値を厳しく（異常や変化と判定しやすく）設定すれば、本当の異常や変化をより早く検知できるようになりますが、その分だけ通常の変動でもアラートが鳴りやすくなり、誤報が増えます。逆に閾値を緩めれば誤報は減りますが、実際に異常が発生してから検知されるまでの遅れが大きくなり、対応が後手に回るリスクが高まります。

この閾値の調整は、CUSUMであれば累積和の閾値、予測ベースの手法であれば予測区間の幅（何%予測区間を使うか）、rupturesのようなライブラリであればpenパラメータというように、手法ごとに異なるパラメータとして現れますが、いずれも本質的には「検知の遅れ」と「誤報率」のバランスを調整しているという点で共通しています。この調整は統計的な最適解が一意に決まるものではなく、見逃した場合の被害の大きさと、誤報1件あたりの対応コストを業務側と擦り合わせながら決めていく必要があります。

誤報が多い監視システムが実務で陥りやすい問題が、アラート疲れ（Alert Fatigue）です。閾値を厳しくしすぎて日常的に大量のアラートが鳴る状態が続くと、現場の担当者は次第にアラートを確認する優先度を下げるようになり、最終的には本当に重要なアラートすら見過ごされるようになります。これは検知精度の問題というより運用設計の問題であり、次のような工夫で緩和できます。

• アラートを重大度で段階分けする（即座に対応が必要な緊急アラートと、日次でまとめて確認すればよい注意アラートを区別する）
• 同じ異常が連続して検知され続ける場合は、最初の1回だけ通知し、以降は状態が続いていることをダッシュボード上で表示するに留める
• 季節性や事業環境の変化に応じて、閾値を定期的（四半期に1回など）に見直すサイクルを運用に組み込む
• 異常の判定結果と、実際にその後何が起きたか（誤報だったか、本物の異常だったか）を記録し、閾値調整の材料として蓄積する

異常検知は最終判断ではなく候補の絞り込みである

ここまで扱ってきた予測ベース・分解ベース・変化点検知のいずれの手法も、統計的な基準に基づいて「異常らしさ」をスコア化し、閾値を超えたものを候補として挙げる仕組みです。どれだけ手法を洗練させても、統計的な閾値だけで異常か正常かを完全に、かつ自動的に判定しきることはできません。閾値の設定自体にトレードオフが伴う以上、一定の誤報と見逃しは構造的に避けられないためです。

したがって実務における異常検知システムの正しい位置づけは、最終判断を下す装置ではなく、人間が確認すべき候補を絞り込む一次スクリーニングの装置です。異常検知モデルが候補としてフラグを立てた時点を、業務のドメイン知識を持つ担当者が確認し、実際に対応が必要な事象なのか、それとも既知のキャンペーンやメンテナンス作業による説明可能な変動なのかを判断する、という運用体制とセットで設計する必要があります。この確認のプロセスをあらかじめ業務フローに組み込んでおくことで、異常検知モデル単体の精度に過度に依存しない、頑健な監視の仕組みを作ることができます。

時系列の異常検知・変化点検知は、統計的な閾値によって「確認すべき候補」を絞り込む一次スクリーニングとして位置づけるのが実務上の基本方針です。最終的な判断は、業務のドメイン知識を持つ担当者による確認プロセスと組み合わせ、異常判定の結果とその後の実際の顛末を記録して閾値調整に生かすサイクルを、運用設計にあらかじめ組み込んでおくことが重要です。

ビジネス適用例の整理

時系列の異常検知・変化点検知が使われる代表的な業務領域ごとに、主に検出したい異常の類型や、適した手法の傾向を整理します。

この整理から分かるとおり、どの領域でも単一の手法だけで完結させるのではなく、複数のアプローチを組み合わせ、かつ最終的には人間による確認プロセスを経て運用するという構図が共通しています。第12章では、こうした異常検知の考え方も踏まえつつ、需要予測を実務のシステムとして設計する際に必要となる、より広い視点での実務設計を扱います。

この章を深めたい方への参考書籍

『入門 機械学習による異常検知 Rによる実践ガイド』（井手剛、コロナ社）：ホテリング理論から始まり、時系列データにおける変化点検知やベイズ的な異常検知の考え方まで、数理的な背景を丁寧に解説した書籍です。本章で直感的に紹介したCUSUMやベイズ的変化点検知の理論的な裏付けを深めたい読者に適しています。実装はRですが、考え方自体はPythonでの実務にもそのまま応用できます。

第12章 需要予測を業務に組み込む実務設計

前章では、異常検知と変化点検知という、予測とは少し異なる角度から時系列データの変化を捉える手法を扱いました。本章は最終章として、第5章から第10章で扱ったARモデル・SARIMA・VARモデル・状態空間モデル・機械学習による予測手法を、実際の需要予測業務にどう組み込むかという実務設計に主眼を置きます。

最初に押さえておきたいのが、「精度の高い予測モデルを作れること」と「需要予測が業務で機能すること」は別の問題だという点です。モデルの当てはまりの良さを競うだけでは、評価・検証・意思決定・運用という設計が伴わない限り、優れたモデルも使われなくなります。本章では評価指標、ベンチマークとバックテスト、階層予測の整合性、間欠需要、予測から意思決定への接続、運用体制の設計という順に整理し、本ガイド全体を締めくくります。

予測モデルが「作れる」ことと「機能する」ことの違い

需要予測のプロジェクトが行き詰まる経緯には共通点があります。高精度なモデルを構築し社内発表では良好な検証結果が示されるものの、実際の発注業務に組み込む段階で「この数字をどう使えばいいのか分からない」といった声が上がり、結局は経験則ベースの見積もりに戻ってしまうというものです。背景には、モデル構築そのものよりも運用の設計が不足しているという事情があります。

需要予測が業務で機能するには、少なくとも4つの設計が必要です。評価指標の設計、ベンチマークとバックテストの設計、複数粒度の予測を整合させる階層設計、点予測を在庫・人員配置に落とし込む接続設計です。本章はこの4つを順に扱ったうえで、最後に運用体制の設計を扱います。

予測精度の指標: MAE・RMSE・MAPEの使い分けと限界

需要予測の精度指標として広く使われるのが、MAE（Mean Absolute Error、平均絶対誤差）、RMSE（Root Mean Squared Error、二乗平均平方根誤差）、MAPE（Mean Absolute Percentage Error、平均絶対パーセント誤差）の3つです。性質を理解せずに使うと、評価そのものが実態を見誤らせます。

この中でもっとも注意が必要なのがMAPEです。誤差を実績値で割るため商品間で比較しやすく、「予測精度85%」のように経営層へ説明しやすい一方、実績値がゼロに近い商品では計算が破綻します。実績2個・予測5個なら誤差率は150%、実績がゼロの月は定義できません（ゼロ除算）。低回転品や新商品を多く抱える需要予測では、こうした少数の低需要品の極端な誤差率にMAPE全体が引きずられ、大部分の商品が当たっていても指標が著しく悪化して見える歪みが生じます。

wMAPE・WAPEとバイアスの監視

この歪みを避けるために実務でよく使われるのが、wMAPE（weighted MAPE）、あるいはWAPE（Weighted Absolute Percentage Error）です。各商品の誤差率を単純平均するのではなく、誤差の絶対値の合計を実績値の合計で割ることで、需要量の大きい商品の誤差により大きな重みを与えます。

3つのSKU、実績100個・予測90個（誤差10個）、実績10個・予測5個（誤差5個）、実績2個・予測6個（誤差4個）を例にとります。MAPEは各SKUの誤差率10%・50%・200%を単純平均するため約87%となり、需要量がわずか2個のSKUの極端な誤差率に評価が引きずられます。wMAPEは誤差の合計19個を実績の合計112個で割るためおよそ17%となり、事業全体の量に対してどれだけ外れているかという、意思決定に近い数字が得られます。

精度と並んで見落とされがちなのが、バイアス（系統的な過大・過小予測の傾向）の監視です。MAE・RMSE・MAPEは誤差の絶対値を扱うため、実績より多い外れと少ない外れを区別せず、一定方向へ偏り続ける傾向を見逃します。バイアスは誤差の符号を残した合計を実績の合計で割ることで計算でき、先の例では予測と実績の差の合計がマイナス11個、これを112個で割るとおよそマイナス9.8%になります。実績を1割弱、系統的に少なく見積もる傾向があると分かり、精度とは別の角度から在庫不足のリスクを示唆する手がかりになります。

精度指標（MAE・RMSE・MAPE・wMAPE）は「どれだけ外れているか」しか教えてくれません。「どちら向きに系統的に外れ続けているか」を知るにはバイアスの監視が別途必要です。精度が良好でもバイアスが一方向に偏っている予測は、欠品や過剰在庫を継続的に生み出している可能性があり、実務では精度とバイアスの両方を毎期並べて確認する運用が望まれます。

次のコードは、wMAPEとバイアスを計算する関数の実装例です。ゼロ除算を避けるため、実績の合計がちょうどゼロになる場合を分けて扱っています。

import numpy as np

def wmape(actual, forecast):
    """重み付けした絶対パーセント誤差(wMAPE)を計算する"""
    actual = np.asarray(actual, dtype=float)
    forecast = np.asarray(forecast, dtype=float)
    denom = np.sum(np.abs(actual))
    if denom == 0:
        return np.nan
    return np.sum(np.abs(actual - forecast)) / denom

def forecast_bias(actual, forecast):
    """系統的な過大・過小予測の傾向(バイアス)を計算する。正なら過大、負なら過小"""
    actual = np.asarray(actual, dtype=float)
    forecast = np.asarray(forecast, dtype=float)
    denom = np.sum(actual)
    if denom == 0:
        return np.nan
    return np.sum(forecast - actual) / denom

actual = np.array([100, 10, 2])
forecast = np.array([90, 5, 6])

print(f"wMAPE: {wmape(actual, forecast):.1%}")
print(f"バイアス: {forecast_bias(actual, forecast):+.1%}")

実行すると、wMAPEはおよそ17.0%、バイアスはおよそマイナス9.8%となり、手計算の数値と一致します。実務では、この2つの関数をSKU単位・カテゴリ単位・全社単位の各粒度で継続的に計算し、精度とバイアスを併記して毎期確認する運用が有効です。

まず問うべきはナイーブ予測に勝てているか

「予測精度85%」のような単独の数字だけでは、予測の良し悪しを判断できません。需要が安定した商品であれば平凡な水準かもしれませんし、変動が激しい商品であれば優れた水準かもしれません。比較対象がなければ、単独の数字は意味を持ちません。

この比較対象として実務でまず設定すべきなのが、ナイーブ予測（素朴な予測）によるベンチマークです。直前期の実績をそのまま使う「前期値予測」、季節性のある商品について1年前の同じ月の実績を使う「前年同月予測」が代表例です。手の込んだモデルを構築する前に、こうした素朴な予測に対しどれだけの改善を示せているかを確認する必要があり、高度なモデルのwMAPEが前期値予測とほとんど変わらない、あるいは劣っているなら、そのモデルは実務上の価値をまだ生み出せていないと判断すべきです。

改善率はMASE（Mean Absolute Scaled Error）で標準化して表すこともできます。これは改善率そのものではなく、ナイーブ予測の誤差を1としたときの倍率で自社モデルの誤差を表す指標で、1を下回れば単純な予測より優れていると読めます。ただし実務では、自社モデルとナイーブ予測のwMAPEを並べて「ナイーブ予測比でどれだけ誤差を削減できたか」を報告する運用でも十分に目的を果たします。この比較を欠いた「精度85%を達成しました」という報告は、実務上ほとんど意味を持たないという認識を関係者全員で共有しておくことが重要です。

バックテストの設計: ローリングオリジン評価

第1章では、時系列データにランダム分割の交差検証を適用してはならないという原則を確認しました。この原則を具体化する標準的な手法が、ローリングオリジン評価（rolling origin evaluation、時間起点を進めながら繰り返し検証する評価）です。学習データの終端（オリジン、起点）を少しずつ未来へずらしながらその都度モデルを学習し直し、それぞれの起点から一定期間先までを予測して誤差を記録することを繰り返し、時期による精度のばらつきを評価します。

設計するうえで欠かせないのが、実際の業務で使うリードタイム（発注から納品までの期間など）に合わせて「何期先の予測か」を分けて評価することです。1期先（来月）の予測は当たりやすくても、3期先や6期先になると精度が大きく落ちるモデルは珍しくありません。リードタイムが3か月であれば業務で使われるのは3期先の予測精度であり、1期先だけを見て「優秀だ」と判断すると、実運用への投入後に想定外の精度低下に直面します。

検証期間の選び方にも注意が必要です。直近の数か月だけでは、たまたま需要が安定していた時期の精度しか確認できません。平常期だけでなくキャンペーン期や季節の谷間、需要が大きく変動した過去の異常期間も含めることで、モデルが平時にしか通用しないのか、変動が大きい局面でも頑健性を持つのかを見極める必要があります。

次のコードは、月次データに対しローリングオリジン評価を行う関数の実装例です。学習期間を1か月ずつ広げながら、指定した期間先（horizon）の予測誤差を、任意の予測関数forecast_funcに対して評価します。

import numpy as np
import pandas as pd

def rolling_origin_backtest(series, forecast_func, min_train_size, horizon, step=1):
    """
    series: 予測対象の時系列(pandas.Series、時間順に並んでいること)
    forecast_func: 学習データとhorizonを受け取り、horizon期先までの予測値を返す関数
    min_train_size: 最初の学習データの最小サイズ
    horizon: 何期先まで予測するか、step: オリジンを何期ずつ進めるか
    """
    records = []
    n = len(series)

    for train_end in range(min_train_size, n - horizon + 1, step):
        train = series.iloc[:train_end]
        actual = series.iloc[train_end:train_end + horizon].values
        pred = forecast_func(train, horizon)

        for h in range(horizon):
            records.append({
                "origin": series.index[train_end - 1],
                "horizon": h + 1,
                "actual": actual[h],
                "forecast": pred[h],
            })

    return pd.DataFrame(records)

def naive_seasonal_forecast(train, horizon, season_length=12):
    """前年同月をそのまま予測値として使うナイーブ予測"""
    last_season = train.iloc[-season_length:].values
    reps = int(np.ceil(horizon / season_length))
    return np.tile(last_season, reps)[:horizon]

# result = rolling_origin_backtest(monthly_sales, naive_seasonal_forecast,
#                                   min_train_size=36, horizon=3, step=1)
# result.groupby("horizon").apply(lambda d: wmape(d["actual"], d["forecast"]))

この関数はオリジンごと・何期先ごとに実績と予測を記録したデータフレームを返すため、horizon（何期先か）ごとにグループ化してwMAPEを計算すれば、リードタイム別の精度劣化を定量的に把握できます。forecast_func部分に、SARIMAや状態空間モデル、機械学習モデルの学習・予測処理を差し込むことで、任意のモデルに対して同じ枠組みでバックテストを行えます。

階層予測の整合性: トップダウン・ボトムアップ・リコンシリエーション

需要予測の多くは、全社合計・カテゴリ別・SKU別といった複数の粒度で必要とされます。各粒度で個別に予測を作ると、SKU単位の予測を合計してもカテゴリ予測と一致せず、カテゴリ予測を合計しても全社予測と一致しないという整合性の問題がしばしば生じます。

対処法は大きく3つです。トップダウン方式は全社合計の予測を各カテゴリ・SKUの過去の構成比で按分する方法で、全体傾向は捉えやすい一方SKU固有の動きに弱く、ボトムアップ方式はSKU単位の予測を積み上げる方法で、個々の動きは捉えやすい一方誤差が積み重なりブレが大きくなりやすい弱点があります。リコンシリエーション（reconciliation、整合化）は、各粒度でいったん作った予測に、階層構造全体を考慮した調整を事後的に加える方法です。

近年はリコンシリエーションを体系化した手法（各粒度の予測誤差の大きさを踏まえて重みを決めるMinTなど）が普及し、Pythonでもhierarchicalforecastのようなライブラリで利用しやすくなっています。重要なのは、複数粒度の予測が矛盾しないことを需要予測の業務組み込み段階であらかじめ設計しておくことです。

間欠需要とCrostonの方法

ここまでの議論は、ある程度まとまった需要が継続的に発生する商品を念頭に置いていました。しかしスペアパーツや保守部品、稀にしか出荷されない特注品のように、大半の期間で需要がゼロで時折まとまった数量の注文が入る間欠需要（intermittent demand）を抱える商品では事情が異なります。

間欠需要には、ARIMAや状態空間モデル、多くの機械学習モデルがそのままでは適しません。需要が連続分布に従うことを前提とする手法をゼロが大半のデータに当てはめると、平坦な予測しか出せなかったり、少数の非ゼロ実績に過剰反応して不安定になったりします。MAPEも実績ゼロの月が頻発すると計算自体が破綻します。

間欠需要向けの代表的な手法がCrostonの方法です。需要を「非ゼロの注文が発生する間隔」と「発生時の数量」という2つの系列に分解し、それぞれに指数平滑法を適用したうえで、間隔と数量の比から単位期間あたりの需要を推定します。需要が発生する頻度自体を明示的にモデル化する点が通常の時系列モデルとの違いです。ゼロが極端に多い場合のバイアスを補正したSBA（Syntetos-Boylan Approximation）など改良版も提案されています。間欠需要を抱える商品群は通常の需要とあらかじめ切り分け、Croston系の手法を用いるといった専用の設計が実務上重要です。

点予測から意思決定へ: サービスレベルと安全在庫

需要予測の目的は、予測値そのものを当てることではなく、それを踏まえて在庫水準や人員配置の意思決定を行うことにあります。必要なのは単一の点予測（来月の需要は100個、という1つの数字）だけではなく、その予測がどの程度ばらつきうるかという分布や区間の情報です。

この考え方が具体化するのが、在庫管理のサービスレベルと安全在庫の設計です。サービスレベルとは欠品を起こさず需要を満たせる確率の目標値であり、安全在庫とは需要や納期のばらつきに備えて上乗せして保持する在庫です。安全在庫は需要のばらつき（標準偏差）、リードタイムの長さ、目標サービスレベルに応じた係数（z値）の3つから計算されます。

たとえば、日次需要の標準偏差が20個、リードタイムが9日間の商品について、サービスレベル95%を目標とする場合、安全在庫はz値の1.65に標準偏差の20個とリードタイムの平方根（9日の平方根である3）を掛け合わせた値、1.65×20×3で99個程度と求められます。目標サービスレベルを99%に引き上げると、z値は2.33に上がり必要な安全在庫はおよそ140個まで増加します。サービスレベルを高く設定するほど在庫コストは非線形的に増えていきます。

ここで重要なのは、この標準偏差が単なる需要の変動の大きさではなく、予測モデルが残す誤差の分布から得られるという点です。予測精度が高いモデルほど誤差の標準偏差は小さくなり、同じサービスレベルを達成するために必要な安全在庫も少なく済みます。予測精度の改善は、在庫コストの削減に直接つながる実務上の効果を持ちます。

予測の外れ方が生む非対称なコスト

多くの指標は誤差の絶対値を扱うため、実績より多く予測した場合と少なく予測した場合を対等に扱います。しかし予測の外れ方が引き起こすコストは多くの場合対称ではありません。予測が実績を下回れば欠品が生じ、販売機会の損失や顧客の信頼低下につながります。予測が実績を上回れば過剰在庫が生じ、保管コストの増加や廃棄・値引き販売の損失につながります。この2つのコストの大きさは商品の性質によって大きく異なります。

鮮度が短い食品や流行に左右されるファッション商品では過剰在庫のコストが欠品のコストを上回ることが多く、予測をやや控えめに立てる方が合理的です。反対に医療関連商品や、欠品が生産ラインの停止に直結する工場向け部品では欠品のコストがはるかに大きく、在庫過剰を許容してでも欠品を避ける方向に設定すべきです。

この非対称性を組み込む簡便な方法は、サービスレベルの目標値を商品カテゴリごとに変え、欠品コストが高い商品には高いサービスレベルを設定することです。より踏み込んだ方法として、非対称な損失関数（分位点回帰のピンボール損失のように上振れと下振れに異なる重みを与える損失関数）でモデルを学習させ、予測値を意図的に偏らせるアプローチもあります。精度指標の改善だけでなく、この非対称性を評価に織り込む姿勢が欠かせません。

運用設計: 再学習・ドリフト監視・判断的調整

需要予測のモデルは構築して終わりではなく、継続的な運用のなかで精度を維持していく必要があります。運用設計で押さえるべき論点は大きく3つあります。

1つ目は再学習の頻度です。需要のパターンは市場環境や競合動向、商品ライフサイクルの変化にともなって変わり続けるため、構築時点のパラメータを使い続けると実態との乖離が広がります。変化の速い商材では月次、比較的安定した商材では四半期ごとといった具合に、商品群の性質に応じて頻度を変える設計が現実的です。

2つ目はモデル劣化（ドリフト）の監視です。第11章の変化点検知の考え方は、予測モデルの誤差の系列を監視する目的にも応用できます。期ごとのwMAPEやバイアスを記録し、CUSUM（累積和）のような手法で傾向的なずれが生じていないかを定期的に確認します。直近1期だけを見ていると、ドリフトの発見が遅れます。

3つ目は、人間による判断的調整（judgmental adjustment）の扱いです。営業担当者などが、競合の新商品投入や大口顧客からの内示といったモデルには見えない情報を踏まえ、予測値を手動で調整することが一般的に行われています。功の面ではモデルが持たない未知の要因を補える一方、罪の面では担当者の楽観・悲観が系統的バイアスとして混入しやすく、根拠が記録されないまま繰り返されると後になって説明できなくなります。

• 調整前の予測値と、人間が調整した後の予測値の両方を必ず記録として残す。
• 調整を行った理由（根拠となった情報）を簡潔にでもログとして残す。
• 期末には調整前後の予測のどちらが実績に近かったかを比較し、判断的調整が付加価値を生んでいるか（Forecast Value Added）を定期的に検証する。

調整前の予測値を残さずモデルの出力を上書きしてしまう運用は、判断的調整が本当に精度向上に寄与しているのかを事後的に検証する手段を失わせます。判断的調整そのものを禁止する必要はありませんが、調整前後の両方を記録しその効果を継続的に検証する仕組みを運用設計に組み込むことが重要です。

予測精度のレビュー会議をどう設計するか

運用を継続する最後の仕組みが、予測精度を定期的に振り返るレビュー会議です。営業・生産・在庫管理といった予測を使う部門を交え月次や四半期ごとに開くことが望ましく、前期の精度とバイアスの確認、誤差が大きかった商品の原因共有、判断的調整の効果検証、次期への変更点の議論という流れで設計します。

重要なのは、会議を「なぜ外れたかの犯人探し」の場にしないことです。予測には本質的に不確実性が伴い、どれほど精緻なモデルでも一定の誤差は避けられません。目的は誤差をゼロに近づけることではなく、誤差の傾向（系統的なバイアスの有無、特定の商品群や期間への偏り、判断的調整の効果）を継続的に把握し、運用ルールを少しずつ改善していくことにあります。

本ガイドの締めくくり: 理解のための時系列分析と予測のための時系列分析

本ガイドは時系列データの性質から出発し、第1章では自己相関がもたらす落とし穴を、第2章ではトレンド・季節性・残差への分解を、第3章では定常性と単位根という概念を、第4章では見せかけの回帰と共和分を扱いました。この4章は、系列そのものの構造を正しく理解するための土台であり、「理解のための時系列分析」に位置づけられます。

続く第5章から第10章では、ARモデル・MAモデルからARIMA・SARIMA、VARモデルとGranger因果、状態空間モデル、機械学習による時系列予測までを扱いました。これらは理解した構造をもとに将来の値を見積もる「予測のための時系列分析」の中核をなす手法群です。第11章の異常検知と変化点検知は、この2つの流れを橋渡しします。構造の変化を捉える技術は、系列を理解する道具であると同時に、予測モデルの劣化を監視する道具でもあるからです。

そして本章では、これらの技術を土台に、需要予測を業務で機能させる実務設計、すなわち評価指標、ベンチマークとバックテスト、階層予測の整合性、間欠需要への対応、意思決定への接続、運用体制の設計を扱いました。系列の構造を正しく理解する力と、その理解を予測という形に落とし込み、組織の意思決定に接続し運用し続ける力の両方が、時系列データをビジネスの成果に結びつける土台になります。

時系列分析は、系列の構造を正しく理解する「理解のための分析」と、その理解を将来の値の見積もりへとつなげる「予測のための分析」という2つの目的の往復から成り立っています。需要予測を業務で機能させるには、精度の高いモデルを作ることに加えて、評価指標・ベンチマーク・階層整合性・意思決定への接続・運用体制という一連の実務設計が欠かせません。

この章を深めたい方への参考書籍

『この1冊ですべてわかる 需要予測の基本』（山口雄大、日本実業出版社）：需要予測を統計的な手法の話に留めず、組織の中でどう機能させるかという実務の視点から丁寧に整理した一冊です。本章で扱った評価指標の使い分けやサービスレベル・安全在庫の考え方、判断的調整の扱いについて、より実務に即した事例とともに理解を深めたい読者に向いています。

終章 時系列分析を経営の武器にするために

12章で積み上げてきたこと

本ガイドでは、時系列データが持つ固有の性質（自己相関・トレンド・季節性・非定常性）の理解から出発し、見せかけの回帰という代表的な罠の見破り方、ARIMAやVARといった古典的なモデル体系、見えない状態を推定する状態空間モデル、機械学習による予測との使い分け、異常検知、そして需要予測を業務に組み込む実務設計までを一続きに見てきました。

序章で掲げた2つの軸を振り返ります。第一に、時間の順序がもたらす前提の崩れを避けること。自己相関のあるデータに独立前提の手法を使わない、非定常な系列同士を安易に回帰しない、検証は必ず時間順を守る。この規律を守るだけで、時系列分析でつまずく原因の大半は防げます。第二に、予測を点ではなく幅で捉え、検証と運用まで設計すること。どんなモデルを使うにせよ、ナイーブなベンチマークとの比較、リードタイムに合わせたバックテスト、導入後の劣化監視までを一続きの作業として設計することが重要です。

手法選びの見取り図

本ガイドで扱った手法群を、実務での選び方として整理します。

• 構造を理解したい・成分に分けたい：STL分解（第2章）、状態空間モデル（第8〜9章）。「季節のせいか実力か」「トレンドは転換したか」に答えます
• 単一〜少数の系列を予測したい：SARIMA・ETS（第6章）を軸に、状態空間モデル（第8章）も候補に。予測区間が自然に得られます
• 大量の系列をまとめて予測したい・特徴量が豊富：GBDTによるグローバルモデル（第10章）。検証設計とリーク対策が生命線です
• 系列同士の関係を知りたい：VARとGranger因果（第7章）。ただし因果の主張には因果推論の設計（続刊）が必要です
• 異常を見つけたい：予測ベース・分解ベースの検知と変化点検知（第11章）。誤報率と検知の速さのバランス設計が本体です

どの道具を選ぶにしても、出発点は「そのデータはどう生成されているか」という問いです。トレンドはあるか、季節性はあるか、ショックは残り続けるか、構造は安定しているか。データの性質の見立てが正しければ、手法選びの大半は自然に決まります。

次に学ぶこと

本ガイドの先には、次の領域が広がっています。

• 因果推論と効果検証編（本シリーズ続刊）：第9章で触れた介入効果の測定を、反実仮想・実験設計の枠組みから本格的に扱います
• ベイズ統計編（計画中）：状態空間モデルをベイズ推定で扱い、階層構造や非ガウスなモデルへ拡張します
• 既刊「ビジネスの要因分析に効く回帰分析と統計モデリングの実践ガイド」：本ガイドの回帰まわりの前提（係数解釈・診断）はこちらで扱っています
• 既刊「機械学習の仕組みと使いどころがわかる実践ガイド」：第10章で使ったGBDTや検証設計の基礎はこちらをご覧ください

Anagraftでは、AIプロジェクトの構想・課題設計から、データ分析・機械学習モデルの開発、AI人材の育成まで一貫したご支援を行っています。ご相談は、以下よりお問い合わせください。

お問い合わせ

本コラムで紹介した書籍一覧

各章末でご紹介した参考書籍の一覧です。学習の段階に合わせてお選びください。

• 『現場ですぐ使える時系列データ分析』（横内大介・青木義充、技術評論社）：ビジネスデータを題材に時系列分析の入口を丁寧に案内する実務寄りの入門書です
• 『基礎からわかる時系列分析』（萩原淳一郎・瓜生真也・牧山幸史、技術評論社）：分解から状態空間モデルまでをRで実践しながら学べる、体系性の高い一冊です
• 『計量経済学』（西山慶彦・新谷元嗣・川口大司・奥井亮、有斐閣）：単位根や時系列回帰を含む計量経済学の現代的な標準テキストです
• 『経済・ファイナンスデータの計量時系列分析』（沖本竜義、朝倉書店）：定常性・単位根・共和分・VARを厳密かつ簡潔に学べる、この分野の定番教科書です
• 『時系列解析 自己回帰型モデル・状態空間モデル・異常検知』（島田直希、共立出版）：ARモデルから状態空間・異常検知までをコンパクトに橋渡しする現代的な入門書です
• 『Forecasting: Principles and Practice（第3版）』（Rob J. Hyndman、George Athanasopoulos、Web公開・英語）：予測実務の世界標準テキスト。Webで無料公開されており、ETSや階層予測の解説が充実しています
• 『実践 時系列解析』（Aileen Nielsen、オライリー・ジャパン）：前処理から機械学習・深層学習まで、時系列実務の幅広い論点をカバーする一冊です
• 『時系列分析と状態空間モデルの基礎 RとStanで学ぶ理論と実装』（馬場真哉、プレアデス出版）：状態空間モデルの考え方を平易な言葉で積み上げる、入門者に評判の良い一冊です
• 『カルマンフィルタ Rを使った時系列予測と状態空間モデル』（野村俊一、共立出版）：カルマンフィルタの仕組みと応用を集中的に学べる統計学One Pointシリーズの一冊です
• 『Kaggleで勝つデータ分析の技術』（門脇大輔ほか、技術評論社）：時系列データの特徴量作成とリークしない検証設計の実践知が詰まっています
• 『入門 機械学習による異常検知 Rによる実践ガイド』（井手剛、コロナ社）：異常検知の理論を体系的に学べる定番書。時系列以外の異常検知にも視野が広がります
• 『この1冊ですべてわかる 需要予測の基本』（山口雄大、日本実業出版社）：需要予測をオペレーションと組織の視点から扱う実務書。モデルの外側の設計を補強できます

参考情報・出典

• statsmodels公式ドキュメント（Time Series Analysis） https://www.statsmodels.org/stable/tsa.html
• pandas公式ドキュメント（Time series） https://pandas.pydata.org/docs/user_guide/timeseries.html
• Rob J. Hyndman, George Athanasopoulos “Forecasting: Principles and Practice （3rd ed.）” https://otexts.com/fpp3/
• scikit-learn公式ドキュメント（Cross validation of time series data） https://scikit-learn.org/stable/modules/cross_validation.html
• LightGBM公式ドキュメント https://lightgbm.readthedocs.io/
• ruptures公式ドキュメント（変化点検知） https://centre-borelli.github.io/ruptures-docs/
• S. Makridakis, E. Spiliotis, V. Assimakopoulos “The M4 Competition: Results, findings, conclusion and way forward” International Journal of Forecasting, 34（4）, 2018, pp. 802-808（第10章の手法比較の根拠）
• S. Makridakis, E. Spiliotis, V. Assimakopoulos “M5 accuracy competition: Results, findings, and conclusions” International Journal of Forecasting, 38（4）, 2022, pp. 1346-1364（第10章の手法比較の根拠）