분산과 공분산¶

import numpy as np

# 공정한 주사위의 분산
values = np.arange(1, 7)          # 나올 수 있는 값 1~6
probs = np.ones(6) / 6            # 각각 확률 1/6

# 기댓값은 "값 x 확률"의 합이다
E_X = np.sum(values * probs)

# E[X^2]은 값을 제곱해서 같은 확률로 가중합한 것.
# 주의: E[X^2] 와 (E[X])^2 은 다르다. 그 차이가 곧 분산이다.
E_X2 = np.sum(values**2 * probs)

# 계산에 편한 공식: Var(X) = E[X^2] - (E[X])^2
# 정의식 E[(X - mu)^2] 를 전개하면 이 형태가 나온다.
var_X = E_X2 - E_X**2
print(f"E[X] = {E_X:.4f}")
print(f"E[X²] = {E_X2:.4f}")
print(f"Var(X) = {var_X:.4f}")
print(f"SD(X) = {np.sqrt(var_X):.4f}")

E[X] = 3.5000
E[X²] = 15.1667
Var(X) = 2.9167
SD(X) = 1.7078

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

def demonstrate_correlation():
    """무상관이 독립을 뜻하지 않는다는 것을 보인다."""
    np.random.seed(42)
    n = 10_000

    X = np.random.randn(n)      # 0을 중심으로 대칭인 표준정규
    Y = X ** 2                  # X만 알면 Y가 완전히 결정된다. 극단적인 종속이다.

    # 그런데 공분산은 0이 나온다. 이유는 대칭성에 있다.
    #   Cov(X, X^2) = E[X^3] - E[X]E[X^2] = 0 - 0*1 = 0
    # X가 0을 중심으로 대칭이면 E[X^3] = 0 이기 때문이다.
    # 양의 X가 만드는 기여와 음의 X가 만드는 기여가 정확히 상쇄된다.
    cov_XY = np.cov(X, Y)[0, 1]
    corr_XY = np.corrcoef(X, Y)[0, 1]

    print(f"Cov(X, X²) = {cov_XY:.4f} (theoretically 0)")
    print(f"Corr(X, X²) = {corr_XY:.4f}")
    print(f"Yet X and X² are clearly dependent!")

    # 표본공분산도 확률변수다. 그 퍼짐을 알아야 -0.0023 을 읽을 수 있다.
    #   Var(표본공분산) ~ E[(X-0)^2 (Y-1)^2]/n = E[X^6 - 2X^4 + X^2]/n = 10/n
    se_cov = np.sqrt(10 / n)
    se_corr = se_cov / np.sqrt(2)          # sd(X)=1, sd(Y)=sqrt(2)
    print(f"\n표본공분산의 SE = {se_cov:.4f},  관측/SE = {cov_XY / se_cov:+.3f}")
    print(f"표본상관계수의 SE = {se_corr:.4f},  관측/SE = {corr_XY / se_corr:+.3f}")

    # 독립이 아님을 사건 하나로 보인다. A = {|X| <= 1}, B = {Y <= 1} 은 같은 사건이다.
    A = np.abs(X) <= 1
    B = Y <= 1
    print(f"\nA = {{|X| <= 1}},  B = {{Y <= 1}}  (실은 같은 사건이다)")
    print(f"  P(A) = {A.mean():.4f}   이론 {2 * stats.norm.cdf(1) - 1:.4f}")
    print(f"  P(A and B) = {(A & B).mean():.4f}")
    print(f"  독립이라면 P(A)P(B) = {A.mean() * B.mean():.4f}")
    print(f"  비 = {(A & B).mean() / (A.mean() * B.mean()):.4f}")

    fig, ax = plt.subplots(figsize=(12, 4))
    ax.scatter(X[:500], Y[:500], alpha=0.3, s=10)
    ax.set_xlabel('X')
    ax.set_ylabel('Y = X²')
    ax.set_title(f'Uncorrelated but Dependent (ρ = {corr_XY:.3f})')
    ax.spines[['top', 'right']].set_visible(False)
    plt.tight_layout()
    plt.show()

