LoginSignup
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

その他技術 Qiita Tech Festa 2025
Qiita Tech Festa20252025年7月15日まで開催中！
@maskot1977(Ikemen Mas Kot)

VARを使わずにグレンジャー因果性検定をするには？OLSベースでラグと変数の因果性を解析

Posted at

user:

ユーザは pandas dataframe 形式の時系列データと、説明変数を表す列名のリスト、目的変数を表す列名、あわせて3個の引数を渡します。すると、その関数は、各説明変数に対して、以下の計算を行います。

・目的変数に対するグレンジャー因果性を計算するためにまず最大ラグを自動で決定します。
・目的変数に対するグレンジャー因果性を計算し、p値や、各ラグに対する係数を返します。

そのようなPython関数を書いてください。

assistant:

以下は、指定された要件を満たすPython関数 granger_causality_analysis です。この関数は、pandas.DataFrame 形式の時系列データ、説明変数の列名リスト、目的変数の列名を引数として受け取り、各説明変数に対してグレンジャー因果性検定を行い、p値およびラグごとの係数を返します。

import pandas as pd
import statsmodels.api as sm
from statsmodels.tsa.stattools import grangercausalitytests
from typing import List, Dict

def granger_causality_analysis(df: pd.DataFrame, 
                                explanatory_vars: List[str], 
                                target_var: str,
                                max_lag_limit: int = 10) -> Dict[str, Dict[str, any]]:
    """
    各説明変数に対するグレンジャー因果性検定を行い、p値と係数を返す。

    Parameters:
    - df: 時系列データを含む pandas DataFrame
    - explanatory_vars: 説明変数の列名リスト
    - target_var: 目的変数の列名
    - max_lag_limit: 最大ラグの上限（デフォルトは10）

    Returns:
    - Dict: 各説明変数に対する結果を格納した辞書
    """
    results = {}

    for var in explanatory_vars:
        # 検定用データ（目的変数 → 説明変数の順）
        test_df = df[[target_var, var]].dropna()

        # 最適ラグ数をAICで自動選択
        best_lag = None
        min_aic = float('inf')
        for lag in range(1, max_lag_limit + 1):
            try:
                model = sm.tsa.VAR(test_df)
                result = model.fit(lag)
                if result.aic < min_aic:
                    min_aic = result.aic
                    best_lag = lag
            except Exception:
                continue

        if best_lag is None:
            results[var] = {
                "error": "モデルフィッティングに失敗しました。"
            }
            continue

        # グレンジャー因果性検定
        try:
            test_result = grangercausalitytests(test_df, maxlag=best_lag, verbose=False)
            p_value = test_result[best_lag][0]['ssr_ftest'][1]

            # 各ラグの係数（目的変数の回帰における説明変数の係数）
            model = sm.tsa.VAR(test_df)
            result = model.fit(best_lag)
            coefs = result.params.filter(like=var).to_dict()

            results[var] = {
                "best_lag": best_lag,
                "p_value": p_value,
                "coefficients": coefs
            }

        except Exception as e:
            results[var] = {
                "error": f"グレンジャー因果性検定に失敗しました: {str(e)}"
            }

    return results

使い方例:

df = pd.read_csv("your_timeseries_data.csv", parse_dates=True, index_col=0)
explanatory_vars = ['x1', 'x2', 'x3']
target_var = 'y'

result = granger_causality_analysis(df, explanatory_vars, target_var)

import pprint
pprint.pprint(result)

出力形式の例:

{
    'x1': {
        'best_lag': 2,
        'p_value': 0.012,
        'coefficients': {
            'L1.x1': 0.45,
            'L2.x1': -0.12
        }
    },
    'x2': {
        'best_lag': 3,
        'p_value': 0.35,
        'coefficients': {
            'L1.x2': 0.10,
            'L2.x2': 0.05,
            'L3.x2': -0.03
        }
    },
    ...
}

必要に応じて、p値の閾値や係数のしきい値によるフィルタ処理なども追加できます。

user:

係数を求めるのはVARで求めるんですね

assistant:

はい、その通りです。
グレンジャー因果性検定（grangercausalitytests）自体は統計的に「ラグ付きの説明変数が目的変数の予測に有意な改善を与えるかどうか」を検定しますが、各ラグにおける係数そのものは提供しません。

そのため、係数を知りたい場合は、以下のように VAR（Vector Auto Regression）モデルを別途フィッティングして、係数行列（params）を取り出す必要があります。

なぜVARを使うのか？

グレンジャー因果性の仮定は「過去の説明変数が目的変数に影響を与えるかどうか」ですが、それをモデリングするためには 多変量の時系列回帰モデルが必要です。

