금융 자료의 Q-Q 그림: 자산 수익률의 비정규성 탐지¶

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

# 넷플릭스 일별 로그수익률을 흉내 낸 자료다. 자유도 8 인 t 를 썼는데,
# 실제 수익률은 이보다도 꼬리가 두껍다.
np.random.seed(42)
returns = np.random.standard_t(df=8, size=1000) * 0.02

# 정규분포와 견주는 Q-Q 그림
fig, ax = plt.subplots(figsize=(8, 6))
stats.probplot(returns, dist="norm", plot=ax)

ax.set_title("Q-Q Plot: Daily Log-Returns vs Normal Distribution", fontsize=12)
ax.set_xlabel("Theoretical Normal Quantiles", fontsize=11)
ax.set_ylabel("Sample Quantiles (Observed Returns)", fontsize=11)

ax.spines[["top", "right"]].set_visible(False)
ax.grid(True, alpha=0.3, linestyle='--')

plt.tight_layout()
plt.show()

osm, osr = stats.probplot(returns, dist="norm", fit=False)
b, a = np.polyfit(osm, osr, 1)
resid = osr - (a + b * osm)

print(f"절편 a = {a:.6f}   표본평균 = {returns.mean():.6f}")
print(f"기울기 b = {b:.6f}   표본표준편차 = {returns.std(ddof=1):.6f}")
print(f"이론 표준편차 0.02*sqrt(4/3) = {0.02 * np.sqrt(4 / 3):.6f}")

for lo, hi, lab in [(0, 1, "|q| < 1"), (2, 99, "|q| > 2")]:
    m = (np.abs(osm) > lo) & (np.abs(osm) < hi)
    print(f"{lab}: {m.sum()}점, 최대 절대잔차 {np.abs(resid[m]).max():.5f}")

print(f"최솟값 {osr[0]:.4f}  적합선 {a + b * osm[0]:.4f}  잔차 {resid[0]:.4f}")
print(f"최댓값 {osr[-1]:.4f}  적합선 {a + b * osm[-1]:.4f}  잔차 {resid[-1]:.4f}")
print(f"g1 = {stats.skew(returns):.4f}   G1 = {stats.skew(returns, bias=False):.4f}")
print(f"g2 = {stats.kurtosis(returns):.4f}   G2 = "
      f"{stats.kurtosis(returns, bias=False):.4f}   이론 = {6 / (8 - 4):.4f}")

z = (returns - returns.mean()) / returns.std(ddof=1)
for k in (3, 4):
    print(f"|z| > {k}: 관측 {np.sum(np.abs(z) > k)}개, "
          f"정규 기대 {1000 * 2 * stats.norm.sf(k):.3f}개")

절편 a = 0.000376   표본평균 = 0.000376
기울기 b = 0.022470   표본표준편차 = 0.022574
이론 표준편차 0.02*sqrt(4/3) = 0.023094
|q| < 1: 682점, 최대 절대잔차 0.00191
|q| > 2: 46점, 최대 절대잔차 0.06371
최솟값 -0.0774  적합선 -0.0715  잔차 -0.0060
최댓값 0.1359  적합선 0.0722  잔차 0.0637
g1 = 0.2257   G1 = 0.2261
g2 = 1.6429   G2 = 1.6572   이론 = 1.5000
|z| > 3: 관측 5개, 정규 기대 2.700개
|z| > 4: 관측 1개, 정규 기대 0.063개

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

# 히스토그램과 Q-Q 그림을 나란히 놓아, 같은 자료에서 둘이 무엇을 보여
# 주는지 견준다. 히스토그램은 가운데를, Q-Q 그림은 꼬리를 말한다.
np.random.seed(42)
df = 5  # 자유도가 작을수록 꼬리가 두껍다
heavy_tailed_returns = stats.t.rvs(df=df, scale=0.02, size=2000)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# 왼쪽: 히스토그램. 가운데가 정규보다 뾰족하지만, 꼬리의 차이는
# 자료가 드물어 거의 보이지 않는다.
ax1.hist(heavy_tailed_returns, bins=50, density=True, alpha=0.6, label='Observed Returns')
x = np.linspace(heavy_tailed_returns.min(), heavy_tailed_returns.max(), 100)
ax1.plot(x, stats.norm.pdf(x, loc=heavy_tailed_returns.mean(),
                            scale=heavy_tailed_returns.std()),
         'r-', lw=2, label='Normal PDF')