demonstrate_correlation()

Cov(X, X²) = -0.0023 (theoretically 0)
Corr(X, X²) = -0.0016
Yet X and X² are clearly dependent!

표본공분산의 SE = 0.0316,  관측/SE = -0.073
표본상관계수의 SE = 0.0224,  관측/SE = -0.072

A = {|X| <= 1},  B = {Y <= 1}  (실은 같은 사건이다)
  P(A) = 0.6820   이론 0.6827
  P(A and B) = 0.6820
  독립이라면 P(A)P(B) = 0.4651
  비 = 1.4663

import numpy as np
import matplotlib.pyplot as plt

def portfolio_variance_demo():
    """분산투자의 이득이 상관계수에 어떻게 달려 있는지 보인다."""
    sigma1, sigma2 = 0.20, 0.30      # 두 자산의 변동성 20%, 30%
    correlations = [-0.5, 0.0, 0.5, 1.0]

    fig, ax = plt.subplots(figsize=(12, 4))
    weights = np.linspace(0, 1, 100)     # 자산1의 비중을 0에서 1까지

    for rho in correlations:
        # 공분산 = 상관계수 x 두 표준편차
        cov_12 = rho * sigma1 * sigma2

        # 포트폴리오 분산: Var(aX + bY) = a^2 Var(X) + b^2 Var(Y) + 2ab Cov(X,Y)
        # 마지막 교차항이 분산투자의 정체다. rho가 작을수록 이 항이 작아지고,
        # 음수이면 아예 위험을 깎아 낸다.
        port_var = (weights**2 * sigma1**2
                    + (1 - weights)**2 * sigma2**2
                    + 2 * weights * (1 - weights) * cov_12)
        port_sd = np.sqrt(port_var)

        # rho = 1 이면 곡선이 직선이 된다. 완전상관이면 섞어도 위험이 줄지 않는다.
        ax.plot(weights, port_sd, label=f'ρ = {rho}')

    ax.set_xlabel('Weight in Asset 1')
    ax.set_ylabel('Portfolio Std Dev')
    ax.set_title('Diversification: Portfolio Risk vs. Allocation')
    ax.legend()
    ax.spines[['top', 'right']].set_visible(False)
    plt.tight_layout()
    plt.show()

portfolio_variance_demo()

import numpy as np
np.set_printoptions(precision=4, suppress=True)

for r12, r13 in [(0.9, 0.9), (0.8, 0.5), (0.5, 0.5)]:
    c = r12 * r13
    h = np.sqrt((1 - r12 ** 2) * (1 - r13 ** 2))
    print(f"  r12={r12}, r13={r13}  ->  r23 는 [{c-h:+.4f}, {c+h:+.4f}] 안에만")

