「時系列分析を用いた回帰モデルの頑健性向上 – 第1部」

「美しさとファッションの専門家による、回帰モデルの頑健性向上のための時系列分析 - 第1部」

シンガポールのHDB再販価格に関するケーススタディ。この記事では、時系列分析を線形回帰に取り入れてモデルの予測力を向上させる方法を示します。

Muhammad Faiz Zulkefleeによる写真、Unsplash

はじめに

自宅から1.5時間の距離に位置するシンガポールはいつも私を魅了しています。より大きな隣国に囲まれているにもかかわらず、この小さな国は逆境を乗り越えてきました。独立後の謙虚な始まりから、現在では重要な金融ハブと国際港湾として台頭しています。それにもかかわらず、シンガポールは土地の不足という大きな課題に直面しています。この問題により、住宅価格が著しく高騰しています。

この問題を予見して、シンガポール政府は1960年に公屋開発庁(HDB)を設立することで、住民に十分な住宅を提供することを試みてきました。この機関を通じて、政府は正確かつ衛生的で手頃な公共住宅を80%以上の人口に提供することに成功しています(HDBの公式 ウェブサイトによると)。

HDBの重要性を認識し、このブログでは、データ駆動のアプローチを使用して、公開されているデータセットを使用してシンガポールのHDB再販価格を理解するために深く掘り下げます。関連する変数からの情報の豊富さから、このデータセットは回帰モデルの構築のための潜在能力を持っているため、興味深いものです。

驚くべきことに、VoAGIやKaggleで見つけることのできる優れた分析の中で、このデータセットの時系列的な性質を深く分析している人物はほとんどいません。この分析の欠落は私にこの問題に取り組むよう余儀なくさせました。

したがって、この記事では、時系列分析と予測を回帰にシームレスに統合することを目標としています。それを念頭に置いて、この視点が既存のモデルの堅牢性を向上させ、新たな洞察を提供できればと思います。

ツール

この分析に使用したプログラミング言語はPythonであり、次のようなライブラリを使用しました(ただし、これらに限定されません)。

import pandas as pd  # 構造化データ操作
import numpy as np # 科学的・数値計算
import scipy # 科学的・数値計算
import matplotlib.pyplot as plt # データ可視化
import seaborn as sns # データ可視化
import statsmodels.api as sm # 統計分析
import datetime as dt # 日付と時刻データの操作

データの抽出、クリーニング、前処理

(この分析のデータの抽出から回帰モデルの構築までのエンドツーエンドのプロセスの完全なコードについては、こちらのGitHubリポジトリをご覧ください)。

公開されていて無料で利用できるHDB再販価格のデータセットは、この ウェブサイトからダウンロードできます。CSVファイルを入手したら、パンダのDataFrameにデータをロードし、次の表が表示されるようにします。

# データセットをデータフレームにロードします
# この場所でデータセットをダウンロードしました:https://data.gov.sg/dataset/resale-flat-prices
# 日付は2017年1月から2023年8月の間としました
data = pd.read_csv('Housing_Resale_Dataset.csv')
data.head()
Raw Dataset — Image by Author

元のデータセットは160,795行(取引数)と11列(変数の数)からなります。HDB再販価格を予測することに重点を置いているため、「resale_price」という列は従属変数またはターゲットとなり、他の列は独立変数または特徴となります。データセットはすでに比較的整然としていましたが、さらなるデータの清掃と前処理を行い、列名の変更、新しい特徴の抽出、重複の削除、および整理を行って、次のような最終的なクリーンなデータセットを取得しました。

cleandata.head()
クリーンなデータセット—著者のイメージ

上記の画像から、’resale_date’列(以前は’month’という名前で表示されず、画像には表示されていません)の情報を抽出して、’resale_year’と’resale_month’という2つの新しい特徴量を作成することができます。このため、クリーンなデータセットには13列があります。同時に、行数はわずかに減少し、160,454行になります。これは重複の削除の結果です。

