LSTMモデルで不動産取引価格の長期予測

###目次
1.検証内容
2.検証動機
3.開発環境
4.開発言語
5.学習用データセット
6.達成目標
7.前処理
8.東京と宮城の各間取りごとの時系列グラフ
9.学習モデルの作成と検証
10.学習モデルを用いた長期予測
11.まとめ
12.今後の展望

###検証内容

LSTMモデルを用いて、不動産取引価格の長期予測をしてみます。
＊学習用のデータは国土交通省の不動産取引価格から宮城と東京分を四半期ごとダウンロード
（2005年〜2021年）
日本を代表する都市である東京と地元の宮城をセレクト

###検証動機
時系列解析で長期予測に取り組んでみたいと思っていましたが、初学者でもあるので、あまり難しくない長期予測解析をしようと思いました。
不動産価格の長期予測をテーマに選んだのは、時系列解析への取組みやすさとして、敷居が高くない題材だと思ったからです。

###開発環境

MacBook Pro（GPUなし）
Google colab(ランタイムタイプ:GPU)

###開発言語
Python3.7.11
tensorflow 2.6.0
numpy 1.19.5
scikit-learn 0.22.2
mojimoji 0.0.11

###学習用データセット
国土交通省不動産取引情報
各都道府県の不動産取引情報が任意の期間分ダウンロードできる

###達成目標
オープンデータセットを使い学習済みのLSTMモデルを生成し、部屋の間取りごとに
取引価格の長期予測をしてみる

###前処理
それでは前処理を始めます。

必要なライブラリをインポート

import pandas as pd
import mojimoji
import datetime
import matplotlib.pyplot as plt

解析に必要な宮城と東京の不動産ファイルの読み込み

df_miyagi = pd.read_csv('/content/drive/MyDrive/RealEstate_Price/Miyagi_20053_20211.csv',engine="python", encoding='cp932') 
df_tokyo = pd.read_csv('/content/drive/MyDrive/RealEstate_Price/Tokyo_20053_20204.csv', engine="python", encoding='cp932')　

間取りごとにグループ分け
row.'列名'　in [セルに格納してある部屋の名前]　で任意のグループを作る
予測対象としては1K〜5LDKを集計・予測対象とする。それ以外はその他として今回の予測には使用しないこととします。

def add_floorplan_group(df):
  floorplan_group_list =[]
  for row in df.itertuples():
    if row.間取り in ['１Ｋ', '１Ｒ']:
      floorplan_group_list.append('1floor')
    elif row.間取り in ['１ＤＫ', '１ＬＤＫ']:
      floorplan_group_list.append('2floors')
    elif row.間取り in ['２ＤＫ', '２ＬＤＫ', '３ＤＫ', '３ＬＤＫ', '４ＤＫ', '４ＬＤＫ', '５ＤＫ', '５ＬＤＫ']:
      floorplan_group_list.append('over_2floors')
    else:
      floorplan_group_list.append('others')
  df['floorplan_group'] = floorplan_group_list
  return df

mojimojiライブラリであらかじめデータセットにある西暦と四半期ごとの数字を全角から半角へ変換し、
さらに、if文で四半期ごとの月を3ヶ月括りに変換

def transform_tradingtime(tt):
  try:
    year = mojimoji.zen_to_han(tt[:4])
    quarter = mojimoji.zen_to_han(tt[6])
    if quarter == '1':
      month = '01'
    elif quarter == '2':
      month = '04'
    elif quarter == '3':
      month = '07'
    elif quarter == '4':
      month = '10'
    else:
      return '0000-00-00'
    return '{}-{}-{}'.format(year, month, '01')
  except:
    return '0000-00-00'

東京と宮城の間取りごとのグループをカウントし昇順に並べかえ
まずは宮城県を間取りごとにまとめ、総数をカウントし、昇順並べかえ。

add_floorplan_group(df_miyagi).groupby("floorplan_group")['floorplan_group'].count().sort_values(
ascending=False)

間取りグループ	合計
その他	9235
3部屋〜５部屋	3461
２部屋	523
１部屋	304

宮城県は３〜５部屋の間取りが一番多いようです。
続いて東京の間取り数はどうなっているかというと、

