あまりに暑いので，140年分の気温をProphetで分析した

# はじめに

統計開始以来最も暑い夏を生きる皆様，お疲れ様です．あまりに暑くてムシャクシャしたので，[気象庁](http://www.data.jma.go.jp/gmd/risk/obsdl/index.php)から約140年分（1872年1月1日-2018年7月21日）の東京の気温データを入手し，[Prophet](https://facebook.github.io/prophet/)で分析しました．

以下は，過去365日分の最高気温の実測値と，2018年7月22日から1ヶ月先までの予測値を表したものです．

<img width="710" alt="2018-07-22-forecast_max.png" src="https://qiita-image-store.s3.amazonaws.com/0/151195/36949132-719d-2ed9-5d03-2fe5858bd89e.png">

また，以下は，最高気温のトレンドと年単位の周期性を表したものです．

<img width="638" alt="2018-07-22-components_max.png" src="https://qiita-image-store.s3.amazonaws.com/0/151195/efde7a04-8f4c-4ea0-58ae-1a0e3b886936.png">

分析の結果，平均・最高・最低気温の全てに関して，1920年付近から上昇し続けており，そのトレンドを考慮してもなお，ここ数日は特に暑いことを確認しました．分析に用いたNotebookは[こちら](https://github.com/haltaro/summer-prophet)です．

> **注**：本記事は，2018年7月22日に[個人サイトに投稿した記事](https://haltaro.github.io/2018/07/22/summer-prophet)を，Qiita向けに再構成したものです．

# 環境

- macOS Sierra, 10.12.6
- Python, 3.6.6
- Prophet, 0.3

# Prophet

## 特徴

[Prophet](https://facebook.github.io/prophet/)は，Facebookが開発したオープンソースの時系列解析ライブラリです．Prophetで想定するモデルは：

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

ここで，$y(t)$は時刻$t$における予測値，$g(t)$は時点$t$におけるトレンド成分，$s(t)$は時点$t$における周期成分，$h(t)$は時点$t$における週末などのイレギュラーな成分，そして$\epsilon(t)$は時点$t$における誤差を表します．各成分のパラメータを[Stan](http://mc-stan.org/)で推定し，モデルを学習します．詳細なモデルは，原論文（[Sean J. Taylor and Benjamin Letham, Forecasting at Scale](https://peerj.com/preprints/3190.pdf)）を参照してください．

Prophetに関する日本語の記事としては，以下のようなものがあります．とても参考になりました．

- [Prophet: forecasting at scale - Facebook](https://facebook.github.io/prophet/)
- [Prophet入門【Python編】Facebookの時系列予測ツール - SlideShare](https://www.slideshare.net/hoxo_m/prophet-facebook-76285278)
- [Facebookの予測ライブラリProphetを用いたトレンド抽出と変化点検知 - Gunosyデータ分析ブログ](https://data.gunosy.io/entry/change-point-detection-prophet)

Prophetの大きな特徴の一つは，時系列解析の前提知識がなくても，誰でも簡単に一定水準の分析をできることです．

## インストール

今回はPythonからProphetを使います．[pip](https://pypi.org/project/pip/)を使う場合は，`pip install fbprophet`で，[conda](https://conda.readthedocs.io/en/latest/)を使う場合は，`conda install -c conda-forge fbprophet`でインストールできます．

# 準備

## データの入手

[過去の気象データ・ダウンロード - 気象庁](https://www.data.jma.go.jp/gmd/risk/obsdl/index.php#)からcsvデータをダウンロードします

- 地点：東京
- 項目：
  + 日別値
  + 気温（平均気温，最高気温，最低気温）
- 期間：1872年1月1日 - 2018年7月21日（53528日）

まず，以下の画面で都道府県および地区を選択します．

<img width="1060" alt="2018-07-22-area.png" src="https://qiita-image-store.s3.amazonaws.com/0/151195/95cd3b7e-538f-b8eb-50e0-8d97cff0ad29.png">

次に，以下の画面で項目を選択します．降水量や湿度も興味深いですが，一度にダウンロードできるデータ量に制限があるため，今回は気温（平均気温，最高気温，最低気温）のみを選択します．

<img width="1037" alt="2018-07-22-item.png" src="https://qiita-image-store.s3.amazonaws.com/0/151195/b81b2450-4ab3-ecb2-eb2e-6539de7f2a8d.png">


最後に，以下の画面で期間を選択します．一度にダウンロードできるデータ量に制限があるため，10年ごとに刻んでダウンロードします．

<img width="1044" alt="2018-07-22-period.png" src="https://qiita-image-store.s3.amazonaws.com/0/151195/4800b36d-0e4c-27f8-80f0-abd898e5e7de.png">


Downloadフォルダの`data.csv`を年代ごとに改名し，`data/raw/`以下に次のように保存します．

- `data_1870.csv`
- `data_1880.csv`
- `data_1890.csv`
- `data_1900.csv`
- `data_1910.csv`
- `data_1920.csv`
- `data_1930.csv`
- `data_1940.csv`
- `data_1950.csv`
- `data_1960.csv`
- `data_1970.csv`
- `data_1980.csv`
- `data_1990.csv`
- `data_2000.csv`
- `data_2010.csv`

## モジュールのインポート

```py
import numpy as np # 数値計算
import pandas as pd # DataFrame
import matplotlib.pyplot as plt # グラフ描画
import seaborn as sns # グラフ描画設定
from tqdm import tqdm # forループの進捗状況確認
import codecs # Shift-JIS読み込み用
from fbprophet import Prophet # prophet

sns.set_style(style='ticks') # グラフスタイルの指定
```

## データの前処理

```py
years = range(1870, 2020, 10)
read_path = 'data/raw/'

all_data = pd.DataFrame()

for year in tqdm(years):
    data_file = 'data_{}.csv'.format(year)
    
    # 補足1：普通にpd.read_csvすると読み込めない
    with codecs.open(read_path+data_file, "r", 
                     "Shift-JIS", "ignore") as f:
        # 補足2：先頭数行のデータは使用しない．
        # また，不要なカラムは除外して利用する．
        data = pd.read_table(
            f, delimiter=",", skiprows=5, index_col=0, 
            usecols=[0, 1, 4, 7])
        
    # datetime形式に変換
    data.index = pd.to_datetime(data.index)
    # カラム名を変更
    data.columns = ['average', 'high', 'low']
    
    # all_dataを更新
    all_data = pd.concat([all_data, data])
    
# 補足3：すべてNanの行は削除
all_data = all_data.dropna(how='all')
```

以下で，詳細を補足します．

### `codecs.open`

`Shift-JIS`で作成されているため，普通に`pandas.read_csv()`を使っても，ファイルを読み込めません．[pandasでread_csv時にUnicodeDecodeErrorが起きた時の対処 (pd.read_table()) - Qiita](https://qiita.com/niwaringo/items/d2a30e04e08da8eaa643)を参考に，`codecs.open()`を利用します．

> **2018年7月25日追記**：@skotaroさんのご指摘の通り，`pandas.read_csv(encoding='Shift-JIS')`で普通に読み込めました．私の不勉強でした．ご指摘ありがとうございました！

### 不要な行・列を除外

元データをそのまま読み込むと，以下のように不要な行・列が含まれてしまいます．

<img width="605" alt="2018-07-22-data.png" src="https://qiita-image-store.s3.amazonaws.com/0/151195/fa3ee962-3026-3654-d43f-072d3bc5731a.png">

そこで，[`pandas.read_table()`](https://pandas.pydata.org/pandas-docs/version/0.23/generated/pandas.read_table.html)のパラメータ`skiprows`と`usecols`を用いて，分析に必要な部分のみ抽出します．

### すべて`np.nan`の行を削除

気象庁のサイト上では1872年1月1日から指定可能でしたが，最初の数年分はデータが含まれていませんでした．`pandas.DataFrame.dropna(how='all')`で，すべての値が`np.nan`の行を削除します．これにより全体の約2%にあたる1272日分のデータが除外されます．

- 地点：東京
- 項目：
  + 日別値
  + 気温（平均気温，最高気温，最低気温）
- 期間：1875年6月10日 - 2018年7月21日（52266日）

# 分析

平均気温，最高気温，および最低気温の分析を行します．

## 平均気温

まず，すべての点を散布図でプロットしてみます．

```py
plt.figure(figsize=(8, 6))
plt.plot(all_data.index, all_data.average, 'o', alpha=.1)
plt.xlabel('Date')
plt.ylabel('Average temparature')
plt.show()
```

<img width="498" alt="2018-07-22-scatter_ave.png" src="https://qiita-image-store.s3.amazonaws.com/0/151195/c82f56a2-498a-0f1a-2b57-a4abccc9e0ef.png">


140年間で平均気温がじわじわ上がっていく様子がわかります．Prophetでより詳細に分析します．

```py
# Prophet用DataFrameを生成
ave_data = pd.DataFrame()
ave_data['ds'] = all_data.index
ave_data['y'] = all_data.average.values

# Prophetモデルの構築
m = Prophet(weekly_seasonality=False)

# 学習
m.fit(ave_data)
```

Prophetで時系列解析を行うためには，以下のカラムを持つ`pandas.DataFrame`を生成する必要があります．

- `ds`：日付カラム．
- `y`：予測対象カラム．

また，Prophetではデフォルトで`yearly_seasonality`と`weekly_seasonality`を仮定して学習します．気温は曜日に依存しないと考えられるため，`weekly_seasonality`を除外します．

以上で学習したモデル`m`を使い，将来30日間の平均気温を予測します．以下のように，`m.make_future_dataframe`で`DataFrame`を生成すると楽です．

```py
# 将来30日間を予測
future = m.make_future_dataframe(periods=30)
forecast = m.predict(future)
```

予測結果を`m.plot()`でプロットします．約140年分プロットするとわけがわからないので，直近1年分のみを表示するようにします．

```python
m.plot(forecast)
plt.xlim(future.ds.iloc[-365-30], future.ds.iloc[-1])
plt.xlabel('Date')
plt.ylabel('Average')
```

<img width="710" alt="2018-07-22-forecast_ave.png" src="https://qiita-image-store.s3.amazonaws.com/0/151195/9854f717-113f-6483-fc1c-e2e5b7e354bb.png">


黒点が実測値，青線の実線が予測値，そして薄い青の塗りつぶしが80%信頼区間[^interval]を示します．ここ数日の平均気温は，信頼区間の上限ギリギリに張り付いていることが確認できます．例年から統計的に予測したとき，いかにここ数日の気温が高いかわかります[^outlier]．また，例年の傾向を考慮すると，気温のピークは8月初旬であるという絶望的な事実がわかります．

[^interval]: デフォルト設定．

[^outlier]: 信頼区間を外れた異常値は他にも結構存在しますが，連日，しかも同じ方向に外れることは珍しいように思います．

次に，`m.plot_components(forecast)`で周期性とトレンドを確認します．

```py
m.plot_components(forecast)
```

<img width="638" alt="2018-07-22-components_ave.png" src="https://qiita-image-store.s3.amazonaws.com/0/151195/bf41f376-6fda-3cc7-d4cb-2033f9525bbb.png">

上の図がトレンド，下の図が年内の周期性を表します．東京の平均気温は1920年あたりから上がり続けていることがわかります．また，年内で見ると，気温のピークは7月下旬から8月上旬にあることがわかります．

## 最高気温

同様の分析を行います．詳細は割愛しますが，平均気温と同様の性質が確認できます．

<img width="710" alt="2018-07-22-forecast_max.png" src="https://qiita-image-store.s3.amazonaws.com/0/151195/d271d934-7501-9730-c46d-a6137c7af4be.png">

<img width="638" alt="2018-07-22-components_max.png" src="https://qiita-image-store.s3.amazonaws.com/0/151195/9649d91b-1a18-d93e-f5a9-ea2f79fa3e7e.png">

## 最低気温

同様の分析を行います．詳細は割愛しますが，平均気温と同様の性質が確認できます．

<img width="710" alt="2018-07-22-forecast_min.png" src="https://qiita-image-store.s3.amazonaws.com/0/151195/7ae07418-4e10-9b47-0132-94b9ffa50180.png">

<img width="638" alt="2018-07-22-components_min.png" src="https://qiita-image-store.s3.amazonaws.com/0/151195/1e7d7425-6348-8c1a-fbdd-18efdf08e68f.png">

# おわりに

地学の知識が中学理科で止まっている私でも，今年の夏はヤバイということがわかりました[^donbiki]．無理は禁物です．ご自愛くださいませ．

あと，この記事を書き終わるころ，東洋経済さんの素晴らしい記事を見つけました（[東京の夏が｢昔より断然暑い｣決定的な裏づけ: 過去140年の最高気温をビジュアル化 - 東洋経済オンライン](https://toyokeizai.net/articles/-/229965)）[^trend]．こんなにわかりやすいヒートマップを初めて見ました．いつかこんなビジュアライゼーションをしてみたいです．

最後まで読んでくださり，ありがとうございました！

[^donbiki]: あまりに明確に温暖化のトレンドが出過ぎてちょっと引きました．

[^trend]: 一応，1920年ごろから気温が上がり始めるという分析結果は，本記事と整合しています．良かった．