時系列データ分析

この記事の狙い・目的

機械学習を取り入れたAIシステムの構築は、
①データセット作成（前処理）→ ②モデルの構築 → ③モデルの適用
というプロセスで行っていきます。
その際「データセット作成（前処理）」の段階では、正しくモデル構築できるよう、事前にデータを整備しておくことが求めらます。
このブログでは、時系列を含むデータに対しての解析方法、および前処理段階で用いられる「時系列データ分析」の手法について解説していきます。

プログラムの実行環境

Python3
MacBook pro（端末）
PyCharm（IDE）
Jupyter Notebook（Chrome）
Google スライド（Chrome）

時系列データとは

時系列データとは、時間の順序にしたがって並べられたデータのことをいいます。
例えば、株価の値動き、気温や降水量などの気象情報、交通量の推移などがこれにあたります。

時系列データがもつ情報

時系列データには大きく分けて「傾向変動」「季節変動」「循環変動」「不規則変動」の４つの情報（要因）を持っています。ひとつひとつ解説していきます。

1.傾向変動（トレンド）

トレンドとは、時系列の⻑期的傾向のことです。時間の経過とともに増加・減少する傾向とも言えます。
トレンドの推定方法には「移動平均」などがあります。詳しくは後述します。

2.季節変動（シーズナル）

季節変動とは、（通常）1年を周期とする規則的な変動のことです。

3.循環変動（サイクル）

傾向変動より短期的で、周期的に繰り返される変動のことです。

4.不規則変動（ノイズ）

トレンド、季節変動、循環変動では説明できない、短期的かつ不規則な変動のことです。

これらを確認するため「季節調整」という手法をまずは簡単に見ていきましょう。

# 成分分解（季節調整）
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
%matplotlib inline

# データ取得
stock_df = pd.read_csv("./Tesla.csv")

# Dataをインデックスに設定
time_series = stock_df.copy()
time_series["Date"] = pd.to_datetime(stock_df['Date'])
time_series = time_series.set_index(['Date']).sort_index(ascending=True)

# 表示サイズ調整
plt.rcParams["figure.figsize"] = [20,10]
plt.rcParams["font.size"] = 15

# 原系列、傾向（トレンド）、季節変動（季節性）、不規則変動（残差）
def plt_seasonal_decompose(df, col_name, period):
    res = sm.tsa.seasonal_decompose(df[col_name], period=period)
    fig = res.plot()
    plt.show()

# 原系列をインプットして確認
plt_seasonal_decompose(time_series, col_name='Close', period=100)

1番目のグラフが原系列、2番目のグラフが傾向変動（Trend）、3番目のグラフが季節変動（Seasonal）、4番目のグラフが不規則変動（Resid（残差））を表しています。

時系列データにおける重要な概念

時系列データ分析を始める前に、まずは時系列データにおける重要な概念を解説していきます。
「定常性」「単位根過程」「自己共分散」の３つを紹介します。

1.定常性

定常性とは、時間経過による影響を（強く）受けない状態のことです。
また、定常過程とは、その確率的な変動の性質が、各時点に依存せず一定であるような確率過程のことを言います。
具体的には、各期の平均が普遍であるとみなせる、上昇（下降）の傾向がないデータのことです。
例えば、携帯アラームの周波数など、平均0の周辺で一定の分散をもってランダムに振動を続けるデータは、（弱）定常性があると言えます。反対に、明らかな上昇傾向が存在し変化し続けている株価の値動きなどは、定常性があるとは言えません。

・定常過程の例

・非定常過程の例

非定常過程のグラフは、グラフの右側へ進むにつれて変動（上昇傾向）があることがわかります。

時系列データについての説明のため、先に基本の4成分（傾向変動等）に触れましたが、（成分分解をする前に）そもそも定常性があるのかを確認する必要があるため、まずは定常性の条件である「平均」と「分散」を確認しておきましょう。

stock_df['Close'].mean()
stock_df['Close'].var()

# 平均： 4270740.898345154
# 分散： 18455369802797.1

（弱）定常性の定義は「平均0の周辺で一定の分散をもってランダムに振動を続けるデータ」でしたが、Close（終値）はかなりばらつきのあるデータと言えるでしょう。後ほど、定常性となりえなかった要因の成分を取り除き、再度確認することにしましょう。

2.単位根過程とは

単位根過程の定義としては、
$ y_t $が非定常過程、また差分系列　$ {\Delta}y_t = y_t – y_{t-1}　$が定常過程である時、$ y_t $は単位根過程である、です。※${\Delta}$は差分を表します。
データの原系列は非定常だけれど、差分を取ると（弱）定常になるようなデータとも言えます。

最も代表的な単位根過程として「ランダムウォーク」があります。
ランダム・ウォーク過程とは、$X_1, X_2, · · ·, X_n $ の系列が、$ Xt = Xt−1 + ut $　と表すことができます。

3.自己共分散

ある時系列 $ y_{1:t} $ に関して、各時点t における期待値を E[$y_t$]$ $= μ_t$ としたとき、
時点差$j$ における $yt$　と $y_$_{t-j}$の共分散を、時差jの自己共分散といいます。
つまり、同一の時系列データにおける異なる時点間での共分散を指します。
（※$Cov$ = Covariance（共分散）、$E$ = Expect（期待値））

Cov(y_t, y_{t-j}) = E[ (y_t - {\mu}_t)(y_{t-j} - {\mu}_{t-j})]

時系列データ分析とは

それでは時系列データ分析について解説していきましょう。
時系列データ分析とは、ある一定間隔の時間に対して観測されているデータの周期性や傾向から将来を予測する手法のことを言います。これにより近い将来を把握し、施策の検討や対策を打つことが可能となります。

時系列データ分析の目的

時系列データ分析には、大きく分けて下記の３つの目的があります。
1.データの可視化
2.情報の抽出
3.将来の予測
それぞれ具体的に見ていきましょう。

1.データの可視化

データの特徴を捉える為、「周期性」および「時間的な相関」を可視化していくことがあります。

①周期性

周期性（周期的変動）とは、時系列が⼀定の間隔で同じような変動を繰り返す成分、またはその状態を指します。
周期性を可視化する方法として「コレログラム」というものがあります。
コレログラムとは、ラグと自己相関を表したグラフです。
ラグとは、元データからどれほど時間をずらしているかを表す指標のことです。
つまり、元データxから時間をずらしたデータyとの相関係数を表すグラフであり、横軸にラグ、縦軸に自己相関をとります。
青で塗られた領域は、データが無相関としたときの95%信頼区間で、この領域の外にある場合この仮定は棄却され、相関があると判断できます。

# 表示サイズ調整
plt.rcParams["figure.figsize"] = [20,2]
plt.rcParams["font.size"] = 10

# コレログラム
def correlogram(cicle, lags, df, feature):
    sm.graphics.tsa.plot_acf(df.resample(cicle)[feature].sum(), title=feature, lags=lags)

# まずは原系列をインプットして確認
correlogram(cicle='Y', lags='7', df=time_series, feature='Close') # 年ベース
correlogram(cicle='M', lags=str(12*6), df=time_series, feature='Close') # 月ベース
correlogram(cicle='D', lags=str(30*12*6), df=time_series, feature='Close') # 日ベース

1番目が年ベースのグラフ、2番目が月ベースのグラフ、3番目が日ベースのグラフです。
これは原系列をコレログラムで表したグラフのため、本来このままでは正確な情報を得ることはできませんが、全体的に緩やかな変化が見て取れ、グラフの後方程、大きな周期変動がないように見られます。

②時間的な相関

時間的な相関を確認する方法として、「相関関数」を用いる方法があります。
これにより離れた時点の時系列との相関を⾒てることができます。
また、この相関にも種類があり、ここでは「自己相関」「偏自己相関」について解説します。

自己相関( ACF : Autocorrelation Function )

自己相関とは、過去の値が現在のデータにどのくらい影響しているかを表しています。
時系列のある時刻と$k$ だけ離れた時刻との相関をみると相関の観点から時系列の特徴を捉えることができるます。
ラグとは、ズラしたデータのステップ数のことです。
自己相関は $|{\rho}(j)|{\leq}1$を満たし、自分自身との自己相関係数は１となり、相関の強いデータほど自己相関係数は１に近づきます。
※ ${\gamma}(0)$ は分散0を表し、${\gamma}(j)$の強度により自己相関の強さを測ることができます。

