Granger因果性の統計的検定をPythonでスクラッチ実装する

Last updated at 2025-04-01Posted at 2025-04-01

user:

statsmodels.tsa.stattools.grangercausalitytests の原理を勉強したいと思います。statsmodels.tsa.stattools.grangercausalitytests と同等の計算を行うPython関数を、 statsmodels.tsa.stattools.grangercausalitytests を使わずに書いてください。

assistant:

Granger因果性検定は、時系列データのある系列が他の系列の予測に役立つかどうかを評価するために使用されます。statsmodels.tsa.stattools.grangercausalitytests を使わずに、同等の計算を行う関数を実装してみます。

import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.regression.linear_model import OLS

def granger_causality_test(data, max_lag, verbose=True):
    """
    Granger causality test without using statsmodels.tsa.stattools.grangercausalitytests.
    
    Parameters:
    data : pd.DataFrame
        A dataframe with two time series columns.
    max_lag : int
        The maximum number of lags to consider.
    verbose : bool, optional
        If True, prints test results.
    
    Returns:
    dict :
        A dictionary with lag values as keys and F-test results as values.
    """
    results = {}
    
    for lag in range(1, max_lag + 1):
        df = data.copy()
        
        # Create lagged variables
        for i in range(1, lag + 1):
            df[f'X_lag{i}'] = df.iloc[:, 0].shift(i)
            df[f'Y_lag{i}'] = df.iloc[:, 1].shift(i)
        
        df.dropna(inplace=True)
        
        # Endogenous variable (Y)
        Y = df.iloc[:, 1]
        
        # Restricted model: Y_t ~ Y_t-1, ..., Y_t-lag
        X_restricted = df[[f'Y_lag{i}' for i in range(1, lag + 1)]]
        X_restricted = sm.add_constant(X_restricted)
        model_restricted = OLS(Y, X_restricted).fit()
        rss_restricted = np.sum(model_restricted.resid ** 2)
        
        # Unrestricted model: Y_t ~ Y_t-1, ..., Y_t-lag + X_t-1, ..., X_t-lag
        X_unrestricted = df[[f'Y_lag{i}' for i in range(1, lag + 1)] + 
                             [f'X_lag{i}' for i in range(1, lag + 1)]]
        X_unrestricted = sm.add_constant(X_unrestricted)
        model_unrestricted = OLS(Y, X_unrestricted).fit()
        rss_unrestricted = np.sum(model_unrestricted.resid ** 2)
        
        # Compute F-statistic
        n = len(Y)
        k = lag
        F_stat = ((rss_restricted - rss_unrestricted) / k) / (rss_unrestricted / (n - 2*k - 1))
        p_value = 1 - sm.stats.f.sf(F_stat, k, n - 2*k - 1)
        
        results[lag] = {'F-stat': F_stat, 'p-value': p_value}
        
        if verbose:
            print(f'Lag {lag}: F-stat={F_stat:.4f}, p-value={p_value:.4f}')
    
    return results

この関数 granger_causality_test は、Granger因果性検定を statsmodels.tsa.stattools.grangercausalitytests を使わずに実装しています。max_lag の範囲で遅延変数を作成し、制約付き・非制約付きモデルの回帰を行い、F検定を実施します。

user:

ありがとうございます。そのコードについて、詳細かつ厳密に説明してください。

assistant:

この granger_causality_test 関数がどのように Granger 因果性検定を行うのか、詳細かつ厳密に説明します。

1. Granger因果性検定の概要

Granger因果性検定は、時系列データにおいて変数 $X$ が変数 $Y$ の将来値を予測するのに有意な情報を提供するかどうかを判定する手法です。

基本的な考え方は次の通りです：

$Y_t$ を時点 $t$ における従属変数とする。
$X_t$ を時点 $t$ における説明変数とする。
$Y_t$ を過去の自身の値のみに基づいて回帰する 制約モデル (restricted model) と、過去の (X) の値も考慮した 非制約モデル (unrestricted model) を構築し、それらの誤差（残差平方和, RSS）を比較する。
2つのモデルの誤差を F 検定を用いて比較し、$X$ の過去の情報が $Y$ の予測に有意に寄与しているかを判断する。

2. `granger_causality_test` の詳細説明

関数は次のステップで実装されています。

(1) 関数の定義

def granger_causality_test(data, max_lag, verbose=True):

data : pandas.DataFrame
- 2列からなる時系列データ (X, Y) を含むデータフレーム。
max_lag : int
- 最大のラグ（遅延）を指定。
verbose : bool
- 結果を出力するかどうかを決定。

(2) 結果を格納する辞書を初期化

results = {}

検定結果（F統計量と p値）をラグごとに格納するための辞書。

(3) ラグごとのループ

for lag in range(1, max_lag + 1):

max_lag の範囲で Granger 因果性を検定する。

(4) 遅延変数を作成

df = data.copy()
for i in range(1, lag + 1):
    df[f'X_lag{i}'] = df.iloc[:, 0].shift(i)
    df[f'Y_lag{i}'] = df.iloc[:, 1].shift(i)
df.dropna(inplace=True)

shift(i) を使い、X および Y の i ステップ前の値を新しいカラムとして作成。
欠損値を除去 (dropna)、これにより最初の lag 行は削除される。

(5) 従属変数 Y_t を設定

Y = df.iloc[:, 1]

Y_t を取り出す（データフレームの2列目）。

(6) 制約モデル (Restricted Model) の構築

X_restricted = df[[f'Y_lag{i}' for i in range(1, lag + 1)]]
X_restricted = sm.add_constant(X_restricted)
model_restricted = OLS(Y, X_restricted).fit()
rss_restricted = np.sum(model_restricted.resid ** 2)

制約モデル（Restricted Model）：
- $ Y_t $ を自身の過去の値 $ Y_{t-1}, Y_{t-2}, \dots, Y_{t-\text{lag}} $ だけで回帰。
- 定数項 (intercept) を追加 (sm.add_constant)。
- OLS（最小二乗法）で回帰。
- 残差平方和 (RSS) を計算。

(7) 非制約モデル (Unrestricted Model) の構築

X_unrestricted = df[[f'Y_lag{i}' for i in range(1, lag + 1)] + 
                     [f'X_lag{i}' for i in range(1, lag + 1)]]
X_unrestricted = sm.add_constant(X_unrestricted)
model_unrestricted = OLS(Y, X_unrestricted).fit()
rss_unrestricted = np.sum(model_unrestricted.resid ** 2)

非制約モデル（Unrestricted Model）：
- $ Y_t $ を自身の過去の値＋ $X$ の過去の値 $X_{t-1}, X_{t-2}, \dots, X_{t-\text{lag}} $ で回帰。
- RSS（残差平方和）を計算。

(8) F検定の計算

