本記事の目的

本記事の目的は、時系列性を伴うデータの分析を行う際に現れる「見せかけの回帰」というとても恐ろしい事象についてご紹介することです。Google Colaboratory を使って実行しながら、その恐ろしさを体感していただければと思っています。

ソースコードは GitHub で公開しておりますので、必要であればご参照ください。

見せかけの回帰とは？

言葉による説明

見せかけの回帰とは本来全く関連性のない２変数の間に関連性があるかのような回帰結果が得られることを意味します。

おそらく言葉だけではイメージがつきにくいと思いますし、まだ懐疑的な方もいらっしゃると思いますので実験をしてみましょう！

実験による確認

問題のないケース

実験に関しては Google Colaboratory で行います。まずは１つ目のセルで必要なパッケージをインポートしましょう。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import statsmodels.api as sm

次に２つ目のセルでシミュレーションを実施します。といってもやることはとてもシンプルで、標準正規分布に従う乱数1000個作成します。これを２セット作ることでそれぞれを説明変数( $u$ )、目的変数( $v$ )とし回帰を実行するだけです。

もちろんこれだけでは $u$ や $v$ に時系列性が一切存在しないので、ここで我々が確認したいケースではないのですが一旦「回帰において有意性が確認されない」ことを確認しておきましょう。

コードは以下のようになります。

# サンプル数とシード値を指定
num_of_sample = 1000
np.random.seed(8888)

# 標準正規分布（もちろん独立）に従って u, v の２系列を作成
u = np.random.randn(num_of_sample)
v = np.random.randn(num_of_sample)

# 回帰分析の定数項を推定するため準備
u_add_const = sm.add_constant(u)

# 線形回帰の推定と推定結果の出力
linear_model = sm.OLS(endog=v, exog=u_add_const)
fitted_linear_model = linear_model.fit()
print(fitted_linear_model.summary())

私の手元では結果は以下のように出力されました。

ということで、p値は 96.7% となり、特段問題ない推定結果が得られています。

問題の発生するケース

次に、乱数に時系列性が存在するケースを考えてみましょう。もっとも典型的なケースはランダムウォークと呼ばれており、以下の数式で表現できます。

\begin{aligned} x_{t} &= x_{t-1} + u_{t},\ t \in {1,2,...,999} \\ y_{t} &= y_{t-1} + v_{t},\ t \in {1,2,...,999} \\ x_0 &= u_{0} \\ y_0 &= v_{0} \end{aligned}

ただし、 $t \in {0,1,...,999}$ はデータが何番目のデータ化を表し、 $u_{x,t}, u_{y,t}$ はそれぞれ互いに独立な標準正規分布に従う変数です。つまり今日のデータの値は、昨日のデータの値に乱数がのっかって出てくるということを言っているわけですね。一方最初の日については参照できる前の日がありませんから、そのまま乱数の値を採用しています。

さて、手元にちょうど1000個ずつ $u$ と $v$ のサンプルがありますから、この乱数を利用して、 $x$ , $y$ による回帰分析を行ってみましょう。元々の実験で推定の有意性が確認されていないのですから、有意性が確認される訳がないですよね？

３つ目のセルに以下のように打ち込み、このコードを実行してみましょう。

# ランダムウォークに従う系列に変更(これでランダムウォークになる理由は以下で説明)
x = np.cumsum(u)
y = np.cumsum(v)

# 回帰分析の定数項を推定するための準備
x_add_const = sm.add_constant(x)

# 線形回帰の推定と推定結果の出力
rw_linear_model = sm.OLS(endog=y, exog=x_add_const)
rw_fitted_linear_model = rw_linear_model.fit()

私の手元では結果は以下のようになりました。

大変恐ろしいことに、それまで無関係であった系列をランダムウォークにした結果、有意な推定係数が得られてしまいました。

なぜ np.cumsum でランダムウォークの系列が得られるのか

$x_2$ について数式を少しいじってみましょう。

\begin{aligned} x_{2} &= \underbrace{x_1}_{=x_0 + u_{1}} + u_{2} \\ &= \underbrace{x_0}_{=u_0} + u_{1} + u_{2} \\ &= u_0 + u_1 + u_2 \end{aligned}

となりますね。つまり、初期値から乱数を足し上げていくとランダムウォークに従う系列 (ここでは $x$ ) を得ることができます。これは $t=3,4,...$ と大きくしていっても同じです。したがって最初の値から任意の順番の値までを足し上げる機能である np.cumsum() を利用することでランダムウォークに従う系列を得ることができます。

納得いかない人のためにさらなる実験