{\rho}(j) = \frac{Cov(y_t, t_{t-j})}{\sqrt{Var(y_t)Var(y_{t-j})}} = \frac{{\gamma}(j)}{{\gamma}(0)}

偏自己相関( PACF : Partial Autocorrelation Function)

ある時点$j$ の偏自己相関とは、時点$t$ と時点$t-j$ の間に存在している$j-1$個のデータ $y_{t-(j-1)}$, $y_{t-(j-2)}$, ..., $y_{t-1}$ の影響を取り除いた後の時点$t$ と時点$t-j$ の間の相関のことです。
例えば、今日と2日前の関係を見たい際、間接的に1日前の影響が含まれる為、偏自己相関を用いることで、1日前の影響を除いて今日と2日前だけの関係を調べる事ができるようになります。
また、時点$t$ と時点$t-2$ の偏自己相関は、「時点$t$ のうち時点$t-1$ で説明できなかった情報$ y_t - {\alpha}y_{t-1} $」と「時点$t-2$ のうち、時点$t-1$ で説明できなかった情報$ y_{t-2}-{\beta}y_{t-1} $」の相関を表しています。通常 ${\alpha}, {\beta}$ は、$E[ (y_t - {\alpha}y_{t-1})^2 ]$ が最小となるように選ようにします。計算式で表すと以下の通りです。

\frac{Cov(y_t - {\alpha}y_{t-1}, y_{t-2} - {\beta}y_{t-1})}{\sqrt{Var(y_t - {\alpha}y_{t-1})Var(y_{t-2} - {\beta}y_{t-1})}}

2.情報（要因）の抽出

次に時系列データの情報抽出の方法として「季節調整」について解説していきます。
季節調整とは、何らかの原因で特定の周期で繰り返す成分を除去して本質的な現象を抽出する手法です。
具体的には、時系列データを「傾向変動」「季節変動」「不規則変動」に３つに分解していきます。これにより時系列データの傾向をはっきり⾒ることができます。

①傾向変動（トレンド）の推定

「移動平均」を取ることで傾向変動を推定する方法があります。
移動平均により、時系列データを平滑化（変動を小さく）し、なめらかな傾向変動（トレンド）を取り出すことができます。
具体的には、時系列データにおいて、ある一定区間ごとの平均値を区間をずらしながら求めていきます。

ある時点$t$ と時点$t-s$ の比の計式は下記の通りです。

Z_t = \frac{y_t}{y_{t-s}}　\\

ある時点$t$ の前後のデータを使って時点$t$ の値を調整し、時点$t$ の前後の値との差を小さくする流れを考えてみましょう。
具体的には、時系列全体に対して、ある時点$t$ のデータを、$t-k$ から $t+k$ までのデータの平均をとっていきます。

Z_t = \frac{1}{k + 1}(y_{t-k} + y_{t-(k+1)} + … + y_t … + y_{t+(k-1)} + y_{t+k})

移動平均を用いてグラフを作成すると、長期的な傾向を表す滑らかな曲線が得られます。
時系列データの変化のブレが細かく、解析への悪影響が懸念される際に使用されることがあります。

# インデックスに['Date']（時間情報）を設定
stock = stock_df.copy()
stock["Date"] = pd.to_datetime(stock['Date'])
stock = stock.set_index(['Date']).sort_index(ascending=True)

# 移動平均
stock['Close_move_mean'] = stock['Close'].rolling(100).mean().round(1)
plt.plot(stock['Close'], label='終値')
plt.plot(stock['Close_move_mean'], label='移動平均')
plt.legend(loc='lower right')
plt.show();

各時点の平均を取ることで、ゆるやかな変化を捉えられるようになったことがわかります。

②季節変動（シーズナル）の推定

季節変動を推定する方法として「季節階差」をとる方法があります。
季節階差とは、トレンド除去後の時系列データから季節変動（季節成分）を求める方法です。
階差とは今期と前期の値の差のことです。そして季節階差は周期pの階差を取ることで、時系列データからゆっくりした傾向変動を除去することができます。
また、季節成分を分離することにより，単なるトレンド推定よりトレンドの微妙な変化を捉えることができるようになります。
周期pの1周期をsとした場合、計算式は下記のように表すことができます。

{\Delta}_sy_t = y_t - y_{t-s}

