''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Copyright Jan Iwanik 2014-2017. All rights reserved.
' Unauthorised copying is prohibited.
' Under no circumstances can the software be considered a deliverable of any
' contract for services between you and Iwanik Ltd. The software is a separate
' proprietary tool developed and used by Iwanik Ltd. You should delete any
' copies of the software if your contract for services with Iwanik Ltd has finished.
' You must not decrypt or remove access protection from this software.
' The software is provided 'as is' and we take not responsibility for the results of
' your work if it is completed using this software.


Option Explicit
 
 
' MIGHT NEED TO BE TESTED....


''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'' Poisson
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

'copied from http://www.actuarialoutpost.com/actuarial_discussion_forum/showthread.php?t=7730
' DELETED

' Returns the smallest integer N such that POISSON(N, Mean, True) >= Prob
Public Function PoissonInvCDF(prob As Double, Mean As Double) As Double
  ' E.g. POISSON(5, 10, TRUE) returns 0.067085962
  ' then PoissonInvCDF(0.067085962, 10) returns 5
  Dim n             As Long     ' number of events
  Dim CDF           As Double   ' cumulative distribution function at N
 
  ' These two variables are used to simplify the probability mass
  ' function summation CDF = CDF + Exp(-Mean) * Mean ^ N / N!
  Dim dExpMean      As Double   ' =Exp(-Mean)
  Dim dK            As Double   ' incremental power & factorial

  ' Return a value error for invalid parameter values
  If prob <= 0 Or prob >= 1 Or Mean <= 0 Then
    PoissonInvCDF = CVErr(xlErrValue)
    
  ' Return a NA error if Mean > 100
  ElseIf Mean > 100 Then
    PoissonInvCDF = CVErr(xlErrNA)

  ' Calculate function
  Else
    dExpMean = Exp(-Mean)
    dK = 1#
    CDF = dExpMean

    Do While CDF < prob
      n = n + 1
      dK = dK * Mean / n
      CDF = CDF + dExpMean * dK
    Loop

    PoissonInvCDF = n
  End If
End Function


' Random variable from Poisson distribution (jii)
Public Function PoissonRand(Lambda As Double) As Long
    Dim this_rand As Double
    Dim p0 As Double, cumm As Double, lambdaKK As Double
    p0 = Exp(-Lambda)
    lambdaKK = 1
    PoissonRand = 0
    cumm = 0
    this_rand = Rnd()
    
    Do While True
        cumm = cumm + p0 * lambdaKK
        If cumm >= this_rand Then Exit Function
        ' next step
        PoissonRand = PoissonRand + 1
        lambdaKK = lambdaKK * Lambda / PoissonRand
    Loop
End Function



''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'' Negative Binomial
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

' Returns probability that X = n for n = 0,1,2,3...
' This function generalizes the NegBinomial distribution to allow for r that is non-integer.
' Inputs
'    mean         the expected value of the distribution = EX
'    q_panjer       = variance / mean
' This function uses the (a, b, 0) formulation of the negative binomial, avoiding the need for the integer R (number of failures)
Function NegBinomPDF(n As Integer, Mean As Double, Q_Panjer As Double) As Double
    Dim r As Double, Beta As Double            ' (a, b, 0) parameters
    Dim PDF            As Double               ' probability distribution
    Dim nn             As Integer
 
    If Q_Panjer = 1 Then                       ' Poisson
        NegBinomPDF = Application.WorksheetFunction.Poisson_Dist(n, Mean, False)
    Else
        Beta = (Q_Panjer - 1)                      ' transform from model paramters
        r = Mean / Beta                            ' transform from model paramters
        
        If n = 0 Then
            If Q_Panjer = 1 Then
                PDF = Exp(-Mean)                        ' for n=0 under Poisson
            Else
                PDF = (1 + Beta) ^ (-r)                 ' for n=0
            End If
        Else
            PDF = r * Beta / (1 + Beta) ^ (r + 1)   ' for n=1
            If n > 1 Then
                For nn = 2 To n                     ' for n=2+
                    PDF = (r + nn - 1) * Beta / (nn * (1 + Beta)) * PDF
                Next
            End If
        End If
    NegBinomPDF = PDF
    End If
End Function