n = len(Y)
k = lag
F_stat = ((rss_restricted - rss_unrestricted) / k) / (rss_unrestricted / (n - 2*k - 1))
p_value = 1 - sm.stats.f.sf(F_stat, k, n - 2*k - 1)

F検定統計量の公式：
$$
F = \frac{(\text{RSS}{\text{restricted}} - \text{RSS}{\text{unrestricted}}) / k}{\text{RSS}_{\text{unrestricted}} / (n - 2k - 1)}
$$
- $ k $ はラグの数。
- $ n $ はサンプル数。
- sm.stats.f.sf を使って p値を求める。

(9) 結果の格納

results[lag] = {'F-stat': F_stat, 'p-value': p_value}
if verbose:
    print(f'Lag {lag}: F-stat={F_stat:.4f}, p-value={p_value:.4f}')

各ラグについて F統計量と p値を results 辞書に格納。
verbose=True の場合、結果を出力。

(10) 結果を返す

return results

results は、キーが lag、値が {'F-stat': F値, 'p-value': p値} の辞書。

3. まとめ

この関数は、statsmodels.tsa.stattools.grangercausalitytests を使用せずに Granger 因果性検定を行います。

max_lag の範囲でラグを設定。
X と Y の遅延変数を作成。
Y_t を従属変数とする 制約モデル（自身の過去のみ）と 非制約モデル（自身の過去＋(X)の過去）を構築。
F 検定を用いて 2 つのモデルの差異を評価。
結果として、各ラグごとの F統計量と p値を出力。

これにより、$X$ の過去の情報が $Y$ の予測に寄与するかどうかを判定できます。

user:

ありがとうございます。次のように使ってみたんですが、

def sum_of_sigmoids(A, x):
    """
    A: M×3の行列 (numpy.ndarray)
       - A[:, 0] -> B (スケール)
       - A[:, 1] -> C (傾き)
       - A[:, 2] -> -D/C (変曲点のx座標)
    x: 計算する点 (float または numpy.ndarray)

    返り値:
        M個のシグモイド関数の和
    """
    A = np.asarray(A)  # 念のため ndarray に変換
    if A.shape[1] != 3:
        raise ValueError("A は M×3 の行列でなければなりません。")

    B = A[:, 0]
    C = A[:, 1]
    x_inflection = A[:, 2]  # 変曲点のx座標 (-D/C)
    D = -C * x_inflection  # Dを計算

    # M個のシグモイド関数を計算し、それらの合計を返す
    return np.sum(B[:, None] / (1 + np.exp(-(C[:, None] * x + D[:, None]))), axis=0)


def sum_of_sin(A, x):
    """
    A: M×3の行列 (numpy.ndarray)
       - A[:, 0] -> 振幅 (Amplitude)
       - A[:, 1] -> 周期 (Period)
       - A[:, 2] -> 位相 (Phase)
    x: 計算する点 (float または numpy.ndarray)

    返り値:
        M個のsin関数の和
    """
    A = np.asarray(A)  # 念のため ndarray に変換
    if A.shape[1] != 3:
        raise ValueError("A は M×3 の行列でなければなりません。")

    amplitudes = A[:, 0]
    periods = A[:, 1]
    phases = A[:, 2]

    # M個の sin 関数を計算し、それらの合計を返す
    return np.sum(amplitudes[:, None] * np.sin(2 * np.pi * x / periods[:, None] + phases[:, None]), axis=0)


def generate_timeline(A, B, x):
    return sum_of_sigmoids(A, x) + sum_of_sin(B, x)

# 使用例
A1 = np.array([[1, 0.1, 100],
              [0.5, -1, 200],
              [2, 0.2, 300]])

A2 = np.array([[1, 0.1, 110],
              [0.5, -1, 210],
              [2, 0.2, 310]])

B1 = np.array([[0.1, 7, 0],
              [0.1, 31, np.pi/4],
              [0.1, 365, np.pi/2]])

B2 = np.array([[0.1, 7, 0],
              [0.1, 31, 0],
              [0.1, 365, 0]])

x = np.linspace(0, 365, 365)
y1 = generate_timeline(A1, B1, x)
y2 = generate_timeline(A2, B2, x)

data = pd.DataFrame({'X': y1, 'Y': y2})
max_lag = 30
granger_causality_test(data, max_lag, verbose=True)

次のようなエラーになりました。

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-7-0ff8b9ad1037> in <cell line: 0>()
     72 data = pd.DataFrame({'X': y1, 'Y': y2})
     73 max_lag = 30
---> 74 granger_causality_test(data, max_lag, verbose=True)

<ipython-input-2-5518c85d0771> in granger_causality_test(data, max_lag, verbose)
     52         k = lag
     53         F_stat = ((rss_restricted - rss_unrestricted) / k) / (rss_unrestricted / (n - 2*k - 1))
---> 54         p_value = 1 - sm.stats.f.sf(F_stat, k, n - 2*k - 1)
     55 
     56         results[lag] = {'F-stat': F_stat, 'p-value': p_value}

AttributeError: module 'statsmodels.stats.api' has no attribute 'f'

assistant:

エラーの原因は、sm.stats.f.sf の部分です。
statsmodels には f という属性は存在しないため、適切な scipy.stats.f.sf を使う必要があります。

修正点

次の部分：

p_value = 1 - sm.stats.f.sf(F_stat, k, n - 2*k - 1)

を、以下のように修正してください：

from scipy.stats import f

p_value = 1 - f.sf(F_stat, k, n - 2*k - 1)

修正後のコード

from scipy.stats import f

def granger_causality_test(data, max_lag, verbose=True):
    results = {}

    for lag in range(1, max_lag + 1):
        df = data.copy()

        # Create lagged variables
        for i in range(1, lag + 1):
            df[f'X_lag{i}'] = df.iloc[:, 0].shift(i)
            df[f'Y_lag{i}'] = df.iloc[:, 1].shift(i)

        df.dropna(inplace=True)

        Y = df.iloc[:, 1]

        # Restricted model: Y_t ~ Y_t-1, ..., Y_t-lag
        X_restricted = df[[f'Y_lag{i}' for i in range(1, lag + 1)]]
        X_restricted = sm.add_constant(X_restricted)
        model_restricted = OLS(Y, X_restricted).fit()
        rss_restricted = np.sum(model_restricted.resid ** 2)

        # Unrestricted model: Y_t ~ Y_t-1, ..., Y_t-lag + X_t-1, ..., X_t-lag
        X_unrestricted = df[[f'Y_lag{i}' for i in range(1, lag + 1)] + 
                             [f'X_lag{i}' for i in range(1, lag + 1)]]
        X_unrestricted = sm.add_constant(X_unrestricted)
        model_unrestricted = OLS(Y, X_unrestricted).fit()
        rss_unrestricted = np.sum(model_unrestricted.resid ** 2)

        # Compute F-statistic
        n = len(Y)
        k = lag
        F_stat = ((rss_restricted - rss_unrestricted) / k) / (rss_unrestricted / (n - 2*k - 1))
        p_value = 1 - f.sf(F_stat, k, n - 2*k - 1)  # 修正

        results[lag] = {'F-stat': F_stat, 'p-value': p_value}

        if verbose:
            print(f'Lag {lag}: F-stat={F_stat:.4f}, p-value={p_value:.4f}')

    return results