VARは、次のような形のモデルです：

$$
\mathbf{y}t = A_1 \mathbf{y}{t-1} + A_2 \mathbf{y}{t-2} + \dots + A_p \mathbf{y}{t-p} + \epsilon_t
$$

ここで $\mathbf{y}_t$ はすべての系列（目的変数と説明変数）のベクトルで、
各 $A_i$ は係数行列です。
この中から、目的変数に対する説明変数のラグ項の係数を抽出することで、**「どのラグがどの程度影響しているか」**を定量的に見ることができます。

補足：単回帰 vs VAR

  • grangercausalitytests() の内部では、回帰モデルは「目的変数だけを従属変数」としたものを使っています（他系列は独立変数として使われるだけ）。

  • 一方、VAR モデルを使うと、全ての系列を同時にモデリングします。

    • これは多変量データの時間的相互依存をより良く捉える方法です。

結論

  • p値の計算 → grangercausalitytests()
  • 係数の取得 → VAR.fit().params

両者を組み合わせるのが、実践上もっとも効果的な方法です。

必要であれば、statsmodels.tsa.stattools.grangercausalitytests() を使わずに、OLS を用いて目的変数だけの回帰モデルを構築し、p値も係数も一括で自前で求めることもできます。それをご希望でしたら、代替バージョンのコードもご用意できます。

user:

ひとつ質問があります。VARで係数を求めるにあたって、全ての説明変数を盛り込んだVARを使うのと、今回のように説明変数を１つだけ入れたVARを使うのと、どう違いますか？

assistant:

とても重要な質問です。
**「すべての説明変数を含むVAR」と、「1つの説明変数だけを含むVAR」**のどちらを使うかで、得られる結果の解釈や信頼性に違いが出てきます。

✅ 1. 「1変数ずつのVAR（ペアごとのVAR）」の特徴

✔ 利点

  • 各説明変数の個別の影響を単純に評価できる。
  • グレンジャー因果性検定の典型的な使い方に合っている（grangercausalitytests はこの前提）。

❌ 限界

  • 他の説明変数の影響を無視している。
  • 共通の原因や交絡因子があると、偽の因果関係を検出するリスクがある（例：x₁とx₂の両方がyに影響を与えている場合、x₁だけ見ても誤解を招く可能性）。

✅ 2. 「すべての変数を含むVAR（多変量VAR）」の特徴

✔ 利点

  • 他の説明変数の影響を同時に制御した上で、各変数の効果を評価できる。
  • 交絡因子をある程度制御できるため、因果性の主張に説得力が増す。
  • 複雑な相互依存をモデリングできる（例：x₁がyに、x₂がx₁に影響など）。

❌ 限界

  • モデルが大きくなり、推定に多くのデータを必要とする。
  • 高次ラグや多変数が入ると過学習のリスクがある。
  • 各係数の解釈が複雑になる（相関構造をすべて考慮する必要あり）。

🎯 まとめ：いつどちらを使うか？

目的 推奨モデル 理由
個別の説明変数が目的変数に影響するかを単純に見たい 1変数ごとのVAR（今の方法） 実装が簡単、直感的、説明しやすい
交絡要因の影響を制御した上で、より厳密な因果性を評価したい 多変量VAR（全変数入り） 現実的な因果解釈、誤検出リスクの低減

🔧 実践的なアプローチ

「両方やって比較」するのが現実的です：

  1. まずは各説明変数で個別にグレンジャー因果性検定（今の方式）

  2. 次に、すべての説明変数を入れたVARで、目的変数に対する係数を見て影響を比較

    • それぞれの係数の有意性（t値、p値）を見る
    • summary() や test_causality() を使うと便利です

ご希望があれば、全変数を含むVARモデルを使って test_causality() で説明変数ごとのグレンジャー因果性を検定する関数もご提供できます。

user:

ありがとうございます。では、話を戻して、「statsmodels.tsa.stattools.grangercausalitytests() を使わずに、OLS を用いて目的変数だけの回帰モデルを構築し、p値も係数も一括で自前で求める」関数を書いてください。

assistant:

こちらこそ、明確なご要望ありがとうございます。
以下は、statsmodels.api.OLS を用いて 目的変数のラグ + 説明変数のラグ を含む回帰モデルを構築し、各説明変数に対するグレンジャー因果性の p 値とラグごとの係数を一括で取得する関数です。

✅ ols_based_granger_causality 関数

import pandas as pd
import numpy as np
import statsmodels.api as sm
from typing import List, Dict

