[R] Discrete Event

Insurance Risk Model

Insurance Risk Model을 적합할 때 필요한 확률변수들은 다음과 같다.

• 보험계약자의 보험금 청구 : Rate가 $\lambda$인 포아송 과정
• 보험 청구 건당 금액 : $F$
• 신규 계약 : Rate가 $\nu$인 포아송 과정
• 보험유지기간 : 평균이 $1/\mu$인 지수분포

계약자는 시간당 $c$ 만큼의 보험금을 납부한다. 시점 $t=0$에서 계약자의 수가 $n_0$, 보험사의 자본금이 $a_0$일 때, 기간 $(0,T)$ 동안 자본금이 마이너스가 되지 않을 확률을 시뮬레이션으로 추정해보자.

System State : $(n,a)$. $n$ : 계약자 수, a는 자본금
Event List : $t_E$ (사건 발생시각)
사건의 유형으로는 세 가지가 존재하는데, 사건의 종류와 각각의 발생확률은 다음과 같다.

• 보험금 청구 : $\frac{n\lambda}{n\lambda+\nu+n\mu}$
• 신규 가입자 발생 : $\frac{\nu}{n\lambda+\nu+n\mu}$
• 기존 계약해지 : $\frac{n\mu}{n\lambda+\nu+n\mu}$
  

Output Variable : $I$ $(0,T)$ 사이에 계속 $a\geq 0$이면 $I=0$, 아니면 $I=1$.

이때, 사건 발생 시간 간격의 분포는 평균이 $1/(n\lambda+\nu+n\mu)$인 지수분포를 따른다.

예제를 보면서 이해한 이론을 점검해보자.

보험금 청구는 $\lambda$가 10인 포아송 과정을 따른다고 가정하자. 청구 금액은 평균이 1000$인 지수분포를 따르는 확률변수이고, 계약자는 보험사에 하루 11000$씩 보험금을 납입한다.
보험사의 초기자본이 25000$인 상태에서 365일 동안 항상 양수일 확률을 추정하라.

위 예제에서는 고객이 1명으로 고정된 상황이기 때문에, 신규가입이나 계약해지는 발생하지 않는다.
따라서 $\nu=\mu=0,\ n=1$으로 변하지 않는다는 것을 생각하며 R로 구현해보자.

> n.sim = 100 # 모의실험 반복횟수
> n0 = 1; a0 = 25000; T=365; c = 11000
> lambda = 10; nu = 0; mu = 0 # 보험금 청구율, 보험 가입율, 해지율
> generate.Y = function() rexp(1, rate = 1/1000)
> I = numeric(n.sim)
> for (i in 1:n.sim){
+   t = 0; a = a0; n = n0
+   total.rate = nu + n*mu + n*lambda
+   tE = rexp(1, total.rate)
+   repeat{
+     if (tE > T) {I[i] = 1; break}
+     if (tE <= T){
+       a = a + n*c*(tE - t)
+       t = tE
+       J = sample(1:3, 1, prob = c(nu, n*mu, n*lambda))
+       if (J == 1) n = n + 1
+       if (J == 2) n = n - 1
+       if (J == 3) {Y = generate.Y(); if (Y > a) {I[i] = 0; break} else a = a - Y}
+       tE = t + rexp(1, total.rate)
+     }
+   }
+ }
> (Ibar = mean(I)); (B = sqrt(var(I)/n.sim))
[1] 0.85
[1] 0.03588703
> c(Ibar-B, Ibar+B)
[1] 0.814113 0.885887

우선 시뮬레이션 횟수만큼의 길이를 갖는 벡터를 만들어 아웃풋을 저장할 수 있게 했다.
한 번의 시뮬레이션 내에는 repeat가 있는데, 자본금이 바닥나거나 365일이 지나면 종료된다.

J는 각 사건의 발생확률에 따라 1,2,3 중 하나의 값을 갖는다.
사건 발생 간격은 rate가 $n\lambda+\nu+n\mu$인 지수분포를 따른다는 것을 기억하자.

100번의 시뮬레이션 결과 초기 자본이 바닥나지 않을 추정확률은 0.85이다.
구간 추정을 하는 경우 표준오차를 같이 구해야한다는 점을 기억하자.

Exercising a Stock Option