この修正を行えば、AttributeError は解消され、granger_causality_test 関数が正常に動作するはずです。

user:

ありがとうございます。次の結果が得られました。

Lag 1: F-stat=40.0606, p-value=1.0000
Lag 2: F-stat=6.2186, p-value=0.9978
Lag 3: F-stat=391.1485, p-value=1.0000
Lag 4: F-stat=80.8124, p-value=1.0000
Lag 5: F-stat=26.9606, p-value=1.0000
Lag 6: F-stat=26.5639, p-value=1.0000
Lag 7: F-stat=17.2936, p-value=1.0000
Lag 8: F-stat=14.9922, p-value=1.0000
Lag 9: F-stat=14.9578, p-value=1.0000
Lag 10: F-stat=16.7254, p-value=1.0000
Lag 11: F-stat=15.8893, p-value=1.0000
Lag 12: F-stat=14.5003, p-value=1.0000
Lag 13: F-stat=67.4585, p-value=1.0000
Lag 14: F-stat=65.0852, p-value=1.0000
Lag 15: F-stat=27792.1321, p-value=1.0000
Lag 16: F-stat=140769238.7171, p-value=1.0000
Lag 17: F-stat=5384829162.1713, p-value=1.0000
Lag 18: F-stat=331749790882.8946, p-value=1.0000
Lag 19: F-stat=1565.4609, p-value=1.0000
Lag 20: F-stat=-5.0478, p-value=0.0000
Lag 21: F-stat=-4.1613, p-value=0.0000
Lag 22: F-stat=323.5229, p-value=1.0000
Lag 23: F-stat=15.3345, p-value=1.0000
Lag 24: F-stat=17.6963, p-value=1.0000
Lag 25: F-stat=12.2656, p-value=1.0000
Lag 26: F-stat=104.4445, p-value=1.0000
Lag 27: F-stat=-3.0060, p-value=0.0000
Lag 28: F-stat=62.1631, p-value=1.0000
Lag 29: F-stat=1.6513, p-value=0.9780
Lag 30: F-stat=2048.3092, p-value=1.0000
Lag 31: F-stat=3.1309, p-value=1.0000
Lag 32: F-stat=47.5816, p-value=1.0000
Lag 33: F-stat=72.8319, p-value=1.0000
Lag 34: F-stat=848.5461, p-value=1.0000
Lag 35: F-stat=-2.8467, p-value=0.0000
Lag 36: F-stat=2.6321, p-value=1.0000
Lag 37: F-stat=1.2057, p-value=0.7964
Lag 38: F-stat=-3.5646, p-value=0.0000
Lag 39: F-stat=56.6989, p-value=1.0000
Lag 40: F-stat=32.6048, p-value=1.0000
Lag 41: F-stat=387.6863, p-value=1.0000
Lag 42: F-stat=4.2368, p-value=1.0000
Lag 43: F-stat=233.2896, p-value=1.0000
Lag 44: F-stat=138.0012, p-value=1.0000
Lag 45: F-stat=336.8925, p-value=1.0000
Lag 46: F-stat=-0.2686, p-value=0.0000
Lag 47: F-stat=2.6543, p-value=1.0000
Lag 48: F-stat=-3.4268, p-value=0.0000
Lag 49: F-stat=5.2814, p-value=1.0000
Lag 50: F-stat=49.9757, p-value=1.0000
<ipython-input-10-b57655b9688e>:16: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  df[f'Y_lag{i}'] = df.iloc[:, 1].shift(i)
{1: {'F-stat': np.float64(40.06058373728911),
  'p-value': np.float64(0.9999999992695425)},
 2: {'F-stat': np.float64(6.218611611429912),
  'p-value': np.float64(0.9977861674306957)},
 3: {'F-stat': np.float64(391.14852505351894), 'p-value': np.float64(1.0)},
 4: {'F-stat': np.float64(80.81244491498757), 'p-value': np.float64(1.0)},
 5: {'F-stat': np.float64(26.960575202299896), 'p-value': np.float64(1.0)},
 6: {'F-stat': np.float64(26.56388258652794), 'p-value': np.float64(1.0)},
 7: {'F-stat': np.float64(17.293553242017474), 'p-value': np.float64(1.0)},
 8: {'F-stat': np.float64(14.992200869018188), 'p-value': np.float64(1.0)},
 9: {'F-stat': np.float64(14.957755827916666), 'p-value': np.float64(1.0)},
 10: {'F-stat': np.float64(16.72536091543979), 'p-value': np.float64(1.0)},
 11: {'F-stat': np.float64(15.889288171018274), 'p-value': np.float64(1.0)},
 12: {'F-stat': np.float64(14.500285656469321), 'p-value': np.float64(1.0)},
 13: {'F-stat': np.float64(67.45850743207275), 'p-value': np.float64(1.0)},
 14: {'F-stat': np.float64(65.08520698708242), 'p-value': np.float64(1.0)},
 15: {'F-stat': np.float64(27792.132131850354), 'p-value': np.float64(1.0)},
 16: {'F-stat': np.float64(140769238.7170694), 'p-value': np.float64(1.0)},
 17: {'F-stat': np.float64(5384829162.17133), 'p-value': np.float64(1.0)},
 18: {'F-stat': np.float64(331749790882.8946), 'p-value': np.float64(1.0)},
 19: {'F-stat': np.float64(1565.4609449227123), 'p-value': np.float64(1.0)},
 20: {'F-stat': np.float64(-5.047849937902715), 'p-value': np.float64(0.0)},
 21: {'F-stat': np.float64(-4.16125584341578), 'p-value': np.float64(0.0)},
 22: {'F-stat': np.float64(323.5228517850055), 'p-value': np.float64(1.0)},
 23: {'F-stat': np.float64(15.334473176133748), 'p-value': np.float64(1.0)},
 24: {'F-stat': np.float64(17.696259285887386), 'p-value': np.float64(1.0)},
 25: {'F-stat': np.float64(12.265558690622997), 'p-value': np.float64(1.0)},
 26: {'F-stat': np.float64(104.4444621773218), 'p-value': np.float64(1.0)},
 27: {'F-stat': np.float64(-3.006033683736338), 'p-value': np.float64(0.0)},
 28: {'F-stat': np.float64(62.163080008426476), 'p-value': np.float64(1.0)},
 29: {'F-stat': np.float64(1.6512764581760477),
  'p-value': np.float64(0.9779727125128722)},
 30: {'F-stat': np.float64(2048.309169135248), 'p-value': np.float64(1.0)},
 31: {'F-stat': np.float64(3.130880523744094),
  'p-value': np.float64(0.9999997012830844)},
 32: {'F-stat': np.float64(47.58158435037054), 'p-value': np.float64(1.0)},
 33: {'F-stat': np.float64(72.83190185946394), 'p-value': np.float64(1.0)},
 34: {'F-stat': np.float64(848.5460924351703), 'p-value': np.float64(1.0)},
 35: {'F-stat': np.float64(-2.846691622944562), 'p-value': np.float64(0.0)},
 36: {'F-stat': np.float64(2.632076403623817),
  'p-value': np.float64(0.9999939299739421)},
 37: {'F-stat': np.float64(1.205671019430491),
  'p-value': np.float64(0.7964032819685993)},
 38: {'F-stat': np.float64(-3.5646055843735653), 'p-value': np.float64(0.0)},
 39: {'F-stat': np.float64(56.698872625205155), 'p-value': np.float64(1.0)},
 40: {'F-stat': np.float64(32.604787157931916), 'p-value': np.float64(1.0)},
 41: {'F-stat': np.float64(387.6862626222316), 'p-value': np.float64(1.0)},
 42: {'F-stat': np.float64(4.236810627155339),
  'p-value': np.float64(0.9999999999994623)},
 43: {'F-stat': np.float64(233.28962364163547), 'p-value': np.float64(1.0)},
 44: {'F-stat': np.float64(138.00120274758103), 'p-value': np.float64(1.0)},
 45: {'F-stat': np.float64(336.89248975126924), 'p-value': np.float64(1.0)},
 46: {'F-stat': np.float64(-0.2685792329381428), 'p-value': np.float64(0.0)},
 47: {'F-stat': np.float64(2.65427598137857),
  'p-value': np.float64(0.999999103472938)},
 48: {'F-stat': np.float64(-3.426819608538715), 'p-value': np.float64(0.0)},
 49: {'F-stat': np.float64(5.281358411022542), 'p-value': np.float64(1.0)},
 50: {'F-stat': np.float64(49.97567633985313), 'p-value': np.float64(1.0)}}

