05_monte_carlo.py

  1"""
  2몬테카를로 시뮬레이션 (Monte Carlo Simulation)
  3Monte Carlo Methods
  4
  5확률적 방법을 사용한 수치 계산 및 시뮬레이션입니다.
  6"""
  7
  8import numpy as np
  9import matplotlib.pyplot as plt
 10from typing import Callable, Tuple
 11
 12
 13# =============================================================================
 14# 1. π 추정 (원의 넓이)
 15# =============================================================================
 16def estimate_pi(n_samples: int) -> Tuple[float, float]:
 17    """
 18    단위 원을 사용한 π 추정
 19
 20    정사각형 [-1, 1] x [-1, 1] 내에 무작위 점을 뿌리고
 21    단위 원 내부에 떨어지는 비율을 계산
 22
 23    원의 넓이 / 정사각형 넓이 = π / 4
 24    → π = 4 * (원 내부 점 / 전체 점)
 25
 26    오차: O(1/√n)
 27    """
 28    # 균일 분포로 점 생성
 29    x = np.random.uniform(-1, 1, n_samples)
 30    y = np.random.uniform(-1, 1, n_samples)
 31
 32    # 원 내부에 있는 점 개수
 33    inside = np.sum(x**2 + y**2 <= 1)
 34
 35    pi_estimate = 4 * inside / n_samples
 36    std_error = 4 * np.sqrt(inside * (n_samples - inside) / n_samples**3)
 37
 38    return pi_estimate, std_error
 39
 40
 41def estimate_pi_convergence(max_samples: int = 100000) -> Tuple[np.ndarray, np.ndarray]:
 42    """π 추정의 수렴 과정"""
 43    x = np.random.uniform(-1, 1, max_samples)
 44    y = np.random.uniform(-1, 1, max_samples)
 45    inside = (x**2 + y**2 <= 1).cumsum()
 46    n = np.arange(1, max_samples + 1)
 47    return n, 4 * inside / n
 48
 49
 50# =============================================================================
 51# 2. 몬테카를로 적분
 52# =============================================================================
 53def monte_carlo_integrate(
 54    f: Callable[[np.ndarray], np.ndarray],
 55    a: float,
 56    b: float,
 57    n_samples: int
 58) -> Tuple[float, float]:
 59    """
 60    1D 몬테카를로 적분
 61
 62    ∫[a,b] f(x)dx ≈ (b-a) * (1/n) * Σf(x_i)
 63
 64    Args:
 65        f: 피적분 함수
 66        a, b: 적분 구간
 67        n_samples: 샘플 수
 68
 69    Returns:
 70        (적분값 추정, 표준 오차)
 71    """
 72    x = np.random.uniform(a, b, n_samples)
 73    fx = f(x)
 74
 75    integral = (b - a) * np.mean(fx)
 76    variance = np.var(fx)
 77    std_error = (b - a) * np.sqrt(variance / n_samples)
 78
 79    return integral, std_error
 80
 81
 82def monte_carlo_integrate_nd(
 83    f: Callable[[np.ndarray], float],
 84    bounds: list,
 85    n_samples: int
 86) -> Tuple[float, float]:
 87    """
 88    다차원 몬테카를로 적분
 89
 90    Args:
 91        f: 다변수 함수 f(x) where x is array
 92        bounds: [(a1,b1), (a2,b2), ...] 각 차원의 범위
 93        n_samples: 샘플 수
 94    """
 95    dim = len(bounds)
 96    volume = np.prod([b - a for a, b in bounds])
 97
 98    # 각 차원에서 균일 샘플링
 99    samples = np.array([
100        np.random.uniform(a, b, n_samples)
101        for a, b in bounds
102    ]).T  # shape: (n_samples, dim)
103
104    values = np.array([f(x) for x in samples])
105
106    integral = volume * np.mean(values)
107    std_error = volume * np.std(values) / np.sqrt(n_samples)
108
109    return integral, std_error
110
111
112# =============================================================================
113# 3. 중요도 샘플링 (Importance Sampling)
114# =============================================================================
115def importance_sampling(
116    f: Callable[[np.ndarray], np.ndarray],
117    g: Callable[[np.ndarray], np.ndarray],
118    sample_g: Callable[[int], np.ndarray],
119    n_samples: int
120) -> Tuple[float, float]:
121    """
122    중요도 샘플링
123
124    목표: ∫f(x)dx 를 더 효율적으로 추정
125
126    g(x)를 중요도 분포로 사용:
127    ∫f(x)dx = ∫(f(x)/g(x))g(x)dx = E_g[f(x)/g(x)]
128
129    분산 감소: f(x)와 비슷한 g(x) 선택
130
131    Args:
132        f: 피적분 함수 * 원래 분포
133        g: 중요도 분포의 PDF
134        sample_g: g에서 샘플 생성 함수
135        n_samples: 샘플 수
136    """
137    x = sample_g(n_samples)
138    weights = f(x) / g(x)
139
140    integral = np.mean(weights)
141    std_error = np.std(weights) / np.sqrt(n_samples)
142
143    return integral, std_error
144
145
146# =============================================================================
147# 4. 랜덤 워크 시뮬레이션
148# =============================================================================
149def random_walk_1d(n_steps: int, n_walks: int = 1) -> np.ndarray:
150    """
151    1D 랜덤 워크
152
153    각 스텝에서 +1 또는 -1로 이동
154    """
155    steps = np.random.choice([-1, 1], size=(n_walks, n_steps))
156    positions = np.cumsum(steps, axis=1)
157    return positions
158
159
160def random_walk_2d(n_steps: int) -> Tuple[np.ndarray, np.ndarray]:
161    """2D 랜덤 워크"""
162    angles = np.random.uniform(0, 2*np.pi, n_steps)
163    dx = np.cos(angles)
164    dy = np.sin(angles)
165    x = np.cumsum(dx)
166    y = np.cumsum(dy)
167    return x, y
168
169
170# =============================================================================
171# 5. 옵션 가격 결정 (Black-Scholes Monte Carlo)
172# =============================================================================
173def black_scholes_mc(
174    S0: float,
175    K: float,
176    T: float,
177    r: float,
178    sigma: float,
179    n_simulations: int,
180    option_type: str = 'call'
181) -> Tuple[float, float]:
182    """
183    유럽형 옵션 가격의 몬테카를로 추정
184
185    기하 브라운 운동:
186    S_T = S_0 * exp((r - σ²/2)T + σ√T * Z)
187
188    Args:
189        S0: 현재 주가
190        K: 행사가
191        T: 만기 (년)
192        r: 무위험 이자율
193        sigma: 변동성
194        n_simulations: 시뮬레이션 횟수
195        option_type: 'call' or 'put'
196    """
197    # 주가 시뮬레이션
198    Z = np.random.standard_normal(n_simulations)
199    ST = S0 * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z)
200
201    # 페이오프 계산
202    if option_type == 'call':
203        payoffs = np.maximum(ST - K, 0)
204    else:  # put
205        payoffs = np.maximum(K - ST, 0)
206
207    # 현재 가치로 할인
208    price = np.exp(-r * T) * np.mean(payoffs)
209    std_error = np.exp(-r * T) * np.std(payoffs) / np.sqrt(n_simulations)
210
211    return price, std_error
212
213
214def black_scholes_analytical(
215    S0: float,
216    K: float,
217    T: float,
218    r: float,
219    sigma: float,
220    option_type: str = 'call'
221) -> float:
222    """Black-Scholes 해석적 해 (비교용)"""
223    from scipy.stats import norm
224
225    d1 = (np.log(S0/K) + (r + 0.5*sigma**2)*T) / (sigma*np.sqrt(T))
226    d2 = d1 - sigma * np.sqrt(T)
227
228    if option_type == 'call':
229        price = S0 * norm.cdf(d1) - K * np.exp(-r*T) * norm.cdf(d2)
230    else:
231        price = K * np.exp(-r*T) * norm.cdf(-d2) - S0 * norm.cdf(-d1)
232
233    return price
234
235
236# =============================================================================
237# 6. 버핏 바늘 문제 (Buffon's Needle)
238# =============================================================================
239def buffon_needle(
240    needle_length: float,
241    line_spacing: float,
242    n_drops: int
243) -> Tuple[float, float]:
244    """
245    버핏 바늘 문제로 π 추정
246
247    평행선 사이에 바늘을 떨어뜨렸을 때
248    선을 교차할 확률 = 2L / (πD)
249
250    → π = 2L * n_drops / (D * crossings)
251
252    Args:
253        needle_length: 바늘 길이 (L)
254        line_spacing: 선 간격 (D), L ≤ D
255        n_drops: 바늘 떨어뜨리기 횟수
256    """
257    if needle_length > line_spacing:
258        raise ValueError("바늘 길이는 선 간격보다 작아야 합니다")
259
260    # 바늘 중심의 위치 (0 ~ D/2)
261    y_center = np.random.uniform(0, line_spacing/2, n_drops)
262
263    # 바늘 각도 (0 ~ π)
264    theta = np.random.uniform(0, np.pi, n_drops)
265
266    # 바늘 끝이 선을 넘는지 확인
267    # 바늘 끝의 y 좌표 변화: (L/2) * sin(θ)
268    crosses = y_center <= (needle_length/2) * np.sin(theta)
269    n_crossings = np.sum(crosses)
270
271    if n_crossings == 0:
272        return float('inf'), float('inf')
273
274    pi_estimate = 2 * needle_length * n_drops / (line_spacing * n_crossings)
275
276    return pi_estimate, 1 / np.sqrt(n_crossings)  # 대략적인 오차
277
278
279# =============================================================================
280# 시각화
281# =============================================================================
282def plot_monte_carlo_examples():
283    """몬테카를로 예제 시각화"""
284    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
285
286    # 1. π 추정
287    ax = axes[0, 0]
288    n = 10000
289    x = np.random.uniform(-1, 1, n)
290    y = np.random.uniform(-1, 1, n)
291    inside = x**2 + y**2 <= 1
292
293    ax.scatter(x[inside], y[inside], c='blue', s=1, alpha=0.5)
294    ax.scatter(x[~inside], y[~inside], c='red', s=1, alpha=0.5)
295    theta = np.linspace(0, 2*np.pi, 100)
296    ax.plot(np.cos(theta), np.sin(theta), 'k-', linewidth=2)
297    ax.set_xlim(-1.1, 1.1)
298    ax.set_ylim(-1.1, 1.1)
299    ax.set_aspect('equal')
300    pi_est = 4 * np.sum(inside) / n
301    ax.set_title(f'π 추정: {pi_est:.4f} (n={n})')
302
303    # 2. π 추정 수렴
304    ax = axes[0, 1]
305    n, estimates = estimate_pi_convergence(50000)
306    ax.semilogx(n, estimates, 'b-', alpha=0.7)
307    ax.axhline(y=np.pi, color='r', linestyle='--', label=f'π = {np.pi:.6f}')
308    ax.set_xlabel('샘플 수')
309    ax.set_ylabel('π 추정값')
310    ax.set_title('π 추정 수렴')
311    ax.legend()
312    ax.grid(True)
313
314    # 3. 랜덤 워크
315    ax = axes[1, 0]
316    for _ in range(5):
317        x, y = random_walk_2d(1000)
318        ax.plot(np.concatenate([[0], x]), np.concatenate([[0], y]), alpha=0.7)
319    ax.plot(0, 0, 'ko', markersize=10)
320    ax.set_xlabel('x')
321    ax.set_ylabel('y')
322    ax.set_title('2D 랜덤 워크 (5개 경로)')
323    ax.axis('equal')
324    ax.grid(True)
325
326    # 4. 옵션 가격 수렴
327    ax = axes[1, 1]
328    S0, K, T, r, sigma = 100, 100, 1, 0.05, 0.2
329    exact = black_scholes_analytical(S0, K, T, r, sigma)
330
331    ns = np.logspace(2, 5, 20).astype(int)
332    estimates = []
333    errors = []
334    for n in ns:
335        est, err = black_scholes_mc(S0, K, T, r, sigma, n)
336        estimates.append(est)
337        errors.append(err)
338
339    ax.errorbar(ns, estimates, yerr=errors, fmt='o-', capsize=3)
340    ax.axhline(y=exact, color='r', linestyle='--', label=f'해석해: {exact:.4f}')
341    ax.set_xscale('log')
342    ax.set_xlabel('시뮬레이션 횟수')
343    ax.set_ylabel('옵션 가격')
344    ax.set_title('콜옵션 가격 수렴')
345    ax.legend()
346    ax.grid(True)
347
348    plt.tight_layout()
349    plt.savefig('/opt/projects/01_Personal/03_Study/Numerical_Simulation/examples/monte_carlo.png', dpi=150)
350    plt.close()
351    print("그래프 저장: monte_carlo.png")
352
353
354# =============================================================================
355# 테스트
356# =============================================================================
357def main():
358    print("=" * 60)
359    print("몬테카를로 시뮬레이션 (Monte Carlo)")
360    print("=" * 60)
361
362    np.random.seed(42)
363
364    # 1. π 추정
365    print("\n[1] π 추정 (원의 넓이)")
366    print("-" * 40)
367    for n in [1000, 10000, 100000, 1000000]:
368        pi_est, std_err = estimate_pi(n)
369        print(f"n={n:>8}: π ≈ {pi_est:.6f} ± {std_err:.6f}, "
370              f"오차: {abs(pi_est - np.pi):.6f}")
371
372    # 2. 몬테카를로 적분
373    print("\n[2] 몬테카를로 적분")
374    print("-" * 40)
375
376    # ∫[0,1] x² dx = 1/3
377    f1 = lambda x: x**2
378    integral, error = monte_carlo_integrate(f1, 0, 1, 100000)
379    print(f"∫x²dx [0,1] = {integral:.6f} ± {error:.6f} (정확: 0.333333)")
380
381    # ∫[0,π] sin(x) dx = 2
382    f2 = lambda x: np.sin(x)
383    integral, error = monte_carlo_integrate(f2, 0, np.pi, 100000)
384    print(f"∫sin(x)dx [0,π] = {integral:.6f} ± {error:.6f} (정확: 2.0)")
385
386    # 3. 다차원 적분
387    print("\n[3] 다차원 적분")
388    print("-" * 40)
389
390    # 단위 구의 부피: 4π/3 ≈ 4.189
391    def sphere_indicator(x):
392        return 1 if np.sum(x**2) <= 1 else 0
393
394    volume, error = monte_carlo_integrate_nd(sphere_indicator, [(-1,1)]*3, 100000)
395    print(f"단위 구 부피 = {volume:.4f} ± {error:.4f} (정확: {4*np.pi/3:.4f})")
396
397    # 4. 옵션 가격
398    print("\n[4] 유럽형 콜옵션 가격")
399    print("-" * 40)
400
401    S0, K, T, r, sigma = 100, 100, 1, 0.05, 0.2
402
403    exact = black_scholes_analytical(S0, K, T, r, sigma)
404    mc_price, mc_error = black_scholes_mc(S0, K, T, r, sigma, 100000)
405
406    print(f"Black-Scholes 해석해: {exact:.4f}")
407    print(f"Monte Carlo (n=100000): {mc_price:.4f} ± {mc_error:.4f}")
408    print(f"오차: {abs(mc_price - exact):.4f}")
409
410    # 5. 버핏 바늘
411    print("\n[5] 버핏 바늘 문제")
412    print("-" * 40)
413
414    pi_est, _ = buffon_needle(1, 2, 100000)
415    print(f"π 추정 (L=1, D=2, n=100000): {pi_est:.6f}")
416
417    # 6. 랜덤 워크 통계
418    print("\n[6] 1D 랜덤 워크 통계")
419    print("-" * 40)
420
421    n_steps = 1000
422    n_walks = 10000
423    positions = random_walk_1d(n_steps, n_walks)
424    final_positions = positions[:, -1]
425
426    print(f"{n_steps}스텝 후 ({n_walks}개 워크):")
427    print(f"  평균 위치: {np.mean(final_positions):.2f} (이론: 0)")
428    print(f"  표준편차: {np.std(final_positions):.2f} (이론: {np.sqrt(n_steps):.2f})")
429
430    # 시각화
431    try:
432        plot_monte_carlo_examples()
433    except Exception as e:
434        print(f"그래프 생성 실패: {e}")
435
436    print("\n" + "=" * 60)
437    print("몬테카를로 방법 정리")
438    print("=" * 60)
439    print("""
440    장점:
441    - 고차원 문제에서도 수렴 속도 유지 (차원의 저주 회피)
442    - 복잡한 영역/조건에 적용 용이
443    - 구현이 단순
444
445    단점:
446    - 수렴 속도가 느림: O(1/√n)
447    - 확률적 → 결과에 불확실성
448    - 많은 샘플 필요
449
450    분산 감소 기법:
451    - 중요도 샘플링 (Importance Sampling)
452    - 대조 변량 (Control Variates)
453    - 안티테틱 변량 (Antithetic Variates)
454    - 층화 샘플링 (Stratified Sampling)
455    - 준난수 (Quasi-random / Low-discrepancy sequences)
456
457    응용:
458    - 금융: 파생상품 가격 결정, 리스크 분석
459    - 물리: 통계역학, 입자 시뮬레이션
460    - 컴퓨터 그래픽: 경로 추적 렌더링
461    - 최적화: 시뮬레이티드 어닐링, MCMC
462    """)
463
464
465if __name__ == "__main__":
466    main()