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

Automatically generated by Colaboratory.

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

# Section 3.2

# Exercise 3.4
"""

import numpy as np

alpha = 3
beta = 2
lambda_ = 10
a = [0,1,4,6]
m = 100000

U = np.random.uniform(size=m)
Z = beta * (np.power(U, -1/alpha)-1)
Zre0 = np.maximum(Z-a[0],0)
Zre1 = np.maximum(Z-a[1],0)
Zre2 = np.maximum(Z-a[2],0)
Zre3 = np.maximum(Z-a[3],0)
print(lambda_*np.mean(Zre0), lambda_*np.mean(Zre1),lambda_*np.mean(Zre2),lambda_*np.mean(Zre3))

"""# Exercise 3.5"""

alpha = 3
beta = 2
lambda_ = 10
b = [2,4,6]
gamma = 0.5
m = 10000

a = 1
U = np.random.uniform(size=m)
Z = beta * (np.power(U, -1/alpha)-1)
Zre1 = np.maximum(Z-a,0)
Zre20 = gamma*np.maximum(Zre1-b[0],0)
Zre21 = gamma*np.maximum(Zre1-b[1],0)
Zre22 = gamma*np.maximum(Zre1-b[2],0)
print(lambda_*np.mean(Zre20),lambda_*np.mean(Zre21),lambda_*np.mean(Zre22))

a = 2
U = np.random.uniform(size=m)
Z = beta * (np.power(U, -1/alpha)-1)
Zre1 = np.maximum(Z-a,0)
Zre20 = gamma*np.maximum(Zre1-b[0],0)
Zre21 = gamma*np.maximum(Zre1-b[1],0)
Zre22 = gamma*np.maximum(Zre1-b[2],0)
print(lambda_*np.mean(Zre20),lambda_*np.mean(Zre21),lambda_*np.mean(Zre22))

"""# Section 3.3

# Exercise 3.6
"""

from scipy.stats import poisson, lognorm

lambda_ = 10
xi = 0
sigma = 1
m = 10000

N = poisson.rvs(lambda_, size=m)
Xcal = np.zeros(m)
for i in range(m):
  Z = lognorm.rvs(sigma, loc=xi, size=N[i]) # remember to check parametrisation
  Xcal[i] = sum(Z)

print(lambda_*np.exp(xi+0.5*np.power(sigma,2)), np.sqrt(lambda_)*np.exp(xi+np.power(sigma, 2))) # theoretical
print(np.mean(Xcal), np.std(Xcal)) # simulated

"""# Exercise 3.7"""

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

# Use Xcal from previous exercise
density = kde.gaussian_kde(Xcal)
x = np.linspace(0, 100, 200)
y = density(x)
plt.plot(x, y)
plt.vlines(np.mean(Xcal), ymin=0, ymax=0.06)
plt.vlines(sorted(Xcal)[int(m*0.95)], ymin=0, ymax=0.06)
plt.vlines(sorted(Xcal)[int(m*0.99)], ymin=0, ymax=0.06)

print("mean ", np.mean(Xcal))
print("q95 ", sorted(Xcal)[int(m*0.95)])
print("q99 ", sorted(Xcal)[int(m*0.99)])

"""# Exercise 3.8"""

from scipy.stats import poisson, lognorm, kde
import numpy as np
import matplotlib.pyplot as plt

lambda_ = 100
xi = 0
sigma = 1
m = 10000

N = poisson.rvs(lambda_, size=m)
Xcal = np.zeros(m)
for i in range(m):
  Z = lognorm.rvs(sigma, loc=xi, size=N[i])
  Xcal[i] = sum(Z)

density = kde.gaussian_kde(Xcal)
x = np.linspace(0, 400, 200)
y = density(x)
plt.plot(x, y)

print(np.mean(Xcal))
print(sorted(Xcal)[int(m*0.95)])
print(sorted(Xcal)[int(m*0.99)])

"""# Exercise 3.9"""

from scipy.stats import gamma

lambda_ = 10
xi = 1
alpha = 5
m = 10000

N = poisson.rvs(lambda_, size=m)
xcal = np.zeros(m)
for i in range(m):
  Z = xi * gamma.rvs(alpha, size=N[i], scale=1/alpha)
  Xcal[i] = sum(Z)

print(lambda_*xi, np.sqrt(lambda_)*xi*np.sqrt(1/alpha+1))
print(np.mean(Xcal), np.std(Xcal))

density = kde.gaussian_kde(Xcal)
x = np.linspace(0, 30, 200)
y = density(x)
plt.plot(x, y)

plt.vlines(np.mean(Xcal), ymin=0, ymax=0.12)
plt.vlines(sorted(Xcal)[int(m*0.95)], ymin=0, ymax=0.12)
plt.vlines(sorted(Xcal)[int(m*0.99)], ymin=0, ymax=0.12)

print(np.mean(Xcal))
print(sorted(Xcal)[int(m*0.95)])
print(sorted(Xcal)[int(m*0.99)])

"""# Exercise 3.10"""

xi = 0
sigma = 1
a = 0
b = 1
m = 100000

Z = lognorm.rvs(sigma, loc=xi, size=m)
H = np.minimum(np.maximum(Z-a, 0), b)

print(np.mean(Z), np.std(Z))
print(np.mean(H), np.std(H))

"""# Exercise 3.11"""

lambda_ = 10
xi = 1
alpha = 5
a = [0.5, 1, 2, 100]
b = 100
m = 100000

N = poisson.rvs(lambda_, size=m)
Xcal0 = np.zeros(m)
Xcal1 = np.zeros(m)
Xcal2 = np.zeros(m)
Xcal3 = np.zeros(m)

for i in range(m):
  Z = xi * gamma.rvs(alpha, size=N[i], scale=1/alpha)
  H0 = np.minimum(np.maximum(Z-a[0], 0), b)
  H1 = np.minimum(np.maximum(Z-a[1], 0), b)
  H2 = np.minimum(np.maximum(Z-a[2], 0), b)
  H3 = np.minimum(np.maximum(Z-a[3], 0), b)
  Xcal0[i] = sum(Z - H0)
  Xcal1[i] = sum(Z - H1)
  Xcal2[i] = sum(Z - H2)
  Xcal3[i] = sum(Z - H3)

print(sorted(Xcal0)[int(m*0.95)], sorted(Xcal0)[int(m*0.99)])
print(sorted(Xcal1)[int(m*0.95)], sorted(Xcal1)[int(m*0.99)])
print(sorted(Xcal2)[int(m*0.95)], sorted(Xcal2)[int(m*0.99)])
print(sorted(Xcal3)[int(m*0.95)], sorted(Xcal3)[int(m*0.99)])

"""# Exercise 3.12"""

lambda_ = 10
xi = 1
alpha = 5
a = [0.5, 1, 2, 100]
b = 100
m = 100000

N = poisson.rvs(lambda_, size=m)
Xcal0 = np.zeros(m)
Xcal1 = np.zeros(m)
Xcal2 = np.zeros(m)
Xcal3 = np.zeros(m)

for i in range(m):
  Z = xi * gamma.rvs(alpha, size=N[i], scale=1/alpha)
  H0 = np.minimum(np.maximum(Z-a[0], 0), b)
  H1 = np.minimum(np.maximum(Z-a[1], 0), b)
  H2 = np.minimum(np.maximum(Z-a[2], 0), b)
  H3 = np.minimum(np.maximum(Z-a[3], 0), b)
  Xcal0[i] = sum(H0)
  Xcal1[i] = sum(H1)
  Xcal2[i] = sum(H2)
  Xcal3[i] = sum(H3)

print(sorted(Xcal0)[int(m*0.95)], sorted(Xcal0)[int(m*0.99)])
print(sorted(Xcal1)[int(m*0.95)], sorted(Xcal1)[int(m*0.99)])
print(sorted(Xcal2)[int(m*0.95)], sorted(Xcal2)[int(m*0.99)])
print(sorted(Xcal3)[int(m*0.95)], sorted(Xcal3)[int(m*0.99)])

print(np.mean(Xcal0))
print(np.mean(Xcal1))
print(np.mean(Xcal2))
print(np.mean(Xcal3))

"""# Exercise 3.13"""

lambda_ = 10
xi = 0
sigma = 1
a = 15
b = 25
m = 10000

N = poisson.rvs(lambda_, size=m)
Xcal = np.zeros(m)
for i in range(m):
  Z = lognorm.rvs(sigma, loc=xi, size=N[i])
  Xcal[i] = sum(Z)
Xcalre = np.minimum(np.maximum(Xcal - a,0), b)
Xcalce = Xcal - Xcalre

print(np.mean(Xcalre))
print(sorted(Xcalce)[int(m*0.95)], sorted(Xcalce)[int(m*0.99)])
print(sorted(Xcalre)[int(m*0.95)], sorted(Xcalre)[int(m*0.99)])
print(sorted(Xcal)[int(m*0.95)], sorted(Xcal)[int(m*0.99)])

"""# Section 3.4

# Exercise 3.17
"""

import matplotlib.pyplot as plt

theta0 = 0.0009
theta1 = 0.000044
theta2 = 0.09076
l0 = 35
K = 50

pl1 = np.exp(-theta0-theta1*np.exp(theta2*np.arange(l0, K -1 +l0)))
pl1 = np.append(1, pl1)
kpl01 = np.cumprod(pl1)

plt.plot(np.arange(K), kpl01)

l0 = 50
pl2 = np.exp(-theta0-theta1*np.exp(theta2*np.arange(l0, K -1 +l0)))
pl2 = np.append(1, pl2)
kpl02 = np.cumprod(pl2)

plt.plot(np.arange(K), kpl02)

l0 = 65
pl3 = np.exp(-theta0-theta1*np.exp(theta2*np.arange(l0, K -1 +l0)))
pl3 = np.append(1, pl3)
kpl03 = np.cumprod(pl3)

plt.plot(np.arange(K), kpl03)

plt.plot(np.arange(K), kpl01)
plt.plot(np.arange(K), kpl02)
plt.plot(np.arange(K), kpl03)

"""# Exercise 3.18"""

theta0 = 0.0009
theta1 = 0.000044
theta2 = 0.09076
s = 1
lr = 65
r = 0.02
d = 1/(1+r)
l0 = [20,35,50,65]
le = 120

pl0 = np.exp(-theta0-theta1*np.exp(theta2*np.arange(l0[0], le)))
pl0 = np.append(1, pl0)
kpl00 = np.cumprod(pl0) 
pl1 = np.exp(-theta0-theta1*np.exp(theta2*np.arange(l0[1],le)))
pl0 = np.append(1, pl1)
kpl01 = np.cumprod(pl1) 
pl2 = np.exp(-theta0-theta1*np.exp(theta2*np.arange(l0[2],le)))
pl0 = np.append(1, pl2)
kpl02 = np.cumprod(pl2) 
pl3 = np.exp(-theta0-theta1*np.exp(theta2*np.arange(l0[3],le)))
pl0 = np.append(1, pl3)
kpl03 = np.cumprod(pl3)

pil00 = s*sum(np.power(d, np.arange(lr - l0[0],le-l0[0]))*kpl00[lr - l0[0]:le-l0[0]])
pil01 = s*sum(np.power(d, np.arange(lr - l0[1],le-l0[1]))*kpl01[lr - l0[1]:le-l0[1]])
pil02 = s*sum(np.power(d, np.arange(lr - l0[2],le-l0[2]))*kpl02[lr - l0[2]:le-l0[2]])
pil03 = s*sum(np.power(d, np.arange(lr - l0[3],le-l0[3]))*kpl03[lr - l0[3]:le-l0[3]])
print(pil00,pil01,pil02,pil03)

r = 0.04
d = 1/(1 + r)
pil00 = s*sum(np.power(d, np.arange(lr - l0[0],le-l0[0]))*kpl00[lr - l0[0]:le-l0[0]])
pil01 = s*sum(np.power(d, np.arange(lr - l0[1],le-l0[1]))*kpl01[lr - l0[1]:le-l0[1]])
pil02 = s*sum(np.power(d, np.arange(lr - l0[2],le-l0[2]))*kpl02[lr - l0[2]:le-l0[2]])
pil03 = s*sum(np.power(d, np.arange(lr - l0[3],le-l0[3]))*kpl03[lr - l0[3]:le-l0[3]])
print(pil00,pil01,pil02,pil03)