add_floorplan_group(df_tokyo).groupby("floorplan_group")['floorplan_group'].count().sort_values(ascending=False)

間取りグループ	合計
その他	255360
3部屋〜５部屋	107607
1部屋	59923
2部屋	26406

やはり東京はデータ数が多い。宮城とは桁違い。3部屋〜５部屋がこちらも多いようです。

続いて間取りごとのグループをdf-miyagi,df_tokyoに格納し、
to_datetimeで取引時点の列をタイムスタンプに変換

df_miyagi = add_floorplan_group(df_miyagi)
df_tokyo = add_floorplan_group(df_tokyo)
df_miyagi['取引時点'] = pd.to_datetime(df_miyagi['取引時点'].apply(transform_tradingtime))
df_tokyo['取引時点'] = pd.to_datetime(df_tokyo['取引時点'].apply(transform_tradingtime))

get_group() の引数内のグループで取引時点での取引価格平均を計算し各price_へ格納させる

price_miyagi_1floor = df_miyagi.groupby('floorplan_group').get_group('1floor').groupby('取引時点')['取引価格（総額）'].mean()
price_miyagi_2floors = df_miyagi.groupby('floorplan_group').get_group('2floors').groupby('取引時点')['取引価格（総額）'].mean()
price_miyagi_over_2floors = df_miyagi.groupby('floorplan_group').get_group('over_2floors').groupby('取引時点')['取引価格（総額）'].mean()
price_tokyo_1floor = df_tokyo.groupby('floorplan_group').get_group('1floor').groupby('取引時点')['取引価格（総額）'].mean()
price_tokyo_2floors = df_tokyo.groupby('floorplan_group').get_group('2floors').groupby('取引時点')['取引価格（総額）'].mean()
price_tokyo_over_2floors = df_tokyo.groupby('floorplan_group').get_group('over_2floors').groupby('取引時点')['取引価格（総額）'].mean()

#東京と宮城の各間取りごとの時系列グラフ

次は、集計した部屋ごとの平均取引価格を時系列順にグラフにしてみます。

plt.figure(figsize=(16, 9))
plt.plot(price_tokyo_1floor.index, price_tokyo_1floor.values, label='tokyo_1floor')
plt.plot(price_tokyo_2floors.index, price_tokyo_2floors.values, label='tokyo_2floors')
plt.plot(price_tokyo_over_2floors.index, price_tokyo_over_2floors.values, label='tokyo_over_2floors')
plt.plot(price_miyagi_1floor.index, price_miyagi_1floor.values, label='miyagi_1floor')
plt.plot(price_miyagi_2floors.index, price_miyagi_2floors.values, label='miyagi_2floors')
plt.plot(price_miyagi_over_2floors.index, price_miyagi_over_2floors.values, label='miyagi_over_2floors')
plt.legend()

東京の1部屋(tokyo_1floor)の取引価格平均は宮城の3〜5部屋(over_2floor)と同じ傾向。物価の違いが確認される。
データ数として、宮城が少ないため、東京の方で長期予測は向いていると思われました。更に時系列グラフの傾向を確認すると、東京の1部屋（tokyo_1floor）が下がったり上がったりと増減傾向が見受けられるため、検証サンプルとして、面白そうと思い検証データセットとして選択してみた。が、しかし、グラフの凹凸がどのグループも酷く予測に影響しそうです。ここは、一手間かけて移動平均を使って再集計しようと思います。

rollingメソッドで４期分移動平均し、dropnaでnanを除去

window_size = 4
price_tokyo_1floor_R = price_tokyo_1floor.rolling(window=window_size).mean().dropna()
price_tokyo_2floors_R = price_tokyo_2floors.rolling(window=window_size).mean().dropna()
price_tokyo_over_2floors_R = price_tokyo_over_2floors.rolling(window=window_size).mean().dropna()
price_miyagi_1floor_R = price_miyagi_1floor.rolling(window=window_size).mean().dropna()
price_miyagi_2floors_R = price_miyagi_2floors.rolling(window=window_size).mean().dropna()
price_miyagi_over_2floors_R = price_miyagi_over_2floors.rolling(window=window_size).mean().dropna()