③不規則変動（ノイズ）の推定

傾向変動と季節変動を除去することで、不規則変動を推定することができます。

3.将来の予測

時系列データ分析により「将来の予測」が可能となります。
その為に様々なモデルが考案されています。代表的なものとして「MAモデル」「ARモデル」「ARMAモデル」「ARIMAモデル」「SARIMA」「状態空間モデル」などがあります。時系列データの予測モデルについては、また別記事で紹介させていただきます。

定常性の確認

時系列データが定常過程を単位根過程かを確認する手法に、ADF（Augmented Dickey-Fuller）検定というものがあります。
多くの時系列モデルでは定常過程を前提としているため、時系列に対してまず最初に単位根を確認する事が多いです。
具体的には、帰無仮説:単位根過程, 対立仮説 : 定常過程、とし、P値（有意水準）が0.05（5%）以下なら帰無仮説が棄却され、定常過程とみなしていきます。実際に実装と結果を見ていきましょう。

def adf_test(series):
    adf_df = pd.DataFrame(
        [
            stattools.adfuller(series)[1]
        ],
        columns=['P値']
    )
    adf_df['P値'] = adf_df['P値'].round(decimals=3).astype(str)
    print(adf_df)

# 原系列のp値を算出
adf_test(series=time_series['Close'])

# 結果：p値 = 0.815

p値が81.5%となり、原系列は非定常性であると判断できます。
定常性を持たない場合、そのままの状態ではデータの特徴や関係性を正しく捉えられないため、定常性を持つようなデータに変換が必要になります。データの変換方法にもいくつかあるため、以下で解説していきます。

時系列データの前処理

検定により時系列データが単位根過程と判断された場合、データの関係性を正しく解析ができない為、定常性を持たせるための変換を行う必要があります。
ここでは原系列に対する、各種の変換方法を解説していきます。
原系列とは、変換前のもともとの時系列データのことです。

1.差分変換

差分変換とは、ある時点$t$ のときの値とその直前の時点$t-1$ のときの値の差分 ${\delta}t = y_t - y_{t-1}$ を計算し、その差を新しいデータとする手法です。
つまり、1時点離れたデータとの差を取る手法です。その結果を差分系列、または階差系列と言います。
例えば、今日の株価と昨日の株価の差を考えたい時などに使用します。
上昇あるいは下降のトレンドを持つ時系列データ$（y_t = a_t + b）$ のとき、差分を計算することで、トレンドを取り除くことができます。計算式は以下の通りです。

Z_t = y_t - y_{t-1} = (a_t + b) - (a(t - 1) + b) = a

# Dataをインデックスに設定
time_series = stock_df.copy()
time_series["Date"] = pd.to_datetime(stock_df['Date'])
time_series = time_series.set_index(['Date']).sort_index(ascending=True)

# 原系列のプロット
plt.plot(time_series['Close'])

# 差分変換後のプロット
time_proc = time_series.copy()
Passengers_diff_c = time_proc["Close"].diff(periods = 1)
Passengers_diff_c.plot(figsize=[20,5]);

Passengers_diff_c.mean()
Passengers_diff_c.var()

# 平均：0.14051448905972796
# 分散：19.280544728402738

上のグラフが原系列（比較用）、下のグラフが差分変換後のグラフです。
差分を取ることで、上昇傾向が取り除かれ、平均が0の正常過程の状態に近づけることができるようになりました。y軸のスケールが変換前後で変わっているのもポイントです。

2.対数変換

対数変換とは、変動が著しく大きな時系列データに対し、対数値変換を施すことで原系列の傾向を保持したまま、値を小さく変換する手法です。
変動の分散を一様にし、複雑な時系列でも変数変換で分析が簡単になることがあります。
対数変換された時系列データは対数系列と呼ばます。数式では${\log}(y_t)$ と表します。
例えば、ビットコインの価格などの極端な動きをするデータに対数系列を用います。

Passengers_log_c = np.log(time_proc["Close"])
Passengers_log_c.plot(figsize=[20,5]);

Passengers_log_c.mean()
Passengers_log_c.var()

# 平均：4.48285833956542
# 分散：0.9894592908040513

