# -*- coding: utf-8 -*-
"""Examples Chapter 3.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1X663eTE6xRxYnD0U-tLbj9dsIdKXLNO2

# Example p. 65-69:
"""

J = 100
mu = 0.1
T = 1
alpha = 3
beta = 2

from scipy.stats import poisson
import numpy as np

m = 1000
N_star = poisson.rvs(J*mu*T, size=m)
Xcal_star = np.zeros(m)
for i in range(m):
  U_star = np.random.uniform(size=N_star[i]) 
  Z_star = beta*(np.power((U_star),(-1/alpha))-1)
  Xcal_star[i] = sum(Z_star)

Xcal_bar_star = np.mean(Xcal_star)
s_star = np.std(Xcal_star)
kappa_star = np.mean(np.power(Xcal_star-np.mean(Xcal_star),4))/(np.power(np.std(Xcal_star), 4))
xi = beta/(alpha-1)
sigma = xi*np.sqrt(alpha/(alpha-2))

print(J*mu*T*xi, Xcal_bar_star, s_star/np.sqrt(m))
print(np.sqrt(J*mu*T*(sigma**2+xi**2)), s_star, s_star/np.sqrt(2*m)*np.sqrt(1+kappa_star/2))

M = 5
m1 = 1000
m2 = 10000
N_star = poisson.rvs(J*mu*T, size=(m2, M))
Xcal1_star = np.zeros((m1,M))
Xcal2_star = np.zeros((m2,M))

for j in range(M):
  for i in range(m2):
    U_star = np.random.uniform(size=N_star[i,j]) 
    Z_star = beta*(np.power(U_star, (-1/alpha))-1)
    if i < m1:
      Xcal1_star[i][j] = sum(Z_star)
    Xcal2_star[i][j] = sum(Z_star)

yps = [0.95,0.99]
print('Xcal_1', yps[0], [np.mean(sorted(Xcal1_star[:, i])[int(yps[0] * m1)]) for i in range(5)])
print('Xcal_1', yps[1], [np.mean(sorted(Xcal1_star[:, i])[int(yps[1] * m1)]) for i in range(5)])
print('Xcal_2', yps[0], [np.mean(sorted(Xcal2_star[:, i])[int(yps[0] * m2)]) for i in range(5)])
print('Xcal_2', yps[1], [np.mean(sorted(Xcal2_star[:, i])[int(yps[1] * m2)]) for i in range(5)])

m = 1000000
a = 0.5
b1 = 4
b2 = 2
b3 = 6

N_star = poisson.rvs(mu*T, size=(m,J))
Xcal1_star = np.zeros(m)
Xcal2_star = np.zeros(m)
Xcal3_star = np.zeros(m)
for i in range(m):
  U_star = np.random.uniform(size=sum(N_star[i,:])) 
  Z_star = beta*(np.power(U_star,-1/alpha)-1)
  H1_star = np.maximum(Z_star,a)-a
  H2_star = np.maximum(Z_star,a)-a
  H1_star = np.minimum(H1_star,b1)
  H2_star = np.minimum(H2_star, np.append(np.ones(sum(N_star[i,0:50]))*b2, np.ones(sum(N_star[i,50:100]))*b3))
  Xcal1_star[i] = sum(Z_star)
  Xcal2_star[i] = sum(H1_star)
  Xcal3_star[i] = sum(H2_star)

yps = [0.9,0.95,0.99,0.9997]
print([sorted(Xcal1_star)[int(m*yps[i])] for i in range(4)])
print([sorted(Xcal2_star)[int(m*yps[i])] for i in range(4)])
print([sorted(Xcal3_star)[int(m*yps[i])] for i in range(4)])

from scipy.stats import kde
import matplotlib.pyplot as plt


density1 = kde.gaussian_kde(Xcal_star)
x = np.linspace(0, 100, 200)
y1 = density1(x)
fig, ax = plt.subplots()
ax.plot(x, y1, label='Unlimited')
ax.set_xlabel('"Million USD')

density2 = kde.gaussian_kde(Xcal2_star)
y2 = density2(x)
ax.plot(x, y2, label='Fixed upper limit')

density3 = kde.gaussian_kde(Xcal3_star)
y3 = density3(x)
ax.plot(x, y3, label='Variable upper limit')
ax.legend()

a = np.arange(0,60,10)
Xcal_star = Xcal1_star
Xcal_re_star = [np.maximum(Xcal_star-a[i],0) for i in range(len(a))]

print([sorted(Xcal_star-Xcal_re_star[i])[int(m*0.99)] for i in range(len(Xcal_re_star))])
print([np.mean(Xcal_re_star[i]) for i in range(len(Xcal_re_star))])

"""# Example p. 72:"""

import numpy as np
import matplotlib.pyplot as plt

