금융 변동성 비교¶

import numpy as np
from scipy import stats

rng = np.random.default_rng(42)

# t(5) 의 표준편차는 sqrt(5/3) 이다. 목표 변동성을 정확히 맞추려면
# 이 값으로 나눠야 한다. 이러면 달라지는 것은 꼬리 두께뿐이다.
scale_correction = np.sqrt(5 / 3)

# 1구간: 일변동성 1.8%
returns_before = rng.standard_t(df=5, size=60) * 0.018 / scale_correction

# 2구간: 일변동성 2.4%. 실제로 분산이 달라진 자료다.
returns_after = rng.standard_t(df=5, size=60) * 0.024 / scale_correction

print(f"sample sd: {returns_before.std(ddof=1):.5f}, "
      f"{returns_after.std(ddof=1):.5f}")

# 중앙값을 중심으로 쓰므로 꼬리가 두꺼운 자료에서도 버틴다.
bf_stat, bf_p = stats.levene(returns_before, returns_after, center='median')
print(f"Brown-Forsythe:  W = {bf_stat:.4f}, p = {bf_p:.4f}")

# 순위만 쓰는 검정이라 셋 중 가장 로버스트하다.
fk_stat, fk_p = stats.fligner(returns_before, returns_after)
print(f"Fligner-Killeen: H = {fk_stat:.4f}, p = {fk_p:.4f}")

# F 검정은 정규성에 매우 민감하다. 수익률처럼 꼬리가 두꺼운 자료에서는
# 분산이 같아도 자주 기각해 버리므로 쓰지 않는 편이 낫다.
f_stat = np.var(returns_before, ddof=1) / np.var(returns_after, ddof=1)
f_p = 2 * min(stats.f.cdf(f_stat, 59, 59), stats.f.sf(f_stat, 59, 59))
print(f"F-test:          F = {f_stat:.4f}, p = {f_p:.4g} "
      f"(unreliable for heavy-tailed data)")

sample sd: 0.01504, 0.02611
Brown-Forsythe:  W = 5.4046, p = 0.0218
Fligner-Killeen: H = 3.2325, p = 0.0722
F-test:          F = 0.3317, p = 3.791e-05 (unreliable for heavy-tailed data)

import numpy as np
from scipy import stats

# (1) 극한 크기를 식으로 먼저 적는다. t_5 의 첨도는 3 + 6/(nu-4) = 9 다.
beta2 = 3 + 6 / (5 - 4)
c = stats.norm.isf(0.025) * np.sqrt(2 / (beta2 - 1))
print(f"t_5 의 첨도 beta2 = {beta2:.1f}   (초과첨도 {beta2 - 3:.1f})")
print(f"산포 배율 sqrt((beta2-1)/2) = {np.sqrt((beta2 - 1) / 2):.3f}")
print(f"1.96*sqrt(2/(beta2-1)) = {c:.4f}   ->  극한 크기 = {2 * stats.norm.sf(c):.4f}")


def bf_p2(A, B):
    """두 표본에 대한 중앙값 중심 Levene(= Brown-Forsythe) p-값을 벡터화.
    scipy.stats.levene(a, b, center='median') 와 같은 값을 준다."""
    Z = [np.abs(X - np.median(X, axis=1, keepdims=True)) for X in (A, B)]
    n1, n2 = A.shape[1], B.shape[1]
    m = [z.mean(axis=1) for z in Z]
    gm = (n1 * m[0] + n2 * m[1]) / (n1 + n2)
    msb = n1 * (m[0] - gm) ** 2 + n2 * (m[1] - gm) ** 2
    msw = (((Z[0] - m[0][:, None]) ** 2).sum(axis=1)
           + ((Z[1] - m[1][:, None]) ** 2).sum(axis=1)) / (n1 + n2 - 2)
    return stats.f.sf(msb / msw, 1, n1 + n2 - 2)


# (2) 표본을 키우며 두 검정의 실제 크기를 잰다. 참 변동성은 두 기간이 같다.
rng = np.random.default_rng(11)
print(f"\n{'n':>6}{'R':>8}{'F-test':>10}{'Brown-Forsythe':>17}")
for n, R in [(60, 20_000), (250, 20_000), (1000, 10_000), (4000, 4_000)]:
    A = rng.standard_t(5, size=(R, n))
    B = rng.standard_t(5, size=(R, n))
    F = A.var(axis=1, ddof=1) / B.var(axis=1, ddof=1)
    d = stats.f(n - 1, n - 1)
    p_f = 2 * np.minimum(d.cdf(F), d.sf(F))
    print(f"{n:>6}{R:>8}{np.mean(p_f < 0.05):>10.4f}{np.mean(bf_p2(A, B) < 0.05):>17.4f}")