"""# Exercise 3.19"""

theta0 = 0.0009
theta1 = 0.000044
theta2 = 0.09076
s = 1
lr1 = 55
lr2 = 70
r = 0.02
d = 1/(1+r)
l0 = 45
le = 120

pl = np.exp(-theta0-theta1*np.exp(theta2*np.arange(l0,le)))
pl = np.append(1, pl)
kpl0 = np.cumprod(pl)

pi = np.zeros(lr2-lr1)
for lr in range(lr1,lr2):
   pi[lr-lr1] = s*sum(np.power(d, np.arange(lr-l0,le-l0))*kpl0[lr-l0:le-l0])

plt.scatter(np.arange(lr1,lr2), pi)

l0 = 35
pl1 = np.exp(-theta0-theta1*np.exp(theta2*np.arange(l0,le)))
pl1 = np.append(1, pl1)
kpl01 = np.cumprod(pl1)
pi1 = np.zeros(lr2-lr1)
for lr in range(lr1,lr2):
   pi1[lr-lr1] = s*sum(np.power(d, np.arange(lr-l0,le-l0))*kpl01[lr-l0:le-l0])

plt.scatter(np.arange(lr1,lr2), pi)
plt.scatter(np.arange(lr1,lr2), pi1)

l0 = 45
r = 0.04
d = 1/(1 + r)
pl2 = np.exp(-theta0-theta1*np.exp(theta2*np.arange(l0,le)))
pl2 = np.append(1, pl2)
kpl02 = np.cumprod(pl2)
pi2 = np.zeros(lr2-lr1)
for lr in range(lr1,lr2):
   pi2[lr-lr1] = s*sum(np.power(d, np.arange(lr-l0,le-l0))*kpl02[lr-l0:le-l0])

l0 = 35
pl3 = np.exp(-theta0-theta1*np.exp(theta2*np.arange(l0,le)))
pl3 = np.append(1, pl3)
kpl03 = np.cumprod(pl3)
pi3 = np.zeros(lr2-lr1)
for lr in range(lr1,lr2):
   pi3[lr-lr1] = s*sum(np.power(d, np.arange(lr-l0,le-l0))*kpl03[lr-l0:le-l0])

plt.scatter(np.arange(lr1,lr2), pi)
plt.scatter(np.arange(lr1,lr2), pi1)
plt.scatter(np.arange(lr1,lr2), pi2)
plt.scatter(np.arange(lr1,lr2), pi3)

"""# Exercise 3.20"""

from scipy.stats import binom

theta0 = 0.0009
theta1 = 0.000044
theta2 = 0.09076
l0 = 30
K = 25
Ncal0 = 1000
m = 50

pl = np.exp(-theta0-theta1*np.exp(theta2*np.arange(l0,K + l0)))

Ncal = np.ones((K+1,m)) * Ncal0
for i in range(m):
  for k in range(K):
    Ncal[k+1,i] = binom.rvs(p=pl[k], n=int(Ncal[k,i]))

for i in range(m):
  plt.plot(np.arange(K+1) ,100*Ncal[:, i]/Ncal0)
# plt.plot(np.arange(K+1) ,100*Ncal/Ncal0)

"""# Exercise 3.21"""

theta0 = 0.0009
theta1 = 0.000044
theta2 = 0.09076
l0 = 65
K = 25
Ncal0 = 1000
s = 1
m = 50

pl = np.exp(-theta0-theta1*np.exp(theta2*np.arange(l0,K + l0)))
pl = np.append(1, pl)
kpl0 = np.cumprod(pl)

Ncal = np.ones((K+1,m)) * Ncal0
for i in range(m):
  for k in range(K):
    Ncal[k+1,i] = binom.rvs(p=kpl0[k+1], n=int(Ncal[k,i]))
Xcal = s*Ncal

for i in range(m):
  plt.plot(np.arange(K+1), Xcal[:, i])

Ncal0 = 10000
Ncal = np.ones((K+1,m)) * Ncal0
for i in range(m):
  for k in range(K):
    Ncal[k+1,i] = binom.rvs(p=kpl0[k+1], n=int(Ncal[k,i]))
Xcal = s*Ncal
for i in range(m):
  plt.plot(np.arange(K+1), Xcal[:, i])