ax1.set_xlabel('Daily Log-Return', fontsize=11)
ax1.set_ylabel('Density', fontsize=11)
ax1.set_title('Distribution Shape: Heavy Tails vs Normal', fontsize=12)
ax1.legend()
ax1.spines[["top", "right"]].set_visible(False)

# 오른쪽: Q-Q 그림. 히스토그램이 감춘 꼬리의 차이가 양끝의 휘어짐으로
# 또렷하게 드러난다.
stats.probplot(heavy_tailed_returns, dist="norm", plot=ax2)
ax2.set_title("Q-Q Plot: Revealing Heavy Tails", fontsize=12)
ax2.set_xlabel('Theoretical Normal Quantiles', fontsize=11)
ax2.set_ylabel('Sample Quantiles', fontsize=11)
ax2.spines[["top", "right"]].set_visible(False)
ax2.grid(True, alpha=0.3, linestyle='--')

plt.tight_layout()
plt.show()

r = heavy_tailed_returns
s = r.std(ddof=1)

# 왼쪽 그림: 히스토그램의 봉우리와 꼬리 칸을 수로 잰다.
dens, edges = np.histogram(r, bins=50, density=True)
print(f"칸 폭 {edges[1] - edges[0]:.5f}   자료 범위 [{r.min():.4f}, {r.max():.4f}]")
print(f"가장 높은 칸의 밀도 {dens.max():.3f}  vs  정규 봉우리 "
      f"{stats.norm.pdf(r.mean(), r.mean(), r.std()):.3f}")
print(f"양 끝 칸의 밀도: 왼쪽 {dens[0]:.3f}  오른쪽 {dens[-1]:.3f}"
      f"  (봉우리의 {dens[0] / dens.max():.4f}배)")
for k in (3, 4):
    n_out = np.sum(np.abs(r - r.mean()) > k * s)
    print(f"|x - 평균| > {k}s: 관측 {n_out}개, 정규 기대 "
          f"{2000 * 2 * stats.norm.sf(k):.2f}개")

# 오른쪽 그림: Q-Q 적합선에서의 잔차.
osm, osr = stats.probplot(r, dist="norm", fit=False)
b, a = np.polyfit(osm, osr, 1)
resid = osr - (a + b * osm)
print(f"\n적합선 절편 {a:.6f}  기울기 {b:.6f}  (표본표준편차 {s:.6f})")
for lo, lab in [(1, "|q| < 1"), (2, "|q| > 2")]:
    m = np.abs(osm) < lo if lo == 1 else np.abs(osm) > lo
    print(f"{lab}: {m.sum()}점, 최대 절대잔차 {np.abs(resid[m]).max():.5f}")
print(f"최솟값 {osr[0]:.4f}  적합선 {a + b * osm[0]:.4f}  잔차 {resid[0]:.4f}")
print(f"최댓값 {osr[-1]:.4f}  적합선 {a + b * osm[-1]:.4f}  잔차 {resid[-1]:.4f}")

print(f"\ng1 = {stats.skew(r):.4f}   G1 = {stats.skew(r, bias=False):.4f}   이론 = 0")
print(f"g2 = {stats.kurtosis(r):.4f}   G2 = {stats.kurtosis(r, bias=False):.4f}"
      f"   이론 = {6 / (5 - 4):.1f}")

칸 폭 0.00820   자료 범위 [-0.1170, 0.2932]
가장 높은 칸의 밀도 18.043  vs  정규 봉우리 15.483
양 끝 칸의 밀도: 왼쪽 0.061  오른쪽 0.061  (봉우리의 0.0034배)
|x - 평균| > 3s: 관측 18개, 정규 기대 5.40개
|x - 평균| > 4s: 관측 7개, 정규 기대 0.13개

적합선 절편 -0.000550  기울기 0.025061  (표본표준편차 0.025772)
|q| < 1: 1366점, 최대 절대잔차 0.00339
|q| > 2: 90점, 최대 절대잔차 0.20869
최솟값 -0.1170  적합선 -0.0856  잔차 -0.0314
최댓값 0.2932  적합선 0.0845  잔차 0.2087

g1 = 0.8632   G1 = 0.8639   이론 = 0
g2 = 10.4205   G2 = 10.4497   이론 = 6.0