def ols_based_granger_causality(df: pd.DataFrame, 
                                 explanatory_vars: List[str], 
                                 target_var: str,
                                 max_lag_limit: int = 10) -> Dict[str, Dict[str, any]]:
    """
    OLS回帰を使って、各説明変数が目的変数にグレンジャー因果性を持つかを検定する。
    各説明変数に対して、最適ラグ数の選定 → ラグ付き回帰 → p値と係数の取得を行う。

    Parameters:
        df: pandas DataFrame（時系列データ）
        explanatory_vars: 説明変数の列名リスト
        target_var: 目的変数の列名
        max_lag_limit: 最大ラグ数（デフォルト10）

    Returns:
        Dict: 各説明変数に対する p値とラグ係数を含む結果辞書
    """
    results = {}

    for ex_var in explanatory_vars:
        min_aic = np.inf
        best_lag = None

        # 最適ラグ選定（AICベース）
        for lag in range(1, max_lag_limit + 1):
            df_lagged = pd.DataFrame()

            # ラグを作成
            for l in range(1, lag + 1):
                df_lagged[f'{target_var}_lag{l}'] = df[target_var].shift(l)
                df_lagged[f'{ex_var}_lag{l}'] = df[ex_var].shift(l)
            df_lagged['y'] = df[target_var]

            df_lagged = df_lagged.dropna()

            # 説明変数と目的変数を分ける
            X = df_lagged.drop(columns='y')
            y = df_lagged['y']

            X = sm.add_constant(X)
            try:
                model = sm.OLS(y, X).fit()
                if model.aic < min_aic:
                    min_aic = model.aic
                    best_lag = lag
            except:
                continue

        # 最適ラグで再構築
        if best_lag is not None:
            df_lagged = pd.DataFrame()
            for l in range(1, best_lag + 1):
                df_lagged[f'{target_var}_lag{l}'] = df[target_var].shift(l)
                df_lagged[f'{ex_var}_lag{l}'] = df[ex_var].shift(l)
            df_lagged['y'] = df[target_var]
            df_lagged = df_lagged.dropna()

            X = df_lagged.drop(columns='y')
            y = df_lagged['y']
            X = sm.add_constant(X)

            model = sm.OLS(y, X).fit()

            # 説明変数のラグ項に対する係数とp値
            coef_dict = {
                k: v for k, v in model.params.items() if ex_var in k
            }
            pval_dict = {
                k: v for k, v in model.pvalues.items() if ex_var in k
            }

            results[ex_var] = {
                'best_lag': best_lag,
                'p_values': pval_dict,
                'coefficients': coef_dict
            }
        else:
            results[ex_var] = {
                'error': '有効なラグでモデル構築できませんでした。'
            }

    return results

✅ 使用例

df = pd.read_csv("your_timeseries.csv", parse_dates=True, index_col=0)
explanatory_vars = ['x1', 'x2', 'x3']
target_var = 'y'

res = ols_based_granger_causality(df, explanatory_vars, target_var)

import pprint
pprint.pprint(res)

✅ 出力例

{
    'x1': {
        'best_lag': 2,
        'p_values': {
            'x1_lag1': 0.03,
            'x1_lag2': 0.10
        },
        'coefficients': {
            'x1_lag1': 0.42,
            'x1_lag2': -0.08
        }
    },
    ...
}

✅ 特徴と利点

  • grangercausalitytests に依存せず、p値と係数を両方取得可能
  • AIC によるラグ自動選択
  • ラグ付き回帰の構造が明示的なのでカスタマイズしやすい

必要であれば、目的変数のラグを含まないモデル（条件付き独立）との F検定や比較 を入れることも可能です。ご希望があれば拡張版もご用意します。

user:

ありがとうございます。もしも、最適ラグが３だったとしましょう。それにもかかわらずラグ１２まで計算した場合、そのp値などは信頼できなくなりますか？

assistant:

とても鋭い質問です。結論からお答えすると：

✅【結論】

「ラグ3が最適であるにもかかわらず、ラグ12まで含めて回帰した場合、係数やp値の信頼性は低下する可能性が高い」です。

🔍 なぜ信頼性が下がるのか？

1. パラメータ数が増えてモデルが過剰になる

  • ラグ12まで使うということは、目的変数のラグ + 説明変数のラグで 合計 24項以上の変数 を回帰モデルに入れることになります（lag12 × 2系列）。
  • サンプルサイズに比してパラメータ数が多くなると、多重共線性や過学習のリスクが高まり、各係数の推定誤差が大きくなります。
  • その結果、p値が不安定になります（偽陽性や偽陰性が増える）。