上のグラフが原系列（比較用）、下のグラフが対数変換後のグラフです。
際立っていた上下運動が平滑かされ、原系列のトレンドに比例した値のばらつきが解消されていることがわかります。また分散を見るとかなり0に近づいており、データのばらつきが減っていることが数値から見てもわかります。

3.対数差分変換

対数差分変換とは、対数変換と差分変換の両方の変換を施す手法です。
時系列データが変化率や成長率などである場合、対数差分 ${\Delta}{\log}(y_t) = {\log}(y_t / y_{t-1}) $ を計算して解析することがあります。

Passengers_log_diff_c = Passengers_log_c.diff(periods = 1)
Passengers_log_diff_c.plot(figsize=[20,5]);

Passengers_log_diff_c.mean()
Passengers_log_diff_c.var()

# 平均：0.0014151238107781069
# 分散：0.001069792377059396

１番目が原系列（比較用）、２番目が対数変換後の時系列（比較用）、３番番目が対数差分変換後の時系列です。
全体を通して平均と分散はのばらつきは解消し、年ごとの周期性がゆるやかになっていることがわかります。
また、平均と分散がかなり0に近づいており、データのばらつきが解消され、定常過程に近づいたことが数値をみてもわかります。

再度ADF検定

変換後のデータで再度ADF検定を実施した結果を見てみましょう。

# 対数変換
Passengers_log_c = np.log(time_proc["Close"])
Passengers_log_c.plot(figsize=[20,5]);

# 季節差分変換
Passengers_log_sdiff_c = Passengers_log_c.diff(periods = 12) # 季節：periods = 12
Passengers_log_sdiff_c.plot(figsize=[20,5]);

# ラグ1の差分を取る
Passengers_sdiff2_c = Passengers_log_sdiff_c.diff()
Passengers_sdiff2_c.plot(figsize=[20,5]);

# 変換済みデータの欠損の削除
Passengers_sdiff2_c2 = Passengers_sdiff2_c.dropna()

# 成分分解（季節調整）
res = sm.tsa.seasonal_decompose(Passengers_sdiff2_c2, period=100)
fig = res.plot()
plt.show()

# 原系列のp値を算出
adf_test(series=Passengers_sdiff2_c2)

# 結果：p値 = 0.534

変換前のp値は「0.815」でしたが、変換後は「0.534」と下り、定常過程に近づいたことがわかります。
ただし、有意水準（0.05以下）には程遠く、まだこのままでは解析が行えません。
その為、定常性となりえなかった要因を、更に取り除いていく必要があります。
そのために、まずはなぜ定常性がないのかを考察する必要があります。
また、今回は特徴量が非常にシンプルなデータセットを用いていますが、本来はそれ以外の特徴量を追加・生成するなどの対応も必要と思われます。
もしくは「定常・非定常」に左右されない分析手法を採用する、などの代替案を検討してもよいでしょう。

時系列データにおける特徴量作成

今回は参考程度ですが、時系列データでは下記の特徴量（生成）が有効となるケースがあります。
月、曜日、週末フラグ、祝日、休日、コロナ患者数、などなど。

まとめ

今回は時系列データ分析について、基本的な事項を解説してきました。
実際に分析を進めてみて、実務ではリリースのタイミングなどで確率分布が途中で変化することが予想される為、今回解説した手法の実践投入は自然科学的なデータセットに限られると思いました。
実際にはパターンマッチング的な考え方(LSTMなどDLモデル)や様々なモデルを組み合わせられる状態空間モデルなどを利用することが多いと思います。
また、コレログラムは途中少し解説しただけとなってしまいましたが、実際はデータ変換時に都度確認に用いるものかと思います。あと、今回用いたデータセットは情報量（≒特徴量）が少なく（考察と取り得る対処の幅に限界があり）、もう少しバリエーションのあるデータで分析や特徴量生成を行った方が良いと思いました。
別途、時系列モデリングとそのアルゴリズムについても解説したいと思います。

最後に

他の記事はこちらでまとめています。是非ご参照ください。

参考文献

東京大学 数理・情報教育研究センター - 時系列データ解析
東京⼤学 - 数理⼿法VII（時系列解析）
学習院大学 - 時系列解析入門
大阪大学 - 時系列分析

解析結果

実装結果：GitHub/chronological_order_analysis.ipynb
データセット：Tesla Stock Price - kaggle

参考資料