探索的データ分析

次に、データの中で最も重要な情報を明らかにするために、可視化を作成します。可視化は主にJupyter NotebookでMatplotlibとSeabornを使用して直接分析しやすく、柔軟性があるため、このブログのこの部分では、見やすく読者にやさしいプレゼンテーションのためにTableauで生成された画像も使用しています。

価格分布

まずはじめに、以下のヒストグラムからHDB再販価格を見てみましょう:

価格分布—著者のイメージ
cleandata['resale_price'].describe()
記述統計の概要—著者のイメージ

ビンのサイズを40,000で設定した場合、’resale_price’の分布は右に偏っており、平均がS$486,519で中央値がS$455,000であることがわかります。

町別の平均価格

棒グラフを使用して、さまざまな町のフラットの平均価格を確認します。紫の棒はHDBフラットの平均価格が最も高い5つの町を示しています。この図からも分かるように、平均価格には町によってかなりのばらつきがあることがわかります。

町別の平均価格—著者のイメージ

月次平均価格

このブログの主な目標は、時系列分析が回帰モデルを向上させる方法を紹介することですので、データの時系列的な側面も可視化することが重要です。以下は2017年1月から2023年8月までのHDBフラットの月次平均価格です。ブルーの破線は実際の月次平均価格を表し、紫の線は6か月移動平均を表しています。

月次平均価格—著者のイメージ

青い線をよく見ると、ジグザグのパターンが見られ、データセットには月次の季節性がある可能性があります。一方、紫の線はデータのトレンドを示しています。紫の線から、データセットには非定常な正のトレンドがあることが分かります。月次平均価格が2018年中ごろに下落し、2019年初めに安定したが、2020年中ごろに平均価格が急上昇し始めました。

月次取引

以下の画像は2017年1月から2023年8月までの月次取引を示しています。グラフからもわかるように、取引データには季節性があり、12月と1月の周辺では通常、月次取引が低下します。また、2020年の第1四半期に取引が急激に減少したという異常事態も見えます。この急激な減少は、Covid-19のロックダウンによるものと考えられます。

月次取引 — 著者による画像

時間系列分析

Jake Hills氏撮影 写真 — Unsplash

前述の通り、このストーリーではHDBフラットの価格設定を理解することに焦点を当てています。そのため、時間系列分析の前提条件として、まずは「resale_date」をインデックスとして使用し、月ごとの平均価格を含むPandasシリーズの時間系列データセットを作成します。

price_time = cleandata.groupby('resale_date')['resale_price'].mean()price_time.head()
月次平均価格 — 著者による画像

分解

月次平均価格の可視化から、データの傾向や季節性を観察しました。次に、時間系列のデータを傾向、季節性、残差の構成要素に分解するために、時間系列分解を使用します。

import statsmodels.api as sm# 季節性は時間とともに増大しないように「加法モデル」を使用# 先頭と末尾の6ヶ月にトレンドを拡張するために6つの予測期間を追加price_decomp = sm.tsa.seasonal_decompose(price_time, model = 'additive', extrapolate_trend=6) price_decomp.plot();
時間系列分解 — 著者による画像

単純な移動平均により、時間系列データの傾向を既に確認しました。一方、分解を行うことで、データが12ヶ月の周期で季節性を持っていることも分かります。上記のグラフからは、季節性は平均価格に対して小さく、-4,000から2,500の範囲内にあります。見やすくするために、季節性のグラフを分離し、折れ線グラフを棒グラフに変更します:

price_season = price_decomp.seasonal # 季節性 price_season_month = price_season.reset_index()price_season_month['resale_month'] = price_season_month['resale_date'].dt.strftime('%b') # 月の情報を抽出months_order = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']price_season_month['resale_month'] = pd.Categorical(price_season_month['resale_month'], categories= months_order, ordered = True) # データ上で月の順序を作成pm = price_season_month.groupby('resale_month')['seasonal'].mean()pm.plot(kind = 'bar')
月次平均価格の季節性 — 著者による画像

