時系列の周期性を発見する：「うちは波があるよね」を仕入計画に変える分析手法

import pandas as pdfrom statsmodels.tsa.seasonal import STL
# 日次売上データを読み込みsales = pd.read_csv("daily_sales.csv", parse_dates=["date"], index_col="date")["amount"]
# STL: period=7 で「曜日周期」を抽出 (月次なら 30, 年次なら 365)stl = STL(sales, period=7, robust=True).fit()
# trend / seasonal / resid の 3 系列に分解されたprint("Trend:    ", stl.trend.tail(3).values)print("Seasonal: ", stl.seasonal.tail(3).values)print("Residual: ", stl.resid.tail(3).values)
# 4 つのチャートで一望stl.plot()

import matplotlib.pyplot as pltfrom statsmodels.graphics.tsaplots import plot_acf, plot_pacf
fig, axes = plt.subplots(2, 1, figsize=(10, 6))plot_acf(sales, lags=60, ax=axes[0])plot_pacf(sales, lags=60, ax=axes[1])axes[0].set_title("ACF: 周期がある場合、ラグ N でピークが立つ")axes[1].set_title("PACF: 単独ラグの寄与だけを抜き出した版")plt.tight_layout()plt.show()
# プログラム的にピークを抽出（青の信頼区間外のラグ）from statsmodels.tsa.stattools import acfimport numpy as np
acf_vals, confint = acf(sales, nlags=60, alpha=0.05, fft=True)significant_lags = [    lag for lag in range(1, len(acf_vals))    if abs(acf_vals[lag]) > (confint[lag, 1] - acf_vals[lag])]print("有意な周期ラグ候補:", significant_lags[:10])# 例) [7, 14, 21, 28, 30, 60]  → 曜日周期 7 と月次周期 30 の両方を検出

import numpy as npimport matplotlib.pyplot as plt
# 平均を引いてから FFT (DC 成分を除去)y = sales.values - sales.mean()N = len(y)freqs = np.fft.rfftfreq(N, d=1.0)        # サンプリング間隔=1日spectrum = np.abs(np.fft.rfft(y))
# 周期 = 1 / 周波数periods = np.where(freqs > 0, 1 / freqs, np.inf)
# パワーが強い順に上位5周期top = np.argsort(spectrum)[::-1][:5]for i in top:    if 2 <= periods[i] <= N / 2:        print(f"周期 {periods[i]:.1f} 日 (強さ {spectrum[i]:.0f})")
# 可視化（横軸を周期に取ると経営層に説明しやすい）plt.figure(figsize=(10, 4))plt.plot(periods[1:N//2], spectrum[1:N//2])plt.xscale("log"); plt.xlabel("周期 (日)"); plt.ylabel("強さ")plt.title("周期スペクトル: 山の位置が「効いている周期」")plt.show()

from scipy.stats import kruskal
df = sales.to_frame("amount")df["weekday"] = df.index.dayofweek  # 0=月 ... 6=日
# 曜日ごとの売上分布が等しい (= 曜日周期なし) という帰無仮説を検定groups = [df[df["weekday"] == d]["amount"].values for d in range(7)]stat, p = kruskal(*groups)print(f"Kruskal-Wallis statistic={stat:.2f}, p={p:.4g}")if p < 0.05:    print("→ 曜日によって売上分布が有意に異なる = 曜日周期あり")else:    print("→ 曜日周期は統計的には有意でない")

import numpy as npimport pandas as pd
np.random.seed(42)days = 365 * 2t = np.arange(days)trend = 1000 + 0.5 * t                                     # 長期成長weekly = 200 * np.sin(2 * np.pi * t / 7)                   # 曜日効果monthly = 150 * np.sin(2 * np.pi * t / 28)                 # 給料日 28 日周期noise = np.random.normal(0, 80, days)
sales = pd.Series(    trend + weekly + monthly + noise,    index=pd.date_range("2025-01-01", periods=days),    name="amount",)sales.head()

# FFT スペクトルの上位 3 周期y = sales.values - sales.mean()freqs = np.fft.rfftfreq(len(y), d=1.0)spectrum = np.abs(np.fft.rfft(y))periods = np.where(freqs > 0, 1 / freqs, np.inf)mask = (periods >= 2) & (periods <= len(y) / 2)top = sorted(np.argsort(spectrum * mask)[::-1][:3], key=lambda i: -spectrum[i])for i in top:    print(f"発見: 周期 {periods[i]:.1f} 日, 強さ {spectrum[i]:.0f}")# 出力例:#   発見: 周期 7.0 日, 強さ 36500#   発見: 周期 28.0 日, 強さ 27300#   発見: 周期 365.0 日, 強さ 18000  ← 長期トレンドが残ったもの
# MSTL で複数周期を同時分解from statsmodels.tsa.seasonal import MSTLmstl = MSTL(sales, periods=[7, 28]).fit()mstl.plot()

import pandas as pdfrom statsmodels.tsa.stattools import acf
# GA4 BigQuery export を SQL で日次 PV にしておく想定# columns: date, page_path, pageviewsdf = pd.read_csv("ga4_daily_pageviews.csv", parse_dates=["date"])
def cyclicality_score(series, period=7):    """指定周期での自己相関の強さをスコア化（-1〜1）"""    if len(series) < period * 2:        return None    s = series - series.mean()    a = acf(s, nlags=period, fft=True)    return float(a[period])
# ページごとに 7日周期のスコアresult = (    df.groupby("page_path")      .apply(lambda g: cyclicality_score(g.set_index("date")["pageviews"], period=7))      .dropna()      .sort_values(ascending=False)      .head(10))print("曜日周期が強いページ TOP10:")print(result)