東京都の降水量データを試行錯誤する①（データの取得編）

今回はSARIMAモデルを試してみたく、季節性のありそうな降水量データを気象庁から取得してデータ分析したいと思います。

先にモデルを作成した結論として、1年間の降水量の上がり下がりをそれなりに追従する形でできました。
一方で、一番降水量の多い台風は異常値として扱われてしまったようで、他の月より多い降水量を予測しているが実態とは乖離が大きかったです。

今回はデータの取得部分と指標部分をノートに記載します。
SARIMAモデルのパラメータを作成したりする部分は後日記載します。
（東京都の降水量データを試行錯誤する②(モデル作成編)）

この記事で実装するプログラムの概要

・気象庁のホームページから東京都降水量データを取得する
・分析対象の期間は直近10年(2010/01～2019/12)
・降水量データの内最初の7年(2010/01～2016/12)を学習用データとし
　残りの3年分をテスト用として分ける
・降水量データを使って「自己相関係数」「編自己相関係数」「ADF検定」「季節調整済系列」を表示する

データの取得

まずは必要なライブラリを読み込みます。
Pandas、NumPy、MatplotLib等お決まりのものを中心に今回はSARIMAXもimportします。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
import itertools
from statsmodels.tsa.statespace.sarimax import SARIMAX

そして、Pandasを使ってデータを読み込み整形する。

##################################################################################################################
# データ取得＆整形

# 気象庁から東京都のデータ取得
df = pd.read_html('https://www.data.jma.go.jp/obd/stats/etrn/view/monthly_s3.php?prec_no=44&block_no=47662&year=1964&month=7&day=&view=p5')

# 降水量に限定
df = df[0]
# 年をindexにする
df = df.set_index('年')
# 初年度と今年は欠損があるため除外
df = df[1:-1]
df = df[-11:-1]
# いらいない列を削除して1月～12月に限定する
df = df.drop('年の値',axis=1)

# 欠損確認
# isnull()で1列でもTrue(欠損あり)の行を抽出する
#df.iat[0,0] = np.NaN
print('欠損数: {}'.format(df[df.isnull().any(axis=1)].size))

# 記号が混ざっており下のfloat変換に失敗するため数値以外を破棄
df = df.replace('[^0-9.]','', regex=True)

# 降水量が文字列として読み込まれているのでdfをグラフ表示できない(floatに変換する)
df = df.astype(float)

# 1次元に変換する
nda = df.to_numpy().reshape(df.size)
index = pd.date_range('{}-01-31'.format(df.index[0]), '{}-12-31'.format(df.index[-1]), freq = "M")
df = pd.DataFrame(data=nda,index=index,columns=['precipitation'])

# 学習用と予測用にデータを分ける
train_num = int(len(df)/12*0.7)*12
train = df[:train_num]
print('降水量データ数：{}か月({}年)　うち学習用 {}か月({}年)'.format(len(df),int(len(df)/12),train_num,train_num/12))

降水量データは下図の通り縦に年、横に月のデータになっています。
※下図は整形前の気象庁で公開しているデータそのままです

指標を見る

４つの指標を表示します。
コードの中にコメントで指標の解釈を記載しています。

自己相関係数

# 自己相関係数
# ⇒ 1から95%信頼区間内にプロットされており、自己相関はほぼない
sm.graphics.tsa.plot_acf(df['precipitation'],lags=60)

編自己相関係数

# 編自己相関係数
# ⇒ 11,12と相関がある
sm.graphics.tsa.plot_pacf(df['precipitation'],lags=60)

ADF検定

# ADF検定
# ⇒ p-value <= 0.05なので定常データと言える
print('ADF検定 p-value: {}'.format(sm.tsa.stattools.adfuller(df['precipitation'], autolag='AIC')[1]))

季節調整済系列の表示

# 季節調整済系列の表示
# ⇒ Trendグラフを見るとトレンドは一定ではない
sm.tsa.seasonal_decompose(df, freq=12).plot()

以上がデータを取得してデータを見てみる部分となります。
次の記事でモデルの作成をしていきます。