2. 情報基準（AIC/BIC）によるラグ選定は「汎化性能の最大化」が目的

  • ラグ3が最適と判定されたのは、AIC（またはBIC）に基づいて「データへの当てはまりとモデルの簡潔さのバランス」が取れたポイントだからです。
  • それより多いラグを入れるのは、ノイズの過剰な取り込みになりやすく、推定の精度や解釈性が悪化します。

✅【例外】

ただし、以下のような場合はラグ12を使うのが理にかなっていることもあります。

🔸 ドメイン知識がある場合：

  • 例えば「四半期データで前年同期の効果がある」という仮説があるなら、lag=4やlag=12を入れる意味があります。

🔸 グレンジャー因果性の「存在」だけ見たい場合：

  • モデルの適合度には関係なく「何かしらのラグで因果性があるか」だけを調べるときは、高めのラグも許容される。

✅【おすすめの実践戦略】

アプローチ 説明
最適ラグ（AIC/BICベース）を使う 汎化性能・安定性が高く、推定結果の信頼性が高い
固定ラグ（例えば12）で全変数一律に使う 比較が簡単、ただし過剰なラグに注意
両方やって比較 最適ラグ vs. 固定ラグでの結果を比較し、頑健性を検証する

✅ 結論まとめ

「最適ラグが3であるにもかかわらず12まで入れると、p値や係数はノイズの影響を受けて不安定になるため、基本的には避けるべき」です。

ただし、目的や仮説、データサイズによっては柔軟に判断することも可能です。

ご希望であれば「最適ラグ vs 固定ラグ」の両方を比較出力する関数にもできます。

user:

ありがとうございます。ここでOLSを使う理由について説明していただけますか？

assistant:

はい、非常に重要なポイントです。
ここで OLS（最小二乗法）を使う理由について、以下の観点から丁寧に解説します：

✅ 1. グレンジャー因果性検定の本質は「線形回帰モデルの比較」

グレンジャー因果性の考え方はこうです：

「ある説明変数 X の過去の値が、目的変数 Y の予測精度を統計的に有意に向上させるなら、XはYのグレンジャー原因である」

この「予測精度の向上が統計的に有意か？」を調べるために、以下の2つの回帰モデルをOLS（線形回帰）で構築して比較します。

モデル名 含まれる説明変数
制限モデル（帰無仮説） Yのラグのみ
拡張モデル（対立仮説） Yのラグ ＋ Xのラグ

→ OLSにより両モデルの残差（二乗和）を比較して、F検定または係数のt検定を行う、というのがグレンジャー因果性検定の基本構造です。

✅ 2. OLSはパラメータ推定・係数・p値の取得が容易

statsmodels.api.OLS は以下を直接提供してくれます：

  • 各ラグ係数（推定値）
  • 各係数のp値
  • モデルのAIC・BIC・決定係数
  • F検定も可能（必要に応じて）

これにより、モデルの解釈がしやすく、自由度の高いカスタマイズも可能です。

✅ 3. VAR（Vector Autoregression）も本質はOLSの拡張

VARモデルは、複数系列の多変量自己回帰モデルですが、内部的には「各系列に対してOLSを適用して推定」しています。

したがって、今回のように：

  • 説明変数を1つずつ調べる
  • 目的変数が1つのみ

という状況では、VARを使わずOLSだけで十分ですし、むしろその方がシンプルです。

✅ 4. 自由度・サンプルサイズに対する柔軟性

OLSは少ないサンプル数でも比較的安定して推定できます。VARのように複数変数を一括で扱うとパラメータ数が爆発しやすいですが、OLSを1対1の関係で使えば過剰パラメータを避けられます。

✅ まとめ：なぜOLSを使うのか？

理由 説明
理論 グレンジャー因果性の検定はOLSモデルの比較に基づいているから
実用 p値や係数を明示的に取り出せて、解釈・可視化がしやすい
柔軟性 カスタマイズ（ラグ選定・回帰式構成など）が容易
効率性 VARよりも少ないパラメータで安定した推定が可能

もしご希望であれば、「OLSを使った2ステップF検定によるグレンジャー因果性検定」を手作業で実装することもできますよ。

user:

ありがとうございます。もう一度、話を戻します。もしも、最適ラグが３だったとしましょう。そうすると、ラグ４やラグ５は因果には影響しないと考えて良いですか？

assistant:

とても重要なご質問です。

✅【結論】

最適ラグが3であるからといって、「ラグ4やラグ5が因果に影響しない」と断言はできません。
ただし、「ラグ3までの情報があればそれ以降のラグを加えてもモデルの性能が大きく向上しない」という意味では、ラグ4や5の影響は相対的に小さいと判断された、ということです。

