必要なモジュールの読み込み¶
In [ ]:
import numpy as np
import pandas as pd
import numpy as np
from statsmodels.tsa.seasonal import STL
from statsmodels.tsa.stattools import adfuller
import matplotlib.pyplot as plt
plt.style.use('ggplot') #グラフスタイル
plt.rcParams['figure.figsize'] = [10, 5] #グラフサイズ
plt.rcParams['font.size'] = 10 #フォントサイズ
import japanize_matplotlib #グラフ内で日本語利用
時系列データのグラフ化¶
In [ ]:
# 月次データの読み込み
df_m = pd.read_csv(
'MMM_ts_manga5_monthly.csv',
index_col='month',
parse_dates=True
)
df_m.index = df_m.index.to_period('M')
print(df_m)
sales traditional internet promotional month 2016-01 4030283000 426400000 666900000 650475000 2016-02 5334600000 465450000 657450000 434325000 2016-03 3356769000 1119000000 0 404475000 2016-04 5004399000 926500000 613725000 197400000 2016-05 11650794000 1523800000 699675000 157162500 ... ... ... ... ... 2023-08 5942020000 1521500000 584475000 193537500 2023-09 3889616000 931600000 469875000 326812500 2023-10 3118571000 0 551700000 182812500 2023-11 3107643000 0 410325000 338625000 2023-12 3942322000 656100000 326250000 437775000 [96 rows x 4 columns]
In [ ]:
# 年と月の変数の作成
df_m['Year'] = df_m.index.year
df_m['Month'] = df_m.index.month
# 年別月間売上の集計
df_yearly_sales = df_m.groupby('Year')['sales'].sum()
# 年別月間売上と年単位の売上の推移のグラフを上下に並べて表示
fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(8, 8))
# 年単位の売上の推移のグラフ
bar = axes[0].bar(
df_yearly_sales.index,
df_yearly_sales.values,
color='skyblue'
)
for rect in bar:
height = rect.get_height()
axes[0].annotate(f'{height/1e8:.0f}億円',
xy=(rect.get_x() + rect.get_width() / 2, height),
xytext=(0, 3), # 3 points vertical offset
textcoords='offset points',
ha='center', va='bottom')
axes[0].set_title('年単位の売上の推移')
axes[0].set_xlabel('年')
axes[0].set_ylabel('売上')
axes[0].set_ylim(0, 1e11)
axes[0].set_xticks(df_yearly_sales.index)
axes[0].tick_params(axis='x', rotation=45)
# 年別月間売上のグラフ
years = df_m['Year'].unique()
colors = plt.cm.viridis(np.linspace(0, 1, len(years)))
for i, year in enumerate(years):
df_year = df_m[df_m['Year'] == year]
axes[1].plot(
df_year['Month'],
df_year['sales'],
label=str(year),
color=colors[i]
)
axes[1].set_title('年別月間売上')
axes[1].set_xlabel('月')
axes[1].set_ylabel('売上')
axes[1].legend()
axes[1].set_xticks(range(1, 13))
axes[1].set_ylim(0, 1.6e10)
axes[1].set_xticklabels([
'1月', '2月', '3月', '4月', '5月', '6月',
'7月', '8月', '9月', '10月', '11月', '12月'
])
plt.tight_layout()
plt.show()
In [ ]:
# 2023年のデータをデータフレームから抽出
df_2023 = df_m.loc['2023']
# 4つのサブプロットを作成
total_rows = 4
fig, axes = plt.subplots(nrows=total_rows, ncols=1, figsize=(8, 5))
# 各列のプロット情報
graph_info = [
('2023年 売上（sales）', 'sales', '売上金額', 1.6e10),
('2023年 マスコミ四媒体広告費（ほぼテレビCM）', 'traditional', 'コスト', 1.6e9),
('2023年 インターネット広告費', 'internet', 'コスト', 1.6e9),
('2023年 プロモーションメディア広告費', 'promotional', 'コスト', 1.6e9)
]
# プロットをループで作成
for ax, (title, column, ylabel, ylim) in zip(axes, graph_info):
ax.set_title(title)
ax.plot(df_2023['Month'], df_2023[column])
ax.set_ylabel(ylabel)
ax.set_xticks(range(1, 13))
ax.set_ylim(0, ylim)
ax.set_xticklabels(['1月', '2月', '3月', '4月', '5月', '6月', '7月', '8月', '9月', '10月', '11月', '12月'])
# レイアウトを調整して表示
plt.tight_layout()
plt.show()
時系列データの基礎分析¶
In [ ]:
# 週次データの読み込み
df = pd.read_csv(
'MMM_ts_manga5.csv',
index_col='week',
parse_dates=True
)
print(df)
sales traditional internet promotional week 2015-12-27 649315000 656100000 0 0 2016-01-03 958266000 0 178050000 188587500 2016-01-10 755973000 0 150450000 162037500 2016-01-17 551603000 0 0 126900000 2016-01-24 520183000 0 175500000 0 ... ... ... ... ... 2023-11-26 839617000 0 154950000 0 2023-12-03 829314000 0 0 249187500 2023-12-10 852512000 0 148200000 0 2023-12-17 963276000 656100000 0 0 2023-12-24 1297220000 0 178050000 188587500 [418 rows x 4 columns]
In [ ]:
# 直近5年間（52.25週×5年）のデータに絞り
# 目的変数yと説明変数Xに分解
# 対象データ期間の設定
data_term = int(52.25 * 5)
# 説明変数Xと目的変数yに分解
X = df.drop(columns=['sales'])[-data_term:]
y = df['sales'][-data_term:]
In [ ]:
from statsmodels.tsa.seasonal import STL
import matplotlib.pyplot as plt
# STL分解の実行
# y: 元の時系列データ (売上データ)
# period: 頻度（ここでは52週=1年間を表す）
stl = STL(y, period=52)
result = stl.fit()
# 分解された成分の取得
# trend: トレンド成分 (全体的な長期的傾向)
# seasonal: 季節成分 (周期的な変動)
# resid: 残差成分 (トレンドと季節成分を除いた部分)
trend = result.trend
seasonal = result.seasonal
resid = result.resid
# 結果の可視化
plt.figure(figsize=(8, 5))
# 元の時系列データをプロット
plt.subplot(411)
plt.plot(y.index, y)
plt.title('元の時系列データ - 売上(sales)')
plt.ylim(bottom=0, top=5e9)
# トレンド成分をプロット
plt.subplot(412)
plt.plot(y.index, trend)
plt.title('トレンド成分')
plt.ylim(bottom=-2e9, top=3e9)
# 季節成分をプロット
plt.subplot(413)
plt.plot(y.index, seasonal)
plt.title('季節成分')
plt.ylim(bottom=-2e9, top=3e9)
# 残差成分をプロット
plt.subplot(414)
plt.plot(y.index, resid)
plt.title('残差成分')
plt.ylim(bottom=-2e9, top=3e9)
# レイアウト調整して表示
plt.tight_layout()
plt.show()
In [ ]:
# STL分解後の残差(resid)に対して拡張ディッキー・フラー検定(Augmented Dickey-Fuller test)を適用
# ADF検定は単位根の存在を確認するためのもので、時系列データが定常かどうかを判定します
result = adfuller(resid)
# ADF検定のテスト統計量を出力
print('テスト統計量:', result[0])
# 検定の有意性を判断するためのp値を出力
print('p値:', result[1])
# 信頼水準（1%、5%、10%）ごとの臨界値を出力
print('臨界値:', result[4])
テスト統計量: -11.58222731294389
p値: 2.924170854469547e-21
臨界値: {'1%': -3.4557539868570775, '5%': -2.8727214497041422, '10%': -2.572728476331361}
In [ ]:
# 定常性検定の実行
# 各説明変数について拡張ディッキー・フラー検定（ADFテスト）を実施し、結果を出力します
# 説明変数ごとに検定を実行
for column in X.columns:
# 各列のデータに対してADF検定を実施
result = adfuller(X[column])
# 結果をコンソールに出力
print(column) # 列名（変数名）を出力
print('- Test Statistic:', result[0]) # テスト統計量を出力
print('- p-value:', result[1]) # p値を出力
print('- Critical Values:', result[4]) # 信頼水準における臨界値を出力
traditional
- Test Statistic: -10.159352311769545
- p-value: 7.579888911588827e-18
- Critical Values: {'1%': -3.4558530692911504, '5%': -2.872764881778665, '10%': -2.572751643088207}
internet
- Test Statistic: -15.681906957312055
- p-value: 1.4691906894548407e-28
- Critical Values: {'1%': -3.4557539868570775, '5%': -2.8727214497041422, '10%': -2.572728476331361}
promotional
- Test Statistic: -5.29778728915373
- p-value: 5.533738398557567e-06
- Critical Values: {'1%': -3.457105309726321, '5%': -2.873313676101283, '10%': -2.5730443824681606}
航空機の乗客数（Passengers）データの例¶
In [ ]:
# 航空機の乗客数（Passengers）データの読み込み
df = pd.read_csv(
'AirPassengers.csv',
index_col='Month',
parse_dates=True
)
In [ ]:
# STL分解の実行
stl = STL(df['Passengers'], period=12)
result = stl.fit()
# 分解された成分の取得
trend = result.trend
seasonal = result.seasonal
resid = result.resid
# トレンド成分を除去したデータの作成
passengers_detrended = df['Passengers'] - trend
# トレンドおよび季節成分を除去したデータの作成
passengers_detrended_deseasonalized = df['Passengers'] - trend - seasonal
# 結果の可視化
plt.figure(figsize=(8, 5))
# 共通のy軸範囲を設定
y_min = min(
df['Passengers'].min(),
passengers_detrended.min(),
passengers_detrended_deseasonalized.min()
)
y_max = max(
df['Passengers'].max(),
passengers_detrended.max(),
passengers_detrended_deseasonalized.max()
)
plt.subplot(3,1,1)
plt.plot(df.index, df['Passengers'])
plt.title('元の時系列データ - 乗客数（Passengers）')
plt.ylim(y_min, y_max) # スケールを統一
plt.subplot(3,1,2)
plt.plot(df.index, passengers_detrended)
plt.title('トレンド成分を調整（除去）した乗客数（Passengers）')
plt.ylim(y_min, y_max) # スケールを統一
plt.subplot(3,1,3)
plt.plot(df.index, passengers_detrended_deseasonalized)
plt.title('トレンドおよび季節成分を調整（除去）した乗客数（Passengers）')
plt.ylim(y_min, y_max) # スケールを統一
plt.tight_layout()
plt.show()