rng = np.random.default_rng(7)
for nu in (5, 12):
    g2 = np.array([stats.kurtosis(stats.t.rvs(df=nu, size=2000, random_state=rng))
                   for _ in range(2000)])
    print(f"t_{nu}  이론 {6 / (nu - 4):.2f}  |  g2 중앙값 {np.median(g2):.2f}"
          f"  평균 {g2.mean():.2f}  5~95% [{np.percentile(g2, 5):.2f}, "
          f"{np.percentile(g2, 95):.2f}]  최대 {g2.max():.1f}")

t_5  이론 6.00  |  g2 중앙값 3.24  평균 4.51  5~95% [1.75, 10.42]  최대 197.5
t_12  이론 0.75  |  g2 중앙값 0.69  평균 0.74  5~95% [0.33, 1.29]  최대 5.0

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

# 1단계: 자료를 읽거나 만든다.
np.random.seed(42)
returns = np.random.standard_t(df=6, size=1500) * 0.025

# 2단계: 요약통계. 왜도와 초과첨도가 정규에서 얼마나 벗어났는지 먼저 본다.
mean_ret = returns.mean()
std_ret = returns.std()
skewness = stats.skew(returns)
kurtosis = stats.kurtosis(returns)  # Excess kurtosis

print(f"Mean:     {mean_ret:.4f}")
print(f"Std Dev:  {std_ret:.4f}")
print(f"Skewness: {skewness:.4f}")
print(f"Ex. Kurtosis: {kurtosis:.4f}")

# 3단계: 검정 셋. KS 는 모수를 자료에서 추정해 넘겼으므로 p-값이
# 실제보다 크게 나온다는 점을 감안해 읽어야 한다.
_, p_ks = stats.kstest(returns, 'norm', args=(mean_ret, std_ret))
_, p_jb = stats.jarque_bera(returns)

# anderson() 은 (통계량, 기각값, 유의수준) 을 돌려준다. p-값은 주지 않는다
ad_result = stats.anderson(returns, dist='norm')

print(f"\nKolmogorov-Smirnov test p-value: {p_ks:.4f}")
print(f"Jarque-Bera test p-value: {p_jb:.4g}")
print(f"Anderson-Darling statistic: {ad_result.statistic:.4f}")
print(f"  critical values (15/10/5/2.5/1%): {ad_result.critical_values}")

# 4단계: 그림으로 마무리. 검정이 기각했다면 어디가 어긋났는지를 여기서 본다.
fig, ax = plt.subplots(figsize=(8, 6))
stats.probplot(returns, dist="norm", plot=ax)
ax.set_title("Diagnostics: Are Returns Normal?", fontsize=12)
ax.spines[["top", "right"]].set_visible(False)
plt.show()

Mean:     -0.0002
Std Dev:  0.0297
Skewness: 0.5411
Ex. Kurtosis: 4.6396

Kolmogorov-Smirnov test p-value: 0.0710
Jarque-Bera test p-value: 9.067e-309
Anderson-Darling statistic: 3.6529
  critical values (15/10/5/2.5/1%): [0.574 0.654 0.785 0.916 1.089]

n = len(returns)
g1 = stats.skew(returns)          # bias=True 가 기본 -> g1
g2 = stats.kurtosis(returns)      # bias=True, fisher=True 가 기본 -> g2
JB = n / 6 * (g1 ** 2 + g2 ** 2 / 4)

print(f"g1 = {g1:.6f}   g2 = {g2:.6f}")
print(f"닫힌 꼴 JB      = {JB:.4f}")
print(f"scipy jarque_bera = {stats.jarque_bera(returns)[0]:.4f}")
print(f"chi2_2 꼬리확률  = {stats.chi2.sf(JB, 2):.4g}")

# 보정판을 넣으면 어긋난다.
G1 = stats.skew(returns, bias=False)
G2 = stats.kurtosis(returns, bias=False)
print(f"보정판으로 계산하면 JB = {n / 6 * (G1 ** 2 + G2 ** 2 / 4):.4f}  (어긋남)")

g1 = 0.541087   g2 = 4.639645
닫힌 꼴 JB      = 1418.5882
scipy jarque_bera = 1418.5882
chi2_2 꼬리확률  = 9.067e-309
보정판으로 계산하면 JB = 1430.0756  (어긋남)

rng = np.random.default_rng(3)
n, B = 1500, 4000
D = np.empty(B)
p = np.empty(B)
for i in range(B):
    x = rng.standard_normal(n)          # 귀무가설이 참인 자료
    D[i], p[i] = stats.kstest(x, 'norm', args=(x.mean(), x.std()))