🔍 解説：ラグ選定と「因果のラグ」の違い

◉ ラグ選定（例えば AIC/BIC で最適ラグ3と判定された）とは…

  • 全体の予測性能が最も良くなるラグ長を意味します。
  • これはモデル全体（自己回帰 + 説明変数）の構造に対しての統計的なバランスを見ています。
  • よって、目的変数や説明変数の個々のラグの有意性を示すものではありません。

◉ 「ラグ4やラグ5は因果に影響しない」とは言い切れない理由

  1. 情報量としては加えられるが、モデルに有意な改善を与えないだけ

    • ラグ4や5が微弱ながら有意に影響する可能性もある。
    • ただし、ラグ3まででほとんどの因果関係が説明されてしまっているため、追加の寄与が小さい（≠ゼロではない）。

  2. p値の検定はラグ個別に確認すべき

    • 仮にラグ12まで使ったとしても、「ラグ4の係数だけ有意」というケースはあります。
    • したがって、「最適ラグ3」という情報だけでラグ4や5の有意性を否定するのは不適切です。

  3. マルチコリニアリティの影響

    • ラグ変数間に強い相関があると、ラグ4や5の影響が他のラグに吸収されて見えなくなることもあります。
    • これは「ラグ4が無意味」ではなく、「モデル内で説明力が重複している」ということ。

✅ 実務での使い分け

意図 推奨アプローチ
モデル全体の予測性能を重視 AIC/BICで選ばれた最適ラグだけ使う
各ラグの因果的影響の有無を知りたい 高めのラグまで含めて、各ラグのp値や係数を確認する

✅ まとめ

「最適ラグ3」という結果は、ラグ4や5が「まったく無意味」という意味ではありません。
あくまで予測性能に対して追加効果が小さい、というだけで、因果的影響の有無は個別に検定して判断すべきです。

必要であれば「ラグ最大12まで入れたモデルを作成し、各ラグに対する係数とp値を出力する関数」も提供できます。ご希望ですか？

user:

そうですね、基本的にはols_based_granger_causalityと同じになりそうですが、「ユーザが指定した最大ラグまで、各ラグに対する係数とp値を出力する関数」を書いてください。

assistant:

以下は、ユーザが指定した最大ラグ数まで、各説明変数の各ラグに対する係数とp値を出力する関数です。ols_based_granger_causality() と似ていますが、最適ラグを自動選択せず、指定された最大ラグまで常に計算します。

✅ 関数の仕様：

  • 引数

    • df: pandas DataFrame（時系列データ）
    • exog_vars: 説明変数の列名リスト
    • endog_var: 目的変数の列名
    • max_lag: ユーザが指定する最大ラグ数

  • 出力

    • 各説明変数について：

      • ラグごとの係数とp値（辞書形式）

✅ コード

import pandas as pd
import numpy as np
import statsmodels.api as sm

def granger_lagwise_ols(df, exog_vars, endog_var, max_lag):
    df = df.copy()
    results = {}

    # ラグ変数を生成
    for lag in range(1, max_lag + 1):
        df[f'{endog_var}_lag{lag}'] = df[endog_var].shift(lag)
        for exog in exog_vars:
            df[f'{exog}_lag{lag}'] = df[exog].shift(lag)

    df.dropna(inplace=True)

    for exog in exog_vars:
        # 説明変数に含めるラグ変数（endogとexogの両方）
        X_cols = [f'{endog_var}_lag{i}' for i in range(1, max_lag + 1)] + \
                 [f'{exog}_lag{i}' for i in range(1, max_lag + 1)]
        X = df[X_cols]
        y = df[endog_var]

        X = sm.add_constant(X)
        model = sm.OLS(y, X).fit()

        # ラグごとの係数とp値
        lag_results = {}
        for i in range(1, max_lag + 1):
            coef = model.params.get(f'{exog}_lag{i}', np.nan)
            pval = model.pvalues.get(f'{exog}_lag{i}', np.nan)
            lag_results[f'lag_{i}'] = {'coef': coef, 'pval': pval}

        results[exog] = lag_results

    return results

✅ 使用例

result = granger_lagwise_ols(df, exog_vars=['x1', 'x2'], endog_var='y', max_lag=12)

# 結果を確認
import pprint
pprint.pprint(result)

✅ 出力形式（例）

{
  'x1': {
    'lag_1': {'coef': 0.213, 'pval': 0.041},
    'lag_2': {'coef': -0.019, 'pval': 0.872},
    ...
    'lag_12': {'coef': 0.093, 'pval': 0.157}
  },
  'x2': {
    ...
  }
}

