Jarque-Bera 검정¶

(1) 맞지 않는다. \(49.9516\)이 나온다. 출력된 값은 bias=False로 구한 보정판 \(G_1 = 0.37074\), \(G_2 = 1.85650\)인데, scipy.stats.jarque_bera는 비보정판 \(g_1\), \(g_2\)로 통계량을 만든다. 둘을 나란히 놓으면

이다. 아래쪽 합이 \(49.95\)로 scipy가 돌려준 \(47.5606\)보다 \(5.0\%\) 크다. \(n = 300\)이나 되는데도 그만큼 어긋난다. 왜도의 보정 계수는

로 \(0.5\%\)에 지나지 않지만, 첨도 쪽 보정은 \(1.85650/1.80578 = 1.02809\)로 더 크고 공식이 이들을 제곱하므로 차이가 두 배로 증폭된다.

그러므로 코드가 찍는 Skewness g1·Excess kurtosis g2라는 이름은 사실 \(G_1\), \(G_2\)이고, 그 값으로는 같은 출력에 찍힌 \(JB\)를 재현할 수 없다. 코드의 주석이 "둘 다 scipy 의 기본값과 다르니 다른 책의 숫자와 견줄 때 확인해야 할 대목이다"라고 경고하는 것이 바로 이 일이다. 손으로 \(JB\)를 확인하려면 bias 인자를 떼고 기본값으로 다시 불러야 한다.

(2) 첨도 항이 \(85.7\%\)다. 비보정 판본으로

이고 합이 정확히 \(47.5606\)이다. 비율은 왜도 \(14.3\%\), 첨도 \(85.7\%\)다. \(g_2 = +1.81\)로 부호가 양수이니 꼬리가 두꺼운 쪽이며, 대수정규 오염이 들여온 몇 개의 큰 값이 네제곱 적률을 지배한 것이다. 치우침(\(g_1 = 0.37\))도 있지만 기각을 끌고 온 주범은 아니다.

import numpy as np
from scipy import stats

# 정규 240개에 로그정규 60개를 섞은 자료다.
rng = np.random.default_rng(0)
x = np.concatenate([rng.normal(0, 1, size=240),
                    rng.lognormal(0, 0.6, size=60)])

jb_stat, p = stats.jarque_bera(x)

# bias=False 는 유한표본 보정을 적용한다는 뜻이다. fisher=True 는 첨도에서
# 3 을 빼 정규를 0 으로 맞춘 초과첨도를 준다. 둘 다 scipy 의 기본값과 다르니
# 다른 책의 숫자와 견줄 때 확인해야 할 대목이다.
g1 = stats.skew(x, bias=False)
g2 = stats.kurtosis(x, fisher=True, bias=False)

print(f"Sample size n = {x.size}")
print(f"Skewness g1 = {g1:.4f}")
print(f"Excess kurtosis g2 = {g2:.4f}")
print(f"Jarque-Bera: JB = {jb_stat:.4f}, p-value = {p:.4g}")
if p < 0.05:
    print("=> Reject normality at alpha = 0.05.")
else:
    print("=> Fail to reject normality at alpha = 0.05.")

출력:

Sample size n = 300
Skewness g1 = 0.3707
Excess kurtosis g2 = 1.8565
Jarque-Bera: JB = 47.5606, p-value = 4.703e-11
=> Reject normality at alpha = 0.05.

\(p\)값은 \(e^{-JB/2} = e^{-23.78} = 4.703\times10^{-11}\)로 사실상 0이다. 같은 자료에 다고스티노 \(K^2\)를 걸면 \(K^2 = 22.21\)(\(p = 1.50\times10^{-5}\))로 훨씬 작은데, 이는 검정력의 차이가 아니라 \(K^2\)가 \(g_2\)를 비선형 변환으로 먼저 길들이기 때문이다. 유한표본에서 \(g_2\)의 분포가 오른쪽으로 크게 치우쳐 있어 \(JB\)가 쉽게 부풀려진다.

출력의 두 수에 속지 말 것. 아래 "해석"에서 다시 쓰는 분해 \(6.80 + 40.76\)은 비보정 \(g_1 = 0.3688\), \(g_2 = 1.8058\)로 계산한 것이고, 그 합만이 \(47.5606\)과 끝자리까지 일치한다. \(\square\)

import numpy as np
from scipy import stats

rng = np.random.default_rng(42)
x = rng.normal(0, 1, size=500)

jb, p = stats.jarque_bera(x)
print(f"JB = {jb:.4f}, p = {p:.4f}")

출력:

JB = 0.7827, p = 0.6762

\(\text{JB} = 0.783\)은 \(\chi^2_2\)의 95백분위수 \(5.99\)보다 훨씬 작고 \(p = 0.676\)이므로 기각하지 않는다. 자료가 참으로 정규이므로 기대한 결과이다.

\(\chi^2_2\)의 평균이 2이므로 \(\text{JB} = 0.783\)은 오히려 평균보다 작은 값이다. 귀무가설 아래에서 흔히 일어나는 일이다. \(\square\)

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
x = np.concatenate([rng.normal(0, 1, 240), rng.lognormal(0, 0.6, 60)])

g1 = stats.skew(x, bias=True)
g2 = stats.kurtosis(x, fisher=True, bias=True)
n = x.size
jb_manual = (n / 6) * (g1**2 + g2**2 / 4)

jb_scipy, _ = stats.jarque_bera(x)

print(f"g1 (biased) = {g1:.4f}, g2 (biased) = {g2:.4f}")
print(f"Manual JB:  {jb_manual:.4f}")
print(f"SciPy JB:   {jb_scipy:.4f}")

출력:

g1 (biased) = 0.3688, g2 (biased) = 1.8058
Manual JB:  47.5606
SciPy JB:   47.5606

수치 정밀도까지 정확히 일치한다.

중요한 함정: stats.jarque_bera는 왜도와 첨도의 편향 추정값(bias=True)을 쓴다. 손으로 계산할 때도 반드시 같게 맞춰야 한다. 편향 보정판(\(g_1 = 0.3707\), \(g_2 = 1.8565\))을 쓰면 \(\text{JB} = 49.95\)가 나와 SciPy 값 \(47.56\)과 어긋난다. 차이가 5%다.

\(n\)이 크면 두 추정값의 차이가 사라지지만, 작은 표본에서는 어느 쪽을 쓰는지가 실질적으로 중요하다. \(\square\)

JB 통계량은 점근 분산 \(\text{Var}(g_1) \approx 6/n\)과 \(\text{Var}(g_2) \approx 24/n\)에서 유도된다. 작은 \(n\)에서는 이 근사가 부정확하다. 참된 유한표본 분산이 다르고, \(g_1\)과 \(g_2\)의 분포가 정규로 잘 근사되지 않는다. 게다가 작은 \(n\)에서 \(g_1\)과 \(g_2\)가 독립이 아닐 수 있다.

특히 \(g_2\)의 유한표본 분포는 오른쪽으로 심하게 치우쳐 있다. 정규성 아래에서도 \(g_2\)는 아래로는 \(-2\) 정도까지밖에 갈 수 없지만 위로는 얼마든지 커질 수 있다. 그 결과 \(g_2^2\)의 분포가 \(\chi^2_1\)과 크게 다르다.

결과적으로 \(\chi^2_2\) 기준분포가 잘 맞지 않는다. JB의 실제 귀무분포는 \(\chi^2_2\)보다 작은 값 쪽에 더 몰려 있어서 명목 \(\alpha = 0.05\) 임계값(\(\approx 5.99\))이 너무 크다. 검정이 과소기각하며(경험적 크기 \(< 0.05\)) 그만큼 검정력이 낮아진다. 연습문제 5에서 수치로 확인한다.

작은 표본에는 D'Agostino \(K^2\) 검정이나 Shapiro-Wilk 검정이 선호된다. \(\square\)

점근 분산에서 \(\text{Var}(g_1) \to 6/n\), \(\text{Var}(g_2) \to 24/n\)이다. 비율은

JB 공식에서 각 적률을 표준화하면

따라서 \(g_2^2\)에 붙은 인자 \(1/4\)은 \(g_2\)의 분산이 \(g_1\)의 네 배라는 사실을 보정하며, 두 표준화 성분이 모두 단위분산을 갖고 그 합이 \(\chi^2_2\)을 따르도록 보장한다.

실무적 함의: 이 가중치 때문에 JB는 같은 크기의 왜도보다 첨도에 덜 민감하다. \(|g_1| = |g_2| = c\)이면 왜도 항이 첨도 항의 네 배를 기여한다. 그러나 실제 자료에서는 비정규성이 흔히 첨도 쪽에서 더 큰 값으로 나타나므로(본문 보기에서 \(g_2 = 1.86\) 대 \(g_1 = 0.37\)), 결국 첨도 항이 지배하는 경우가 많다. \(\square\)

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

rng = np.random.default_rng(0)
reps, alpha = 10000, 0.05
ns = [30, 100, 500, 2000]
sizes = []

for n in ns:
    rej = sum(1 for _ in range(reps)
              if stats.jarque_bera(rng.normal(0, 1, n))[1] < alpha)
    sizes.append(rej / reps)
    print(f"n = {n:>4}: empirical size = {rej/reps:.4f}")

fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(ns, sizes, marker="o")
ax.axhline(0.05, color="red", linestyle="--", label="Nominal 0.05")
ax.set_xscale("log")
ax.set_xlabel("Sample size n (log scale)")
ax.set_ylabel("Empirical rejection rate")
ax.set_title("JB Test: Empirical Size vs n")
ax.legend()
plt.tight_layout()
plt.show()

출력:

n =   30: empirical size = 0.0319
n =  100: empirical size = 0.0470
n =  500: empirical size = 0.0467
n = 2000: empirical size = 0.0461

\(n = 30\)에서 경험적 크기 \(0.0319\)는 명목값의 64%에 불과하다. 연습문제 3에서 설명한 과소기각을 확인해 준다. 실제 유의수준이 0.032라면 검정력도 그만큼 떨어진다.

\(n = 100\)에서 \(0.047\)로 크게 개선되며 사실상 여기서 안정된다.

다만 눈여겨볼 점이 있다. \(n = 100\)에서 \(n = 2000\)으로 가도 더 나아지지 않는다(\(0.0470 \to 0.0467 \to 0.0461\)). 여전히 명목값보다 살짝 낮은 수준에 머문다. 몬테카를로 표준오차가 \(\sqrt{0.05 \times 0.95/10000} = 0.0022\)이므로 \(0.046\)과 \(0.05\)의 차이는 약 1.8 표준오차로 경계선상이지만, 세 표본크기에서 일관되게 같은 방향으로 나타난다는 점이 우연이 아님을 시사한다.

이는 JB 검정의 잘 알려진 성질이다. \(\chi^2_2\) 근사로의 수렴이 매우 느리다. \(g_2\)의 유한표본 분포가 치우쳐 있어 아주 큰 표본에서도 잔여 왜곡이 남는다. 실무적 결론: \(n \geq 100\)이면 JB를 쓸 만하지만, 크기 정확도가 결정적이면 모의실험 임계값을 쓰라. \(\square\)