assistant:

結果を確認すると、いくつかの点が気になります。

気になる点

異常に大きな F-stat 値
- 例えば Lag 16: F-stat=140769238.7171 や Lag 18: F-stat=331749790882.8946 など、F統計量が極端に大きくなっています。
- これは何らかの数値計算の問題（スケーリングの問題、データの相関が異常に高いなど）が影響している可能性があります。
p-value が 1.0000 または 0.0000 に集中
- F-stat の計算結果が極端すぎるためか、ほとんどの p-value が 1.0 または 0.0 になっています。
- これは統計的に不自然であり、何らかのバグやデータの特性の影響を受けている可能性があります。
負の F-stat
- Lag 20: F-stat=-5.0478 など、F-stat が負になっているケースがあります。
- F統計量は通常、分散比を表すため負になることは理論的にありえません。
- これは rss_unrestricted や rss_restricted の計算に誤りがある可能性があります。

原因と考えられる対策

(1) データの確認

まず、入力データ data の統計情報を確認してみてください：

print(data.describe())

特に、X と Y の標準偏差や最小値・最大値が極端になっていないか確認してください。

(2) ラグの多重共線性

ラグを多く取りすぎると、多重共線性の影響で数値が不安定になることがあります。

max_lag=50 では多すぎる可能性があるため、 max_lag=10 など適切な範囲で試してみてください。

(3) RSS (残差平方和) の計算をチェック

負の F-stat の問題を解決するため、RSS の計算をデバッグすることをおすすめします。

以下のように、rss_restricted と rss_unrestricted の値をログ出力して、異常な値になっていないか確認してください：

print(f'Lag {lag}: RSS Restricted={rss_restricted}, RSS Unrestricted={rss_unrestricted}')

もし rss_unrestricted > rss_restricted の関係が崩れているなら、モデルのフィッティングに問題がある可能性があります。

次に試すべきこと

data.describe() でデータのスケールを確認
max_lag=10 くらいに減らして再試行
rss_restricted と rss_unrestricted のログ出力で異常値を調査
必要なら F_stat = max(0, F_stat) で負の値を防ぐ (ただし根本解決にはならない)

この辺りを試してみて、結果を教えてもらえますか？

user:

ありがとうございます。

y1 = (y1 - np.mean(y1)) / np.std(y1)
y2 = (y2 - np.mean(y2)) / np.std(y2)

を入れると、結果が次のように変わりました。