必要があれば、これを DataFrame に変換して表形式で出力する関数も提供できます。ご希望ですか？

user:

はい、DataFrame に変換して表形式で出力する関数にしてください。

assistant:

以下は、前回の関数を拡張して、結果を pandas DataFrame 形式で出力するようにしたバージョンです。

✅ 出力形式（DataFrame）

exog_var lag coef pval
x1 1 0.213 0.041
x1 2 -0.019 0.872
... ... ... ...
x2 1 0.091 0.012
x2 2 -0.053 0.453

✅ 関数コード

import pandas as pd
import numpy as np
import statsmodels.api as sm

def granger_lagwise_ols_table(df, exog_vars, endog_var, max_lag):
    df = df.copy()
    rows = []

    # ラグ変数を生成
    for lag in range(1, max_lag + 1):
        df[f'{endog_var}_lag{lag}'] = df[endog_var].shift(lag)
        for exog in exog_vars:
            df[f'{exog}_lag{lag}'] = df[exog].shift(lag)

    df.dropna(inplace=True)

    for exog in exog_vars:
        # 説明変数に含めるラグ変数（endogとexogの両方）
        X_cols = [f'{endog_var}_lag{i}' for i in range(1, max_lag + 1)] + \
                 [f'{exog}_lag{i}' for i in range(1, max_lag + 1)]
        X = df[X_cols]
        y = df[endog_var]

        X = sm.add_constant(X)
        model = sm.OLS(y, X).fit()

        # ラグごとの係数とp値を保存
        for i in range(1, max_lag + 1):
            coef = model.params.get(f'{exog}_lag{i}', np.nan)
            pval = model.pvalues.get(f'{exog}_lag{i}', np.nan)
            rows.append({'exog_var': exog, 'lag': i, 'coef': coef, 'pval': pval})

    result_df = pd.DataFrame(rows)
    return result_df

✅ 使用例

result_df = granger_lagwise_ols_table(df, exog_vars=['x1', 'x2'], endog_var='y', max_lag=12)
print(result_df)

✅ オプション（必要なら追加）

  • 指定した p 値未満の結果だけフィルタする
  • プロット（棒グラフやヒートマップ）で可視化
  • CSV/Excel に出力する

必要であればこれらもご用意できます。

user:

ありがとうございます。ラグの影響に関するp値は上記の方法で算出できましたが、ラグ変数ではなく元の変数間の因果のp値は計算できますか？

assistant:

とても良い質問です。

✅ 要点

  • ラグ変数に対する p 値は、各ラグの係数が「ゼロでないか」を検定しています。
  • しかし、元の変数（たとえば x1）が目的変数 y にグレンジャー因果性を持つかどうかは、「全ラグの係数が同時にゼロではない」ことを検定する必要があります。
  • これは F検定（多重制約の仮説検定） であり、statsmodels で可能です。

✅ どうやって「元の変数が因果性を持つか」の p 値を求めるか？

1. フルモデル（目的変数＋説明変数のラグ全て）を構築

2. 制約付きモデル（説明変数のラグなし）を構築

3. F検定で「制約が妥当か」を比較 → p 値が低ければ因果あり

✅ 実装例（変数単位の因果性検定）

def granger_variable_level_pval(df, exog_vars, endog_var, max_lag):
    df = df.copy()

    # ラグ変数の作成
    for lag in range(1, max_lag + 1):
        df[f'{endog_var}_lag{lag}'] = df[endog_var].shift(lag)
        for var in [endog_var] + exog_vars:
            df[f'{var}_lag{lag}'] = df[var].shift(lag)

    df.dropna(inplace=True)
    results = []

    # フルモデルのベース（目的変数自身のラグ）
    y = df[endog_var]
    base_X = [f'{endog_var}_lag{i}' for i in range(1, max_lag + 1)]

    for exog in exog_vars:
        # モデル1（フルモデル）
        X_full = base_X + [f'{exog}_lag{i}' for i in range(1, max_lag + 1)]
        X1 = sm.add_constant(df[X_full])
        model_full = sm.OLS(y, X1).fit()

        # モデル2（制約付きモデル）
        X2 = sm.add_constant(df[base_X])
        model_restricted = sm.OLS(y, X2).fit()

        # F検定
        f_test = model_full.compare_f_test(model_restricted)
        f_stat, p_value, _ = f_test

        results.append({
            'exog_var': exog,
            'f_stat': f_stat,
            'pval': p_value
        })

    return pd.DataFrame(results)

