Retail sales forecastのデータ分析

## はじめに
Aidemyでデータ分析コースを受講しました。その成果をまとめていきたいと思います。
## この記事の概要
この記事では、Kaggleで提供されているRetailsalesForecastを使用してデータ分析を行なっていきます。
## 開発環境
Mac OS Sonoma 14.3 Apple M1
Python 3.9.13
Goolge colab
## 動機
前職はアパレル業界で働いていた為、小売業界のデータ分析をしたいと思いました。そこで、Kaggleで提供されていたRetailsaleForecastを利用して、成果物を作ろうと思います。

## 目的
小売データのモデリングの1つの課題は、限られた履歴に基づいて意思決定を行う必要があるということです。休日や特定の主要イベントは年に一度しか発生せず、戦略的な決定が業績にどのような影響を与えたかを見るチャンスも一年に一度しかありません。また、値引きは売上に影響を与えることが知られています - 課題は、どの部門が影響を受け、どの程度影響を受けるかを予測することです。
## 内容
異なる地域にある45店舗の過去の販売データが提供されています。各店舗には複数の部門が含まれています。会社はまた、1年を通じていくつかのプロモーション価格改定イベントを実施しています。これらの価格改定は、スーパーボウル、労働者の日、感謝祭、クリスマスの4つの主要な休日の前に行われます。これらの休日を含む週は、非休日の週よりも評価が5倍高くなります。

Excelシートには、Stores、Features、Salesの3つのタブがあります。
## 方向性
以下の年の各店舗の各部門の売上を予測する
休日週における価格改定の効果をモデル化する
得られた洞察に基づいて推奨アクションを提供する。最大のビジネス影響を優先して考慮する
## 分析方法
小売りデータ分析を使用して、販売量を予測するために販売のみと販売＋外部情報を使用した単純な自己回帰（AR）モデルをテストします。

## 目次
1.データの収集
2.データの整理
3.モデルの構築
4.データの予測・可視化
5.考察
## 1.データの収集
まずは、ライブラリをインポートします。
```python: 
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
```
Date型に変換
```python: 
features['Date'] = pd.to_datetime(features['Date'])
sales['Date'] = pd.to_datetime(sales['Date'])
```