' Inverse CDF of negoative binomial
' This function uses the (a, b, 0) formulation of the negative binomial, avoiding the need for the integer R (number of failures)
' panjer_q = variance / mean
Public Function NegBinomInvCDF(prob As Double, Mean As Double, Q_Panjer As Double) As Double
    Dim r              As Double           ' (a, b, 0) parameters
    Dim Beta           As Double           ' (a, b, 0) parameters
    Dim CDF            As Double           ' cumulative probability
    Dim PDF            As Double           ' probability distribution
    Dim n              As Integer          ' observed frequency
 
    If Q_Panjer = 1 Then                       ' Poisson
        NegBinomInvCDF = PoissonInvCDF(prob, Mean)
    Else
        Beta = (Q_Panjer - 1)                  ' transform from model paramters
        r = Mean / Beta                        ' transform from model paramters
        n = 0
        PDF = (1 + Beta) ^ (-r)                   ' for n=0
        CDF = PDF
     
        If CDF >= prob Then
            NegBinomInvCDF = n
        Else
            n = 1
            PDF = r * Beta / (1 + Beta) ^ (r + n)  ' for n=1
            CDF = CDF + PDF
        
            If CDF >= prob Then
                NegBinomInvCDF = n
            Else
                Do While CDF < prob
                    n = n + 1
                    PDF = (r + n - 1) * Beta / (n * (1 + Beta)) * PDF
                    CDF = CDF + PDF
                    If CDF >= prob Then
                        NegBinomInvCDF = n
                        Exit Do
                    End If
                Loop
            End If
        End If
    End If
End Function


' Random variable from Negative Binomial distribution - jii
' Inputs:
'   mean - the expected value of this distribution
'   p    - basic p parameter, relates to mean and variance relationship 1 - p = mean / variance
Public Function NegBinomRand(Mean As Double, panj_q As Double) As Long
    Dim rand_x As Double
    rand_x = Rnd()
    NegBinomRand = NegBinomInvCDF(rand_x, Mean, panj_q)
End Function






''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'' Normal, Gaussian
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

' Alternatives:
' based on the following article: http://www.johndcook.com/blog/normal_cdf_inverse/
' a c# version can be found here: http://www.johndcook.com/blog/csharp_phi_inverse/
' based on the following article: http://www.johndcook.com/blog/python_phi/
' a c# version can be found here: http://www.johndcook.com/blog/csharp_phi/

Public Function StandardNormalCDF(x As Double) As Double
    StandardNormalCDF = Application.WorksheetFunction.Norm_S_Dist(x, True)
End Function

Public Function StandardNormalInvCDF(prob As Double) As Double
    StandardNormalInvCDF = Application.WorksheetFunction.Norm_S_Inv(prob)
End Function




''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'' Log-Normal
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

' Probability F(x) for log-normal
Public Function LognormCDF(x As Double, mu As Double, sigma As Double) As Double
    LognormCDF = StandardNormalCDF((Log(x) - mu) / sigma)
End Function

' This function returns the inverse CDF for a lognormal distribtion given Prob
Public Function LognormInvCDF(prob As Double, mu As Double, sigma As Double)
    ' Return a value error for an invalid parameter values
    If prob <= 0 Or prob >= 1 Or mu <= 0 Or sigma <= 0 Then
        LognormInvCDF = CVErr(xlErrValue)
    ' Calculate function
    Else
        LognormInvCDF = Exp(mu + sigma * StandardNormalInvCDF(prob))
    End If
End Function

' Capped Expected Value of log-normal
Public Function LognormCappedEX(cap As Double, mu As Double, sigma As Double) As Double
    LognormCappedEX = Exp(mu + sigma * sigma / 2) * StandardNormalCDF((Log(cap) - mu - sigma * sigma) / sigma) + cap - cap * StandardNormalCDF((Log(cap) - mu) / sigma)
End Function




''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'' Generalized Pareto Distribution (GPD)
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

' Inverted F() for Generalised Pareto Distribution, for random variable generation
' From Wikipedia (Which matches Parodi book), https://en.wikipedia.org/wiki/Generalized_Pareto_distribution
'           mu      = Location
'           sigma   = Scale
'           xi      = Shape
Public Function GeneralizedParetoInvCDF(prob As Double, mu As Double, _
                                 sigma As Double, xi As Double) As Double
    ' Return a value error for an invalid parameter values
    If prob <= 0 Or prob >= 1 Or mu <= 0 Or sigma <= 0 Then
        GeneralizedParetoInvCDF = CVErr(xlErrValue)
    ' Calculate function
    Else
        ' Note that form for generating GPD's in Wikipedia inverts the Prob term. No harm in isolation but we needed to adjust formula below, from Prob to 1-Prob.
        GeneralizedParetoInvCDF = mu + sigma * (((1 - prob) ^ (-xi)) - 1) / xi
    End If
