SARIMAモデルによる米国航空機の遅延時分に関する時系列解析

■はじめに

本noteをご覧いただき、ありがとうございます。
本noteは、Aidemy Premium Plan「データ分析コース」を受講した成果を記すために開設しました。

■自己紹介

まず、はじめに簡単な自己紹介をさせていただきます。

□経歴

大学卒業後、鉄道会社に就職。

□受講のきっかけ

学生時代に、少しだけSPSSをいじったこともあり、従前から「データ分析」に興味を持っていました。
そして、入社後、「いまの会社でどんな仕事を将来的にしたいのか？」と自問したところ、「数値」を扱って、多くのお客さまにとって利便性を高い運行計画を策定する業務に就きたいと思うに至りました。
そのためには、データ分析のスキルが必須だと考え、今回、受講に至りました。

■本記事の概要

□想定する読者

本記事は、Aidemyのカリキュラムの一環である「成果物作成」を目的に執筆しているものです。そのため、「Aidemyではどんなことを学ぶことができるだろう」と思っている方に手を取っていただけるとうれしいと思います。

□本記事で記述する内容

本記事では、米国航空機による到着の遅れ時間をSARIMAモデルを用いて時系列解析を行った結果について記述します。

■分析結果

それでは、さっそく、分析をしていきましょう。

□実行環境

Google Colaboratoryを用いました。

□データセットの取得

kaggleから取得しました。具体的には、"Airline Delay"と検索し、"Airline Delay Analysis"というページからデータセットを取得しました。
ここには、2009年～2020年までの米国国内におけるフライトのキャリアから始まり、発着地、到着時刻、到着にどの程度遅れたかに至るまでデータが各年ごとにCSVファイルにて格納されています。
そこで、まずは各CSVファイルをダウンロードしますが、それに先立ち、各ライブラリをインポートします。

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

続いて、ダウンロードした各CSVファイルを読み込ませます。

#データの読み込み
df_2009 = pd.read_csv("/content/drive/MyDrive/kaggle/2009.csv", header=0)
df_2010 = pd.read_csv("/content/drive/MyDrive/kaggle/2010.csv", header=0)
df_2011 = pd.read_csv("/content/drive/MyDrive/kaggle/2011.csv", header=0)
df_2012 = pd.read_csv("/content/drive/MyDrive/kaggle/2012.csv", header=0)
df_2013 = pd.read_csv("/content/drive/MyDrive/kaggle/2013.csv", header=0)
df_2014 = pd.read_csv("/content/drive/MyDrive/kaggle/2014.csv", header=0)
df_2015 = pd.read_csv("/content/drive/MyDrive/kaggle/2015.csv", header=0)
df_2016 = pd.read_csv("/content/drive/MyDrive/kaggle/2016.csv", header=0)
df_2017 = pd.read_csv("/content/drive/MyDrive/kaggle/2017.csv", header=0)
df_2018 = pd.read_csv("/content/drive/MyDrive/kaggle/2018.csv", header=0)
df_2019 = pd.read_csv("/content/drive/MyDrive/kaggle/2019.csv", header=0)
df_2020 = pd.read_csv("/content/drive/MyDrive/kaggle/20.csv", header=0)

ここで、問題が発生しました。
2009年～2020年までのデータを読み込む途中で、notebookがクラッシュしメモリ不足となってしまいました。
CSVファイルの中を確認すると、1年間だけで100万件近い数のフライトのデータが格納されてケースもあり、データ数が大きすぎるのが原因だと考えられます。
そこで、今回は2017年及び2018年の2年間に限定して分析を行うこととします。
これは、2019年以降のデータを用いると、新型コロナウイルスの影響が反映されてしまうことを考慮したためです。

#データの読み込み
df_2017 = pd.read_csv("/content/drive/MyDrive/kaggle/2017.csv", header=0)
df_2018 = pd.read_csv("/content/drive/MyDrive/kaggle/2018.csv", header=0)