✅ 出力形式

| exog_var | f_stat | pval  |
|----------|--------|-------|
| x1       | 4.89   | 0.003 |
| x2       | 0.57   | 0.841 |

✅ 解釈

  • pval < 0.05 なら、その変数は目的変数にグレンジャー因果性ありと判断。
  • この関数で「変数レベルの因果性の有無（統計的に有意かどうか）」をチェックできます。

ご希望があれば、この出力と先ほどのラグ別出力を 1つの表に統合することも可能です。ご希望ですか？

user:

pval が None になることがあります。

assistant:

はい、それは重要な観察です。pval が None になるのは、以下のような状況で compare_f_test() がF検定を正しく計算できなかった場合に発生します。

✅ pval が None になる主な原因

  1. 自由度の不足（自由度がゼロ以下）

    • ラグ数が多すぎて、サンプル数が足りないと発生しやすいです。
    • 特に小さいデータセットに max_lag=12 などを指定した場合。

  2. モデルが完全に同一になる（差がない）

    • exog のラグ系列がすべてゼロや定数になっている、または endog と完全に同じ値を取っている。

  3. 数値的な問題（線形依存、多重共線性）

    • 説明変数に多重共線性があり、行列が特異（Singular）になっている。

✅ 対処方法

方法 1：ラグ数とデータ数をチェック

print(f"データ数: {len(df)}、推定に使う行数: {len(df.dropna())}、自由度: {len(df.dropna()) - パラメータ数}")
  • 一般に「推定に使える行数 ＞ パラメータ数（定数項 + ラグ × 説明変数数）」である必要があります。

方法 2：compare_f_test() の失敗を補足し、None を回避

try:
    f_stat, p_value, _ = model_full.compare_f_test(model_restricted)
except Exception as e:
    f_stat, p_value = np.nan, np.nan

✅ 修正済み関数（安全なバージョン）

def granger_variable_level_pval(df, exog_vars, endog_var, max_lag):
    df = df.copy()

    # ラグ変数の作成
    for lag in range(1, max_lag + 1):
        df[f'{endog_var}_lag{lag}'] = df[endog_var].shift(lag)
        for var in [endog_var] + exog_vars:
            df[f'{var}_lag{lag}'] = df[var].shift(lag)

    df.dropna(inplace=True)
    results = []

    y = df[endog_var]
    base_X = [f'{endog_var}_lag{i}' for i in range(1, max_lag + 1)]

    for exog in exog_vars:
        try:
            X_full = base_X + [f'{exog}_lag{i}' for i in range(1, max_lag + 1)]
            X1 = sm.add_constant(df[X_full])
            model_full = sm.OLS(y, X1).fit()

            X2 = sm.add_constant(df[base_X])
            model_restricted = sm.OLS(y, X2).fit()

            f_stat, p_value, _ = model_full.compare_f_test(model_restricted)

        except Exception as e:
            f_stat, p_value = np.nan, np.nan  # または None, None

        results.append({
            'exog_var': exog,
            'f_stat': f_stat,
            'pval': p_value
        })

    return pd.DataFrame(results)

✅ アドバイス

  • 時系列データの行数が 数百行以上 あるのが理想です（ラグ数が多い場合）。
  • max_lag を抑える（例：5以下）とエラーが出にくくなります。

もしご希望であれば、「データ数とラグ数から自由度を自動チェックし、警告を出す仕組み」も追加可能です。どうされますか？

user:

ありがとうございます。今日はこの辺で大丈夫です。ここまでの議論を踏まえて、ブログに記録を残したいのですが、タイトル案をいくつか考えていただけますか？

assistant:

こちらこそ、非常に本質的で実践的なご質問をありがとうございました。以下は、今回の議論（グレンジャー因果性・OLSによる検定・ラグ選択・変数レベルでのF検定など）を踏まえたブログ記事のタイトル案です：

🔍 技術的・実装重視のタイトル案

  1. 「グレンジャー因果性をOLSで実装してみた：p値とラグ係数を自前で出力するPython関数」
  2. 「VARを使わずにグレンジャー因果性検定をするには？OLSベースでラグと変数の因果性を解析」
  3. 「グレンジャー因果性を1行ずつ理解する：回帰式・ラグ構造・F検定をPythonで解剖」
  4. 「statsmodelsの裏側に迫る：グレンジャー因果性をF検定で手計算する方法」
  5. 「ラグごとのp値、変数ごとの因果性：Pythonでグレンジャー因果を分解可視化」

