自己相関関数と相互相関関数

自己相関関数と相互相関関数#

データの特徴を表す指標として、平均や分散、最大値最小値などがあることは「基本的な要約統計量の計算」で説明しました。
　時系列データに特有の指標として、自己相関関数・相互相関関数があげられます。これらは、データ間の時間方向に対する類似性を評価可能で、モデル選択や特徴量設計に有用です。

自己相関関数#

1つの時系列データを対象として、少しずつ時間をずらして相関を求めます。自分自身の過去のデータとどのくらい類似してそうかを調べられます。自分自身のラグだけから相関を求めるため自己相関と呼ばれます。
ある時系列データ\(x(t)\) に対して、自己相関は次式で求めることができます。

\[\frac{E[(x(t)-\mu)(x(t-k)-\mu)]}{\sigma^2}\]

ここで、\(k\) はラグで、\(x(t-k)\) は\(k\)だけ時間をずらした\(x(t)\)、 \(\mu, \sigma^2\) はそれぞれ \(x(t)\) の平均と分散です。分子は自己共分散で、分母は自己共分散を\([-1, 1]\)に正規化する役割です。なお、ここで\(x(k)\)は定常であるとします（平均と分散が時間依存ではない）。定常性の詳細については「定常性の確認」を参照してください。

statsmodels には、自動で計算してくれるモジュールがあるのでこれを使います。

# Warningsの消去
import warnings

warnings.filterwarnings("ignore")

# モジュールの読み込み
from datetime import datetime, timedelta
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm

%matplotlib inline
# matplotlibでは日本語にデフォルトで対応していないので、日本語対応フォントを読み込む
plt.rcParams["font.family"] = "BIZ UDGothic"  # Windows
# plt.rcParams['font.family'] = 'Hiragino Maru Gothic Pro' # Mac

# その他使えるフォント一覧は以下のコードで取得できます。
# import matplotlib.font_manager
# [f.name for f in matplotlib.font_manager.fontManager.ttflist]

# データの読み込み
fname = "../data/2020413.csv"
df = pd.read_csv(
    fname, index_col=0, parse_dates=True, na_values="None", encoding="shift-jis"
)

# 例として1つ抜き出してみる
col_name = "気化器温度_PV"

acf = sm.tsa.stattools.acf(df[col_name])
print(acf)

[1.         0.62974482 0.62350951 0.60453868 0.61395385 0.61202015
 0.60929824 0.61134841 0.61230394 0.61473471 0.60889061 0.60857538
 0.61025949 0.61053135 0.60902582 0.60661423 0.60890996 0.60705146
 0.60689684 0.60717655 0.60692475 0.60628056 0.60740425 0.60671663
 0.60975835 0.60938884 0.60623177 0.6064455  0.60623949 0.60781058
 0.60603641 0.60813034 0.6048722  0.60773175 0.60724384 0.60875073
 0.60552821 0.60416543 0.6068039  0.60355536 0.60653161]

ラグがゼロのときは当然のことながら、

\[\frac{E[(x(t)-\mu)(x(t-0)-\mu)]}{\sigma^2}=\frac{\sigma^2}{\sigma^2}=1\]

になります。

numpyで実装することもできます。

nlags = 40
x = df[col_name].values
x_mean = np.mean(x)

acf_np = [np.sum((x - x_mean) * (x - x_mean) / np.sum((x - x_mean) ** 2))]
for i in range(1, nlags + 1):
    _acf = np.sum((x[i:] - x_mean) * (x[:-i] - x_mean)) / np.sum((x - x_mean) ** 2)
    acf_np.append(_acf)

acf_np = np.array(acf_np)
print(acf_np)  # 一致した結果が出てくる

[1.         0.62974482 0.62350951 0.60453868 0.61395385 0.61202015
 0.60929824 0.61134841 0.61230394 0.61473471 0.60889061 0.60857538
 0.61025949 0.61053135 0.60902582 0.60661423 0.60890996 0.60705146
 0.60689684 0.60717655 0.60692475 0.60628056 0.60740425 0.60671663
 0.60975835 0.60938884 0.60623177 0.6064455  0.60623949 0.60781058
 0.60603641 0.60813034 0.6048722  0.60773175 0.60724384 0.60875073
 0.60552821 0.60416543 0.6068039  0.60355536 0.60653161]

statsmodels は自動でグラフを作成してくれるので、それを使って可視化してみます。このデータは、1秒ごとに記録されており長周期の挙動をとらえることが難しいので、間隔が60秒ごとになるように間引きをします。本来はしっかりとしたデシメーションを行うべきですが、今回は雑に間引きます。