□データの整理

次に、データフレームの連結を行った後、不要なカラムの削除を行うことで、データの整理を行います。
今回は、日付及び到着地への遅れ時分を必要としているため、以下のコードを用いました。

#データフレームの連結
df = pd.concat([df_2017, df_2018], axis=0)
#不要なカラムの削除
df_drop = df.drop(df.columns[range(1, 14)], axis=1)
df_drop = df_drop.drop(df.columns[range(15, 28)], axis=1)

ここで、一度、それぞれのカラムのデータの形を確認しておきます。

#データの形を表示
print(f'df_drop_shape : {df_drop.shape}')

出力結果は以下のようになり、2年分だけでも「12,888,067」件となり、どれだけ大きなデータであるかがわかりました。

df_drop_shape : (12888067, 2)

次に、各カラムのデータの型を確認します。

#各columnのデータ型を確認
print(f'{df_drop.dtypes}')

出力結果は次のようになりました。

FL_DATE       object
ARR_DELAY    float64
dtype: object

"FL_DATE"はフライト日、"ARR_DELAY"は到着地への遅れ時分を表します。
時系列解析を行うにあたり、FL＿DATEをobject型→datetime型へと変換します。

df_drop['FL_DATE'] = pd.to_datetime(df_drop['FL_DATE'])

再度、データの型の確認です。

#各columnのデータ型を確認
print(f'{df_drop.dtypes}')

無事に変換できました。

FL_DATE      datetime64[ns]
ARR_DELAY           float64
dtype: object

このままでは、データ数が多すぎて分析に適さないと考え、日付単位で遅れ時間の合計値を算出しました。この際、groupbyメソッドを用いました。

#'FL_DATE'ごとに'ARR_DELAY'の合計値sumを算出
sum = df_drop.groupby('FL_DATE').sum()

#sumを出力
print(sum)

"sum"を表示した結果は以下の通りです。

            ARR_DELAY
FL_DATE              
2017-01-01    90473.0
2017-01-02   372600.0
2017-01-03   212466.0
2017-01-04   137457.0
2017-01-05   217315.0
...               ...
2018-12-27   414647.0
2018-12-28   425751.0
2018-12-29   174401.0
2018-12-30    74027.0
2018-12-31    97465.0

[730 rows x 1 columns]

最終行の"[730 rows x 1 columns]"という出力、及び、"730 = 365 × 2"であることから、無事2年分の遅れ時分を日付単位でまとめることが出来ました。

最後に、データフレームのインデックスを"FL_DATE"に設定します。

df_drop = df_drop.set_index(df_drop['FL_DATE'])

□データの可視化

次に、ここまでまとめてきたデータをグラフという形で出力し可視化することを試みます。

#データを折れ線グラフで表す
#インデックスの設定
df_drop.index = df_drop['FL_DATE']
#グラフのタイトルを設定
plt.title("Daily-ArrDelay")
#グラフのx軸とy軸の名前設定
plt.xlabel('FL_DATE')
plt.ylabel('ARR_DELAY')

#グラフの出力
plt.plot(sum)
plt.show()

出力された図は次の通り。

周期性があるようなないような感じなので、自己相関係数を算出した上で、コレログラムを得ます。

#自己相関係数の算出
df_drop_acf = sm.tsa.stattools.acf(sum, nlags=40)
print(df_drop_acf)

print('----------\n')
#自己相関係数のグラフの出力
sm.graphics.tsa.plot_acf(sum,lags=40)
plt.show()

加えて、「トレンド」「季節性」「残差」にも分けてみました。

7from statsmodels.tsa.seasonal import seasonal_decompose

stl=seasonal_decompose(sum)

original = sum # オリジナルデータ
trend = stl.trend # トレンドデータ
seasonal = stl.seasonal # 季節性データ
residual = stl.resid # 残差データ

plt.figure(figsize=(8, 8)) # グラフ描画枠作成、サイズ指定