R = np.array([[1, .9, .9], [.9, 1, -.5], [.9, -.5, 1]])   # 범위 밖
ev, V = np.linalg.eigh(R)
a = V[:, 0]
print(f"\n  고유값 {np.round(ev, 4)}  <- 음수가 있다")
print(f"  a = {np.round(a, 3)} 일 때 Var(a'X) = {a @ R @ a:.4f}  <- 음의 분산")

  r12=0.9, r13=0.9  ->  r23 는 [+0.6200, +1.0000] 안에만
  r12=0.8, r13=0.5  ->  r23 는 [-0.1196, +0.9196] 안에만
  r12=0.5, r13=0.5  ->  r23 는 [-0.5000, +1.0000] 안에만

  고유값 [-0.5471  1.5     2.0471]  <- 음수가 있다
  a = [-0.635  0.546  0.546] 일 때 Var(a'X) = -0.5471  <- 음의 분산

import numpy as np

s2, cbar = 0.04, 0.012            # 평균분산 0.04 (sd 20%), 평균공분산 0.012
print(f"{'n':>6}{'포트폴리오 분산':>18}{'표준편차':>12}")
for n in (1, 2, 5, 10, 50, 500):
    v = s2 / n + (1 - 1 / n) * cbar
    print(f"{n:>6}{v:>18.5f}{np.sqrt(v):>12.4f}")
print(f"{'무한':>6}{cbar:>18.5f}{np.sqrt(cbar):>12.4f}  <- 넘을 수 없는 바닥")

     n          포트폴리오 분산        표준편차
     1           0.04000      0.2000
     2           0.02600      0.1612
     5           0.01760      0.1327
    10           0.01480      0.1217
    50           0.01256      0.1121
   500           0.01206      0.1098
    무한           0.01200      0.1095  <- 넘을 수 없는 바닥

import numpy as np

rng = np.random.default_rng(0)
print("참 공분산 = 단위행렬 (모든 고유값이 정확히 1이어야 한다)")
print(f"{'p':>5}{'n':>7}{'p/n':>7}{'최소 고유값':>14}{'최대 고유값':>14}{'조건수':>12}")
for p, n in [(5, 500), (50, 500), (200, 500), (400, 500), (600, 500)]:
    S = np.cov(rng.normal(size=(n, p)), rowvar=False)
    ev = np.linalg.eigvalsh(S)
    cond = ev[-1] / ev[0] if ev[0] > 1e-12 else float("inf")
    print(f"{p:>5}{n:>7}{p / n:>7.2f}{ev[0]:>14.4f}{ev[-1]:>14.4f}{cond:>12.1f}")

print("\n마르첸코-파스투르 법칙의 예측 구간  [(1-r)^2, (1+r)^2],  r = sqrt(p/n)")
for p, n in [(50, 500), (200, 500), (400, 500)]:
    r = np.sqrt(p / n)
    print(f"  p={p:>3}: [{(1 - r) ** 2:.4f}, {(1 + r) ** 2:.4f}]")

참 공분산 = 단위행렬 (모든 고유값이 정확히 1이어야 한다)
    p      n    p/n        최소 고유값        최대 고유값         조건수
    5    500   0.01        0.8270        1.1561         1.4
   50    500   0.10        0.4915        1.6853         3.4
  200    500   0.40        0.1531        2.5739        16.8
  400    500   0.80        0.0123        3.5212       286.6
  600    500   1.20       -0.0000        4.2999         inf

마르첸코-파스투르 법칙의 예측 구간  [(1-r)^2, (1+r)^2],  r = sqrt(p/n)
  p= 50: [0.4675, 1.7325]
  p=200: [0.1351, 2.6649]
  p=400: [0.0111, 3.5889]

import numpy as np

rng = np.random.default_rng(1)
p, n = 100, 150
beta = rng.uniform(0.5, 1.5, p)                       # 1-요인 구조
Sigma = np.outer(beta, beta) * 0.04 + np.diag(rng.uniform(0.01, 0.05, p))
one = np.ones(p)

def gmv(C):                                           # 최소분산 포트폴리오
    w = np.linalg.solve(C, one)
    return w / w.sum()

best = 1 / (one @ np.linalg.solve(Sigma, one))         # 참 Sigma 를 알 때의 하한
X = rng.standard_normal((n, p)) @ np.linalg.cholesky(Sigma).T
S = np.cov(X, rowvar=False)
target = np.trace(S) / p * np.eye(p)

print(f"참 공분산을 알 때 달성 가능한 최소 분산 = {best:.6f}\n")
print(f"{'축소계수':>10}{'조건수':>12}{'실현 분산':>13}{'하한 대비':>12}")
for a in (0.0, 0.1, 0.3, 0.5, 0.9):
    w = gmv((1 - a) * S + a * target)
    realized = w @ Sigma @ w
    print(f"{a:>10.1f}{np.linalg.cond((1 - a) * S + a * target):>12.1f}"
          f"{realized:>13.6f}{realized / best:>11.2f}배")

참 공분산을 알 때 달성 가능한 최소 분산 = 0.003391

      축소계수         조건수        실현 분산       하한 대비
       0.0      5734.8     0.010957       3.23배
       0.1       509.9     0.004808       1.42배
       0.3       142.2     0.004359       1.29배
       0.5        62.4     0.005170       1.52배
       0.9         7.9     0.019602       5.78배