### Kaggleからデータ収集
[RetailsaleForecast](https://www.kaggle.com/datasets/manjeetsingh/retaildataset)
### 提供されているCSVデータを読み込み
```python:
features=pd.read_csv('../input/Features data set.csv')
sales=pd.read_csv('../input/sales data-set.csv')
stores=pd.read_csv('../input/stores data-set.csv')
```
## データの確認
### 特徴量
Store 店舗：店舗番号
Date 日付：週
Temperature 気温：地域の平均気温
Fuel_Price 燃料価格：地域の燃料費
MarkDown1-5：プロモーション価格改定に関連する匿名化されたデータ。MarkDownデータは2011年11月以降のみ利用可能であり、すべての店舗が常に利用できるわけではありません。欠損値はNAでマークされています。
消費者物価指数（CPI）：消費者物価指数
Unemployment 失業率：失業率
IsHoliday：その週が特別な休日週であるかどうか
### データの詳細を確認
```python:
print(features.shape)
print(sales.shape)
print(stores.shape)

print(sales[0:1].Date, sales[-1:].Date)

print(features[0:1].Date, features[-1:].Date)
```
```
(8190, 12)
(421570, 5)
(45, 3)
0   2010-05-02
Name: Date, dtype: datetime64[ns] 421569   2012-10-26
Name: Date, dtype: datetime64[ns]
0   2010-05-02
Name: Date, dtype: datetime64[ns] 8189   2013-07-26
Name: Date, dtype: datetime64[ns]
```

'sales' DataFrameと 'features' DataFrameを左結合して新しいDataFrameを作成します。これにより、販売データと特徴量データが結合されます。作成したDataFrame（df）と 'stores' DataFrameを左結合して、新しいDataFrameを作成します。これにより、店舗情報が追加されます。新しいDataFrameの欠損値を0で埋めます。変更後のDataFrameの先頭行を表示します。

```python:
df=pd.merge(sales,features, on=['Store','Date', 'IsHoliday'], how='left')
df=pd.merge(df,stores, on=['Store'], how='left')

df=df.fillna(0)
df['Temperature'] = (df['Temperature']- 32) * 5./9.

types_encoded, types =df['Type'].factorize()
df['Type'] = types_encoded

df.head()
```

```
	Store	Dept	Date	Weekly_Sales	IsHoliday	Temperature	Fuel_Price	MarkDown1	MarkDown2	MarkDown3	MarkDown4	MarkDown5	CPI　Unemployment	Type	Size
0	1	1	2010-05-02	24924.50	False	5.727778	2.572	0.0	0.0	0.0	0.0	0.0	211.096358	8.106	0	151315
1	1	1	2010-12-02	46039.49	True	3.616667	2.548	0.0	0.0	0.0	0.0	0.0	211.242170	8.106	0	151315
2	1	1	2010-02-19	41595.55	False	4.405556	2.514	0.0	0.0	0.0	0.0	0.0	211.289143	8.106	0	151315
3	1	1	2010-02-26	19403.54	False	8.127778	2.561	0.0	0.0	0.0	0.0	0.0	211.319643	8.106	0	151315
4	1	1	2010-05-03	21827.90	False	8.055556	2.625	0.0	0.0	0.0	0.0	0.0	211.350143	8.106	0	151315
```

## グラフで視覚化

データフレーム内に存在する変数のうちいくつかをプロットします。例えば、温度、燃料価格、消費者物価指数、失業率などです。

```python:
df[['Date', 'Temperature', 'Fuel_Price', 'CPI', 'Unemployment', 
    'MarkDown1', 'MarkDown2', 'MarkDown3', 'MarkDown4', 'MarkDown5']].plot(x='Date', subplots=True, figsize=(20,15))
plt.show()
```
![スクリーンショット 2024-02-26 18.11.07.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3633845/60441286-61d8-2ea8-ffdb-9674baceecf6.png)

小売チェーンの週次売上高を得るために、すべての店舗からの週次売上高を合計する。'Date' 列と週次売上合計の2つの列を持つ新しいDataFrameであり、df_average_sales_weekに格納されます。週次売上合計に基づいてDataFrameを降順でソートします。つまり、週次売上が最も高い週から順に並び替えます。

```python:
df_average_sales_week = df.groupby(by=['Date'], as_index=False)['Weekly_Sales'].sum()
df_average_sales = df_average_sales_week.sort_values('Weekly_Sales', ascending=False)

plt.figure(figsize=(20,5))
plt.plot(df_average_sales_week.Date, df_average_sales_week.Weekly_Sales)
plt.show()
```

## 総売上高の予測
```python:
from statsmodels.graphics.tsaplots import acf, pacf, plot_acf, plot_pacf

fig, axes = plt.subplots(1,2, figsize=(20,5))
plot_acf(ts, lags=100, ax=axes[0])
plot_pacf(ts, lags=100, ax=axes[1])
plt.show()
```
![スクリーンショット 2024-02-26 18.27.07.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3633845/fd3f3920-318c-fcf3-b975-8a52a287ba67.png)

## モデルの定義
statsmodelsライブラリで実装されているARモデルは、連続した季節性項目のみを含めることができます。そのため、相関分析から観察されるように、週1、6、52などの連続でない季節性項目を使用できるように、モデルを手動で実装します。

```python:
from sklearn.linear_model import LinearRegression

def fit_ar_model(ts, orders):
    
    X=np.array([ ts.values[(i-orders)].squeeze() if i >= np.max(orders) else np.array(len(orders) * [np.nan]) for i in range(len(ts))])
    
    mask = ~np.isnan(X[:,:1]).squeeze()
    
    Y= ts.values
    
    lin_reg=LinearRegression()
    
    lin_reg.fit(X[mask],Y[mask])
    
    print(lin_reg.coef_, lin_reg.intercept_)

    print('Score factor: %.2f' % lin_reg.score(X[mask],Y[mask]))
    
    return lin_reg.coef_, lin_reg.intercept_
    
def predict_ar_model(ts, orders, coef, intercept):
    return np.array([np.sum(np.dot(coef, ts.values[(i-orders)].squeeze())) + intercept  if i >= np.max(orders) else np.nan for i in range(len(ts))])
```

自己相関関数（ACF）と偏自己相関関数（PACF）から観察される季節性に基づいて、週（1、6、52）からの季節性を含むARモデルが実装されます。

```python:
orders=np.array([1,6,52])
coef, intercept = fit_ar_model(ts,orders)
pred=pd.DataFrame(index=ts.index, data=predict_ar_model(ts, orders, coef, intercept))
plt.figure(figsize=(20,5))
plt.plot(ts, 'o')
plt.plot(pred)
plt.show()
```
![スクリーンショット 2024-02-26 18.57.34.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3633845/96ea175a-24b1-773e-5b25-7e4a917a3ddf.png)

```python:
diff=(ts['Weekly_Sales']-pred[0])/ts['Weekly_Sales']

print('AR Residuals: avg %.2f, std %.2f' % (diff.mean(), diff.std()))
 
plt.figure(figsize=(20,5))
plt.plot(diff, c='orange')
plt.grid()
plt.show()
```
![スクリーンショット 2024-02-26 19.00.52.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3633845/77207c97-370c-4fd1-16ac-cd604f2ae16c.png)

予測モデルの決定係数（R^2）は0.41であり（完全な予測の場合の最大スコアは1です）、残差の分布は0を中心にしており、標準偏差は7％です。

## 店舗別売上高の予測

最も高い売上高を示す店舗番号20のための予測モデルを開発します。
'Store' 列が20の行だけを残し、それ以外の行をNaNに置き換えます。
NaN値を持つ行を削除します。
Date' 列でデータをグループ化し、各週の売上の合計を計算します。
'Date' 列をDataFrameのインデックスに設定します。

```python:
df20=df.where( df['Store'] == 20)
df20=df20.dropna()
df20=df20.groupby(by=['Date'], as_index=False)['Weekly_Sales'].sum()
df20 = df20.set_index('Date')
df20.head()
```
![スクリーンショット 2024-02-26 19.05.43.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3633845/4e8b7cfc-f3d9-0e33-7fd0-c71a8e398342.png)


```python:
plt.figure(figsize=(20,5))
plt.plot(df20.index, df20.values)
plt.show()
```
![スクリーンショット 2024-02-26 19.07.17.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3633845/38e5d4bb-41d0-e0b3-7e85-cd86a1f0107a.png)

```python:
fig, axes = plt.subplots(1,2, figsize=(20,5))
plot_acf(df20.values, lags=100, alpha=0.05, ax=axes[0])
plot_pacf(df20.values, lags=100, alpha=0.05, ax=axes[1])
plt.show()
```

![スクリーンショット 2024-02-26 19.12.24.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3633845/2d3f7706-18fd-d0bd-a427-05c5bd16edd1.png)

この特定の店舗では、週29と週46から追加の季節性が観察されます。

```python:
orders=np.array([1,6,29,46,52])
coef, intercept = fit_ar_model(df20,orders)
pred=pd.DataFrame(index=df20.index, data=predict_ar_model(df20, orders, coef, intercept))
plt.figure(figsize=(20,5))
plt.plot(df20, 'o')
plt.plot(pred)
plt.show()
```
![スクリーンショット 2024-02-26 19.17.19.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3633845/2ea4920d-9e17-57e9-ecd8-2f42e310bdc9.png)

```python:
diff=(df20['Weekly_Sales']-pred[0])/df20['Weekly_Sales']

print('AR Residuals: avg %.2f, std %.2f' % (diff.mean(), diff.std()))
 
plt.figure(figsize=(20,5))
plt.plot(diff, c='orange')
plt.grid()
plt.show()
```

![スクリーンショット 2024-02-26 19.19.30.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3633845/5055c0a6-83b6-621e-a3f1-acd62fe0765f.png)
予測モデルの決定係数（R^2）は0.34であり（完全な予測の場合の最大スコアは1です）、残差の分布は0を中心にしており、標準偏差は8％です。

## 外部変数から予測

```python:
import seaborn as sns
corr = dfext.corr()
plt.figure(figsize=(10,10))
sns.heatmap(corr, 
            annot=True, fmt=".3f",
            xticklabels=corr.columns.values,
            yticklabels=corr.columns.values)
plt.show()
corr['shifted_sales'].sort_values(ascending=False)
```
![__results___46_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3633845/c276ac0e-bd13-55d2-f0cb-20f3befa2d6f.png)

![スクリーンショット 2024-02-26 19.29.09.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3633845/396e2f16-85c8-9221-710d-9cab39275c2f.png)


利用可能な外部変数は、1日の遅延売上時系列といくらか相関しています。これは、1日間に予測力があることを意味し、モデルを改善するために使用できます。 'MarkDown' と 'Temperature' が最も相関しており、それぞれ相関があります。


```python:
def fit_ar_model_ext(ts, orders, ext, fitter=LinearRegression()):
    
    X=np.array([ ts.values[(i-orders)].squeeze() if i >= np.max(orders) else np.array(len(orders) * [np.nan]) for i in range(len(ts))])
    
    X = np.append(X, ext.values, axis=1)
    
    mask = ~np.isnan(X[:,:1]).squeeze()
    
    Y= ts.values
    
    fitter.fit(X[mask],Y[mask].ravel())
    
    print(fitter.coef_, fitter.intercept_)

    print('Score factor: %.2f' % fitter.score(X[mask],Y[mask]))
    
    return fitter.coef_, fitter.intercept_
    
def predict_ar_model_ext(ts, orders, ext, coef, intercept):

    X=np.array([ ts.values[(i-orders)].squeeze() if i >= np.max(orders) else np.array(len(orders) * [np.nan]) for i in range(len(ts))])
    
    X = np.append(X, ext.values, axis=1)
    
    return np.array( np.dot(X, coef.T) + intercept)

dfexte=dfext[['Unemployment','Fuel_Price','CPI','Temperature',
              'MarkDown1', 'MarkDown2', 'MarkDown3', 'MarkDown4', 'MarkDown5']]

orders=np.array([1,6,29,46,52])
coef, intercept = fit_ar_model_ext(df20,orders,dfexte)
pred_ext=pd.DataFrame(index=df20.index, data=predict_ar_model_ext(df20, orders, dfexte, coef, intercept))
plt.figure(figsize=(20,5))
plt.plot(df20, 'o')
plt.plot(pred)
plt.plot(pred_ext)
plt.show()
```
![スクリーンショット 2024-02-26 19.37.15.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3633845/9f9185c9-dbb7-ed00-f689-6a236e218b71.png)

```python:
diff=(df20['Weekly_Sales']-pred[0])/df20['Weekly_Sales']
diff_ext=(df20['Weekly_Sales']-pred_ext[0])/df20['Weekly_Sales']

print('AR Residuals: avg %.2f, std %.2f' % (diff.mean(), diff.std()))
print('AR wiht Ext Residuals: avg %.2f, std %.2f' % (diff_ext.mean(), diff_ext.std()))
 
plt.figure(figsize=(20,5))
plt.plot(diff, c='orange', label='w/o external variables')
plt.plot(diff_ext, c='green', label='w/ external variables')
plt.legend()
plt.grid()
plt.show()
```
![スクリーンショット 2024-02-26 19.38.42.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3633845/8cbbbf1d-2202-0744-6cd7-edb22f3acf90.png)

## 結論

外部変数を含むモデルは、予測の精度を40％以上向上させました（R^2スコア：0.58、対して0.34）。残差の標準偏差も約30％改善しました（7％、対して8％）

## 最後に
データ分析を勉強する事でデータの取捨選択、統合、相関関係の有無やグラフの視覚化で、どの要素が結果に影響するのかを分析する知識・技術を身につける事ができました。
今回、勉強した内容を基に日々の業務にデータを活用して、効率的に結果を導き出していきたいです。