移動平均のグラフを表示します。

plt.figure(figsize=(16, 9))
plt.plot(price_tokyo_1floor_R.index, price_tokyo_1floor_R.values, label='tokyo_1floor')
plt.plot(price_tokyo_2floors_R.index, price_tokyo_2floors_R.values, label='tokyo_2floors')
plt.plot(price_tokyo_over_2floors_R.index, price_tokyo_over_2floors_R.values, label='tokyo_over_2floors')
plt.plot(price_miyagi_1floor_R.index, price_miyagi_1floor_R.values, label='miyagi_1floor')
plt.plot(price_miyagi_2floors_R.index, price_miyagi_2floors_R.values, label='miyagi_2floors')
plt.plot(price_miyagi_over_2floors_R.index, price_miyagi_over_2floors_R.values, label='miyagi_over_2floors')
plt.legend()

だいぶ穏やかなグラフになりました。
傾向は移動平均をとっても変わりません。
これを元に長期予測をしたいと思います。

###学習モデルの作成と検証

続いてLSTMモデルを学習させ、モデルの検証をします。

import pandas as pd
import numpy as np
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error
import tensorflow as tf
from tensorflow import keras
from tensorflow.keras import layers
import math
from keras.callbacks import EarlyStopping
import random
import os

#乱数を固定させる
seed = 0
tf.random.set_seed(seed)
np.random.seed(seed)
random.seed(seed)
os.environ["PYTHONHASHSEED"] = str(seed)

# 入力データ・正解ラベルを作成する関数を定義
def create_dataset(dataset, look_back):
    data_X, data_Y = [], []
    for i in range(look_back, len(dataset)):
        data_X.append(dataset[i-look_back:i, 0])
        data_Y.append(dataset[i, 0])
    return np.array(data_X), np.array(data_Y)


# データセットの読み込み
dataset = price_tokyo_1floor_R.values
dataset = dataset.astype('float32') #kerasの仕様によりfloatへ変換する必要あり

# trainデータとtestデータの割合を決める
train_size = int(len(dataset) * 0.67)
train, test = dataset[0:train_size], dataset[train_size:len(dataset)]

# データのスケーリング（正規化）
scaler = MinMaxScaler(feature_range=(0, 1))
scaler_train = scaler.fit(train.reshape(-1, 1))
train = scaler_train.transform(train.reshape(-1, 1))
test = scaler_train.transform(test.reshape(-1, 1))

# 入力データと正解ラベルを作成
look_back = 18
train_X, train_Y = create_dataset(train, look_back)
test_X, test_Y = create_dataset(test, look_back)

# データの整形
train_X = train_X.reshape(train_X.shape[0], train_X.shape[1], 1)
test_X = test_X.reshape(test_X.shape[0], test_X.shape[1], 1)
test_Y = test_Y.reshape(-1, 1)

# LSTMモデルの作成と訓練
model = keras.Sequential()
model.add(layers.LSTM(16, return_sequences=False, input_shape=(look_back, 1)))
model.add(layers.Dense(1))
model.compile(loss='mean_squared_error', optimizer='adam')
es = EarlyStopping(monitor = 'val_loss', min_delta = 0.0000, patience = 5)
print(test_X.shape)
print(test_Y.shape)
model.fit(train_X, train_Y, epochs=500, validation_data=(test_X, test_Y), batch_size=1, callbacks = [es], verbose=1)

# 予測データを作成
train_predict = model.predict(train_X)

# Future prediction（未来予測）
future_test = train_X[-1, :, :]    
future_result = []

for i in range(len(dataset)-train_size):
    test_X_tmp = np.reshape(future_test, (1, look_back, 1))
    batch_predict = model.predict(test_X_tmp)
    future_test = np.delete(future_test, 0)
    future_test = np.append(future_test, batch_predict)
    future_result = np.append(future_result, batch_predict)
print(future_result)

# スケールしたデータを元に戻す
train_predict = scaler_train.inverse_transform(train_predict)
train_Y = scaler_train.inverse_transform([train_Y])
test_Y = test_Y.reshape(-1,1)
test_Y = scaler_train.inverse_transform(test_Y)
future_result = scaler_train.inverse_transform(future_result.reshape(-1, 1))