End Function

Public Function GeneralizedParetoCDF(x As Double, mu As Double, sigma As Double, xi As Double) As Double
    ' Return a value error for an invalid parameter values
    If mu <= 0 Or sigma <= 0 Then
        GeneralizedParetoCDF = CVErr(xlErrValue)
    ' Calculate function
    Else
        If x <= mu Then GeneralizedParetoCDF = 0: Exit Function
        If xi <> 0 Then GeneralizedParetoCDF = 1 - (1 + xi * (x - mu) / sigma) ^ (-1 / xi)
        If xi = 0 Then GeneralizedParetoCDF = 1 - Exp(-(x - mu) / sigma)
    End If
End Function





''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'' Log-Normal with GPD tail
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

' jii: new function
' Returns the CDF of combined lognormal - GPD distribution function
Public Function LognormalWithGPDTailCDF(x As Double, _
                                 Location_LN As Double, Scale_LN As Double, _
                                 weight_tail As Double, _
                                 Location_GPD As Double, Scale_GPD As Double, Shape_GPD As Double) As Double
    ' Return a value error for an invalid parameter values
    If Location_LN <= 0 Or Scale_LN <= 0 Or Location_GPD <= 0 Or Scale_GPD <= 0 Or weight_tail < 0 Or weight_tail > 1 Then _
        LognormalWithGPDTailCDF = CVErr(xlErrValue): Exit Function
    If x <= 0 Then LognormalWithGPDTailCDF = 0: Exit Function
        
    If x <= Location_GPD Then
        ' log-normal area, need to rescale probablity to relect the conditional lognormal distribution under the Location_GPD threshold
        LognormalWithGPDTailCDF = StandardNormalCDF((Log(x) - Location_LN) / Scale_LN)
        LognormalWithGPDTailCDF = LognormalWithGPDTailCDF / StandardNormalCDF((Log(Location_GPD) - Location_LN) / Scale_LN) * (1 - weight_tail)
    Else
        If Shape_GPD <> 0 Then LognormalWithGPDTailCDF = 1 - (1 + Shape_GPD * (x - Location_GPD) / Scale_GPD) ^ (-1 / Shape_GPD)
        If Shape_GPD = 0 Then LognormalWithGPDTailCDF = 1 - Exp(-(x - Location_GPD) / Scale_GPD)
        LognormalWithGPDTailCDF = (1 - weight_tail) + LognormalWithGPDTailCDF * weight_tail
    End If
End Function


' jii: change this function
' This function returns the severity of a combined distribution with LogNormal body and Generalised Pareto tail
' Extra parameter weight_tail. If weight_tail = 10% then one in 10 losses will be from the tail.
Public Function LognormWithGPDTailInvCDF(prob As Double, _
                                     Location_LN As Double, Scale_LN As Double, _
                                     weight_tail As Double, _
                                     Location_GPD As Double, Scale_GPD As Double, Shape_GPD As Double) As Double
    ' This function returns the severity of a combined distribution with (shifted) LogNormal body and Generalised Pareto tail
    ' Extra parameter weight_tail. If weight_tail = 10% then one in 10 losses will be from the tail.

    ' Return a value error for an invalid parameter values
    If prob < 0 Or prob >= 1 Or Location_LN <= 0 Or Scale_LN <= 0 Or Location_GPD <= 0 Or Scale_GPD <= 0 Or weight_tail < 0 Or weight_tail > 1 Then _
        LognormWithGPDTailInvCDF = CVErr(xlErrValue): Exit Function
    
    ' Calculate function
    If prob <= 1 - weight_tail Then
        ' log-normal area, need to rescale probablity to relect the conditional lognormal distribution under the Location_GPD threshold
        prob = prob / (1 - weight_tail) * StandardNormalCDF((Log(Location_GPD) - Location_LN) / Scale_LN)
        LognormWithGPDTailInvCDF = LognormInvCDF(prob, 0, Location_LN, Scale_LN) ' Return Lognormal for body of distribution
    Else
        ' tail of the distribution area
        prob = (prob - (1 - weight_tail)) / weight_tail ' Rebase probability once after GPD threshold
        LognormWithGPDTailInvCDF = GeneralizedParetoInvCDF(prob, Location_GPD, Scale_GPD, Shape_GPD) ' Return GPD for tail of distribution
    End If