Lag 1: F-stat=40.0606, p-value=1.0000
Lag 2: F-stat=6.2186, p-value=0.9978
Lag 3: F-stat=391.1485, p-value=1.0000
Lag 4: F-stat=80.8124, p-value=1.0000
Lag 5: F-stat=26.9606, p-value=1.0000
Lag 6: F-stat=26.5639, p-value=1.0000
Lag 7: F-stat=17.2936, p-value=1.0000
Lag 8: F-stat=14.9922, p-value=1.0000
Lag 9: F-stat=14.9578, p-value=1.0000
Lag 10: F-stat=16.7254, p-value=1.0000
Lag 11: F-stat=15.8893, p-value=1.0000
Lag 12: F-stat=14.5003, p-value=1.0000
Lag 13: F-stat=67.4585, p-value=1.0000
Lag 14: F-stat=65.0852, p-value=1.0000
Lag 15: F-stat=27792.1321, p-value=1.0000
Lag 16: F-stat=140769238.7038, p-value=1.0000
Lag 17: F-stat=5386075220.8519, p-value=1.0000
Lag 18: F-stat=1026325558562962.1250, p-value=1.0000
Lag 19: F-stat=142619.3983, p-value=1.0000
Lag 20: F-stat=58956.5418, p-value=1.0000
Lag 21: F-stat=10550.2958, p-value=1.0000
Lag 22: F-stat=26.5232, p-value=1.0000
Lag 23: F-stat=2807.4149, p-value=1.0000
Lag 24: F-stat=840.0745, p-value=1.0000
Lag 25: F-stat=47.0092, p-value=1.0000
Lag 26: F-stat=3805.1295, p-value=1.0000
Lag 27: F-stat=121.9484, p-value=1.0000
Lag 28: F-stat=326.0130, p-value=1.0000
Lag 29: F-stat=404.3329, p-value=1.0000
Lag 30: F-stat=974.9442, p-value=1.0000
Lag 31: F-stat=68.0484, p-value=1.0000
Lag 32: F-stat=383.3183, p-value=1.0000
Lag 33: F-stat=168.0869, p-value=1.0000
Lag 34: F-stat=937.8194, p-value=1.0000
Lag 35: F-stat=57.5015, p-value=1.0000
Lag 36: F-stat=897.2447, p-value=1.0000
Lag 37: F-stat=352.1889, p-value=1.0000
Lag 38: F-stat=79.6217, p-value=1.0000
Lag 39: F-stat=35.9442, p-value=1.0000
Lag 40: F-stat=50.6164, p-value=1.0000
Lag 41: F-stat=580.4667, p-value=1.0000
Lag 42: F-stat=7.7125, p-value=1.0000
Lag 43: F-stat=10.3796, p-value=1.0000
Lag 44: F-stat=82.8213, p-value=1.0000
Lag 45: F-stat=414.3063, p-value=1.0000
Lag 46: F-stat=438.0053, p-value=1.0000
Lag 47: F-stat=642.6906, p-value=1.0000
Lag 48: F-stat=10.2997, p-value=1.0000
Lag 49: F-stat=191.0983, p-value=1.0000
Lag 50: F-stat=11.4681, p-value=1.0000
<ipython-input-10-b57655b9688e>:16: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  df[f'Y_lag{i}'] = df.iloc[:, 1].shift(i)
{1: {'F-stat': np.float64(40.06058373728919),
  'p-value': np.float64(0.9999999992695425)},
 2: {'F-stat': np.float64(6.21861161142987),
  'p-value': np.float64(0.9977861674306956)},
 3: {'F-stat': np.float64(391.14852505351405), 'p-value': np.float64(1.0)},
 4: {'F-stat': np.float64(80.81244491498437), 'p-value': np.float64(1.0)},
 5: {'F-stat': np.float64(26.960575202295477), 'p-value': np.float64(1.0)},
 6: {'F-stat': np.float64(26.563882586524734), 'p-value': np.float64(1.0)},
 7: {'F-stat': np.float64(17.29355324200711), 'p-value': np.float64(1.0)},
 8: {'F-stat': np.float64(14.992200869032736), 'p-value': np.float64(1.0)},
 9: {'F-stat': np.float64(14.957755827899712), 'p-value': np.float64(1.0)},
 10: {'F-stat': np.float64(16.725360915428617), 'p-value': np.float64(1.0)},
 11: {'F-stat': np.float64(15.889288171008136), 'p-value': np.float64(1.0)},
 12: {'F-stat': np.float64(14.500285656472457), 'p-value': np.float64(1.0)},
 13: {'F-stat': np.float64(67.45850743207374), 'p-value': np.float64(1.0)},
 14: {'F-stat': np.float64(65.08520698710566), 'p-value': np.float64(1.0)},
 15: {'F-stat': np.float64(27792.132131912324), 'p-value': np.float64(1.0)},
 16: {'F-stat': np.float64(140769238.7037663), 'p-value': np.float64(1.0)},
 17: {'F-stat': np.float64(5386075220.851897), 'p-value': np.float64(1.0)},
 18: {'F-stat': np.float64(1026325558562962.1), 'p-value': np.float64(1.0)},
 19: {'F-stat': np.float64(142619.39826526796), 'p-value': np.float64(1.0)},
 20: {'F-stat': np.float64(58956.54180730601), 'p-value': np.float64(1.0)},
 21: {'F-stat': np.float64(10550.29577200691), 'p-value': np.float64(1.0)},
 22: {'F-stat': np.float64(26.52324697317801), 'p-value': np.float64(1.0)},
 23: {'F-stat': np.float64(2807.4149383694566), 'p-value': np.float64(1.0)},
 24: {'F-stat': np.float64(840.0744980036021), 'p-value': np.float64(1.0)},
 25: {'F-stat': np.float64(47.009206214161786), 'p-value': np.float64(1.0)},
 26: {'F-stat': np.float64(3805.129477610578), 'p-value': np.float64(1.0)},
 27: {'F-stat': np.float64(121.9483671748811), 'p-value': np.float64(1.0)},
 28: {'F-stat': np.float64(326.01296173362203), 'p-value': np.float64(1.0)},
 29: {'F-stat': np.float64(404.3328901298167), 'p-value': np.float64(1.0)},
 30: {'F-stat': np.float64(974.944226101577), 'p-value': np.float64(1.0)},
 31: {'F-stat': np.float64(68.04838394628801), 'p-value': np.float64(1.0)},
 32: {'F-stat': np.float64(383.3183085186771), 'p-value': np.float64(1.0)},
 33: {'F-stat': np.float64(168.08687308205353), 'p-value': np.float64(1.0)},
 34: {'F-stat': np.float64(937.8194436244197), 'p-value': np.float64(1.0)},
 35: {'F-stat': np.float64(57.50146962704562), 'p-value': np.float64(1.0)},
 36: {'F-stat': np.float64(897.244652721118), 'p-value': np.float64(1.0)},
 37: {'F-stat': np.float64(352.18892655351505), 'p-value': np.float64(1.0)},
 38: {'F-stat': np.float64(79.62166561492785), 'p-value': np.float64(1.0)},
 39: {'F-stat': np.float64(35.944176873745064), 'p-value': np.float64(1.0)},
 40: {'F-stat': np.float64(50.61638758633128), 'p-value': np.float64(1.0)},
 41: {'F-stat': np.float64(580.4666576909693), 'p-value': np.float64(1.0)},
 42: {'F-stat': np.float64(7.7125052835917725), 'p-value': np.float64(1.0)},
 43: {'F-stat': np.float64(10.379567053237169), 'p-value': np.float64(1.0)},
 44: {'F-stat': np.float64(82.82134600526793), 'p-value': np.float64(1.0)},
 45: {'F-stat': np.float64(414.30631541356394), 'p-value': np.float64(1.0)},
 46: {'F-stat': np.float64(438.0053483072181), 'p-value': np.float64(1.0)},
 47: {'F-stat': np.float64(642.6905550697393), 'p-value': np.float64(1.0)},
 48: {'F-stat': np.float64(10.299716963825787), 'p-value': np.float64(1.0)},
 49: {'F-stat': np.float64(191.0983181457845), 'p-value': np.float64(1.0)},
 50: {'F-stat': np.float64(11.468071628778453), 'p-value': np.float64(1.0)}}

assistant:

標準化（Zスコア変換）を行った結果、F-statの値が極端に大きくなり、p-valueがすべて1.0になっています。このような結果になる原因として、以下の点が考えられます。

考えられる原因

標準化の影響
- 標準化後のデータのスケールが変わり、分散の変化が統計計算に影響を与えている可能性があります。
- 特にF統計量が極端に大きくなることは、分母（誤差分散）が非常に小さくなったためと考えられます。
分散の縮小（スケーリングエラー）
- 標準化によって分散が1に固定されるため、グレンジャー因果テストの計算に影響を及ぼしている可能性があります。
データの構造の変化
- 標準化前後でデータの自己相関構造が変わってしまい、適切な結果が得られなくなった可能性があります。