# 予測精度の計算
train_score = math.sqrt(mean_squared_error(train_Y[0], train_predict[:, 0]))
print('Train Score: %.2f RMSE' % (train_score))

tmp = future_result[:, 0]
future_score = math.sqrt(mean_squared_error(dataset[train_size:], future_result[:, 0]))
print('Future  Score: %.2f RMSE' % (future_score))



# プロットのためのデータ整形
train_predict_plot = np.empty_like(dataset)
train_predict_plot[:] = np.nan #初期化したいが、numpy配列にはデータなしは格納できないためnanを格納
train_predict_plot[look_back:len(train_predict)+look_back] = train_predict[:, 0]
test_predict_plot = np.empty_like(dataset)
test_predict_plot[:] = np.nan
test_predict_plot[train_size:] = future_result[:, 0]

# データのプロット
plt.figure(figsize=(12, 8), dpi=100)
plt.title("Real_estate_Price_Prediction")
plt.xlabel("Quarter")
plt.ylabel("Mean_TransactionPrice")
# 読み込んだままのデータをプロット
plt.plot(dataset, label='dataset')
# トレーニングデータから予測した値をプロット
plt.plot(train_predict_plot, label='train_predict')
# テストデータから予測した値をプロット
plt.plot(test_predict_plot, label='test_predict')

plt.legend()
plt.show()

試しに東京１部屋（price_tokyo_1floor_R）の学習モデルを作成し、予測してみました。
中間層がない方が精度としては良く、look_backは高ければ高いほど良い傾向にありました。
グラフにすると、若干オーバーシュート気味ですが、そこそこ妥当な予測ができました。
testデータの予測（test_predict）結果はdatasetの少し前のデータを反映しているので、その傾向を引き継いだと思われます。

平均２乗誤差結果
Train Score: 250458円
Future Socre: 704364円

###学習モデルを用いた長期予測

次に宮城の3部屋〜5部屋（price_miyagi_over_2floors_R）を用いて長期予測をしてみます。
またバリデーションデータも取り入れて、実際のデータとの差異も検証確認してみます。

平均２乗誤差結果
Train Score: 852031円
Future Score: 1488320円

バリデーションデータはデータ範囲の都合上、２点しか取れなかったのですが、まあ、、、良しとしたいです。過学習はなっていない！と思いたいです。また、長期予測にしては、現状の傾向を掴んだような結果となりました。

import pandas as pd
import numpy as np
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error
import tensorflow as tf
from tensorflow import keras
from tensorflow.keras import layers
import math
from keras.callbacks import EarlyStopping
import random
import os

#乱数を固定させる
seed = 0
tf.random.set_seed(seed)
np.random.seed(seed)
random.seed(seed)
os.environ["PYTHONHASHSEED"] = str(seed)

# 入力データ・正解ラベルを作成する関数を定義
def create_dataset(dataset, look_back):
    data_X, data_Y = [], []
    for i in range(look_back, len(dataset)):
        data_X.append(dataset[i-look_back:i, 0])
        data_Y.append(dataset[i, 0])
    return np.array(data_X), np.array(data_Y)


# データセットの読み込み
dataset = price_miyagi_over_2floors_R.values
dataset = dataset.astype('float32') #kerasの仕様によりfloatへ変換する必要あり
# 
train_size = int(len(dataset) * 0.67)
train, test = dataset[0:train_size], dataset[train_size:len(dataset)]

# データのスケーリング（正規化）
scaler = MinMaxScaler(feature_range=(0, 1))
scaler_train = scaler.fit(train.reshape(-1, 1))
train = scaler_train.transform(train.reshape(-1, 1))
test = scaler_train.transform(test.reshape(-1, 1))

# 入力データと正解ラベルを作成
look_back = 18
train_X, train_Y = create_dataset(train, look_back)
test_X, test_Y = create_dataset(test, look_back)

# データの整形
train_X = train_X.reshape(train_X.shape[0], train_X.shape[1], 1)
test_X = test_X.reshape(test_X.shape[0], test_X.shape[1], 1)
test_Y = test_Y.reshape(-1, 1)