Ncal0 = 100000
Ncal = np.ones((K+1,m)) * Ncal0
for i in range(m):
  for k in range(K):
    Ncal[k+1,i] = binom.rvs(p=kpl0[k+1], n=int(Ncal[k,i]))
Xcal = s*Ncal
for i in range(m):
  plt.plot(np.arange(K+1), Xcal[:, i])

"""# Exercise 3.22"""

theta0 = 0.0009
theta1 = 0.000044
theta2 = 0.09076
l0 = 30
K = 25
Ncal0 = 1000
m = 50

pl = np.exp(-theta0-theta1*np.exp(theta2*np.arange(l0,K + l0)))

Ncal = np.ones((K+1,m)) * Ncal0
for i in range(m):
  for k in range(K):
    Ncal[k+1,i] = binom.rvs(p=pl[k], n=int(Ncal[k,i]))

Xcal = s*(Ncal[:K, :] - Ncal[1:K+1,:])

for i in range(m):
  plt.plot(np.arange(K), Xcal[:, i])

Ncal0 = 10000
m = 50
pl = np.exp(-theta0-theta1*np.exp(theta2*np.arange(l0,K + l0)))

Ncal = np.ones((K+1,m)) * Ncal0
for i in range(m):
  for k in range(K):
    Ncal[k+1,i] = binom.rvs(p=pl[k], n=int(Ncal[k,i]))

Xcal = s*(Ncal[:K, :] - Ncal[1:K+1,:])

for i in range(m):
  plt.plot(np.arange(K), Xcal[:, i])

Ncal0 = 100000
m = 50
pl = np.exp(-theta0-theta1*np.exp(theta2*np.arange(l0,K + l0)))

Ncal = np.ones((K+1,m)) * Ncal0
for i in range(m):
  for k in range(K):
    Ncal[k+1,i] = binom.rvs(p=pl[k], n=int(Ncal[k,i]))

Xcal = s*(Ncal[:K, :] - Ncal[1:K+1,:])

for i in range(m):
  plt.plot(np.arange(K), Xcal[:, i])

"""# Section 3.5

# Exercise 3.23
"""

from scipy.stats import norm

nu0 = 1
sigma = 0.25
r = 0.03
T = 1
rg = np.arange(-0.02,0.10,0.02)

a = (np.log(1+rg)-r*T+0.5*sigma**2*T)/(sigma*np.sqrt(T))
pi_nu0 = ((1+rg)*np.exp(-r*T)*norm.cdf(a)-norm.cdf(a-sigma*np.sqrt(T)))*nu0
print(pi_nu0)

"""# Exercise 3.24"""

nu0 = 1
sigma = 0.25
r = 0.03
T = 1
rg = 0.04

m1 = 10000
m2 = 100000
m3 = 1000000
eps = norm.rvs(size=m3) 
R = np.exp(r*T-0.5*sigma**2*T+sigma*np.sqrt(T)*eps)-1
X1 = np.maximum(rg-R[:m1],0)
X2 = np.maximum(rg-R[:m2],0)
X3 = np.maximum(rg-R[:m3],0)
pi_nu01 = np.exp(-r*T)*np.mean(X1)*nu0
pi_nu02 = np.exp(-r*T)*np.mean(X2)*nu0
pi_nu03 = np.exp(-r*T)*np.mean(X3)*nu0

print(pi_nu01,pi_nu02,pi_nu03)

"""# Exercise 3.29"""

nu0 = 1
sigma = 0.25
r = 0.03
T = 1
rg = 0.04
rc = np.array([0.09,0.12,0.15,0.5])

ag = (np.log(1+rg)-r*T+0.5*sigma**2*T)/(sigma*np.sqrt(T))
ac = (np.log(1+rc)-r*T+0.5*sigma**2*T)/(sigma*np.sqrt(T))
pi_nu0 = ((1+rg)*np.exp(-r*T)*norm.cdf(ag)-norm.cdf(ag-sigma*np.sqrt(T)))*nu0-((1+rc)*np.exp(-r*T)*norm.cdf(ac)-norm.cdf(ac-sigma*np.sqrt(T)))*nu0-(1-np.exp(-r*T)*(1+rc))*nu0

print(pi_nu0)

"""# Section 3.6

# Exercise 3.30
"""

from scipy.stats import poisson, lognorm
import matplotlib.pyplot as plt

lambda_ = 10
xi = 0
sigma = 1
Pi = 18
nu0 = 35
K = 20
m = 200

N = poisson.rvs(lambda_,size=K*m)
Xcal = np.zeros(K*m)
for i in range(int(K*m)):
  Z = lognorm.rvs(sigma, size=N[i],loc=xi)
  Xcal[i] = np.sum(Z)
Xcal = np.reshape(Xcal, (-1, m))

Ycal = np.expand_dims(np.arange(1, K+1)*Pi +nu0, axis=1) - np.cumsum(Xcal, axis=0)
Ycal = np.append(np.ones((1,m))*nu0,Ycal,axis=0)

for i in range(m):
  plt.plot(np.arange(K+1), Ycal[:, i])
  plt.ylim((-200,200))

Pi = 16.4
Ycal = np.expand_dims(np.arange(1, K+1)*Pi +nu0, axis=1) - np.cumsum(Xcal, axis=0)
Ycal = np.append(np.ones((1,m))*nu0,Ycal,axis=0)

for i in range(m):
  plt.plot(np.arange(K+1), Ycal[:, i])
  plt.ylim((-200,200))

"""# Exercise 3.31"""

from scipy.stats import poisson, lognorm
import numpy as np
import matplotlib.pyplot as plt

lambda_ = 1000
xi = 0
sigma = 1
Pi = 1800
nu0 = 3500
K = 20
m = 200

N = poisson.rvs(lambda_, size=m*K)
Xcal = np.zeros(K*m)
for i in range(m*K):
  Z = lognorm.rvs(sigma, size=N[i], loc=xi)
  Xcal[i] = sum(Z)
Xcal = np.reshape(Xcal, (-1, m))

Ycal = np.tile(np.reshape(np.arange(1,K+1)*Pi + nu0, (K,-1)), (1,m)) - np.cumsum(Xcal, axis=0)
Ycal = np.insert(Ycal, 0, np.ones(m)*nu0, axis=0)

for i in range(m):
  plt.plot(np.arange(K+1), Ycal[:, i])

"""# Exercise 3.32"""

lambda_ = 10
xi = 0.455
sigma = 0.3
Pi = 18
nu0 = 35
K = 20
m = 200

N = poisson.rvs(lambda_, size=m*K)
Xcal = np.zeros(K*m)
for i in range(m*K):
  Z = lognorm.rvs(sigma, size=N[i], loc=xi)
  Xcal[i] = sum(Z)
Xcal = np.reshape(Xcal, (-1, m))

Ycal = np.tile(np.reshape(np.arange(1,K+1)*Pi + nu0, (K,-1)), (1,m)) - np.cumsum(Xcal, axis=0)
Ycal = np.insert(Ycal, 0, np.ones(m)*nu0, axis=0)

for i in range(m):
  plt.plot(np.arange(K+1), Ycal[:, i])
  plt.ylim((-100, 200))

"""# Exercise 3.33"""

lambda_ = 10
xi = 0
sigma = 1
Pi = 18
nu0 = 0
K = 20
yps = np.array([0.05,0.01])
m = 10000

N = poisson.rvs(lambda_, size=m*K)
Xcal = np.zeros(K*m)
for i in range(m*K):
  Z = lognorm.rvs(sigma, size=N[i], loc=xi)
  Xcal[i] = sum(Z)
Xcal = np.reshape(Xcal, (-1, m))
Ycal = np.tile(np.reshape(np.arange(1,K+1)*Pi + nu0, (K,-1)), (1,m)) - np.cumsum(Xcal, axis=0)
Ycal = np.insert(Ycal, 0, np.ones(m)*nu0, axis=0)

Ymin = np.min(Ycal, axis=0)
print(-np.sort(Ymin)[int(m*yps[0])], -np.sort(Ymin)[int(m*yps[1])])

Ymin = np.min(Ycal[:10], axis=0)
print(-np.sort(Ymin)[int(m*yps[0])], -np.sort(Ymin)[int(m*yps[1])])