l0 = 60
le = 93
s = 0.1
r = 0.03
J1 = 100
J2 = 1000
K = 20

def copy_elements(arr, nr_times):
  assert len(arr) == len(nr_times)
  new_arr = np.array([])
  for i in range(len(arr)):
    new_arr = np.append(new_arr, np.tile(arr[i], int(nr_times[i])))
  return new_arr

pl = np.exp(-0.0009-0.000044*np.exp(0.09076*np.arange(l0,(le+K)+1)))
c = sum(pl[1:1+(le-l0)])
Nl1 = np.round(J1/c*pl[1:1+(le-l0)])
Nl2 = np.round(J2/c*pl[1:1+(le-l0)])
m = 50
PV01_star = np.zeros(shape=(K+1, m))
PV02_star = np.zeros(shape=(K+1, m))
for i in range(m):
   d = 1
   alive1 = np.ones(int(sum(Nl1)))
   alive2 = np.ones(int(sum(Nl2)))
   pj1 = copy_elements(pl[:(le-l0)], Nl1)
   pj2 = copy_elements(pl[:(le-l0)], Nl2)
   # pj1 = np.tile(pl[1:1+(le-l0)], Nl1)
   # pj2 = np.tile(pl[1:1+(le-l0)], Nl2)
   PV01_star[0][i] = sum(alive1)*s*d
   PV02_star[0][i] = sum(alive2)*s*d
   U1 = np.random.uniform(size=int(sum(Nl1)))
   alive1[U1 > pj1] = 0
   U2 = np.random.uniform(size=int(sum(Nl2)))
   alive2[U2 > pj2] = 0

   pj1 = copy_elements(pl[1:1+(le-l0)], Nl1)
   pj2 = copy_elements(pl[1:1+(le-l0)], Nl2)
   d = d/(1+r)
   for k in range(K):
      PV01_star[k+1][i] = PV01_star[k][i]+sum(alive1)*s*d
      PV02_star[k+1][i] = PV02_star[k][i]+sum(alive2)*s*d
      U1 = np.random.uniform(size=int(sum(Nl1)))
      alive1[U1 > pj1] = 0
      U2 = np.random.uniform(size=int(sum(Nl2)))
      alive2[U2 > pj2] = 0
      pj1 = copy_elements(pl[k+1:k+1+(le-l0)], Nl1)
      pj2 = copy_elements(pl[k+1:k+1+(le-l0)], Nl2)
      d = d/(1+r)

plt.plot(np.arange(K+1), PV01_star)
plt.xlabel('Years ahead')
plt.ylabel('Present value')
plt.title('100 policies')
plt.show()

plt.plot(np.arange(K+1), PV02_star)
plt.xlabel('Years ahead')
plt.ylabel('Present value')
plt.title('1000 policies')
plt.show()

"""# Example p. 76 - 77"""

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

J = 4
T = 1
r = 0.04
rg = 0.07
rc1 = 0.1
rc2 = 0.15
rho = 0.5
w = 0.25
m = 100000
M = 10
sigma = np.arange(0.05,0.41,0.01)

pi11_star = np.zeros((len(sigma), M))
pi12_star = np.zeros((len(sigma), M))
pi21_star = np.zeros((len(sigma), M))
pi22_star = np.zeros((len(sigma), M))
for i in range(M):
  eta = norm.rvs(size=m)
  eta = np.tile(eta, (J, 1))
  eps1 = norm.rvs(size=(J, m))
  eps2 = np.sqrt(rho)*eta+np.sqrt(1-rho)*eps1

  Rcal1_star = np.zeros((len(sigma), m))
  Rcal2_star = np.zeros((len(sigma), m))
  for s in range(len(sigma)):
    for j in range(J):
      Rcal1_star[s, :] = Rcal1_star[s, :]+w*np.exp(r*T-0.5*sigma[s]**2*T+sigma[s]*np.sqrt(T)*eps1[j, :])-w
      Rcal2_star[s, :] = Rcal2_star[s, :]+w*np.exp(r*T-0.5*sigma[s]**2*T+sigma[s]*np.sqrt(T)*eps2[j, :])-w

  X11_star = np.maximum(rg-Rcal1_star,0)
  X12_star = np.maximum(rg-Rcal2_star,0)
  X21_star = np.minimum(np.maximum(rg-Rcal2_star,0),rc1-Rcal2_star)
  X22_star = np.minimum(np.maximum(rg-Rcal2_star,0),rc2-Rcal2_star)
  pi11_star[:,i] = np.exp(-r*T)*np.mean(X11_star, axis=1)
  pi12_star[:,i] = np.exp(-r*T)*np.mean(X12_star, axis=1)
  pi21_star[:,i] = np.exp(-r*T)*np.mean(X21_star, axis=1)
  pi22_star[:,i] = np.exp(-r*T)*np.mean(X22_star, axis=1)

