はじめに

初めまして。初投稿です。今後のデータ駆動社会を考えてただの営業からデータサイエンティストを目指すべく勉強を始めた40代のサラリーマンです。
これまで、Pythonの本をみたりして独学してましたが、何を書いているのかさっぱりであり、わからなくても質問もできないので基礎から勉強できチューターとの会話もできるAidemyを選択しました。

自己紹介

某IT企業に勤める営業です。最近はマーケティングの仕事に着いています。外勤営業からマーケティングに配置転換し、様々なデータに触れる機会が増えてきています。
PCは毎日触れるが、これまで一度もプログラムは勉強経験なし
データサイエンティストを目指すべく、統計検定2級、G検定は開始時には取得済み（ただし、2年前に取得したので忘れてる・・・・汗）

目的

私の住んでいる岡山県は通称’晴れの国おかやま’と呼ばれています。
ちなみに岡山県の説明をみると下記のように記載されています。
ダントツではないが、確実に晴れの日が多い。
・降水量1mm未満の日：276.7日で全国第1位
（気象庁「全国気候表（1991～2020年（30年間）の降水量1ミリ以上の日数の平年値をもとに算出）」）
https://miryoku-harenokuni-okayama.jp
今回は、過去のデータを用いて降水量の予測をしてみたいと思います。

データについて

利用したデータは下記です。
データは気象庁のホームページからダウンロードしました。
気象庁：各種データ・資料
【地点】岡山県岡山市
【項目】降水量の月合計
【期間】2011年1月〜2021年12月（11年132ヶ月分）
月毎のデータを11年分学習させて予測してみることにしました。
図は一部表記ですが、年月日、月の降水量の合計(mm)で出力しています。

実行環境

Google Colaboratory
python 3.7.12
利用モジュール

import warnings
import itertools
import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt

プログラム

1.必要なモジュールをインポート

import warnings
import itertools
import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
%matplotlib inline

2.データの読み込み

weather = pd.read_csv(
    filepath_or_buffer="/content/drive/MyDrive/monthly-weather.csv",
    dtype={'x02': 'float64'}, encoding='utf_8')

マウントの仕方もわからず苦戦しました。また、iloc等で不要なデータを削除しようとしましたが、チューターとの会話の中で実務では加工してから読み込む場合もあるとのことでしたので、あらかじめcsvを加工してから取り込んでいます。

3.データの整理

index = pd.date_range("2011-1-31", "2021-12-31",freq = "M")

期間を、2011-1-31から2021-12-31で取り込み、MonthのMで月を指定

4.インデックスにindexを代入し"Month"カラムを削除

weather.index = index

del weather["Month"]

わかりやすいようにweatherへ変換しています。

5.orderの最適化関数

ここでパラメータを選択するために、時系列データ:DATA, パラメータs(周期):sを入力すると、最も良いパラメーターとそのBICを出力するselectparameteという関数を定義していきます。

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)]

6.モデルの構築

今回降水量予測は時系列データなのでSARIMAモデルを用いて時系列解析をしていきます。

SARIMAモデルは、季節周期を持つ時系列データにも拡張できるようにしたモデルの為今回の予測に適したモデルとなっています。

best_params = selectparameter(weather, 12) 
SARIMA_weather = sm.tsa.statespace.SARIMAX(weather,order=best_params[0],seasonal_order=best_params[1]).fit()

周期は月ごとのデータであることも考慮して1年は12カ月なので、s=12としました。

orderはselectparameter関数の0インデックス, seasonal_orderは1インデックスに格納していきます。

7.予測

まずは、精度を確認するために同期間で記載

pred = SARIMA_weather.predict("2011-1-31", "2021-12-31")

8.グラフを可視化。予測値は赤色でプロット

plt.xlabel("time(month)")
plt.ylabel("Precipitation amount")
plt.plot(weather)
plt.plot(pred,color="r")
plt.show()

せっかくならったラベルも貼り付け

精度確認

出力結果が以下となります。

なんかずれてないか？(特に2011年〜2012年)と考え周期を24ヶ月にするも同じような結果となりました。
そこで、周期を36ヶ月へ変更すると以下の結果となりました。

このグラフを見ると、赤の部分の予測が正解(青の部分)に比べてある程度正しく予測することができていることが分かりました。

そこで、期間を2023-12-31で予測し直します。
したがって、6.モデルの構築と7.予測を以下に修正しています。

6.モデルの構築

best_params = selectparameter(weather, 36) 
SARIMA_weather = sm.tsa.statespace.SARIMAX(weather,order=best_params[0],seasonal_order=best_params[1]).fit()

7.予測

pred = SARIMA_weather.predict("2011-1-31", "2023-12-31")

結果

予測周期を変えるだけで、結果は予測らしいものに仕上がりました。
これだけでは、岡山県が降水量が少ない理由にはならないので、比較データ(全国平均)を入れたりすると信頼性があがるのかなと思いました。

まとめ

ある程度教材で学んだプログラムを元にデータを入れ替えたり、期間を変更したりしてみました。
結果としては、動きましたが月毎、日毎のデータ取得により全く異なる動き方をしたり、大文字・小文字を間違えるだけで動かなくなったりと完成まで何度も同じエラーを出したり、いい経験ができました。
ビジネスとなると、様々な外的要因が多数ありこのようなイメージでは使えない場面が多いと思いますが、一つの通過点として、残せたのは大きな一歩と考えてます。
今後は、元々やりたかったGoogleAnalyticsを用いて外的要因（ホームページ更新頻度、技術ブログ投稿頻度、メルマガ送信頻度等）を加味したアクセス予測及び今の広告の効果を測定する最適な広告の出稿量予測とかができればいいなと考えてます。

参考文献

時系列解析手法のSARIMAモデルを試してみた
https://qiita.com/DS27/items/b70aad44721dfa7ac5f7
気象庁　平年値（年・月ごとの値）
https://www.data.jma.go.jp/obd/stats/etrn/view/nml_sfc_ym.php?prec_no=66&block_no=47768&year=&month=&day=&view=p1