# LSTMモデルの作成と訓練
model = keras.Sequential()
model.add(layers.LSTM(16, return_sequences=False, input_shape=(look_back, 1)))
model.add(layers.Dense(1))
model.compile(loss='mean_squared_error', optimizer='adam')
es = EarlyStopping(monitor = 'val_loss', min_delta = 0.0000, patience = 5)
print(test_X.shape)
print(test_Y.shape)
model.fit(train_X, train_Y, epochs=500, validation_data=(test_X, test_Y), batch_size=1, callbacks = [es], verbose=1)

# 予測データと検証データを作成
train_predict = model.predict(train_X)
test_predict = model.predict(test_X)

# Future prediction（未来予測）
future_test = train_X[-1, :, :]    
future_result = []
future_predict_quarters = 40 #len(dataset)-train_size　← 長期予測をしなくてもいい場合はこちらを使用

for i in range(future_predict_quarters):
    test_X_tmp = np.reshape(future_test, (1, look_back, 1))
    batch_predict = model.predict(test_X_tmp)
    future_test = np.delete(future_test, 0)
    future_test = np.append(future_test, batch_predict)
    future_result = np.append(future_result, batch_predict)
print(future_result)

# スケールしたデータを元に戻す
train_predict = scaler_train.inverse_transform(train_predict)
test_predict = scaler_train.inverse_transform(test_predict)
train_Y = scaler_train.inverse_transform([train_Y])
test_Y = test_Y.reshape(-1,1)
test_Y = scaler_train.inverse_transform(test_Y)
future_result = scaler_train.inverse_transform(future_result.reshape(-1, 1))

# 予測精度の計算
train_score = math.sqrt(mean_squared_error(train_Y[0], train_predict[:, 0]))
print('Train Score: %.2f RMSE' % (train_score))
tmp = future_result[:, 0]
future_score = math.sqrt(mean_squared_error(dataset[train_size:len(dataset)], future_result[:len(dataset)-train_size, 0]))
print('Future  Score: %.2f RMSE' % (future_score))

print(test_X.shape)
print(test_Y.shape)

# プロットのためのデータ整形
train_predict_plot = np.empty(train_size+future_predict_quarters) #長期予測の期間分train_sizeにプラスする
train_predict_plot[:] = np.nan #初期化したいが、numpy配列にはデータなしは格納できないためnanを格納
train_predict_plot[look_back:len(train_predict)+look_back] = train_predict[:, 0]
test_predict_plot = np.empty(train_size+future_predict_quarters)
test_predict_plot[:] = np.nan
test_predict_plot[train_size+look_back:len(dataset)] = test_predict[:, 0]
future_predict_plot = np.empty(train_size+future_predict_quarters)
future_predict_plot[:] = np.nan
future_predict_plot[train_size:] = future_result[:, 0]

# データのプロット
plt.figure(figsize=(12, 8), dpi=100)
plt.title("Real_estate_Price_Prediction")
plt.xlabel("Quarter")
plt.ylabel("Mean_TransactionPrice")
# 読み込んだままのデータをプロット
plt.plot(dataset, label='dataset')
# トレーニングデータから予測した値をプロット
plt.plot(train_predict_plot, label='train_predict')
# テストデータから予測した値をプロット
plt.plot(test_predict_plot, label='test_predict')
# トレーニングデータから長期予測した値をプロット
plt.plot(future_predict_plot, label='future_predict')

plt.legend()
plt.show()

###まとめ
長期予測においては、宮城と東京それぞれ単独でのデータ数が不足しており、与えらたデータ範囲が限定されている状況でした。もっとデータ数があればより良い学習モデルやそれに伴う適切な長期予測ができたのではないかと思いました。バリデーションデータも２期分しか計算できず、申し訳ないレベルです。look_backをこれ以上増やすとできなくなり、減らすと性能が悪くなるのでこれが限界なのかなと思いました。

###今後の展望
今回、時系列解析をやってみましたが、データの内容としては、東京と宮城のそれぞれのデータからの予測です。ですので、今後は全県のデータを学習モデルとし、それを用いて必要な長期予測をしていきたいなと思いました。