fig, ax = plt.subplots()
ax.plot(sigma, pi11_star, label='Unlimited')
ax.set_xlabel('Volatility')
ax.set_ylabel('Price')

ax.plot(sigma, pi12_star, label='Unlimited')
ax.set_xlabel('Volatility')
ax.set_ylabel('Price')
ax.set_title('Put options')
plt.show()

fig, ax = plt.subplots()
ax.plot(sigma, pi21_star, label='Unlimited')
ax.set_xlabel('Volatility')
ax.set_ylabel('Price')

ax.plot(sigma, pi22_star, label='Unlimited')
ax.set_xlabel('Volatility')
ax.set_ylabel('Price')
ax.set_title('Cliquet options')
plt.show()

"""# Example p. 80-82:"""

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

J = 100
mu = 0.1
T = 1
K = 20
alpha = 3
beta = 2
a = 0.5
b = 4
Pi_0 = 6
Rcal1 = 0.04
xi = 0.07531
sigma = 0.2
rho = 0.5 
w = 0.25
yps = np.arange(0.005,0.251,0.001)

m = 10000
Xcal_star = np.zeros((K,m))
Scal1_star = np.zeros((K+1,m))
Scal2_star = np.zeros((K+1,m))
Rcal2_star = np.zeros(m)
for k in range(K):
  N_star = poisson.rvs(J*mu*T, size=m)
  for i in range(m):
    U_star = np.random.uniform(size=N_star[i]) 
    Z_star = np.minimum(np.maximum(beta*((U_star)**(-1/alpha)-1)-a,0),b)
    Xcal_star[k,i] = sum(Z_star)
  eta = norm.rvs(size=(m,1))
  eps1 = norm.rvs(size=(m,4))
  eps2 = np.sqrt(rho)*eta+np.sqrt(1-rho)*eps1
  Rcal2_star = (1+Rcal2_star)*np.sum(w*np.exp(xi+sigma*eps2), axis=1)-1

  Scal1_star[k+1,:] = Scal1_star[k,:]+(Xcal_star[k,:]-Pi_O)/(1+Rcal1)**k  
  Scal2_star[k+1,:] = Scal2_star[k,:]+(Xcal_star[k,:]-Pi_O)/(1+Rcal2_star)

Ycal_star = np.tile(np.reshape(np.arange(1,K+1)*Pi_0, (K,-1)), (1,m)) - np.cumsum(Xcal_star, axis=0)
Ycal_star = np.insert(Ycal_star, 0, np.zeros(m), axis=0)
Ycal_bar_star1 = np.min(Ycal_star[:2,], axis=0)
Ycal_bar_star2 = np.min(Ycal_star[:6,], axis=0)
Ycal_bar_star3 = np.min(Ycal_star[:11,], axis=0)
Ycal_bar_star4 = np.min(Ycal_star[:16,], axis=0)
Ycal_bar_star5 = np.min(Ycal_star[:21,], axis=0)

nu0yps11 = [-np.sort(Ycal_bar_star1)[int(yps[i]*m)] for i in range(len(yps))]
nu0yps12 = [-np.sort(Ycal_bar_star2)[int(yps[i]*m)] for i in range(len(yps))]
nu0yps13 = [-np.sort(Ycal_bar_star3)[int(yps[i]*m)] for i in range(len(yps))]
nu0yps14 = [-np.sort(Ycal_bar_star4)[int(yps[i]*m)] for i in range(len(yps))]
nu0yps15 = [-np.sort(Ycal_bar_star5)[int(yps[i]*m)] for i in range(len(yps))]

Scal_bar_star1 = np.max(Scal1_star[:11,:], axis=0)
Scal_bar_star2 = np.max(Scal2_star[:11,:], axis=0)

nu0yps21 = [np.sort(Scal_bar_star1)[int((1-yps[i])*m)] for i in range(len(yps))]
nu0yps22 = [np.sort(Scal_bar_star2)[int((1-yps[i])*m)] for i in range(len(yps))]

import matplotlib.pyplot as plt

for i in range(100):
  plt.plot(np.arange(11), 15.6+Ycal_star[:11,i], zorder=0)
plt.xlabel('Years ahead')
plt.ylabel('Account')
plt.title('Simulating underwriting')
plt.hlines(0, xmin=0, xmax=10, zorder=1)

plt.plot(yps, nu0yps13, label='Underwriting')
plt.xlabel('Ruin probability')
plt.ylabel('Initial capital')
plt.title('Financial income included')
plt.plot(yps, nu0yps21, label='Risk-free return')
plt.plot(yps, nu0yps22, label='Equity return')
plt.legend()