2022年夏の気温予測～東京都心部編～

2022年3月28日 19:49

1. 自己紹介

　はじめまして。第二子の育休中にプログラミングの勉強を始めたai-itoと
申します。

　文系学部卒。新卒で入社した銀行の総合職11年目。これまでプログラミングは一切やったことがありませんでしたが、近年バズワードとなっている「DX（デジタルトランスフォーメーション）」の大波は当然銀行業界にも
押し寄せています。
　近い将来、関連スキルを習得する必要に迫られることは明白なので、ある程度まとまった時間を取れる育休中に少しでも教養を身に付けるべく、勉強を始めました。

　とはいえ、「転職したい！」「自分でこういうモノを作ってみたい！」といった強い動機があるわけではないため、何から始めるべきか散々迷いました。そして迷った末、数多あまたあるプログラミング言語の中でも初心者が学びやすいと評判のPythonを選び、Python学習に特化したAidemyさんで学ぶことに決めました。

2. 目的

　乳幼児を抱えていると、これからの季節、気になるのが熱中症。
　そこで、東京都心部における過去約10年間の日別最高気温を基に、2022年夏の気温を予測し、子ども達との外遊び計画の参考にしたいと思います。

3. 使用するデータ

　東京都心部における2011年1月1日～2022年3月18日の日別最高気温
　※ 気象庁のホームページからダウンロードしました。以下のキャプチャの
　　とおり、取得したい地点・項目・期間・表示オプションを選びます。

　何かに使えるかなと思い、表示オプションで「日付に曜日を表示」を選択しましたが、今回の予測では曜日情報は使いませんでした。

4. 実行環境

　Google Colaboratory
　Python 3.7.13
　NumPy 1.21.5
　Pandas 1.3.5
　Matplotlib 3.2.2

5. 流れ

　① データの読み込み
　② データの整理
　③ パラメータの決定
　④ モデルの構築
　⑤ 予測
　⑥ グラフの可視化

① データの読み込み

　気象庁HPからダウンロードしたCSVファイルは以下のとおりです。

まずは、このファイルをpandasを用いて読み込みます。

import pandas as pd
df = pd.read_csv("/content/drive/MyDrive/Colab Notebooks/blog/highest_temp_tokyo_20110101_20220318.csv", header=2, encoding="shift-jis")
print(df.head(10))

※ ' header=2 '：上から1・2行目は不要なので、3行目からデータが読み込ま
　　　　　　　れるよう ' header ' を指定します。
※ ' encoding="shift-jis" '：当該CSVファイル内に日本語文字が含まれるため
　　　　　　　　　　　　こちらを指定しなければエラーとなります。

② データの整理

　さらに不要な行列を削除します。

df=df.drop([0,1], axis=0)
df=df.drop(columns=['曜日', '最高気温(℃).1', '最高気温(℃).2'], axis=1)
print(df.head(10))

カラム名を年月日→'Date'、最高気温(℃)→'Highest_temp'に変更します。
また、今回は時系列データの分析をするので、時間情報（日付データ＝'Date'）をdatetimeに変換してindexとすることで、データを扱いやすく
します。

df = df.rename(columns={'年月日' : 'Date', '最高気温(℃)' : 'Highest_temp'})
df['Date'] = pd.to_datetime(df['Date'])
df = df.set_index('Date')
print(df.head())

データの形が整いました。
ここで、必要なモジュールをインポートしておきます。

import warnings
import itertools
import numpy as np
import statsmodels.api as sm
import datetime
import matplotlib.pyplot as plt
warnings.simplefilter('ignore')

③ パラメータの決定

　今回の最高気温予測は時系列データの解析なので、季節周期のある時系列
データを扱うことができるSARIMAモデルを用います。
　なお、PythonにはSARIMAモデルのパラメータ（p, d, q）（sp, sd, sq, s）を
自動的に最適化してくれる機能はありません。
　そのため、情報量規準（今回の場合はBIC（ベイズ情報量基準））によって最適な値を求めるプログラムを書く必要があります。
　以下のコードで時系列データ：DATA, パラメータs（周期）：s を入力すると、最適なパラメータとそのBICを出力します。