• American Call Option
  특정 주식 한 주를 임의의 시점 $t=1,...,N$에서 가격 $K$에 살 수 있는 권한을 Call Option이라고 한다. 이 권한을 어느 시점에 행사하는 것이 가장 좋은 선택일까?

확률변수 $X_i\sim iid \ N(\mu,\sigma^2)$에 대해서 주가(Stock Price) 모형을 다음과 같이 정의한다.

$\alpha=\mu+\frac{\sigma^2}{2}$라고 할 때, $\alpha \geq0$ 이면 권한을 행사하는 최적의 전략이 알려져 있다.
바로 마지막 시점 $N$까지 기다렸다가 그 시점에서의 주가가 $K$를 넘기면 행사하고, 그렇지 않으면 행사를 포기하는 전략이다.

하지만 $\alpha \leq0$일 때의 최적의 전략을 찾기는 어렵다. 그래서 모의실험을 통해 최적의 시점을 찾아보려고 하는데, 알려져 있는 것들 중 하나를 소개해보려고 한다.

• $\alpha \leq0$일 때 권장하는 전략
  $P_m(=S_{N-m})$을 만기를 $m$ 앞둔 시점에서의 주가라고 할 때 다음 조건이 동시에 만족되면 만기를 $m$ 앞둔 시점에서 권한을 행사한다.

1. $P_m>K$
2. 모든 $i=1,...,m$ 에 대하여

이때 $\Phi$는 표준정규분포의 cdf이다. 이 식은 Black - Scholes Option Pricing Formula를 참조하자.
위 두 가지 조건을 만족하면 이익은 $P_m - K$이고, 만족하지 않으면 이익은 0이 된다.

• 위 전략에 따라 옵션을 행사했을 때 얻는 기대이익을 추정하는 모의실험 절차
  Step 1 $\ N, K,\mu, \sigma, S_0$의 초기값을 설정한다. 단, $\alpha = \mu+\frac{\sigma^2}{2}<0$ 이어야 한다.
  Step 2 $\ P_m \leftarrow P_{m+1}e^{X}, \quad X \sim N(\mu,\sigma^2)$
  Step 3 앞에서 제시한 전략의 두 조건을 만족한다면 이익은 $P_m - K$ 이며 실행이 완료된다. 그렇지 않으면 $m \leftarrow m-1$ 으로 설정 후 $m\geq0$ 이면 Step 2 반복. $m<0$ 이면 이익은 0이되고 실행 완료.

이 절차를 이용하여 예제를 하나 해결해보자.

20일 동안 보유하는 100$ Call Option의 기대가치를 시뮬레이션으로 추정하라. 현재 주식의 가격이 100$이고, $\mu=-0.05, \ \sigma=0.3$이라고 가정한다.

Step 1 과 같이 변수들의 초기값을 지정한 뒤, 각 모의실험에서의 이익을 저장할 벡터를 만들어준다.
R에서는 0번 인덱스가 존재하지 않으므로 $P_0$과 같은 값을 P[0]에 저장할 수 없다. 따라서 P의 길이를 1만큼 늘려준 다음 1번 인덱스부터 P[0]을 저장하자.

call_option = function(N = 20, K = 100, Szero = 100, mu = -0.05, sig = 0.3){
  alpha = mu + 0.5 * sig^2
  P = numeric(N + 1) # 이익을 저장할 벡터 생성
  P[N+1] = Szero
  m = N - 1
  flag = FALSE
  for (i in N:1) P[i] = P[i+1] * exp(rnorm(1, mean = mu, sd = sig))
  repeat{
    m.plus = m + 1
    if (P[m.plus] > K) flag = TRUE
    if (flag & m > 0){
      bi = ((1:m)*mu - log(K/P[m.plus]))/(sig*sqrt(1:m))
      op = P[m.plus]*exp((1:m)*alpha)*pnorm(sig*sqrt(1:m)+bi) - K*pnorm(bi)
      flag = all(P[m.plus] > K + op)
    }
    if (flag) break else m = m - 1
    if (m < 0) break
  }
  if (flag) E = P[m.plus] - K else E = 0
  return (E)
}

n.sim = 1000
E = replicate(n.sim, call_option()) # 1000회 시뮬레이션
c(mean(E),sd(E)/sqrt(n.sim),mean(E > 0)) # 기대이익, 표준오차, 권한 행사 확률