print(f"명목 0.05 인데 실제 기각률 = {np.mean(p < 0.05):.5f}"
      f"  ({np.sum(p < 0.05)}/{B})")
print(f"명목 0.10 인데 실제 기각률 = {np.mean(p < 0.10):.5f}")
print(f"D 의 참 95 백분위수   = {np.percentile(D, 95):.5f}")
print(f"표준 KS 5% 임계값      = {1.358 / np.sqrt(n):.5f}")
print(f"릴리에포르 임계 근사    = {0.886 / np.sqrt(n):.5f}")
print(f"관측 D = 0.033246 의 백분위수 = {100 * np.mean(D < 0.033246):.3f}")

명목 0.05 인데 실제 기각률 = 0.00025  (1/4000)
명목 0.10 인데 실제 기각률 = 0.00100
D 의 참 95 백분위수   = 0.02317
표준 KS 5% 임계값      = 0.03506
릴리에포르 임계 근사    = 0.02288
관측 D = 0.033246 의 백분위수 = 99.925

from scipy.stats import t as student_t

# 정규가 아니라면 대안 분포를 찾는 것이 다음 걸음이다. t 를 적합하면
# 자유도가 꼬리 두께를 재는 값이 된다. 작을수록 두껍다.
df, loc, scale = student_t.fit(returns)
print(f"Fitted df: {df:.2f} (lower df → heavier tails)")

Fitted df: 6.68 (lower df → heavier tails)

print(f"적합 df {df:.4f}  loc {loc:.6f}  scale {scale:.6f}")
print(f"참값   df 6       loc 0        scale 0.025")
print(f"환산한 표준편차 scale*sqrt(df/(df-2)) = "
      f"{scale * np.sqrt(df / (df - 2)):.6f}")
print(f"표본표준편차                          = {returns.std(ddof=1):.6f}")
print(f"참 표준편차 0.025*sqrt(6/4)            = {0.025 * np.sqrt(6 / 4):.6f}")

# 참 nu = 6 인 자료를 되풀이 적합해 nu-hat 의 표본분포를 본다.
rng = np.random.default_rng(11)
fits = np.array([student_t.fit(student_t.rvs(df=6, scale=0.025, size=1500,
                                             random_state=rng))[:3]
                 for _ in range(400)])
nu_hat, sc_hat = fits[:, 0], fits[:, 2]
print(f"\nnu-hat 표본분포 (nu=6, n=1500, 400회)")
print(f"  중앙값 {np.median(nu_hat):.2f}  평균 {nu_hat.mean():.2f}"
      f"  5~95% [{np.percentile(nu_hat, 5):.2f}, "
      f"{np.percentile(nu_hat, 95):.2f}]  최대 {nu_hat.max():.2f}")
print(f"  P(nu-hat >= 6.68) = {np.mean(nu_hat >= 6.68):.3f}")
print(f"scale-hat 표본분포: 중앙값 {np.median(sc_hat):.5f}"
      f"  5~95% [{np.percentile(sc_hat, 5):.5f}, "
      f"{np.percentile(sc_hat, 95):.5f}]")

# nu 를 5 로 보든 8 로 보든 99% 분위수가 얼마나 달라지는가.
# 표준편차를 1 로 맞춰 견준다.
for nu in (5, 8):
    print(f"t_{nu} 의 표준화 99% 분위수 = "
          f"{stats.t.ppf(0.99, nu) / np.sqrt(nu / (nu - 2)):.3f}")
print(f"정규의 99% 분위수            = {stats.norm.ppf(0.99):.3f}")

적합 df 6.6758  loc -0.000468  scale 0.024562
참값   df 6       loc 0        scale 0.025
환산한 표준편차 scale*sqrt(df/(df-2)) = 0.029349
표본표준편차                          = 0.029671
참 표준편차 0.025*sqrt(6/4)            = 0.030619

nu-hat 표본분포 (nu=6, n=1500, 400회)
  중앙값 6.12  평균 6.26  5~95% [4.96, 8.10]  최대 11.68
  P(nu-hat >= 6.68) = 0.275
scale-hat 표본분포: 중앙값 0.02500  5~95% [0.02395, 0.02653]
t_5 의 표준화 99% 분위수 = 2.606
t_8 의 표준화 99% 분위수 = 2.508
정규의 99% 분위수            = 2.326