def selectparameter(DATA, s):
   p = d = q = range(0, 2)
   pdq = list(itertools.product(p, d, q))
   seasonal_pdq = [(x[0], x[1], x[2], s)
                   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(DATA,
                                               order=param,
                                               seasonal_order=param_seasonal)
               results = mod.fit()
               parameters.append([param, param_seasonal, results.bic])
               BICs = np.append(BICs, results.bic)
           except:
               continue
   return parameters[np.argmin(BICs)]

selectparameter(df,12)

扱うデータが気温なので、月毎のデータに季節性があると考え、パラメータs（周期）は ' 12 ' とします。
以下の通り、最適な（p, d, q）,（sp, sd, sq, s）が求められました。
※周期を春夏秋冬と考え ' s=4 ' でも試してみたところ、結果は同じでした。

[(1, 1, 1), (0, 0, 0, 12), 20565.644832040496]

④ モデルの構築

　モデルの構築には sm.tsa.statespace.SARIMAX(DATA,order=(p, d, q),
seasonal_order=(sp, sd, sq, s)).fit() を用います。

SARIMA_model = sm.tsa.statespace.SARIMAX(df.Highest_temp, order = (1, 1, 1), seasonal_order=(0, 0, 0, 12)).fit()

⑤ 予測

　predに予測期間の予測値を代入します。

pred = SARIMA_model.predict('2020-01-01', '2022-09-30')

⑥ グラフの可視化

　データを折れ線グラフで示し、予測値は赤色でプロットします。

#グラフのサイズを指定
plt.figure(figsize=(16, 9))

#グラフのタイトルを設定
plt.title("Highest_temperature_Tokyo2022")

#グラフのx軸とy軸の名前設定
plt.xlabel("date")
plt.ylabel("Highest_temperature_Tokyo")

#グラフにグリッドを表示
plt.grid(True)

#データのプロット
plt.plot(df, color="b", label="actual")
plt.plot(pred, color="r", label="predict")

#系列ラベルを表示
plt.legend()
plt.show()

約10年分の日別データなので、データ量が多すぎて分かりづらいです。
そこで、表示範囲を2022年1月1日～2022年9月30日に絞ってみます。

# グラフのサイズを指定
plt.figure(figsize=(16, 9))

# グラフのタイトルを設定
plt.title("Highest_temperature_Tokyo2022")

# グラフのx軸とy軸の名前設定
plt.xlabel("date")
plt.ylabel("Highest_temperature_Tokyo")

# x軸の表示範囲を[2022-01-01, 2022-09-30]に指定
sxmin='2022-01-01'
sxmax='2022-09-30'
xmin = datetime.datetime.strptime(sxmin, '%Y-%m-%d')
xmax = datetime.datetime.strptime(sxmax, '%Y-%m-%d')
plt.xlim([xmin,xmax])

# グラフにグリッドを表示
plt.grid(True)

# データのプロット
plt.plot(df, color="b", label="actual")
plt.plot(pred, color="r", label="predict")

# 系列ラベルを表示
plt.legend()
plt.show()

　
　かなり見やすくはなりましたが、、、
　青線（実際の気温）と赤線（モデルによる予測気温）が重なっている部分については、ある程度正しく予測できたものの、未来日の部分がうまく予測できていません。
　Aidemyのチューターさんにアドバイスをいただき、上記③パラメータの
決定の部分で ' p = d = q = range(0, 2) ' の範囲を ' range(0, 3) '、
' range(0, 4) '、' range(0, 5)' と広げてみましたが、膨大な実行時間がかかった挙句完了に至らず断念。データ量を減らすために日別データを月平均データに変更した上で試みても、同様の結果でした。

　パラメータの決定方法は他にも色々あるため、未来日部分の予測は今後の課題とし、引き続き試行錯誤したいと思います。
　また、過去の気温以外にも、湿度、気圧、日照時間、海水温等といった
説明変数を加えると、より正確な予測ができるのではと期待しています。

6. まとめ

　育休中とはいえ、家事育児の傍らプログラミングという未知の分野を学習するには予想以上に時間が足りず、自分がどの程度成長できたのかは疑問
です。しかし、不明点があればチューターの皆さんがいつも迅速・丁寧に
教えてくださり、なんとか完走することができました。

　今回の記事では触れていませんが、講座の中で学んだ " パラメータチューニングの重要性 " や、" 一見役立ちそうにないデータでも、切り口次第では予測性能の向上に貢献しうる " といった視点は、私にとって大きな気づきでした。これは、Pythonの使用未使用に関わらず、仕事であらゆるデータ分析をする際に即実践しようと思います。

　この3ヵ月でPythonの基礎の基礎は身に付けられたつもりでいますが、Python関連の様々なWEBページや書籍を読んでいると「？？？」となることはいまだに多く、まだまだ学ぶべきことは山ほどあります。今回踏み出した「はじめの一歩」を無駄にしないためにも、今後も学習を継続していきたいと思います。

2022年夏の気温予測 ～東京都心部編～