#demo of excess severity for pareto-exponential/pareto-gamma import pandas as pd import matplotlib.pyplot as plt import numpy as np from scipy.integrate import trapz # define functions of integrand for numerical integration computation def integrand(x, a, b): return a*x/((x+b)*(a+np.log(x/b+1))**2) # define CDF and survival functions def ce(x,a,b): return 1-a/(a+np.log(x/b+1)) def se(x,a,b): return a/(a+np.log(x/b+1)) # similarly, define functions for pareto-gamma (implement by remove '#') #def se(x,a,b,c): # return a**c/(a+np.log(x/b+1))**c #def integrand(x, a, b, c): # return (c*x*se(x,a,b,c)**(1+1/c))/(a*x+a*b) #def ce(x,a,b,c): # return 1-se(x,a,b,c) # numerical experiments by different inputs M = 5500 a = 5 aa = 10 aaa = 15 b = 10 bb = 25 bbb = 50 p = np.zeros(200) q = np.zeros(200) r = np.zeros(200) s = np.linspace(35, 215, 200) for i in range(0,200): t = np.linspace(35+i,M, 200) I = (trapz(integrand(t,a,b),t) - se(i+35,a,b)*(i+35)+se(M,a,b)*M)/(1-ce(i+35,a,b)) h = (trapz(integrand(t,aa,b),t) - se(i+35,a,bb)*(i+35)+se(M,a,bb)*M)/(1-ce(i+35,a,bb)) k = (trapz(integrand(t,aaa,b),t) - se(i+35,a,bbb)*(i+35)+se(M,a,bbb)*M)/(1-ce(i+35,a,bbb)) p[i] = I q[i] = h r[i] = k plt.plot(s,p,'b-') plt.plot(s,q,'r-') plt.plot(s,r,'g-') plt.xlabel("Thresholds (a)") plt.ylabel("Mean Excess Severity") plt.title("Excess Severity") plt.legend(["lambda =5", "lambda= 10", "lambda = 15"],loc ="lower right") plt.show()