データ分割

信頼性を評価できるモデルを構築する前に、データ分割と呼ばれる方法を使用します。このプロセスでは、データをトレーニングデータセットとテストデータセットの2つに分割します。トレーニングデータセットを使用して時間系列モデルを構築し、テストデータセットを使用してモデルのパフォーマンスを評価します。したがって、時間系列モデリングの文脈では、トレーニングデータセットが「過去のデータ」を、テストデータセットが「未来のデータ」を表します。

# モデルのパフォーマンスを向上させるために対数変換を使用する
transform = np.log(price_time)
# Train-Testスプリット
price_test = price_time['Sep 2022':] # テストデータセット(直近12ヶ月)
price_train = price_time[:'Aug 2022'] # トレーニングデータセット(直前のデータ)

データセットは2017年1月から2023年8月までの期間にわたります。テストデータセットとしては、最後の12ヶ月のデータを、トレーニングデータセットとしては直前のデータを選択します。12ヶ月を選ぶ理由について疑問に思うかもしれませんが、これは実際的な目的に沿った選択です。12ヶ月の期間は、予測の振る舞いを観察するのに十分な長さを提供します。一方、その期間は実際のデータとあまりにも大きく乖離することはありません。

タイムシリーズモデリング
事前に得た情報に基づいて、タイムシリーズモデルを構築することができます。最終的な目標は、モデルの理解を活用して、将来のデータの振る舞いを予測することです。このプロセスを予測と呼びます。
様々なタイムシリーズモデリング方法があります。古典的な方法には、移動平均、指数平滑法、ARIMA(自己回帰和分移動平均)などがあります。より高度な手法には、Facebook ProphetおよびLSTMネットワークがあります。この特定のコンテキストでは、古典的な方法がかなりうまく機能し、その方法の中でもARIMAの拡張であるSARIMA(季節自己回帰和分移動平均)が最も優れたパフォーマンスを発揮します。

SARIMAモデルを構築するためには、ARIMAモデルのp、d、qと季節性のP、D、Q、sのパラメータを決定する必要があります。パラメータの説明はこちらとこちらで読むことができます。パラメータを決定するためのより厳密な方法はありますが、この特定のインスタンスでは、グリッドサーチの実用的なアプローチを使用し、結果のAIC(赤池情報量基準)スコアが最も低いモデルを選択します。

from statsmodels.tsa.statespace.sarimax import SARIMAX # statsmodelsからSARIMA関数をインポート
s = 12 # 季節性は12ヶ月の周期からなる
for p in range(1, 4):
    for d in range(0, 3):
        for q in range(0, 4):
            for P in range(1, 4):
                for D in range(0, 3):
                    for Q in range(0, 4):
                        model = SARIMAX(price_train, order=(p, d, q), seasonal_order=(P, D, Q, s))
                        results = model.fit()
                        print(f'ARIMA{p},{d},{q}-SARIMA{P},{D},{Q}-{s} - AIC: {results.aic}')

このデータに対して、最適なパラメータはARIMAパラメータの3、1、3、および季節性パラメータの1、1、3、12です。以下のコードブロックでは、ライブラリStatsmodelsを使用してトレーニングデータセットをモデルにフィットさせる方法が示されています。

p,d,q = 3,1,3
P,D,Q,s = 1,1,3,12
sarima_model = SARIMAX(price_train, order=(p, d, q), seasonal_order=(P, D, Q, s))
sarima_result = sarima_model.fit()

予測
モデルを取得した後、次にするべきことは、残差プロットをチェックすることです。以下の画像のように、モデルのパラメータをどのように調整しても、最初と13番目のインデックスの残差値が常にゼロから非常に大きく乖離している理由はまだ説明できません。この問題は試みても解決されません。

sarima_result.resid.plot()