End Function



' jii: new function
' This function returns the severity (inverse CDF) of a combined distribution with LogNormal body and Generalised Pareto tail
' conditional that loss exceedes a certain amount - useful in practical applications
' Input:
'         x_prob                  - the point of CDF inversion
'         cond_above_this_only    - losses below this number are not considered, this number must be less than GPD threshold.
'                                   the log-normal parameters are from-ground-up but the inverted CFD being returned is of a
'                                   conditional distribution conditional on X > cond_above_this_only, so this function will always
'                                   return numbers which are greater than cond_above_this_only
'         mu_LN                   - main parameter of the standard from-ground-up log-normal distribution
'         sigma_LN                - second parameter of the standard from-ground-up log-normal distribution
'         weight_tail             - required probability that a number returnd by this function will be from the GPD tail and not from the
'                                   conditional log-norm
'         mu_GPD                  - main parameter of the GPD tail
'         sigma_GPD               - second parameter of the GPD tail
'         xi_GPD                  - shape parameter of the GPD tail
'
'Private Function LognormWithGPDTailCondAboveInvCDF_DEPRECIATED(x_prob As Double, _
'                                     cond_above_this_only As Double, _
'                                     mu_ln As Double, sigma_ln As Double, _
'                                     weight_tail As Double, _
'                                     mu_gpd As Double, sigma_gpd As Double, xi_gpd As Double) As Double
'    Dim prob_above_this_only As Double
'    prob_above_this_only = 1 - LognormalWithGPDTailCDF(cond_above_this_only, mu_ln, sigma_ln, weight_tail, mu_gpd, sigma_gpd, xi_gpd)
'
'    ' Return a value error for an invalid parameter values
'     If x_prob < 0 Or x_prob >= 1 Or mu_ln <= 0 Or sigma_ln <= 0 Or mu_gpd <= 0 Or sigma_gpd <= 0 Or weight_tail < 0 Or weight_tail > 1 Or cond_above_this_only < 0 Then _
'        LognormWithGPDTailCondAboveInvCDF_DEPRECIATED = CVErr(xlErrValue): Exit Function
'
'    ' un-conditionalise
'    x_prob = x_prob * prob_above_this_only + (1 - prob_above_this_only)
'    LognormWithGPDTailCondAboveInvCDF_DEPRECIATED = LognormWithGPDTailInvCDF(x_prob, mu_ln, sigma_ln, prob_above_this_only * weight_tail, mu_gpd, sigma_gpd, xi_gpd)
'End Function


' (c) Jan Iwanik 2014-2019, part of proInsurance package
' Used this for severity distributions.
' Combined 1) conditional LogNormal body + 2) Generalised Pareto tail
' Conditional losses must exceeded a certain amount - useful in practical applications
' Input:
'         x_prob                  - the point of CDF inversion
'         cond_above_this_only    - losses below this number are not considered, this number must be less than GPD threshold.
'                                   the log-normal parameters are from-ground-up but the inverted CFD being returned is of a
'                                   conditional distribution conditional on X > cond_above_this_only, so this function will always
'                                   return numbers which are greater than cond_above_this_only
'         mu_LN                   - main parameter of the standard from-ground-up log-normal distribution
'         sigma_LN                - second parameter of the standard from-ground-up log-normal distribution
'         weight_tail             - required probability that a number returned by this function will be from the GPD tail and not from the
'                                   conditional log-norm
'         mu_GPD                  - main parameter of the GPD tail
'         sigma_GPD               - second parameter of the GPD tail
'         xi_GPD                  - shape parameter of the GPD tail
'
Public Function LognormWithGPDTailCondAboveInvCDF(x_prob As Double, _
                                  cond_above_this_only As Double, _
                                  mu_ln As Double, sigma_ln As Double, _
                                  weight_tail As Double, _
                                  mu_gpd As Double, sigma_gpd As Double, xi_gpd As Double) As Double
    If cond_above_this_only > mu_gpd Then LognormWithGPDTailCondAboveInvCDF = -1: Exit Function
    Dim p1 As Double, p2 As Double
    
    ' un-conditionalise
    If x_prob <= 1 - weight_tail Then ' body log normal
        p1 = Application.WorksheetFunction.Norm_Dist((Log(cond_above_this_only) - mu_ln) / sigma_ln, 0, 1, True)
        p2 = Application.WorksheetFunction.Norm_Dist((Log(mu_gpd) - mu_ln) / sigma_ln, 0, 1, True)
        x_prob = x_prob / (1 - weight_tail) * (p2 - p1) + p1
        LognormWithGPDTailCondAboveInvCDF = Exp(mu_ln + sigma_ln * Application.WorksheetFunction.Norm_S_Inv(x_prob))
    Else ' tail GPD
        x_prob = (x_prob - 1 + weight_tail) / weight_tail
        LognormWithGPDTailCondAboveInvCDF = mu_gpd + sigma_gpd * (((1 - x_prob) ^ (0 - xi_gpd)) - 1) / xi_gpd
    End If