たまたまでしょ？という風に思われる方もいらっしゃるかもしれませんので、ここで行った実験を1000回繰り返してみましょう。4つ目のセルに以下を打ち込み、実行してください(for文を使って繰り返し実験を行っているだけです。)。

num_of_simulation = 1000

# 推定係数を格納するための箱
estimated_coefficients = []
p_values = []

for i in np.arange(num_of_simulation):

    # 乱数を変更する
    np.random.seed(i)
    # 標準正規分布（もちろん独立）に従って u, v の２系列を作成
    u = np.random.randn(num_of_sample)
    v = np.random.randn(num_of_sample)
    # ランダムウォークに従う系列に変更
    x = np.cumsum(u)
    y = np.cumsum(v)

    # 回帰分析のための定数項を準備
    x_add_const = sm.add_constant(x)

    # 線形回帰の推定と推定結果の出力
    rw_linear_model = sm.OLS(endog=y, exog=x_add_const)
    rw_fitted_linear_model = rw_linear_model.fit()

    # 推定係数とp値を取得する
    estimated_coefficient = rw_fitted_linear_model.params[1]
    p_value = rw_fitted_linear_model.pvalues[1]
    
    # それぞれを格納する
    estimated_coefficients.append(estimated_coefficient)
    p_values.append(p_value)

コードの実行が終わりましたら、次のセルに以下を打ち込み、有意水準5%で帰無仮説が棄却された（すなわち有意な推定係数が得られた）数を確認してみましょう。

# p値が5%未満だった（有意水準5%で帰無仮説が棄却されるケース）割合を計算
null_hypothesis_testing_results = [p_value < 0.05 for p_value in p_values]
print(np.mean(null_hypothesis_testing_results))

私の手元では 91.9% (すなわち、1000回の実験中919回で) 有意な推定係数が得られたことになりました。

最後に以下コードにより、係数の推定値をヒストグラムにしてみましょう。

# 年のため推定係数を確認
fig, ax = plt.subplots(figsize=(15,6))
sns.histplot(
    x = estimated_coefficients,
    ax = ax,
    kde = False
)
ax.grid()
plt.show()

私の手元では以下のグラフが得られました。

つまり、推定係数は別に常に正に有意になったり負に有意になったりするわけではなく、正にも負にも有意な結果をもたらしてしまっているわけですね。めちゃくちゃですね。

対処法

さて、この実験が示すように、時系列性が疑われる（ここでは、今日のある値が昨日の値に影響を受けて出現する）ケースにおいてはそのまま回帰分析を行うと「本来全く関係がなかった２つの系列の間に関係を認めてしまう」というとんでもない誤謬を生むことになります。これは計量経済学でいうところの「系列相関が疑われるケース」であり、分析から頑健な示唆を得るためには絶対に対処せねばなりません。

さてどうすればこの問題に対処できるでしょうか。

解決法１. モデリングする前に差分を取る

ある種当たり前の解決法ですが、足し上げて回帰を行った結果有意な推定係数が得られてしまったのですから、差分を取ってから回帰分析を行うことができれば対処できるのは自明です。

!

ただし今回は実験でしたので、ランダムウォークの系列が得られることがあらかじめ分かっていたということに注意してください。実データの分析では時系列性があるとしても、本当にランダムウォークに従っているかは分かりません。例えば

x_t = 0.8 \times x_{t-1} + u_t

という関係性が実データにおいて成立していた場合、この $x_t$ は時系列性を持ちますが、単に差分を取ることはベストではありません。そのため計量経済学では、「まず変数がどういった時系列性を持つかといったことを検定や別の分析を利用して明らかにする」というステップを挟み、その後どういった差分を取るかを適切に検証することがあります。詳しくはこちらなどを参考にしてください。

解決法２．時系列性を含んだモデリングを行う

今回の分析は、「ある時点の目的変数を同一時点の説明変数により説明する」モデリングとなっており、それ以上のモデリングを何も行っていません。「 $x$ , $y$ のそれぞれが時系列性を持つ」ことを前提にモデリングができれば今回の問題は解決することができます。

具体的な選択肢の一つに「状態空間モデル」を構築するというアプローチがありますが、こちらについてはまたいずれ記事にしたいと思います。こちらのページが非常に参考になりますので、気になる方はご参照ください。

終わりに

今回は知らないでやってしまいがち（かつ本当に恐ろしい誤謬をもたらす）見せかけの回帰について紹介しました。もし少しでも参考になった方などいらっしゃれば、是非「いいね」やシェアなどいただけますと幸いです。

最後まで読んでくださり有難うございました！

Discussion