fig = plt.figure(figsize=(8, 5))
ax1 = fig.add_subplot(111)
sm.graphics.tsa.plot_acf(df.loc[::60, col_name], lags=40, ax=ax1)  # グラフを自動作成
ax1.set_xlabel("ラグ数")
ax1.set_ylabel("自己相関")
plt.show()
plt.close("all")
../_images/914dbec75f99a8160603b29f52a55ae8cb27005a96cadf8b3f265bbad44c2bfc.svg

グラフの横軸はラグ数で、縦軸は自己相関の値です。今回、データ点は60秒ごとになるように間引かれているので、1ラグあたり60秒の遅れです。また、水色の領域は95%信頼区間で、これの範囲外にあるものは統計的に有意差がある値（≒意味がある値）とされています。
　このグラフをみると、周期性があるような様子が見られます。ラグが32のところに次のプラスのピークが現れてるので、周期が32分ぐらいのようです。

実際に、生データを時間ずらして並べてみます。

fig = plt.figure(figsize=(8, 5))
ax = fig.add_subplot(111)
ax.plot(df.index, df[col_name], label="ラグなし")
ax.plot(df.index - timedelta(minutes=32), df[col_name], label="ラグあり")
ax.legend()
ax.set_xlim([datetime(2020, 4, 13, 0)], datetime(2020, 4, 13, 1))
plt.show()
../_images/8d542583e2bf086623e751b8c37ffec2b8ea81681e40e35d6f89cf436add9f5a.svg

確かに、似たような傾向を示しています。

もう一つ例を見てみましょう

col_name = "蒸留塔第5トレイ温度_PV"
fig = plt.figure(figsize=(8, 10))
fig.tight_layout()
ax1 = fig.add_subplot(211)
ax1.plot(df.index[::60], df.loc[::60, col_name])
ax1.set_title(col_name)
ax2 = fig.add_subplot(212)
sm.graphics.tsa.plot_acf(df.loc[::60, col_name], lags=60, ax=ax2)

plt.show()
../_images/ea9cd8a4bb55cfeeefd52b99785e1ca4fde3861e592dc259424ade905fec9356.svg

こちらの場合だと、自己相関がほとんどありません。このことがわかると時系列モデルの際に有用です。例えば、時系列モデリングでよく用いられる自己回帰（AR）モデルは過去の自身の値から、未来を予測するようなモデリング手法です。しかし、上記のデータだと自分の過去の値には関係がほとんどないので、ARモデルを用いた予測は難しいことが推察されます。

自己相関がわかるといいこと#

  • 過去の値がどれだけ現在の値に影響しそうか？がわかる

  • 周期性がわかる

    • これらの情報から、どのようなモデルを設計すべきかがわかる

  • 逆に相関がない（ホワイトノイズみたい）場合だと、単変数から回帰することは難しい

    • 解析対象から外したほうがいいかもしれない。あるいは別のモデルを選択すべきかもしれない。という指標になる。

偏自己相関関数#

偏自己相関は、ある時刻間だけの相関を計算します。例えば、\(x(t)\)と\(x(t-2)\)は、\(x(t-1)\)を通じて相関している可能性があるので、この\(x(t-1)\)の影響を取り除いた上で\(x(t)\)と\(x(t-2)\)の関係を評価します。
statsmodels の graphics.tsa.plot_pacf() でグラフ作成できます。自己相関関数の項目で解析した 気化器温度_PV の偏自己相関を例示します。

col_name = "気化器温度_PV"

fig = plt.figure(figsize=(8, 5))
ax1 = fig.add_subplot(111)
sm.graphics.tsa.plot_pacf(df[col_name], lags=40, ax=ax1)  # グラフを自動作成
ax1.set_xlabel("ラグ数")
ax1.set_ylabel("偏自己相関")
plt.show()
plt.close("all")
../_images/5a14dc96ffface8ca5d9de3dcdd71b76431ba1e5046760a1b4fba9e9a83e1aaf.svg

ラグが大きくなるにつれて、指数的に0に収束していく様子が見られます。このとき、MA（移動平均）項が存在していると考えられます。MA項についての詳細は、「モデリング」の項で解説予定です。

相互相関関数#

前述の自己相関関数と同様のアイデアで、2種類の時系列データ間の相関を計るのが相互相関関数です。自己相関関数の \(x(t)\) の片方を別の時系列データ \(y(t)\) で置き換えて次式でかけます。

\[\frac{E[(x(t)-\mu_x)(y(t-k)-\mu_y)]}{\sqrt{\sigma_x^2}\sqrt{\sigma_y^2}}\]