End Function


Public Function LognormWithGPDTailCondAboveCDF(x As Double, _
                                  cond_above_this_only As Double, _
                                  mu_ln As Double, sigma_ln As Double, _
                                  weight_tail As Double, _
                                  mu_gpd As Double, sigma_gpd As Double, xi_gpd As Double) As Double
    Dim p1 As Double, p2 As Double
    If x <= cond_above_this_only Then LognormWithGPDTailCondAboveCDF = 0: Exit Function
    If x <= mu_gpd Then
        p1 = LognormCDF(cond_above_this_only, mu_ln, sigma_ln)
        p2 = LognormCDF(mu_gpd, mu_ln, sigma_ln)
        LognormWithGPDTailCondAboveCDF = (LognormCDF(x, mu_ln, sigma_ln) - p1) / (p2 - p1) * (1 - weight_tail)
        Exit Function
    End If
    'tail
    LognormWithGPDTailCondAboveCDF = (1 - weight_tail) + weight_tail * GeneralizedParetoCDF(x, mu_gpd, sigma_gpd, xi_gpd)
End Function



' Capped expected value for the above distribution being a combined distribution with LogNormal body and Generalised Pareto tail
'
Public Function LognormWithGPDTailCondAboveCappedEX(cap As Double, _
                                     cond_above_this_only As Double, _
                                     mu_ln As Double, sigma_ln As Double, _
                                     weight_tail As Double, _
                                     mu_gpd As Double, sigma_gpd As Double, xi_gpd As Double) As Double
    Dim lgnorm_cap As Double, lgnorm_contrib As Double, lgnorm_exp_above2cap As Double, lgnorm_exp_above2cap_cond As Double, lgnorm_prob_undercap_cond As Double
    
    ' Error handling
    If mu_ln <= 0 Or sigma_ln <= 0 Or mu_gpd <= 0 Or sigma_gpd <= 0 Or weight_tail < 0 Or weight_tail > 1 Or cond_above_this_only < 0 Then LognormWithGPDTailCondAboveCappedEX = CVErr(xlErrValue): Exit Function
        
    ' simple case
    If cap <= cond_above_this_only Then
        LognormWithGPDTailCondAboveCappedEX = cap
        Exit Function
    End If
    
    ' capping within log-norm
    lgnorm_cap = WorksheetFunction.min(cap, mu_gpd)
    lgnorm_exp_above2cap = LognormCappedEX(lgnorm_cap, mu_ln, sigma_ln) - LognormCappedEX(cond_above_this_only, mu_ln, sigma_ln) _
                           - (lgnorm_cap - cond_above_this_only) * (1 - LognormCDF(lgnorm_cap, mu_ln, sigma_ln))
    If cond_above_this_only = lgnorm_cap Then
        lgnorm_exp_above2cap_cond = cond_above_this_only
        lgnorm_prob_undercap_cond = 0
    Else
        lgnorm_exp_above2cap_cond = cond_above_this_only + lgnorm_exp_above2cap / (LognormCDF(lgnorm_cap, mu_ln, sigma_ln) - LognormCDF(cond_above_this_only, mu_ln, sigma_ln))
        lgnorm_prob_undercap_cond = (LognormCDF(lgnorm_cap, mu_ln, sigma_ln) - LognormCDF(cond_above_this_only, mu_ln, sigma_ln)) / (LognormCDF(mu_gpd, mu_ln, sigma_ln) - LognormCDF(cond_above_this_only, mu_ln, sigma_ln))
    End If
    lgnorm_contrib = lgnorm_prob_undercap_cond * lgnorm_exp_above2cap_cond + (1 - lgnorm_prob_undercap_cond) * lgnorm_cap
    If cap < mu_gpd Then
        LognormWithGPDTailCondAboveCappedEX = lgnorm_contrib * (1 - weight_tail) + lgnorm_cap * weight_tail
        Exit Function
    End If
    
    ' capping within gpd
     If cap >= mu_gpd Then
        LognormWithGPDTailCondAboveCappedEX = (1 - weight_tail) * lgnorm_contrib + _
                                weight_tail * (mu_gpd + sigma_gpd / (1 - xi_gpd) - _
                                sigma_gpd / (1 - xi_gpd) * Application.WorksheetFunction.Power(1 + xi_gpd * (cap - mu_gpd) / sigma_gpd, 1 - 1 / xi_gpd))
     End If
