# -*- coding: utf-8 -*-
"""Examples chapter 6.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1kkwOdrmODOFgNiyBRG4R1FmQG84ivFpG

# Example p. 186:
"""

import numpy as np
from scipy.stats import gamma, poisson

xi = 0.05
alpha = 100
T = 1
J1 = 10000
J2 = 1000000
K = 25
m = 100000

G = gamma.rvs(alpha, size=(K,m), scale=1/alpha)
mu = xi*G
Ncal1 = np.zeros((K,m))
Ncal2 = np.zeros((K,m))
for k in range(K):
  Ncal1[k,:] = poisson.rvs(J1*mu[k,:]*T, size=m)
  Ncal2[k,:] = poisson.rvs(J2*mu[k,:]*T, size=m)

import matplotlib.pyplot as plt


plt.plot(np.arange(K), Ncal1[:,1:20])
plt.xlabel("Year")
plt.ylabel("Number of claims")
plt.title("10,000 policies")
plt.hlines(y=J1*xi, xmin=0, xmax=24)

plt.plot(np.arange(K), Ncal2[:,1:20])
plt.xlabel("Year")
plt.ylabel("Number of claims")
plt.title("1,000,000 policies")
plt.hlines(y=J2*xi, xmin=0, xmax=24)

"""# Example p. 200:"""

import numpy as np
import matplotlib.pyplot as plt


alpha0 = 0.0009
alpha1 = 0.000044
alpha2 = 0.09076

Ncal0 = 10000
plia = 0.007
plai = 0.001

l0 = 30
le = 75

pl = np.exp(-alpha0 -alpha1 * np.exp(0.09076 * np.arange(l0, le)))
kpl0 = np.cumprod(pl)
kpl0 = np.insert(kpl0, 0, 1.0)

ql = 1 - pl

plt.plot(np.arange(l0, le), ql)

from scipy.stats import bernoulli

# states:
# 0 - active
# 1 - disabled
# 2 - dead
policies = np.ones((le-l0,Ncal0)) * 2
for i in range(Ncal0):
  policies[0][i] = 0
  for k in range(1,le-l0):
    if policies[k-1][i] != 2: # if not dead
      Bde = bernoulli.rvs(pl[k])
      if Bde == 0: # if died
        break
      else: # if survived
        if policies[k-1][i] == 0: # if active
          Bdi = bernoulli.rvs(plia)
          if Bdi == 1: # if got disabled
            policies[k][i] = 1
          else:
            policies[k][i] = 0
        else:
          if policies[k-1][i] == 1: # if disabled
            Ba = bernoulli.rvs(plai)
            if Ba == 1: # if got rehabilitated
              policies[k][i] = 0
            else:
              policies[k][i] = 1

Nde = np.count_nonzero(policies == 2, axis=1)
Na = np.count_nonzero(policies == 0, axis=1)
Ndi = np.count_nonzero(policies == 1, axis=1)

import matplotlib.pyplot as plt

plt.plot(np.arange(le-l0), Na, label='active')
plt.plot(np.arange(le-l0), Nde, label='dead')
plt.plot(np.arange(le-l0), Ndi, label='disabled')
plt.legend()
plt.show()