Python初心者 コロナの影響は？ 箱根のパティシエが需要予測してみた
はじめに
はじめまして。Python初心者…
どころかパソコンもほとんど触ったことのないパティシエです…！
この度はAidemy講座のデータ分析コースの成果物としてこちらに掲載することとなりました。
「なぜパティシエがデータ分析を…？」と思いますよね。笑
コロナ渦に箱根にOPENしたお菓子屋で働いているのですが、
現在２０２３年１０月はコロナが５類感染症になり、売り上げも安定していてむしろ上昇傾向にあります。
嬉しい反面、もしやコロナの反動で一時的に需要が増えただけなのでは…？
と考えるようになったのがきっかけです。
そこで、過去のデータからコロナが無かった場合で売り上げを予測し、現在の状況と比べることで、コロナが無かった状態に近づいているのかそうでないかを確認していこうと考えました。
OPENしたばかりお店のため過去のデータがないので（あったとしても公開できませんが…）今回は国土交通省観光庁が公開している「国内旅行の平均単価」の統計をサンプルにして時系列解析をすすめていきます。
時系列解析とは
時間の経過とともに変化するデータ のことをいいます。具体的には東京の毎時間の気温の記録、ある会社の毎月の売上高の記録、ある国のGDPの記録や日経平均株価の記録などは時系列データといえます。
時系列分析は、会社の売り上げや商品の売り上げ予測、さらには来店者数の予測など、ビジネスにおいて非常に重要な分析技術となります。
時系列…普段あまり使わない言葉で正直ピンとこなかった私が説明すると、
・９月１日 ケーキ 売上げ１００個
・９月２日 ケーキ 売上げ１２０個
・９月３日 ケーキ 売上げ５０個
︙
以上のように日付（時間の経過）とともに売上げ（データ）が変化していくことをさします。
まずは現状をグラフで可視化していき、その後パターンや傾向を加味した今後の予想を可視化したものが時系列解析です。
以下の手法を用いて分析していきます。
・ARモデル(自己回帰モデル)
ARモデルは過去の値から回帰的にある時点の値を推定するモデルです。
以前の値に影響されるモデルであり、直前p個の値と相関をもちます。
例)今日のケーキの売上げを先週７日間（ｐ個＝ｐ日間）の売り上げを用いて推定する（単回帰）。
・MAモデル(移動平均モデル)
MAモデルは単に時系列が異なる式との間に共通部分を持つために相関性が存在するという性質を持つモデルです。
以前の誤差に影響され、直前q個の値の影響を受けます。
例)今年は８月限定でケーキのクーポンを有効とした。クーポン有効のため、以前と比べて売り上げがあがる誤差ある。また、その１か月間は自己相関ができている。
・ARMAモデル
ARモデルとMAモデルは互いに競合する性質がありません。そのためこの2つを組み合わせて定式化することが可能です。AR(p)のARモデルとMA(q)のMAモデルを組み合わせたときARMA(p,q)のように表すことができます。
ARMAモデルは定常過程にしか適応できません。
・ARIMAモデル
先程のARMAモデルへ原系列を階差系列に変換し適応したものをARIMAモデルと言います。
d時点前との差分をとった場合のARMA(p,q)で構築したARIMAモデルをARIMA(p,d,q)と表します。
☆pを自己相関度（モデルが直前 p 個の値を用いて予測されるのかを表す）
☆dを誘導（時系列データを定常にするために d 次の階差が必要だったことを表す）
☆qを 移動平均（モデルが直前 q 個の値に影響を受けることを表す）
・SARIMAモデル
ARIMAモデルをさらに季節周期を持つ時系列データにも拡張できるようにしたモデルです。SARIMAモデルは(p, d, q)のパラメーターに加えて(sp, sd, sq, s)というパラメーターも持ちます。
sp,sd,sq はそれぞれ 季節性自己相関、 季節性導出、 季節性移動平均といいます。
例えば12ヶ月周期の季節変動を持つデータの場合、 s のパラメーターは周期を表すため、s=12 となります。
sq=1 ならちょうど12ヶ月(1周期前)のデータ, sq=2 なら12ヶ月前のデータと24ヶ月前のデータの影響をモデルが受けることを表します。
分かりにくければ単純に q を sq に置き換えてみてください。
馴染みのない英語がたくさんでたきましたが…
超簡単にいうと、
((AR+MA)=ARMA)<ARMA<ARIMA<SARIMA
SARIMAモデルが様々な分析要素を持ちあわせているため、必要事項を入力するだけで、季節性周期を加味しながら過去のデータを反映した予測ができる便利なモデルということです。
では、SARIMAモデルを使用して時系列解析を進めましょう♪
環境、サンプルデータ
環境
Windows 11 Pro
Chrome
Python3
google colaboratory
サンプルデータ
サンプルデータとして、国土交通省観光庁が公開している「旅行・観光消費動向調査」の日本国内居住者であり、住民基本台帳をもとに無作為に抽出した約2万6000人を対象とした統計（※平成30年1-3月期調査までは約2万5000人統計）で、国内旅行平均単価をサンプルデータとして使用しております。
URL：旅行・観光消費動向調査 | 統計情報 | 統計情報・白書 | 観光庁 (mlit.go.jp)
SARIMAモデルの実装
以下の手順で分析を進めていきます。
1.データの読み込み
2.データの整理
3.データの可視化
4.データの周期の把握 (パラメーターsの決定)
5.s以外のパラメーターの決定
6.モデルの構築
7.データとの予測とその可視化
1.データの読み込み
今回は、google driveの「MyDrive」にExcelで「旅行平均単価 1」という名前で保存しています。
また、実装に必要なライブラリをインポートしていきます。
from google.colab import drive
drive.mount('/content/drive')
from google.colab import drive
drive.mount('/content/drive')
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.discrete.discrete_model import L1MultinomialResultsWrapper
import itertools #日本語表示をサポートするためのライブラリ
!pip install japanize-matplotlib
import japanize_matplotlib #Matplotlibライブラリ内のサブモジュール
import matplotlib.dates as mdates
df=pd.read_excel("/content/drive/MyDrive/旅行平均単価 1.xlsx")
df=pd.read_excel("/content/drive/MyDrive/旅行平均単価 1.xlsx")
2.データの整理
ｄｆの中身を確認しましょう。
２０１３年１月から２０２３年６月までの旅行平均単価（１人あたり）の消費額が記載されています。
日付 平均単価
0 2013-01 47365.353017
1 2013-02 47680.069905
2 2013-03 45832.642869
3 2013-04 47358.971187
4 2013-05 46304.197595
.. ... ...
121 2023-02 56112.000000
122 2023-03 62131.000000
123 2023-04 60574.000000
124 2023-05 62208.000000
125 2023-06 58862.000000dfのタイプをdf.dtypesで確認していきます。
# データ型の確認
df.dtypes
この後の作業で、「年だけを取得したい、月だけを取得したい、日付だけを取得したい」などをしたいと思った時、dateカラムのデータ型がobject型だとうまく認識されないため、datetime型に変換しておきます。
変数"df"の 日付カラム のデータ型を変換したいので、loc日付の列を抽出し、applyでpandasのto_datetime()関数を使い、データ型の変換を行います。
また、不要なインデックスをdf.dropで削除し、日付をdateの名称に変えます。
df.loc[:, 'date'] = df.loc[:, '日付'].apply(pd.to_datetime)
df = df.set_index("date")
df = df.drop('日付', axis=1)
print(df)きちんと変換されているか、再度確認しておきます。
現時点のdfはコロナ前から現在までの公開されているすべてのデータですが、コロナが蔓延しなかった世界線もプロットしたいので、コロナ前だけのデータを抽出しておきます。
#コロナ前のデータ
before_COVID_data = df['2013-01':'2019-11']
print(before_COVID_data)3.データの可視化
statsmodels（Pythonの統計分析とモデリングのためのライブラリ）のsm.tsa.seasonal_decompose で、時系列データの季節性成分、トレンド成分、および残差成分を分解していきます。
res = sm.tsa.seasonal_decompose(before_COVID_data)
trend = res.trend # トレンドデータ
seasonal = res.seasonal # 季節性データ
resid=res.resid
plt.figure(figsize=(12, 8)) # グラフ描画枠作成、サイズ指定
# 原型列のプロット
plt.subplot(511)
plt.plot(df)
plt.ylabel('Original')
# トレンドデータのプロット
plt.subplot(512)
plt.plot(trend)
plt.ylabel('Trend')
# 季節性データのプロット
plt.subplot(513)
plt.plot(seasonal)
plt.ylabel('Seasonality')
# 残差データ のプロット
plt.subplot(514)
plt.plot(resid)
plt.ylabel('resid')plt.figure(figsize=(12, 8))
以上のコードでグラフ描画枠作成、サイズ指定をしましたので、plt.subplot()でグラフの位置を指定します。
#trendデータplt .subplot(512)のプロットを例にすると、
「５」サブプロットの行数を５個に設定
「１」サブプロットの列数を１個に指定
「２」サブプロットの位置を２番目に指定
#原型列データ （Original）で元のデータを確認します。
#トレンドデータ （Trend）で時間の経過とともに変化するデータの傾向やパターンを示すデータを確認します。
→上グラフは全体の傾向でいうと右肩あがりなので平均消費単価が上昇していると推測できます。
#季節性データ （Seasonality）を分析することで、特定の季節や周期に関連する傾向やパターンを確認します。
→上グラフは１年（１２か月）周期のパターンの傾向があると推測できます。
#残差データ （resid)は原型列データからトレンドや季節性のデータを除いて残ったデータ
→２０１８年６月頃に降下８月頃上昇しているが、推測は今回できませんでした。
4.データの周期の把握 (パラメーターsの決定)
SARIMAモデルには、7つのパラメーター(p, d, q) (sp, sd, sq, s)があります。
p：自己相関度といい、モデルが直前 p 個の値を用いて予測されるのかを表しています。
d：誘導といい、 時系列データを定常にするために d 次の階差が必要だったことを表しています。
q：移動平均といい、モデルが直前 q 個の値に影響を受けることを表しています。
sp,sd,sq：それぞれ季節性自己相関、季節性導出、季節性移動平均といい、上記3つのパラメーターに季節性の要素を持たせたものを表しています。
s：季節周期を表します。例えば12ヶ月周期のデータであれば s = 12 とすればよいです。
ここまで、データの整理・可視化を行い、1年単位の周期がありそうなことが分かりました。月毎のデータにしたので、季節周期を表すパラメーター s = 12 としていきます。
5.s以外のパラメーターの決定
PythonにはSARIMAモデルのパラメーター(p, d, q) (sp, sd, sq, s)を自動で最も適切にしてくれる機能はありません。そのため情報量規準 (今回は BIC)によってどの値が最も適切なのか調べるコードを書かなければなりません。BICの場合は、値が低ければ低いほどパラメーターが適切であると言えます。
このコードの書き方はまだ理解できていないので、今回はaidemy教材の中から引用していきます。selectparameter()関数を定義し、 p,d,q を0から1の範囲の中から最適な値を探していきます。
返却されたベストパラメータには、季節データから周期が12ヵ月と推定されたため、ｓ＝12に設定する。
def selectparameter(DATA,s):
# p,d,qの範囲を0～1に設定
# p = d = q = range(0, 1)
p = d = q = range(0, 2)
pdq = list(itertools.product(p, d, q))
seasonal_pdq = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]
parameters = []
BICs = np.array([])
for param in pdq:
for param_seasonal in seasonal_pdq:
try:
mod = sm.tsa.statespace.SARIMAX(resid,
order=param,
seasonal_order=param_seasonal)
results = mod.fit()
parameters.append([param, param_seasonal, results.bic])
BICs = np.append(BICs,results.bic)
except:
continue
# bestパラメータの返却
print(parameters[np.argmin(BICs)])
return parameters[np.argmin(BICs)]
# 返却されたベストパラメータを受け取るよう変更
best_parameter = selectparameter(before_COVID_data,s=12)6.モデルの構築
ベストパラメーターを使用して、SARIMAモデルに当てはめていきます。
# モデルの当てはめ
SARIMA_before_COVID_data = sm.tsa.statespace.SARIMAX(
before_COVID_data,
# order=(0, 0, 0),
order=best_parameter[0],
# seasonal_order=(0, 0, 0, 12)
seasonal_order=best_parameter[1]
).fit()予測したい時期入力し、.predict 予測データをpredに代入します。
# predに予測データを代入する
pred = SARIMA_before_COVID_data.predict("2019-10", "2027-11")
7.データとの予測とその可視化
predを可視化します。
fig,ax = plt.subplots(1,1, figsize=(30,10))でpltのサイズを指定し、ax.setでタイトルやラベル等の詳細設定をおこないます。
# predデータともとの時系列データの可視化
plt.subplot(515)
plt.plot(df)
plt.plot(before_COVID_data)
plt.plot(pred, color="r")
plt.ylabel('comparison')
plt.show()
fig,ax = plt.subplots(1,1, figsize=(30,10))
ax.set_title("旅行にかける平均消費量",size=20)
ax.set_xlabel('日付',size=15) # x軸のラベルを設定
ax.set_ylabel("平\n均\n単\n価", labelpad=15, size=15,rotation=0,va='center') # y軸のラベルを縦に表示
ax.plot(df)
ax.plot(before_COVID_data)
ax.plot(pred, color="r")
months = mdates.MonthLocator(bymonth=[5,8,11])#５月、８月、１１月の繁忙期にメモリ設定
years = mdates.YearLocator()
majorFmt = mdates.DateFormatter('%Y_%m') # 大きい軸用フォーマット
minorFmt = mdates.DateFormatter('_%m') # 小さい軸用フォーマット
ax.grid(which='major')
ax.grid(which='minor', linestyle=':')
ax.xaxis.set_major_formatter(majorFmt) # majorFmtを適用
ax.xaxis.set_minor_formatter(minorFmt) # minorFmtを適用
ax.xaxis.set_major_locator(years)
ax.xaxis.set_minor_locator(months)
ax.tick_params(axis='x', which='minor', rotation=45, labelsize='small')
plt.xlim(df.index.min(), pred.index.max())
plt.xticks(rotation=45);
plt.show()
考察
グラフから、コロナ渦の下降している部分は緊急事態宣言（1回目～3回目）がでている時期でした。２０２２年１０月に観光業界活性化のための全国旅行支援がスタートし、平均単価が急上昇しているのがわかります。
その後も旅行支援が延長したことをきっかけに予測データを上回る消費量の月がありましたが、２０２３年からは消費量が予測データに近づいてきたので、コロナの観光業界への影響が収束にむかっているのではないかと考えました。
今回は国内の旅行にかける平均消費量から今後の予測データをだしましたが、訪日外国人が年々増えていますので（JNTO（日本政府観光局）による統計データより）、そちら併せても予測していく必要があると考えました。
おわりに
今後実用化に向けてのデータ分析の課題としては、現地（箱根）のイベント、天気、気温を加味した精度の高い予測データをだすことで、製造計画に役立てたいと思います。
はじめてのPython…難しい！と思うことが沢山ありました。笑
けれど、難しいことにチャレンジしていくことがとても楽しかったです。
去年まで全く興味がなかった分野ですが、真剣にお店の売り上げを考えたときに、やみくもに走るよりデータをもとに今後の予測をし、根拠のもとに予測を実現に変える…叶えたい未来に近づけるのではないかと考え、一歩を踏み出しました。
まだまだこれからも難しい場面に直面と思いますが、今後も難しい課題にこそやりがいを感じていけたらと思います。
最後までご覧いただきありがとうございました。
この記事が気に入ったらサポートをしてみませんか？