# オリジナルデータのプロット
plt.subplot(411) # グラフ4行1列の1番目の位置（一番上）
plt.plot(original)
plt.ylabel('Original')

# trend データのプロット
plt.subplot(412) # グラフ4行1列の2番目の位置
plt.plot(trend)
plt.ylabel('Trend')

# seasonalデータ のプロット
plt.subplot(413) # グラフ4行1列の3番目の位置
plt.plot(seasonal)
plt.ylabel('Seasonality')

# residual データのプロット
plt.subplot(414) # グラフ4行1列の4番目の位置（一番下）
plt.plot(residual)
plt.ylabel('Residuals')

plt.tight_layout() # グラフの間隔を自動調整

ここからも、Trendも、Seasonalityも、周期性を有することを確認できませんでした。
ですが、ここでは一旦、"s=7"として解析を進めることにしました。

□パラメータの探索からSARIMAモデルによる解析

# orderの最適化関数
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)]

# selectparameter(sum, 7)

# 5.モデルの構築

# SARIMAモデルを用いて時系列解析
# orderはselectparameter関数の0インデックス, seasonal_orderは1インデックスに格納
best_params = selectparameter(sum, 7)
SARIMA_sum =  sm.tsa.statespace.SARIMAX(sum, order=(0, 0, 0), seasonal_order=(0, 1, 1, 4)).fit()

# 予測
# predに予測期間での予測値
pred = SARIMA_sum.predict("2018-12-01", "2018-12-31")

# グラフを可視化してください。予測値は赤色でプロット
plt.plot(sum)
plt.plot(pred, color="r")
plt.show()

上記のコードにより実行します。ここでは、予測期間を2018年12月の1か月間としました。
出力結果は次の通りです。

赤線がSARIMAモデルにより予測された各日のフライトの遅れ時分の合計値です。
この結果を見る限り、十分に予測できていないことがうかがえます。

■考察と今後に向けて

□考察及び今後の課題

ここでは、この分析が「なぜ十分に行えなかったのか」についての考察していきます。

①対象となる期間の短さ
冒頭で述べたことと矛盾しているかもしれませんが、対象とした原データの少なさが今回の結果に影響を与えていることが考えられます。
時系列解析では原データの期間が長ければ長いほど良いと言われることがあるようですが、今回のようなデータの読み込み方法では十分な期間のデータを読み込むことが行えませんでした。
そのため、十分な予測を行えなかったと考えられます。
これを解決するために、pickle関数などを用いることが考えられそうなので、他のデータの読み込み手法を修得することが今後の課題の1つとなりそうです。

②適用したモデルの妥当性
今回はSARIMAモデルを用いて分析を行いましたが、可視化したコレログラムなどで確認したように、今回扱ったデータには明らかな周期性が認められませんでした。
そもそもSARIMAモデルはトレンドや季節性といった周期性が確認できるデータに対して適しているようなので、適用するモデルを再考する必要性があります。
ARIMAモデルなど他のモデルでも分析を行えるよう、種々の分析手法の考え方やコードの書き方を学び、目の前にあるデータや分析の目的に対して適切なモデルを適用できるように学習していく必要があることを実感しました。

□今後にむけて

今回の分析はうまくいきませんでした。
きっとこの記事をお読みの読者のみなさんには「こうすればよかったのに…」「なんで筆者はこんなことしているんだろう…」という思いを抱いていることと思います。
しかし、私自身は今回の成果物作成を通して、Pythonでのデータの収集からデータ整理方法について基礎的なスキルを身に着けることが出来たのではないかと考えています。
また、データ分析を行うにあたっては、適切なモデルを選択できるように、自分自身の知識とスキルを磨き続けなければいけないことを実感しました。
受講期間はまだまだ残っていますので、今後はこれまでの復習を通してPythonの使い方をさらにマスターするとともに、データ分析に関する知識及び技術を身に着けていきたいと考えています。
ここまでお読みいただきありがとうございました。