🧠 理論と実装のバランスを取ったタイトル案

  1. 「ラグと因果を分けて考える：グレンジャー因果性をOLSとF検定で実装する」
  2. 「時系列の因果をどう捉える？ラグ係数 vs. 変数レベルのグレンジャー因果性」
  3. 「グレンジャー因果性を再定義する：ラグ単位の検定と変数単位のF検定」
  4. 「『この変数は因果あり？』をp値で判断する：Pythonで実践するグレンジャー検定」

🗂 少しライトな（読者を惹きつける）タイトル案

  1. 「そのp値、本当に意味ある？グレンジャー因果性とラグの罠」
  2. 「グレンジャー因果性を疑え：ラグと変数、因果の真相をPythonで探る」
  3. 「時系列データに『因果』はあるのか？ラグと変数をOLSで検証してみた」

必要に応じて、本文の構成案やサブタイトル案もご提案できます。
気に入ったタイトルがあれば、それに合わせた冒頭文や構成テンプレートもお手伝いできますので、遠慮なくお知らせください！

user:

「VARを使わずにグレンジャー因果性検定をするには？OLSベースでラグと変数の因果性を解析」でいきたいと思います。最後に、詳細なまとめを書いてください。

assistant:

素晴らしいタイトルのご選択です！
以下に、**「VARを使わずにグレンジャー因果性検定をするには？OLSベースでラグと変数の因果性を解析」という記事にふさわしい詳細なまとめ（まとめ＋学びの振り返り）**をお届けします。

📝 まとめ：VARを使わずにグレンジャー因果性検定をするには？

本記事では、**時系列における因果性（グレンジャー因果性）**を検証するために、一般的な VAR モデルを使わず、OLS（最小二乗法）による回帰分析をベースとした方法でアプローチしました。

その中で、以下のポイントを体系的に整理しました：

✅ OLSベースのグレンジャー因果性検定の枠組み

  • 従来のアプローチでは、statsmodels.tsa.stattools.grangercausalitytests() が内部的に VAR モデルを構築してグレンジャー因果性の検定を行います。
  • しかし、今回のアプローチでは、目的変数に対してラグ付き説明変数をOLSで逐次追加していき、そのラグに対する回帰係数とp値を個別に算出しました。
  • 加えて、変数単位での因果性有無の検定には、**制約付きモデル（目的変数の自己ラグのみ）と拡張モデル（説明変数のラグも含む）**を構築し、**F検定（モデル間の有意差検定）**を適用しました。

✅ 本手法の意義と活用シーン

  • ラグごとの影響を可視化できるため、どの時間差（遅れ）が因果に関与しているのかを丁寧に分析可能。
  • 変数レベルの因果性を明示的にテストでき、従来のブラックボックス的な p 値解釈よりも、構造的な理解が深まります。
  • 自前実装により、p値の検定手順やラグ選択の根拠を透明化できる点も大きなメリット。

⚠️ 注意点と実務上の留意事項

  • ラグ数が増えると自由度が減少し、小さなデータセットではモデルが不安定になります。
  • ラグごとのp値が有意でなくても、変数全体としては有意な因果関係がある場合があります（F検定の活用が有効）。
  • p値が None になる場合は、**自由度不足や数値的な特異性（多重共線性など）**が背景にあるため、ラグ数やデータ数に注意が必要です。

📌 今回の手法で得られるアウトプット

ラグ 回帰係数（β） p値
lag1 0.123 0.045
lag2 -0.032 0.610
lag3 0.210 0.002
... ... ...

加えて：

説明変数 F値 p値（変数レベル）
varA 5.43 0.008
varB 1.15 0.318

🎯 最後に：因果性の「構造」を自分の手で理解するということ

VAR を使わずに OLS ベースでグレンジャー因果性を検定する方法は、「モデルの構造を理解しながら因果性を調べたい」というニーズに非常に適しています。
ラグの影響・変数の寄与・モデルの自由度などを意識しながら検証できるため、単なる「Yes/No判定」ではなく、解釈可能な因果探索が可能になります。

今後は、この手法を土台にして、

  • ラグ選択の自動化（AIC/BICベース）
  • 多変量系（複数の目的変数）
  • 階層的因果分析
    などに展開していくことも視野に入るでしょう。

ご希望があれば、記事用のイントロ・セクション構成案・図表の挿入提案などもお手伝いできます。お気軽にどうぞ！

user:

ありがとうございます。これで充分です。また次回もよろしくお願いします！

assistant:

こちらこそ、深く本質的な議論をありがとうございました！
またいつでもお力になれたら嬉しいです。次回もどうぞよろしくお願いいたします。
良い執筆＆分析ライフを！📝📊

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?