対策

標準化の代わりに差分系列を使用
- 時系列データのスケールの影響を除去する方法として「対数差分」や「一次差分」を試すのが有効です。
- 例: y1 = np.diff(y1) (または y1 = np.log(y1).diff())
標準化後にスケールを再調整
- 標準化を行った後に、元の分散を考慮した形に戻す。
別の因果探索手法を検討
- グレンジャー因果以外に、PCMCI+（時系列の因果探索）やLiNGAMといった手法も有効です。

一度、標準化せずに np.diff(y1) などを試して、再度グレンジャー因果テストを行ってみてはいかがでしょうか？

user:

ありがとうございます。こうなりました。

Lag 1: F-stat=0.4419, p-value=0.4934
Lag 2: F-stat=58.1738, p-value=1.0000
Lag 3: F-stat=113.9499, p-value=1.0000
Lag 4: F-stat=29.7000, p-value=1.0000
Lag 5: F-stat=30.5699, p-value=1.0000
Lag 6: F-stat=16.1384, p-value=1.0000
Lag 7: F-stat=15.8381, p-value=1.0000
Lag 8: F-stat=14.9387, p-value=1.0000
Lag 9: F-stat=17.2086, p-value=1.0000
Lag 10: F-stat=16.5551, p-value=1.0000
Lag 11: F-stat=15.0460, p-value=1.0000
Lag 12: F-stat=58.1217, p-value=1.0000
Lag 13: F-stat=56.9289, p-value=1.0000
Lag 14: F-stat=27925.7748, p-value=1.0000
Lag 15: F-stat=58788223.9535, p-value=1.0000
Lag 16: F-stat=66428845.6745, p-value=1.0000
Lag 17: F-stat=2175195376.7232, p-value=1.0000
Lag 18: F-stat=247709117482074944.0000, p-value=1.0000
Lag 19: F-stat=70180245.0597, p-value=1.0000
Lag 20: F-stat=55753291.5825, p-value=1.0000
Lag 21: F-stat=5990341.9511, p-value=1.0000
Lag 22: F-stat=369478.5841, p-value=1.0000
Lag 23: F-stat=47231.4921, p-value=1.0000
Lag 24: F-stat=78768.8270, p-value=1.0000
Lag 25: F-stat=34594.9203, p-value=1.0000
Lag 26: F-stat=12009.9203, p-value=1.0000
Lag 27: F-stat=6446.5397, p-value=1.0000
Lag 28: F-stat=72336.3307, p-value=1.0000
Lag 29: F-stat=4009.4216, p-value=1.0000
Lag 30: F-stat=13579.9389, p-value=1.0000
Lag 31: F-stat=4228.7320, p-value=1.0000
Lag 32: F-stat=8434.9794, p-value=1.0000
Lag 33: F-stat=12964.7404, p-value=1.0000
Lag 34: F-stat=6902.7875, p-value=1.0000
Lag 35: F-stat=4754.2714, p-value=1.0000
Lag 36: F-stat=2121.8779, p-value=1.0000
Lag 37: F-stat=1586.8502, p-value=1.0000
Lag 38: F-stat=3980.7686, p-value=1.0000
Lag 39: F-stat=2073.2140, p-value=1.0000
Lag 40: F-stat=1222.8925, p-value=1.0000
Lag 41: F-stat=4781.1033, p-value=1.0000
Lag 42: F-stat=1161.4193, p-value=1.0000
Lag 43: F-stat=247.6235, p-value=1.0000
Lag 44: F-stat=5836.7575, p-value=1.0000
Lag 45: F-stat=902.9757, p-value=1.0000
Lag 46: F-stat=1808.2229, p-value=1.0000
Lag 47: F-stat=2608.4918, p-value=1.0000
Lag 48: F-stat=1752.3705, p-value=1.0000
Lag 49: F-stat=887.0917, p-value=1.0000
Lag 50: F-stat=540.1859, p-value=1.0000
<ipython-input-10-b57655b9688e>:16: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead. To get a de-fragmented frame, use `newframe = frame.copy()`
  df[f'Y_lag{i}'] = df.iloc[:, 1].shift(i)
{1: {'F-stat': np.float64(0.44186724377802167),
  'p-value': np.float64(0.49335183277295636)},
 2: {'F-stat': np.float64(58.1738058960903), 'p-value': np.float64(1.0)},
 3: {'F-stat': np.float64(113.94994285576401), 'p-value': np.float64(1.0)},
 4: {'F-stat': np.float64(29.700000433659074), 'p-value': np.float64(1.0)},
 5: {'F-stat': np.float64(30.569927128563105), 'p-value': np.float64(1.0)},
 6: {'F-stat': np.float64(16.13836327340939),
  'p-value': np.float64(0.9999999999999998)},
 7: {'F-stat': np.float64(15.838106212469997), 'p-value': np.float64(1.0)},
 8: {'F-stat': np.float64(14.938705485070003), 'p-value': np.float64(1.0)},
 9: {'F-stat': np.float64(17.208638877621663), 'p-value': np.float64(1.0)},
 10: {'F-stat': np.float64(16.55507975986224), 'p-value': np.float64(1.0)},
 11: {'F-stat': np.float64(15.04600963919553), 'p-value': np.float64(1.0)},
 12: {'F-stat': np.float64(58.121725329016776), 'p-value': np.float64(1.0)},
 13: {'F-stat': np.float64(56.92894414363909), 'p-value': np.float64(1.0)},
 14: {'F-stat': np.float64(27925.77484135017), 'p-value': np.float64(1.0)},
 15: {'F-stat': np.float64(58788223.95349363), 'p-value': np.float64(1.0)},
 16: {'F-stat': np.float64(66428845.67453396), 'p-value': np.float64(1.0)},
 17: {'F-stat': np.float64(2175195376.7231984), 'p-value': np.float64(1.0)},
 18: {'F-stat': np.float64(2.4770911748207494e+17),
  'p-value': np.float64(1.0)},
 19: {'F-stat': np.float64(70180245.05971923), 'p-value': np.float64(1.0)},
 20: {'F-stat': np.float64(55753291.582456574), 'p-value': np.float64(1.0)},
 21: {'F-stat': np.float64(5990341.951070447), 'p-value': np.float64(1.0)},
 22: {'F-stat': np.float64(369478.58414511255), 'p-value': np.float64(1.0)},
 23: {'F-stat': np.float64(47231.492088521176), 'p-value': np.float64(1.0)},
 24: {'F-stat': np.float64(78768.82700606743), 'p-value': np.float64(1.0)},
 25: {'F-stat': np.float64(34594.920291934875), 'p-value': np.float64(1.0)},
 26: {'F-stat': np.float64(12009.920303466248), 'p-value': np.float64(1.0)},
 27: {'F-stat': np.float64(6446.539668614562), 'p-value': np.float64(1.0)},
 28: {'F-stat': np.float64(72336.33069007496), 'p-value': np.float64(1.0)},
 29: {'F-stat': np.float64(4009.421561729301), 'p-value': np.float64(1.0)},
 30: {'F-stat': np.float64(13579.938907517808), 'p-value': np.float64(1.0)},
 31: {'F-stat': np.float64(4228.731965296063), 'p-value': np.float64(1.0)},
 32: {'F-stat': np.float64(8434.97941463604), 'p-value': np.float64(1.0)},
 33: {'F-stat': np.float64(12964.740421739989), 'p-value': np.float64(1.0)},
 34: {'F-stat': np.float64(6902.787529768494), 'p-value': np.float64(1.0)},
 35: {'F-stat': np.float64(4754.271439535096), 'p-value': np.float64(1.0)},
 36: {'F-stat': np.float64(2121.8778776553804), 'p-value': np.float64(1.0)},
 37: {'F-stat': np.float64(1586.8502165807395), 'p-value': np.float64(1.0)},
 38: {'F-stat': np.float64(3980.76857906592), 'p-value': np.float64(1.0)},
 39: {'F-stat': np.float64(2073.213964993784), 'p-value': np.float64(1.0)},
 40: {'F-stat': np.float64(1222.89251008512), 'p-value': np.float64(1.0)},
 41: {'F-stat': np.float64(4781.103255417053), 'p-value': np.float64(1.0)},
 42: {'F-stat': np.float64(1161.4193005859959), 'p-value': np.float64(1.0)},
 43: {'F-stat': np.float64(247.62352649449224), 'p-value': np.float64(1.0)},
 44: {'F-stat': np.float64(5836.757538222318), 'p-value': np.float64(1.0)},
 45: {'F-stat': np.float64(902.9757063030826), 'p-value': np.float64(1.0)},
 46: {'F-stat': np.float64(1808.2228995625492), 'p-value': np.float64(1.0)},
 47: {'F-stat': np.float64(2608.491768909889), 'p-value': np.float64(1.0)},
 48: {'F-stat': np.float64(1752.370537710653), 'p-value': np.float64(1.0)},
 49: {'F-stat': np.float64(887.0916985260111), 'p-value': np.float64(1.0)},
 50: {'F-stat': np.float64(540.1859056393564), 'p-value': np.float64(1.0)}}

