モモノキ＆ナノネと一緒にPythonで時系列分析を覚えよう（7）
時系列分析の実践練習（自己相関）
今回は時系列分析で重要な自己相関の分析方法について練習するよ。
自己相関分析って？
簡単に言うと、自己相関分析は過去の値が現在のデータにどれくらい影響しているか？その関係性を調べることだよ。
普通の相関分析だと異なるデータ同士の関係性を調べるよね、例えばX、Yの値の相関みたいに。自己相関の場合はキホン的に一種類のデータを使い、自分自身のデータをズラして関係性を分析するんだよ。
自分自身のデータを使って関係性を調べるから、相関の前に自己って付くんだ。
日次データだったら、一つズラして自己相関を確認すれば、一日前の値が今日にどれくらい影響しているかが分かるよ。ズラしたデータのステップ数をラグ（lag)と呼ぶから覚えておいてね。
タイムラグのラグ（遅れ）と一緒？
うん、たぶん仲間だね。
相関の強さは相関係数みたいな値で出るの？
うん、自己相関係数といって値の範囲は-1～1、正負の関係や絶対値が1に近いほど相関が強いなども同じだよ。
自己相関係数はどうやって求めるの？
普通の相関係数を求めるみたいに、自己相関係数も（自己）共分散をもとにして計算するだよ。式は下の通りでラグはkで表記しているよ。自己相関係数は自動で簡単に計算できるから心配しなくても大丈夫だよ。
自動計算でよかった。コードの書き方は？
ナノネ、前回使ったコードを呼び出してみて。データ準備まででOKだよ。
# 飛行機乗客数の時系列データ（前回使用コード抜粋）
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use('ggplot') # グラフのデザイン変更（お好みで利用）
plt.xkcd() # グラフのデザイン変更（お好みで利用、日本語は未対応）
import statsmodels.api as sm # version 0.8.0以上
# CSVファイル読み込み（事前にダウンロードしたCSVファイルを利用）
df = pd.read_csv('AirPassengers.csv')
# Pandas.Seriesにデータを格納（データに乗客数、インデックスは日付）
passengers = pd.Series(df['#Passengers'], dtype='float') # ①
passengers.index = pd.to_datetime(df['Month']) # ②
飛行機乗客数のデータ準備OK。
自己相関係数（ACF: Autocorrelation Function）は、StatsModelsのtsa.stattools.acf()で求めらるよ。第一引数はデータ、オプション引数nlagsでラグ数を指定してね（デフォルト40）。
# 自己相関係数の出力
passengers_acf = sm.tsa.stattools.acf(passengers, nlags=40) #ラグ=40(デフォルト値)、省略可能
passengers_acf
41個の自己相関係数が表示された。一番最初の相関係数は1.0？
ラグ40を指定したから、ラグ0～40まで計41個の自己相関係数が出力されたんだよ。一番最初はラグ0の結果で常に1.0だよ。ラグ0はデータをズラさない、同じ値同士の相関だから当然だね。
自己相関係数を自力計算した試しコードを載せておくね。上の式をそのまま再現してみたけど、結果は一致するみたい。
# 自己相関係数を試しに計算してみた（Numpy利用）
# import numpy as np インポート済み
LAG = 40 # 計算ラグ数
rk = np.empty(LAG+1) # 自己相関係数の計算結果を保持用
y = np.array(passengers) # 乗客数データのndarray作成（計算用）
y_mean = np.mean(y) # 乗客数の平均値
# ラグ0の自己相関係数：1.0
rk[0] = np.sum((y - y_mean)**2) / np.sum((y - y_mean)**2)
# ラグ1〜40の自己相関係数：-1.0〜1.0
for k in np.arange(1, LAG+1):
rk[k] = np.sum((y[k:] - y_mean)*(y[:-k] - y_mean)) / np.sum((y - y_mean)**2)
print(rk) # -> acf()結果と同じ
自己相関係数の結果も、いい感じにグラフ化できる？
StatsModelsのgraphics.tsa.plot_acf()で、計算とグラフ作成まで自動でできるよ。第一引数にデータ、ラグ数はオプション引数lagsで指定できるよ。
# 自己相関(ACF)のグラフ自動作成
fig = plt.figure(figsize=(12, 4))
ax1 = fig.add_subplot(111)
sm.graphics.tsa.plot_acf(passengers, lags=40, ax=ax1) #飛行機乗客数データ、ラグ40、グラフaxes
plt.show() # インライン表示で同じグラフが2個出でる？ので回避
簡単にできた。横軸がラグ数で、縦軸が自己相関係数。波みたいな形でラグ数が増えると自己相関係数の値が小さくなってる。薄い赤で塗り潰された領域はなに？
95%信頼区間だよ。信頼区間の領域を超えてプロットされているデータは、統計的に有意差がある値とみれるよ。
ラグと自己相関でプロットしたグラフを『コレログラム』と言うから覚えておいてね。
plt.plot(passengers) # 飛行機乗客数の元データのグラフ（再掲）
飛行機乗客数のデータは右肩上がりのトレンドで季節性もあったから、いかにも自己相関ありがそうな結果がでたけど。もし、自己相関が全くない場合はどんな結果になるの？
それなら、疑似的にランダムなデータを作って自己相関係数を確認してみよう。ランダムはデータはNumpyで簡単につくれるよ。np.random.normal(平均, 分散, 出力数)とすればOK。
# import numpy as np # インポート済み
np.random.seed(0) # 乱数の種固定
random_data = np.random.normal(0, 1, 144) # 平均0、分散1、144個（乗客数データと同じ）のランダムデータ作成
plt.plot(random_data) # 作成したランダムデータをそのままプロット
# acf＆グラフの自動作成は前のコードをコピペ
fig = plt.figure(figsize=(8, 4))
ax1 = fig.add_subplot(111)
sm.graphics.tsa.plot_acf(random_data, lags=40, ax=ax1) # ランダムデータを指定
plt.show()
なるほど。意味のないラグ0は除いて、自己相関係数の値はどれも小さい。95%信頼区間にほぼ全部のデータが入っている。
ランダムデータは過去の値と全く関係ないから、分析で予測とかできないけどね。自己相関は過去の影響を調べる重要な分析で、予測モデルをつくるときにも必須情報だよ。
自己相関の分析は重要ね。
あと、自己相関の仲間で偏自己相関（PACF:Pertial Autocorrelation Function）も覚えておいてね。偏自己相関はある時点同士だけの関係性だよ。よくある例で、今日と二日前の関係には間接的に一日前の影響が含まれるよね。偏自己相関を使うと、一日前の影響を除いて今日と二日前だけの関係を調べられるよ。
偏自己相関もStatsModelsのgraphics.tsa.plot_pacf()でグラフが作成できるよ。自己相関（ACF）と偏自己相関（PACF）のグラフを縦に並べて描いてみよう。
fig = plt.figure(figsize=(8, 8))
# 自己相関(ACF)のグラフ
ax1 = fig.add_subplot(211)
sm.graphics.tsa.plot_acf(passengers, lags=40, ax=ax1) #ACF計算とグラフ自動作成
# 偏自己相関(PACF)のグラフ
ax2 = fig.add_subplot(212)
sm.graphics.tsa.plot_pacf(passengers, lags=40, ax=ax2) #PACF計算とグラフ自動作成
plt.tight_layout() # グラフ間スキマ調整
データの呼び名も覚えておいてね。時系列分析をするとき、元のデータを一旦変換してから利用するケースが多いから。
なにも変換処理していない元データは『原型列』と呼ぶよ。時間順に採取した値をそのまま並べただけのデータ。さっき使った飛行機乗客数のデータも原型列だよ。
対数変換した場合は『対数系列』と呼ぶよ。データに大きな値のバラつきがあるときは、対数変換すると数値が扱いやすくなるよ。
階差に変換した場合は『階差系列』と呼ぶよ。時系列で隣り合うデータとの差をとったものだよ（今日の値－昨日の値みたに）。階差をとる回数で1次の階差、2次の階差、n次の階差とも表現されるよ（通常は1次階差が多く使われる）。階差をとるとトレンド成分の除去につながるよ。あと、見せかけの回帰というのを回避するのに有効みたい。
季節変動を取り除いたデータは『季節調整済み系列』と呼ぶよ。他にも平方根をとったり、変換の組み合わせで色々な呼び名があるよ。
名前がいっぱい出てきて覚えきれないなぁ。。。
時系列分析の重要キーワードは、まだまだあって『定常化』もその一つだよ。時系列分析は定常化を前提にしたデータが求められるので、定常・非定常状態なども調べてながら追々覚えて行こうね。
さっき挙がった階差系列も定常化する手法の一つだよ。今回は階差系列でコレログラムを作成して、原型列のときと違いをみて終わりにしよう。
OK、階差系列でコレログラムを作るんだね。階差系列はどうやって作るの？
Pandasのsift()でデータを一つ分シフトできるから、差分をとれば階差系列が作れるよ。データ － データ.shift() みたいな感じでOKだよ。shit()は引数にシフト量も指定できるよ。もし7個シフトしたい場合はshit(7)で。
passengers_diff = passengers - passengers.shift() # 階差系列データの作成
passengers.head() # 原型列（元データ参考）
print(len(passengers_diff))
passengers_diff.head() # 階差系列(1次階差)
階差系列を作れた。でも1行目がNaNデータになっちゃうけど。
今回のNaNデータはPandasのdropna()で捨てちゃおう。
passengers_diff = passengers_diff.dropna() # Nanデータ削除
print(len(passengers_diff))
passengers_diff.head()
次は原型列と階差系列のグラフを横に並べて比べてみよう。
fig = plt.figure(figsize=(12, 4))
# 原型列のグラフ
ax1 = fig.add_subplot(121)
plt.plot(passengers)
# 階差系列のグラフ
ax2 = fig.add_subplot(122)
plt.plot(passengers_diff)
plt.tight_layout() # グラフ間スキマ調整
左側のグラフが原型列で、右側が階差系列でプロットした。階差系列のほうは、右肩上がりの傾向が消えたみたい。周期的なギザギザは残ったままで、振れ幅がだんだん大きくなるのは原型列と似ている。
うん、階差をとったことで原型列に見られた右肩上がりのトレンド成分が取り除かれたんだよ。
ナノネ、最後に階差系列でコレログラムを描いてみて。
# 階差系列データでコレログラム作成
fig = plt.figure(figsize=(8, 8))
# 自己相関(ACF)のグラフ
ax1 = fig.add_subplot(211)
sm.graphics.tsa.plot_acf(passengers_diff, lags=40, ax=ax1) #ACF計算とグラフ自動作成
# 偏自己相関(PACF)のグラフ
ax2 = fig.add_subplot(212)
sm.graphics.tsa.plot_pacf(passengers_diff, lags=40, ax=ax2) #PACF計算とグラフ自動作成
plt.tight_layout() # グラフ間スキマ調整
データ名だけ変えてコピペでグラフが描けた。コレログラムの結果も、かなり変わったみたい。
今使っている飛行機乗客数は月次データで、夏場の乗客数が多かったよね。1年（12カ月）周期の変動傾向があるから、ラグ12のところで強い相関が確認できるんだよ。ACFグラフだと、ラグ12、24、36で95%信頼区間を大きく飛び出した強い相関。周期も明確に分かるね。
グラフで周辺と比べて突出したところを『スパイク』とも呼ぶから覚えておいてね。
次回以降は、モデルの扱い方を練習していくよ。
モデルね、ラジャー。
あけましておめでとうございます。とてもわかりやすく、初心者ですが楽しくフォローしています！
返信削除偏自己相関(PACF)のグラフのところなのですが、ローカル環境でプロットしてみるとlag=39のところで３２くらいに突出していて、この記事のグラフと異なります。なぜでしょうか？