はじめに

普段気象学を中心に学んでいる身ですが、金銭的なモチベーションが欲しいなと思うことがありまして、そういえば気象要素で日経平均予測できれば目をつけている人も少なそうだし儲かりそうだな（単純）と思ったのでやってみることにしました。

データの取得

気象データは気象庁のHP からダウンロードできます。今回利用したデータは平均海面気圧、最高気温、平均気温、平均相対湿度、降水量です。
日経平均のデータは macrotrends からダウンロードしました。

前処理

まずライブラリを読み込みます

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import signal

from sklearn import linear_model
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error
from sklearn.metrics import r2_score

気象データの読み込み及び前処理です。（主に）降水量の欠損値を0埋めし、降水量を基に降水の有無のラベルを作成しました。


tokyo1 = pd.read_csv("data/tokyo_1961-1980.csv", header=2,skiprows=[4],encoding="shift-jis",parse_dates=["date"])
tokyo2 = pd.read_csv("data/tokyo_1981-2000.csv", header=2,skiprows=[4],encoding="shift-jis",parse_dates=["date"])
tokyo3 = pd.read_csv("data/tokyo_2001-.csv", header=2,skiprows=[4],encoding="shift-jis",parse_dates=["date"])
df = pd.concat([tokyo1, tokyo2, tokyo3], ignore_index=True)
df = df.sort_values("date", ignore_index=True)

df=df.drop(["wdr_pr"],axis=1)
df["year"] = df["date"].dt.year
df["month"] = df["date"].dt.month
df["day"] = df["date"].dt.day

df = df.fillna(0)
df["rain_yn"] = 0
df.loc[df["smpre"] >= 1, "rain_yn"] = 1

日経平均も読み込んでグラフ化してみます。ここではバブル崩壊が多少落ち着き始めたと思われる1995年以降のデータを用います。

nikkei = pd.read_csv("data/nikkei-225-index-historical-chart-data.csv", header=8, parse_dates=["date"],names=["date","value"])
nikkei["year"]=nikkei["date"].dt.year
nikkei=nikkei[nikkei["year"]>=1995]
plt.plot(nikkei["date"],nikkei["value"])
plt.xlabel("year")
plt.ylabel("nikkei")
plt.savefig("transition")

だいぶ長期的な変動があることがわかります。この長期的な変動はおそらく気象要因によるものではなさそうなので~~ホントはそれを予想できないと意味ないのですが、~~ハイパスフィルタを適用します。


samplerate = 6270  # 波形のサンプリングレート
x = np.arange(0, round(samplerate/2)) / samplerate  # 波形生成のための時間軸の作成

fp = 34  # 通過域端周波数[Hz]
fs = 17  # 阻止域端周波数[Hz]
gpass = 3  # 通過域端最大損失[dB]
gstop = 40  # 阻止域端最小損失[dB]

def highpass(x, samplerate, fp, fs, gpass, gstop):
    fn = samplerate / 2                           #ナイキスト周波数
    wp = fp / fn                                  #ナイキスト周波数で通過域端周波数を正規化
    ws = fs / fn                                  #ナイキスト周波数で阻止域端周波数を正規化
    N, Wn = signal.buttord(wp, ws, gpass, gstop)  #オーダーとバターワースの正規化周波数を計算
    b, a = signal.butter(N, Wn, "high")           #フィルタ伝達関数の分子と分母を計算
    y = signal.filtfilt(b, a, x)                  #信号に対してフィルタをかける
    return y

nikkei["high"] = highpass(nikkei["value"], samplerate, fp, fs, gpass, gstop)

# グラフを出力
plt.plot(np.abs(np.fft.fft(nikkei["high"])[:3000]))
plt.xlabel("frequency[Hz]")
plt.ylabel("amplitude")
plt.yscale('log')
plt.savefig("high")

このサイトを参考にさせていただきました。その結果FFTした後のグラフは以下のようになりました。

rawデータ	ハイパスフィルタ適用後

スケールが見にくいかもしれませんがだいぶ低周波の成分を取り除けました。またこの結果日経平均のグラフは以下のようになりました。

これぐらいなら気象要素で説明できるかもしれません。

データを結合し2017年までのものをtrainに、残りをtestに分割します


df=df.fillna(0)
df_merge=pd.merge(df, nikkei, on="date")
df_merge

train = df_merge[df_merge["year_y"] <= 2017]
x_train = train.drop(["high","value", "date"], axis=1)
y_train = train["high"]
test = df_merge[df_merge["year_y"] > 2017]
x_test = test.drop(["high","value", "date"], axis=1)
y_test = test["high"]

予測にはランダムフォレストを用いました。


model_fore=RandomForestRegressor(n_estimators=50,max_depth=5).fit(x_train, y_train)
y_pred=model_fore.predict(x_test)
print("test")
print("RMSE:", np.sqrt(mean_squared_error(y_test, y_pred)))
print("r^2 :", r2_score(y_test, y_pred), "\n")

print("train")
y_pred=model_fore.predict(x_train)
print("RMSE:", np.sqrt(mean_squared_error(y_train, y_pred)))
print("r^2 :", r2_score(y_train, y_pred), "\n")


dic={}
for key, value in zip(df.columns, model_fore.feature_importances_):
    dic[key]=value
sorted(dic.items(), key=lambda x: x[1], reverse=True)

結果としてはこんな感じでした。

test
RMSE: 938
r^2 : -0.031
train
RMSE: 513
r^2 : 0.124

全く使い物にならないですね。線形回帰やニューラルネットなどでも試したのですがほぼ同じ結果でした。やはり日経平均と気象要素は関係なさそうです。残念です。
ちなみに変数重要度はこんな感じでした。

変数	重要度
年	0.376
月	0.197
平均気温	0.123
日	0.088
湿度	0.067
気圧	0.066
最高気温	0.063
降水量	0.019
降水の有無	0.001

予測に使った変数も年月など、特に気象と関係ないものが重視されていることがわかりまます…
今回は全く持って予測できなかったので、今度は逆に統計学的に相関があるかどうか検証する形でリベンジしたいところです。

参考文献

PythonのSciPyでローパスフィルタをかける！