End Function




''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'' MBBEFD
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

''''''''''''''''''''''''''''''''''''''''''''''''''''''
' Copyright Jan Iwanik 2014-2017. All rights reserved.
' Unauthorised copying is prohibited.
' Under no circumstances can the software be considered a deliverable of any
' contract for services between you and Iwanik Ltd. The software is a separate
' proprietary tool developed and used by Iwanik Ltd. You should delete any
' copies of the software if your contract for services with Iwanik Ltd has finished.
' You must not decrypt or remove access protection from this software.
' The software is provided 'as is' and we take not responsibility for the results of
' your work if it is completed using this software.


'' https://www.casact.org/library/astin/vol27no1/99.pdf
Public Function MbbefdCDF(x As Double, b As Double, g As Double) As Double
    If x = 1 Then MbbefdCDF = 1: Exit Function
    If g = 1 Or b = 0 Then MbbefdCDF = 0: Exit Function
    If b = 1 And g > 1 Then MbbefdCDF = 1 - 1 / (1 + (g - 1) * x): Exit Function
    If b * g = 1 And g > 1 Then MbbefdCDF = 1 - b ^ x: Exit Function
    If b > 0 And b <> 1 And b * g <> 1 And g > 1 Then _
        MbbefdCDF = 1 - (1 - b) / ((g - 1) * b ^ (1 - x) + (1 - b * g))
End Function

Public Function MbbefdInvCDF(prob As Double, b As Double, g As Double) As Double
    If prob = 1 Then MbbefdInvCDF = 1: Exit Function
    If prob >= 1 - 1 / g Then MbbefdInvCDF = 1: Exit Function
    If g = 1 Or b = 0 Then MbbefdInvCDF = 0: Exit Function
    If b = 1 And g > 1 Then MbbefdInvCDF = 1 / (g - 1) * (1 / (1 - prob) - 1)
    If b * g = 1 And g > 1 Then MbbefdInvCDF = Log(1 - prob) / Log(b): Exit Function
    If b > 0 And b <> 1 And b * g <> 1 And g > 1 Then _
        MbbefdInvCDF = 1 - Log((1 - b) / ((1 - prob) * (g - 1)) - (1 - b * g) / (g - 1)) / Log(b)
End Function



' Two parameter Swiss Re excess curve
Public Function MbbefdG(x As Double, b As Double, g As Double) As Double
    If g = 1 Or b = 0 Then
        MbbefdG = x
    ElseIf b = 1 And g > 1 Then
        MbbefdG = Log(1 + (g - 1) * x) / Log(g)
    ElseIf b * g = 1 And g > 1 Then
        MbbefdG = (1 - b ^ x) / (1 - b)
    ElseIf b > 0 And b <> 1 And b * g <> 1 And g > 1 Then
        MbbefdG = Log(((g - 1) * b + (1 - g * b) * (b ^ x)) / (1 - b)) / Log(g * b)
    End If
End Function


' One parameter Swiss Re excess curve
Public Function MbbefdGc(x As Double, c As Double) As Double
    MbbefdGc = MbbefdG(x, Exp(3.1 - 0.15 * (1 + c) * c), Exp((0.78 + 0.12 * c) * c))
End Function




' Two parameter Swiss Re excess curve
Public Function MbbefdEX(b As Double, g As Double) As Double
    If g = 1 Or b = 0 Then
        MbbefdEX = 1
    ElseIf b = 1 And g > 1 Then
        MbbefdEX = Log(g) / (g - 1)
    ElseIf b * g = 1 And g > 1 Then
        MbbefdEX = (b - 1) / Log(b)
    ElseIf b > 0 And b <> 1 And b * g <> 1 And g > 1 Then
        MbbefdEX = Log(b * g) * (1 - b) / (Log(b) * (1 - g * b))
    End If
End Function