assistant:

F検定の結果を見ると、すべてのp値が1.0000になっていますね。これは統計的にかなり異常な結果です。

考えられる原因

データのスケーリングや前処理の問題
- 変数のスケールが極端に大きい、または小さい可能性がある。
- 標準化や対数変換が必要かもしれません。
多重共線性の問題
- ラグ変数間で強い相関があり、回帰分析がうまく機能していない可能性。
- Variance Inflation Factor（VIF）を計算して、多重共線性をチェックすると良い。
データのサンプルサイズ
- サンプルサイズが小さいと、統計的な信頼性が低くなる。
- 時系列データの場合、十分な観測数があるか確認。
計算の不安定性
- 例えば、データがすべて同じ値であるなど、分散がゼロに近いとF値が異常になることがある。

次のステップ

データの要約統計を確認（平均・標準偏差・最小値・最大値）
VIFを計算して多重共線性をチェック
データのスケーリングを試す（標準化、対数変換など）

どの変数を使ってF検定を行ったか、データのサンプルサイズや基本統計を教えてもらえれば、もう少し具体的なアドバイスができると思います！

user:

ありがとうございます。Variance Inflation Factor（VIF）を計算して、多重共線性をチェックする方法について教えてください。

assistant:

多重共線性をチェックするために、**VIF（Variance Inflation Factor, 分散膨張係数）**を計算する方法を説明します。

VIF（分散膨張係数）とは？

VIFは、多重共線性がどれくらいあるかを定量的に示す指標です。

( VIF = \frac{1}{1 - R^2} )
- ( R^2 ) は、対象の変数を他の説明変数で回帰したときの決定係数
VIFが10以上なら多重共線性の問題がある可能性が高い
VIFが5以上なら注意が必要

VIFを計算する方法（Python）

以下の手順でVIFを計算できます。

1. 必要なライブラリをインストール

import pandas as pd
import numpy as np
from statsmodels.stats.outliers_influence import variance_inflation_factor

2. 説明変数のデータを準備

# 説明変数のデータフレームを用意（目的変数は含めない）
X = df[['説明変数1', '説明変数2', '説明変数3', '説明変数4']]

3. VIFの計算

# 各変数のVIFを計算
vif_data = pd.DataFrame()
vif_data["変数名"] = X.columns
vif_data["VIF"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]

# 結果を表示
print(vif_data)

VIFの結果の解釈

VIF < 5：多重共線性の心配なし
5 ≦ VIF < 10：注意が必要
VIF ≥ 10：多重共線性の問題あり（変数の削除や工夫が必要）

多重共線性を解決する方法

もしVIFが高い変数がある場合、以下の対策を考えましょう。

相関の高い変数を削除
- 例えば、VIFが非常に高い変数を1つ削除することで問題を軽減できることがある。
主成分分析（PCA）を適用
- 似たような情報を持つ変数をまとめて、次元を減らす。
標準化（正規化）
- 変数のスケールを統一すると、多重共線性が改善される場合がある。
正則化回帰（Ridge回帰など）
- 多重共線性があっても回帰分析を安定させるための方法。

次のステップ

まず、あなたのデータに対してVIFを計算してみてください。
もしVIFが高い変数があれば、その変数の削除や変換を検討しましょう。
結果を教えてくれたら、次の対策を一緒に考えます！

user:

なるほど、ありがとうございます。時系列の変数2つに対して、Granger因果性を計算しようと思ったら、多重共線性をチェックしようという流れなのですが、それについては次回にしますね。今回の議論をブログにまとめたいのですが、タイトル案を複数提案していただけますか？

assistant:

それはいいですね！今回の議論を踏まえて、以下のようなタイトル案を考えました。

シンプル＆分かりやすい系

Granger因果性と多重共線性の関係を考える
時系列分析における多重共線性チェックの重要性
Granger因果性を検証する前にやるべきこと：VIFで多重共線性を確認しよう
多重共線性がGranger因果性の検証に与える影響とは？

キャッチー＆興味を引く系

時系列データで「因果関係」を調べる前に、落とし穴をチェックしよう！
VIFは要チェック！Granger因果性を正しく検証するための前準備
多重共線性を無視すると危険？Granger因果性分析の前に知っておくべきこと

専門性を感じさせる系