t_5 의 첨도 beta2 = 9.0   (초과첨도 6.0)
산포 배율 sqrt((beta2-1)/2) = 2.000
1.96*sqrt(2/(beta2-1)) = 0.9800   ->  극한 크기 = 0.3271

     n       R    F-test   Brown-Forsythe
    60   20000    0.2028           0.0474
   250   20000    0.2399           0.0481
  1000   10000    0.2696           0.0499
  4000    4000    0.2755           0.0493

import numpy as np
from scipy import stats

rng = np.random.default_rng(3)
R = 5000
counts = [0, 0, 0]

for _ in range(R):
    a = rng.standard_t(5, 60)
    b = rng.standard_t(5, 60)
    f = a.var(ddof=1) / b.var(ddof=1)
    counts[0] += 2 * min(stats.f.cdf(f, 59, 59),
                         stats.f.sf(f, 59, 59)) < 0.05
    counts[1] += stats.levene(a, b, center='median')[1] < 0.05
    counts[2] += stats.fligner(a, b)[1] < 0.05

print(f"F-test:          {counts[0] / R:.4f}")
print(f"Brown-Forsythe:  {counts[1] / R:.4f}")
print(f"Fligner-Killeen: {counts[2] / R:.4f}")

F-test:          0.1948
Brown-Forsythe:  0.0438
Fligner-Killeen: 0.0448

import numpy as np
from scipy import stats

def garch11(n, rng, omega=1e-5, alpha=0.1, beta=0.85):
    """GARCH(1,1) 수익률을 만든다. 무조건분산은 omega/(1-alpha-beta) 다.

        변동성이 뭉쳐 다니는 실제 수익률의 성질을 흉내 내기 위한 모형이다.
        """
    h = omega / (1 - alpha - beta)
    r = np.zeros(n)
    for t in range(n):
        r[t] = np.sqrt(h) * rng.normal()
        h = omega + alpha * r[t] ** 2 + beta * h
    return r

rng = np.random.default_rng(3)
R = 5000
counts = [0, 0, 0]

for _ in range(R):
    x = garch11(120, rng)
    a, b = x[:60], x[60:]           # two adjacent periods, same true variance
    f = a.var(ddof=1) / b.var(ddof=1)
    counts[0] += 2 * min(stats.f.cdf(f, 59, 59),
                         stats.f.sf(f, 59, 59)) < 0.05
    counts[1] += stats.levene(a, b, center='median')[1] < 0.05
    counts[2] += stats.fligner(a, b)[1] < 0.05

print(f"F-test:          {counts[0] / R:.4f}")
print(f"Brown-Forsythe:  {counts[1] / R:.4f}")
print(f"Fligner-Killeen: {counts[2] / R:.4f}")

F-test:          0.3312
Brown-Forsythe:  0.2968
Fligner-Killeen: 0.2868

import numpy as np

print(f"{'nu':>5} {'sd(t_nu)':>10} {'sd from *0.018':>15}")
for nu in [3, 4, 5, 6, 10, 30, 100]:
    sd = np.sqrt(nu / (nu - 2))
    print(f"{nu:>5} {sd:>10.4f} {0.018 * sd:>15.5f}")

   nu   sd(t_nu)  sd from *0.018
    3     1.7321         0.03118
    4     1.4142         0.02546
    5     1.2910         0.02324
    6     1.2247         0.02205
   10     1.1180         0.02012
   30     1.0351         0.01863
  100     1.0102         0.01818

rng = np.random.default_rng(0)
nu, n, sigma = 5, 100_000, 0.018

naive = rng.standard_t(nu, n) * sigma
correct = rng.standard_t(nu, n) * sigma / np.sqrt(nu / (nu - 2))

print(f"target sd:  {sigma:.5f}")
print(f"naive sd:   {naive.std(ddof=1):.5f}")
print(f"correct sd: {correct.std(ddof=1):.5f}")

target sd:  0.01800
naive sd:   0.02316
correct sd: 0.01805