January 01st, 2024

In [2]:

import numpy as np
from scipy.stats import lognorm
from scipy.integrate import quad
import pandas as pd

from math import exp
from math import log
from math import sqrt

In [6]:

def integrand(x, excess, limit):
    return min(limit, max(x - excess, 0)) * lognorm.pdf(x, s=sigma, scale=np.exp(mu))


means = []
cvs = []
outputs = []

for mean in [500,1000,1500,2000,2500,3000]:
    for cv in [0.05,0.1,0.15,0.2,0.25]:
        
        means.append(mean)
        cvs.append(cv)
        
        stddev = mean*cv
        mu = log(mean/(sqrt(1+stddev **2/mean**2)))
        sigma = sqrt(log(1+stddev**2/mean**2))
        excess = mean
        limit = lognorm.ppf(0.995, s=sigma, scale=np.exp(mu)) - mean
        result, _ = quad(integrand, excess, excess+limit, args=(excess, limit))
        outputs.append(result/limit)

# Create DataFrame
df = pd.DataFrame({
    'Mean': means,
    'Coefficient of Variation': cvs,
    'Output': outputs
})

df_pivot = df.pivot(index='Mean', columns='Coefficient of Variation', values='Output')

# Print the DataFrame
print(df_pivot)

Coefficient of Variation      0.05      0.10      0.15      0.20      0.25
Mean                                                                      
500                       0.140955  0.133088  0.125675  0.118723  0.112229
1000                      0.140955  0.133088  0.125675  0.118723  0.112229
1500                      0.140955  0.133088  0.125675  0.118723  0.112229
2000                      0.140955  0.133088  0.125675  0.118723  0.112229
2500                      0.140955  0.133088  0.125675  0.118723  0.112229
3000                      0.140955  0.133088  0.125675  0.118723  0.112229

In [ ]: