ブロッコリーの値段分析にProphetを使ってみた

## はじめに
ブロッコリーの値段の変動を分析・予測するために、[Facebook Prophet](https://facebook.github.io/prophet/)を利用しました。データは月次データであり、内部的に日付形式へ変換しています（例: "2023年4月" → "2023-04-01"）。

データは、[政府統計の総合窓口（e-Stat）](https://www.e-stat.go.jp/)から取得しました。具体的には、[小売物価統計調査](https://www.e-stat.go.jp/dbview?sid=0003421913)のデータを加工して使用しています。

出典：政府統計の総合窓口(e-Stat)（https://www.e-stat.go.jp/ ）

今回使用したコードはこちらです。

https://github.com/makinzm/hajimeteno-prophet?tab=readme-ov-file


## Prophetのモデル

### 全体モデル
時系列データは次のようなモデルで表されます。

```math
y(t) = g(t) + s(t) + h(t) + \epsilon_t
```

- $y(t)$: 時刻 $t$ における観測値
- $g(t)$: トレンド
- $s(t)$: 季節性
- $h(t)$: 休日やイベントの影響
- $\epsilon_t$: 独立したノイズ（平均0の正規分布）


Prophetは加法モデル（Additive Model）を基本としますが、データの性質によっては乗法モデル（Multiplicative Model）も使用可能です。乗法モデルでは、季節性やイベントの影響がトレンドに比例します。

- Additive（加法モデル）:
  ```math
  y(t) = g(t) + s(t) + h(t) + \epsilon_t
  ```

- Multiplicative（乗法モデル）:
  ```math
  \log(y(t)) = \log(g(t)) + s(t) + h(t) + \epsilon_t
  ```

:::note info
[Forecasting at scale [PeerJ Preprints]](https://peerj.com/preprints/3190/)のpp.7-8より。
:::

### トレンドモデル
Prophetは、以下の2種類のトレンドモデルをサポートしています。

1. **非線形な飽和成長モデル**（logistic growth）
2. **区分線形モデル**（piecewise linear）

デフォルトで区分線形モデルが使用されます。このモデルは次のように表されます。

```math
g(t) = (k + a(t)^T \delta)t + (m + a(t)^T \gamma)
```

- $k$: 基本的な成長率
- $\delta$: 各変化点での成長率の調整ベクトル
- $m$: オフセットパラメータ
- $\gamma_j$: 変化点の影響を補正するためのパラメータ

```math
\gamma_j = -s_j \delta_j
```

ここで $s_j$ は $j$ 番目の変化点の位置を表します。

#### a(t) の定義
$a(t)$ は各時刻 $t$ における変化点の影響を表すベクトルで、次のように定義されます：

```math
a_j(t) =
\begin{cases} 
1 & \text{if } t \geq s_j \\
0 & \text{otherwise}
\end{cases}
```


:::note info
[Forecasting at scale [PeerJ Preprints]](https://peerj.com/preprints/3190/)のpp.9-10より。
:::

デフォルトで区分線形モデルが使用される理由は、GitHubの[Prophetコード](https://github.com/facebook/prophet)における`growth='linear'`が標準設定であることから判断しています。

### changepointの自動選択と事前分布
changepoint（成長率が変化する点）は、次のスパースなラプラス分布を使用して自動で選択されます。

```math
\delta_j \sim Laplace(0, \tau)
```

- $\tau$: changepointの柔軟性を調整するパラメータ

:::note info
[Forecasting at scale [PeerJ Preprints]](https://peerj.com/preprints/3190/)のpp.10より。
:::

### 季節性モデル
季節性はフーリエ展開を用いて次のように表されます。

```math
s(t) = \sum_{j=1}^{n} \left[ \beta_{j,\text{cos}} \cos\left( \frac{2\pi jt}{P} \right) + \beta_{j,\text{sin}} \sin\left( \frac{2\pi jt}{P} \right) \right]
```

- $P$: 季節性の周期（例: 年次データの場合、$P = 365.25$ 日）
- $\beta_{j,\text{cos}}, \beta_{j,\text{sin}}$: フーリエ係数

これらのフーリエ係数は次の正規分布に従います。

```math
\beta_j \sim \mathcal{N}(0, \sigma^2)
```

:::note info
[Forecasting at scale [PeerJ Preprints]](https://peerj.com/preprints/3190/)のpp.11より。
:::

### 休日やイベントの影響（h(t)）
$h(t)$は特定の日付における値への影響をモデル化します。次の式で表されます。

```math
h(t) = Z(t)\kappa
```

- $Z(t)$: イベントまたは休日の時点を示すワンホットベクトル
- $\kappa$: イベントや休日による影響の強度

イベントリストを指定することで、特定の休日やイベントの影響を考慮できます。影響は正規分布に従います。

```math
\kappa \sim \mathcal{N}(0, \nu^2)
```

:::note info
[Forecasting at scale [PeerJ Preprints]](https://peerj.com/preprints/3190/)のpp.12より。
:::

## 使用したコードとデフォルト値

以下のコードでモデルを構築しました。

```python
model = Prophet(
    seasonality_mode="multiplicative",
    yearly_seasonality=True,
    weekly_seasonality=False,
    daily_seasonality=False,
    mcmc_samples=300
)
model.fit(df_city)
```

- **`seasonality_mode="multiplicative"`**: 季節性が値段に比例すると仮定
- **`yearly_seasonality=True`**: 年次の季節性を有効化
- **`weekly_seasonality=False`**: 週次の季節性を無効化
- **`daily_seasonality=False`**: 日次の季節性を無効化
- **`mcmc_samples=300`**: MCMCサンプリングを300回実施

## 結果と考察

実際の分析結果のコードはこちらです。

https://github.com/makinzm/hajimeteno-prophet?tab=readme-ov-file

地域ごとに予測結果を可視化しました。以下は、札幌市のデータを例にしたプロットです。

![札幌市](https://github.com/makinzm/hajimeteno-prophet/blob/main/notebooks/output/%E6%9C%AD%E5%B9%8C%E5%B8%82_row_plot.png?raw=true)

横一列に並んだ3つのグラフで構成しています。それぞれのグラフの内容と解釈のポイントを以下に示します。

### 1. 左列: **予測結果**
このプロットは、過去の観測値と未来の予測値を可視化しています。

#### 内容
- **青い点**: 過去の観測データ（実測値）
- **濃い青の線**: Prophetモデルがフィッティングしたトレンド（予測値）
- **薄い青の帯（信頼区間）**: 予測の不確実性を示します（[80](https://github.com/facebook/prophet/blob/dc1df4cb23a150e14858afb34c9442401c0eb2fc/python/prophet/forecaster.py#L97)%の信頼区間）

#### 解釈のポイント
1. 観測値にモデルがどれだけうまくフィットしているかを確認します。
2. 信頼区間が広がっている箇所は、予測が不確実であることを示しています。
3. トレンドが上昇・下降している部分が長期的な価格の変動傾向を表しています。

### 2. 中央列: **コンポーネントプロット**
このプロットは、予測結果を構成する要素（コンポーネント）を分解して表示しています。

#### 内容
- **トレンド（Trend）**: 時系列全体での長期的な変化を表します。
- **年次の季節性（Yearly Seasonality）**: 1年間を通じた周期的な価格変動を表します（例えば、夏に価格が上昇するなど）。
    - こちらの縦軸は％で表記されています。
    - こちらは、乗法モデルを使用しているため、季節性の影響がどれだけ大きいかを示しています。（[該当コード](https://github.com/facebook/prophet/blob/dc1df4cb23a150e14858afb34c9442401c0eb2fc/python/prophet/plot.py#L435C1-L440C14)）

#### 解釈のポイント
1. トレンドコンポーネントを確認して、全体的な価格の上昇や下降を把握します。
2. 季節性コンポーネントでは、特定の月や時期に価格が上昇または下降する特徴が明確になります。
   - 例えば、冬季に価格が安定する場合、年次季節性コンポーネントにその傾向が現れます。

### 3. 右列: **各年ごとの月別データプロット**
このプロットは、各年の月ごとのデータを個別に可視化したものです。Prophetとは直接関係ありませんが、データの特徴を把握するために使用します。

#### 内容
- 各年のデータが異なる色の線でプロットされ、**年次ごとの変化**を比較できるようになっています。
- x軸: 月（1月から12月）
- y軸: ブロッコリーの価格

#### 解釈のポイント
1. 年次ごとの月別価格の差を比較します。
   - 例えば、ある年に特定の月の価格が大きく異なる場合、その年に特異な要因（異常気象、需給バランスの変化など）が存在した可能性があります。
2. 季節ごとのパターンが年をまたいで一貫しているかを確認します。
   - 例えば、毎年夏に価格が上昇する傾向が一貫していれば、強い季節性の特徴を示していると考えられます。

---

正直なところ、地域差がバラバラであるため、一般的な結論を得るのは難しいと感じました。可視化をする前には、ブロッコリーの旬が冬から春にかけてであるため、その時期の値段が低くなると予想していました。しかし、実際には地域によって異なるため、一概には言えませんでした。

何個かをピックアップして説明をしたいと思います。

### 1. 札幌市

![札幌市](https://github.com/makinzm/hajimeteno-prophet/blob/main/notebooks/output/%E6%9C%AD%E5%B9%8C%E5%B8%82_row_plot.png?raw=true)

最も右の図を見ると7月は値段が安いように見えます。実際真ん中の図を見ると、少し見づらいですが、7月（左から数えて膨らんでいない箇所の7個目のy軸に注目）に値段が下がっていることがわかります。

妄想でしかありませんが、夏野菜が旬を迎えるため、ブロッコリーの需要が減少し価格が下がっている可能性があります。

また、全体的なトレンドとしては、昔から値段の変動はあまりないように見えます。

### 2. 津市

![津市](https://raw.githubusercontent.com/makinzm/hajimeteno-prophet/refs/heads/main/notebooks/output/%E6%B4%A5%E5%B8%82_row_plot.png)

真ん中の図から9月に値段が上昇していることがわかります。また、右の図からも9月に値段が上昇していることがわかります。

これは札幌市とは異なります。地域によって特性があることが伺えます。

また、真ん中の全体のトレンドを見ると2016年ごろに価格のピークを迎えていて、その前後は値段が下がっているようにみえます。

この傾向は他の地域で見られるため、全国的になにかしらの要因がある可能性があります。

(他の地域について説明はしませんがGitHubに画像はあるため、気になる方は見てみてください。)

## まとめ
Prophetを使うと、トレンドや季節性を考慮した時系列データの予測が簡単に行えました。

## 感想

Prophetは非常に便利なツールだと感じました。

一方で時系列データの分析は地域差が多くデータを見るだけでは、因果関係を特定するのは難しいと感じました。そのため、データ分析を行う際は、関連の有りそうな他のデータも併せて分析することが重要だと感じました。