Ymin = np.min(Ycal[:5], axis=0)
print(-np.sort(Ymin)[int(m*yps[0])], -np.sort(Ymin)[int(m*yps[1])])

"""# Exercise 3.34"""

lambda_ = 10
xi = 0
sigma = 1
Pi = 18
nu0 = 35
K = 20
r = 0.04
m = 200

N = poisson.rvs(lambda_, size=m*K)
Xcal = np.zeros(K*m)
for i in range(m*K):
  Z = lognorm.rvs(sigma, size=N[i], loc=xi)
  Xcal[i] = sum(Z)
Xcal = np.reshape(Xcal, (-1, m))

for i in range(m):
  Xcal[:, i] = (Xcal[:,i] - Pi)/(1+r)**np.arange(1, K+1)
Scal = np.cumsum(Xcal, axis=0)

for i in range(m):
  Scal[:, i] = (1+r)**np.arange(1, K+1)*(nu0-Scal[:, i])
Ycal = np.insert(Scal, 0, np.ones(m)*nu0, axis=0)

for i in range(m):
  plt.plot(np.arange(K+1), Ycal[:, i])
  # plt.ylim((-100, 200))

"""# Exercise 3.35"""

lambda_ = 10
xi = 0.455
sigma = 0.3
Pi = 18
nu0 = 35
K = 20
r = 0.04
m = 200

N = poisson.rvs(lambda_, size=m*K)
Xcal = np.zeros(K*m)
for i in range(m*K):
  Z = lognorm.rvs(sigma, size=N[i], loc=xi)
  Xcal[i] = sum(Z)
Xcal = np.reshape(Xcal, (-1, m))

for i in range(m):
  Xcal[:, i] = (Xcal[:,i] - Pi)/(1+r)**np.arange(1, K+1)
Scal = np.cumsum(Xcal, axis=0)

for i in range(m):
  Scal[:, i] = (1+r)**np.arange(1, K+1)*(nu0-Scal[:, i])
Ycal = np.insert(Scal, 0, np.ones(m)*nu0, axis=0)

for i in range(m):
  plt.plot(np.arange(K+1), Ycal[:, i])

"""# Exercise 3.36"""

lambda_ = 10
xi = 0
sigma = 1
Pi = 18
nu0 = 0
K = 20
r = 0.04
yps = [0.05,0.01]
m = 10000

N = poisson.rvs(lambda_, size=m*K)
Xcal = np.zeros(K*m)
for i in range(m*K):
  Z = lognorm.rvs(sigma, size=N[i], loc=xi)
  Xcal[i] = sum(Z)
Xcal = np.reshape(Xcal, (-1, m))

Scal = np.zeros_like(Scal)
for i in range(m):
  Scal[:, i] = (Xcal[:,i] - Pi)/(1+r)**np.arange(1, K+1)
Scal = np.cumsum(Scal, axis=0)
Scalmax = np.max(Scal, axis=0)

print(np.sort(Scalmax)[int(m*(1-yps[0]))], np.sort(Scalmax)[int(m*(1-yps[1]))])

Scalmax = np.max(Scal[:10, :], axis=0)
print(np.sort(Scalmax)[int(m*(1-yps[0]))], np.sort(Scalmax)[int(m*(1-yps[1]))])

Scalmax = np.max(Scal[:5, :], axis=0)
print(np.sort(Scalmax)[int(m*(1-yps[0]))], np.sort(Scalmax)[int(m*(1-yps[1]))])

Ycal = np.tile(np.reshape(np.arange(1,K+1)*Pi + nu0, (K,-1)), (1,m)) - np.cumsum(Xcal, axis=0)
Ycal = np.insert(Ycal, 0, np.ones(m)*nu0, axis=0)
Ymin = np.min(Ycal, axis=0)
print(-np.sort(Ymin)[int(m*yps[0])], -np.sort(Ymin)[int(m*yps[1])])

Ymin = np.min(Ycal[:10, :], axis=0)
print(-np.sort(Ymin)[int(m*yps[0])], -np.sort(Ymin)[int(m*yps[1])])

Ymin = np.min(Ycal[:5, :], axis=0)
print(-np.sort(Ymin)[int(m*yps[0])], -np.sort(Ymin)[int(m*yps[1])])