目的変数 リサイクルガス組成（O2）_PV と、それに対するいくつかの説明変数との相互相関を計算してみます。statsmodels に相互相関の実装がされているのでそれを使います。なぜか、自動でグラフを作ってくれるモジュールはありません…
　上式からわかるように、\(x\)と\(y\)の順番で相互相関の値は異なります（ccf(x, y) ≠ ccf(y, x)）。

col_obj = "リサイクルガス組成（O2）_PV"
cols = [
    "熱交換器入口流量（製品ガス）_PV",
    "反応器流入組成（O2）_PV",
    "リサイクルガス流量_PV",
]
plot_size = 60  # プロットするラグの数

y = df.loc[::60, col_obj]  # .values

for col in cols:
    x = df.loc[::60, col]
    ccf_xy = sm.tsa.ccf(x, y)[1 : plot_size + 1]
    ccf_yx = sm.tsa.ccf(y, x)[
        :plot_size
    ]  # (x, y)を基準にしているので、(y, x)はマイナスのラグ
    ccf = np.concatenate((ccf_yx[::-1], ccf_xy))

    x_axis = np.arange(-plot_size, plot_size)
    fig = plt.figure(figsize=(6, 3), dpi=120)
    ax = fig.add_subplot(111)
    ax.stem(x_axis, ccf)
    ax.set_xlim([-30, 30])
    ax.set_ylim([-1, 1])
    ax.set_title("{} vs {}".format(col, col_obj))
    ax.set_xlabel("ラグ数")
    ax.set_ylabel("相互相関")
    plt.show()
../_images/5abbec755609b32b9684912ee31cab2bc70a76db784e2db631337fe288f94ed7.svg../_images/22741607277a53bbf89857c3c9bfaf7a96ba09260bc6fc19e9404ff5d38e82a3.svg../_images/e8360cdc9fdb74438c7c901998731fb49d046b576237e108c504acaa080a0909.svg

こうしてみると 熱交換器入口流量（製品ガス）_PV の相互相関はかなり小さく、目的変数にあまり影響を及ぼしていなさそうな気がします。パラメータを少なくしたいという要請があるのであれば、説明変数から省いてしまってもいいかもしれません。
逆に、 反応器流入組成（O2）_PV の相互相関がかなり高く、予測に大きく影響を及ぼしていそうです。むしろ影響が強すぎて、 反応器流入組成（O2）_PV を何倍かしたら 目的変数の値になってしまい、予測の意味をなさない恐れもあります。
リサイクルガス流量_PV は（わずかながら）振動しており、何かしら時間遅れを持つような系であることが推察されます。

相互相関関数をみると変数間の類似度をみることができ、さまざまな考察が可能です。また、説明変数間の類似度も同様に計算可能です。しかし、これをすべての変数について行うことは大きな労力がかかるため、ヒートマップで見ることもあります。pandasのcorr()メソッドで計算できます。ここで計算する相互相関はラグ0です。

# ヒートマップで表してみる。
# ここではラグなし（k=0）の時をしめす
def plot_corr_heatmap(df, figsize=(16, 16)):
    """
    相互相関のヒートマップをプロット
    Parameters
    ----------
    df: pandas Dataframe
    figsize: int
        figure size
    Returns
    -------

    """

    df_corr = df.corr()  # 相互相関を計算
    fig, ax = plt.subplots(figsize=figsize)
    fig.tight_layout()
    heatmap = ax.pcolor(df_corr, cmap=plt.cm.jet)

    ax.set_xticks(np.arange(df_corr.shape[0]) + 0.5, minor=False)
    ax.set_yticks(np.arange(df_corr.shape[1]) + 0.5, minor=False)

    ax.invert_yaxis()
    ax.xaxis.tick_top()

    ax.set_xticklabels(df_corr.columns.values, minor=False, rotation=90)
    ax.set_yticklabels(df_corr.index.values, minor=False)
    pp = fig.colorbar(heatmap, aspect=20)
    plt.show()


# スペースがないので一部だけ表示
plot_corr_heatmap(df.iloc[:60:, :40])
../_images/5675515577706e875341fc04637c6b0e03e729a9da8665d05e2c556340f2617a.svg

説明変数間でも強い相関を持っているものが多数存在していることがわかります。説明変数間の相関が強いと、線形なモデル（線形回帰など）で分析する際に悪い影響を及ぼすことが知られています（多重共線性）。多重共線性を回避し、適切な特徴量選択をするためにも、相互相関をみることが重要になります。

小括#

時系列データの特徴を表す統計量の一つである、自己相関関数・相互相関関数について概説しました。これらの統計量から、データ間や過去との類似度を知ることができ、時系列モデルの選択や特徴量選択のために有用です。

目次