Granger因果性の検証における多重共線性の影響と対策
時系列データ分析におけるVIF計算の意義：Granger因果性を正しく評価するために
多重共線性とGranger因果性の関係：時系列データ分析のベストプラクティス

どれか気に入ったものはありますか？
あるいは、もっと特定のニュアンスを出したい場合は調整できますよ！

user:

いえ、多重共線性は最後の方に問題になっただけで、メインテーマではないです。Granger因果性の計算を、statsmodels.tsa.stattools.grangercausalitytests を使わずに実装するというのが今回のメインテーマです。

assistant:

なるほど！では、「Granger因果性の計算を statsmodels.tsa.stattools.grangercausalitytests を使わずに実装する」 ことをメインテーマにしたタイトル案を考えますね。

シンプル＆分かりやすい系

Granger因果性を自作する：ライブラリなしで実装する方法
statsmodels なしでGranger因果性を計算する方法
Granger因果性の手動計算に挑戦！統計モデルを自作してみよう
時系列データでGranger因果性を検証する：ライブラリに頼らないアプローチ

キャッチー＆興味を引く系

ブラックボックスに頼らない！Granger因果性をゼロから実装する方法
statsmodels の裏側を探る：Granger因果性をPythonで手動計算してみた
時系列データ分析の裏側：Granger因果性をスクラッチ実装して理解する

専門性を感じさせる系

Granger因果性の統計的検定をPythonでスクラッチ実装する
Granger因果性のF検定を手動で実装する：statsmodels の内部ロジックを理解する
Granger因果性の数理的背景と実装：統計ライブラリを使わずに検証する方法

どれか気に入るものはありますか？また、もう少し調整したいニュアンスがあれば教えてください！

user:

ありがとうございます。「Granger因果性の統計的検定をPythonでスクラッチ実装する」にしますね。最後に、今回のまとめを書いてください。

assistant:

了解しました！以下に、今回の議論をまとめます。

Granger因果性の統計的検定をPythonでスクラッチ実装する

1. Granger因果性とは？

Granger因果性とは、過去のある変数 $X_t$ の情報が、もう一つの変数 $Y_t$ の予測精度を向上させるかどうかを統計的に検証する手法です。単なる相関とは異なり、時間的な因果関係（予測可能性）を測る点が特徴です。

2. `statsmodels` を使わずにGranger因果性を検定する理由

通常、Pythonでは statsmodels.tsa.stattools.grangercausalitytests を使えば簡単にGranger因果性を検定できます。しかし、今回は以下の理由でスクラッチ実装を行いました。

統計的なロジックを理解するため
ブラックボックス化を避けるため
カスタマイズ性を高めるため

3. Granger因果性検定の実装手順

(1) データの前処理

時系列データを整え、適切なラグ（遅れ）を考慮する。
各ラグのデータを作成する（例：$Y_t$ を $Y_{t-1}, Y_{t-2}, ..., Y_{t-k}$ のように展開）。

(2) 回帰モデルの構築

帰無仮説（$X_t$ は $Y_t$ の予測に寄与しない）を検証する
2つの回帰モデルを構築
1. 制約付きモデル（$Y_t$ を過去の $Y$ のみで回帰）
2. フルモデル（$Y_t$ を過去の $Y$ と $X$ で回帰）

(3) F検定を実施

F統計量の計算
$$
F = \frac{(RSS_r - RSS_u) / k}{RSS_u / (n - 2k - 1)}
$$
- $RSS_r$：制約付きモデルの残差平方和（Restricted Residual Sum of Squares）
- $RSS_u$：フルモデルの残差平方和（Unrestricted Residual Sum of Squares）
- $k$：ラグの数
- $n$：データのサンプル数
p値を計算し、統計的有意性を判断する
- p値が小さい場合（一般的に 0.05 未満）、帰無仮説を棄却し、「$X_t$ が $Y_t$ のGranger因果性を持つ」と結論付ける。

4. 実験結果と考察

実際にPythonで実装し、複数のラグを試した結果を確認。
多重共線性の問題が発生しやすい点も考慮し、今後の解析ではVIF（分散膨張係数）なども検討すべき。

5. まとめと今後の課題

Granger因果性の検定をスクラッチ実装することで、単なるライブラリの使用では見えにくい統計モデルの裏側を理解することができました。特に、F検定の役割や、回帰モデルの設定が因果関係の検定にどのように関わるかを深掘りできました。

今後の課題として、

ラグの選択方法の最適化（AIC/BICを用いたモデル選択など）
多重共線性の影響の評価（VIFや主成分回帰の導入）
非線形因果性の検討（ニューラルネットやVARモデルの拡張）
が挙げられます。

今回の実装はGranger因果性の基本的な理解を深める良い機会となりましたが、より実用的なモデルへ発展させる余地は大いにあります。今後も改良を重ね、より精度の高い因果推論を目指していきましょう！

こんな感じでどうでしょうか？
もし補足したい点や強調したいポイントがあれば教えてください！

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up

Granger因果性の統計的検定をPythonでスクラッチ実装する

user:

assistant:

user:

assistant:

1. Granger因果性検定の概要

2. granger_causality_test の詳細説明

(1) 関数の定義

(2) 結果を格納する辞書を初期化

(3) ラグごとのループ

(4) 遅延変数を作成

(5) 従属変数 Y_t を設定

(6) 制約モデル (Restricted Model) の構築

(7) 非制約モデル (Unrestricted Model) の構築

(8) F検定の計算

(9) 結果の格納

(10) 結果を返す

3. まとめ

user:

assistant:

修正点

修正後のコード

user:

assistant:

気になる点

原因と考えられる対策

(1) データの確認

(2) ラグの多重共線性

(3) RSS (残差平方和) の計算をチェック

次に試すべきこと

user:

assistant:

考えられる原因

対策

user:

assistant:

考えられる原因

次のステップ

user:

assistant:

VIF（分散膨張係数）とは？

VIFを計算する方法（Python）

1. 必要なライブラリをインストール

2. 説明変数のデータを準備

3. VIFの計算

VIFの結果の解釈

多重共線性を解決する方法

次のステップ

user:

assistant:

シンプル＆分かりやすい系

キャッチー＆興味を引く系

専門性を感じさせる系

user:

assistant:

シンプル＆分かりやすい系

キャッチー＆興味を引く系

専門性を感じさせる系

user:

assistant:

Granger因果性の統計的検定をPythonでスクラッチ実装する

1. Granger因果性とは？

2. statsmodels を使わずにGranger因果性を検定する理由

3. Granger因果性検定の実装手順

(1) データの前処理

(2) 回帰モデルの構築

(3) F検定を実施

4. 実験結果と考察

5. まとめと今後の課題

2. `granger_causality_test` の詳細説明

2. `statsmodels` を使わずにGranger因果性を検定する理由