それにもかかわらず、データの残りの部分は通常の振る舞いを示しています。したがって、この問題に取り組むために、私は13番目のインデックスを超えるデータを使用してモデルを分析することにしました。この実用的なアプローチにより、データセットのより安定した部分で作業することができ、極端な外れ値の影響を回避することができました。
{
  "html": "
sarima_result.resid[13:].plot()
"モデルの残差プロット、外れ値なし
モデルの残差プロット、外れ値なし ― 著者によるイメージ

それでは、SARIMAモデルの構築が完了したので、次に訓練データセットとテストデータセットの両方で予測を作成し、実際のデータセットと一緒にプロットして視覚的に評価します。

"実際の値
実際の値 vs 予測値、時系列分析 ― 著者によるイメージ

モデル評価

上記のグラフを見ると、訓練データセットとテストデータセットから得られた予測値と実際の値の間には顕著な類似性があります。予測値はパターンと傾向についても実際の値に密接に従っています。したがって、視覚的な評価に基づいて、SARIMAモデルは訓練データセット内のデータをうまく予測することができます。さらに、この精度はテストデータセットにも適用されます。つまり、モデルは未知のデータに対してもかなり良い汎化能力を持っていると言えます(少なくとも、12ヶ月の期間での話です)。

視覚的な評価に加えて、時系列モデルの正確さを評価するために特定の指標を使用することもできます。そのうちの1つが平均絶対誤差(MAE)です。MAEは、データセット内の予測値と真の値の平均の差の絶対値を取ることによって計算します。要するに、MAEは、モデルが真の値にどれだけ近いかを平均して測定します。MAEが低いほど、モデルの予測が実際のデータに近いことを意味し、MAEが高いほど、モデルの誤差が大きいことを意味します。MAEを計算するためには、Scikit-learnライブラリで提供されている関数を使用することができます。

from sklearn.metrics import mean_absolute_error as MAE# do not forget to reverse our log transformation MAE_train = MAE(price_train[13:].apply(reverse_transform),train_predict[13:].apply(reverse_transform))MAE_test = MAE(price_test.apply(reverse_transform),test_predict.apply(reverse_transform))print(f'MAE of training dataset = {MAE_train:.2f}')print(f'MAE of testing dataset = {MAE_test:.2f}')
"平均絶対誤差
平均絶対誤差 ― 著者によるイメージ

上記の計算から、訓練データセットとテストデータセットのMAEはそれぞれS$5748とS$7549です。MAEの値を比較すると、時間系列データの平均値であるS$481,070と比較して、これらの値はかなり小さいことがわかります。これは、モデルが真の値からわずかにずれているだけであることを示しています。

さらに、テストデータセットのMAEは、訓練データセットのMAEとほとんど変わらないです。この微小な差は、私たちの時系列モデルが、再度、少なくとも短期間では、未知のデータにもうまく汎化できることをさらに確認しています。

次回に続く...

最後に、私たちの時系列モデルが完成したので、この物語の第二部に入る準備が整いました。次のセグメントでは、時系列分析から得た洞察を活用し、回帰モデルの信頼性を向上させる方法をお見せします。

さらなる読み物:

[1] 住宅開発庁、HDBの歴史と町、2023年。〔オンライン〕利用可能:https://www.hdb.gov.sg/cs/infoweb/about-us/history [アクセス日:2023年10月30日]

"
}

[2] J. Brownlee、Pythonにおける時系列予測のためのARIMAモデルの作成方法、2020年。[オンライン]。利用可能:https://machinelearningmastery.com/arima-for-time-series-forecasting-with-python/ [アクセス日:2023年10月30日]

[3] J. Brownlee、Pythonにおける時系列予測のためのSARIMAの基礎的な紹介、2019年。[オンライン]。利用可能:https://machinelearningmastery.com/sarima-for-time-series-forecasting-in-python/ [アクセス日:2023年10月30日]

We will continue to update VoAGI; if you have any questions or suggestions, please contact us!

Share:

Was this article helpful?

93 out of 132 found this helpful

Discover more