From 648ae127219c75dbe72d22238d63a4756a1073c5 Mon Sep 17 00:00:00 2001
From: Alec Loudenback <alecloudenback@gmail.com>
Date: Tue, 16 Apr 2024 23:16:34 -0500
Subject: [PATCH] WIP Risk Measure Revamp

---
 Project.toml            |   2 +-
 src/ActuaryUtilities.jl |   1 +
 src/risk_measures.jl    | 177 ++++++++++++++++++++++++++++++----------
 test/risk_measures.jl   |  34 +++++---
 4 files changed, 158 insertions(+), 56 deletions(-)

diff --git a/Project.toml b/Project.toml
index 21a2a4d..22140e1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -16,7 +16,7 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
 [compat]
-Distributions = "0.24,0.25"
+Distributions = "0.25"
 FinanceCore = "^2"
 FinanceModels = "^4"
 ForwardDiff = "^0.10"
diff --git a/src/ActuaryUtilities.jl b/src/ActuaryUtilities.jl
index 73e0d35..d5a2064 100644
--- a/src/ActuaryUtilities.jl
+++ b/src/ActuaryUtilities.jl
@@ -10,6 +10,7 @@ using MuladdMacro
 import FinanceModels
 import StatsBase
 using PrecompileTools
+import Distributions
 
 include("financial_math.jl")
 include("risk_measures.jl")
diff --git a/src/risk_measures.jl b/src/risk_measures.jl
index 472b0bb..7c8f45c 100644
--- a/src/risk_measures.jl
+++ b/src/risk_measures.jl
@@ -1,49 +1,136 @@
-"""
-    VaR(v::AbstractArray,p::Real;rev::Bool=false)
-
-The `p`th quantile of the vector `v` is the Value at Risk. Assumes more positive values are higher risk measures, so a higher p will return a more positive number, but this can be reversed if `rev` is `true`.
-
-Also can be called with `ValueAtRisk(...)`.
-"""
-function VaR(v, p; rev=false)
-    if rev
-        return StatsBase.quantile(v, 1 - p)
-    else
-        return StatsBase.quantile(v, p)
-    end
+abstract type DistortionFunction end
+
+
+struct Expectation <: DistortionFunction end
+(g::Expectation)(x) = x
+
+struct VaR{T} <: DistortionFunction
+    α::T
+end
+(g::VaR)(x) = x < (1 - g.α) ? 0 : 1
+
+struct CTE{T} <: DistortionFunction
+    α::T
+end
+(g::CTE)(x) = x < (1 - g.α) ? x / (1 - g.α) : 1
+
+struct WangTransform{T} <: DistortionFunction
+    α::T
+end
+function (g::WangTransform)(x)
+    Φ_inv(x) = Distributions.quantile(Distributions.Normal(), x)
+    Distributions.cdf(Distributions.Normal(), Φ_inv(x) + g.α)
 end
-"""
-[`VaR`](@ref)
-"""
-ValueAtRisk = VaR
-
-"""
-    CTE(v::AbstractArray,p::Real;rev::Bool=false)
-
-The average of the values ≥ the `p`th percentile of the vector `v` is the Conditiona Tail Expectation. Assumes more positive values are higher risk measures, so a higher p will return a more positive number, but this can be reversed if `rev` is `true`.
-
-May also be called with `ConditionalTailExpectation(...)`.
-
-Also known as Tail Value at Risk (TVaR), or Tail Conditional Expectation (TCE)
-"""
-function CTE(v, p; rev=false)
-    # filter has the "or approximately equalt to quantile" because
-    # of floating point path might make the quantile slightly off from the right indexing 
-    # e.g. if values should capture <= q, where q should be 10 but is calculated to be 
-    # 9.99999...
-    if rev
-        q = StatsBase.quantile(v, 1 - p)
-        filter = (v .<= q) .| (v .≈ q)
-    else
-        q = StatsBase.quantile(v, p)
-        filter = (v .>= q) .| (v .≈ q)
-    end
-
-    return sum(v[filter]) / sum(filter)
 
+struct DualPower{T} <: DistortionFunction
+    v::T
 end
+(g::DualPower)(x) = 1 - (1 - x)^g.v
+
+struct ProportionalHazard{T} <: DistortionFunction
+    y::T
+end
+(g::ProportionalHazard)(x) = x^(1 / g.y)
+
+
+function ρ(g::DistortionFunction, risk)
+
+    # integral from 0 to infinity of g(S(x))dx
+    # where S(x) is 1-cdf(risk,x)
+    F(x) = cdf_func(risk)(x)
+    S(x) = 1 - F(x)
+    H(x) = 1 - g(1 - x)
+    integral1, _ = quadgk(x -> 1 - H(F(x)), 0, Inf)
+    integral2, _ = quadgk(x -> H(F(x)), -Inf, 0)
+    return integral1 - integral2
+end
+
+
+cdf_func(S::AbstractArray{<:Real}) = StatsBase.ecdf(S)
+cdf_func(S::Distributions.UnivariateDistribution) = x -> Distributions.cdf(S, x)
+
+# """
+#     VaR(v::AbstractArray,p::Real;rev::Bool=false)
+
+# The `p`th quantile of the vector `v` is the Value at Risk. Assumes more positive values are higher risk measures, so a higher p will return a more positive number, but this can be reversed if `rev` is `true`.
+
+# Also can be called with `ValueAtRisk(...)`.
+# """
+# function VaR(v::T, p; sorted=false) where {T<:AbstractArray}
+#     if sorted
+#         _VaR_sorted(v, p)
+#     else
+#         _VaR_sorted(sort(v), p)
+#     end
+# end
+
+# # Core VaR assumes v is sorted
+# function _VaR_sorted(v, p)
+#     i = 1
+#     n = length(v)
+#     q_prior = 0.0
+#     x_prior = first(v)
+#     for (i, x) in enumerate(v)
+#         q = i / n
+#         if q >= p
+#             # return weighted between two points
+#             return x * (p - q_prior) / (q - q_prior) + x_prior * (q - p) / (q - q_prior)
+
+#         end
+#         x_prior = x
+#         q_prior = q
+#     end
+
+#     return last(v)
+# end
+
+
+# """
+# [`VaR`](@ref)
+# """
+# ValueAtRisk = VaR
+
+# """
+#     CTE(v::AbstractArray,p::Real;rev::Bool=false)
+
+# The average of the values ≥ the `p`th percentile of the vector `v` is the Conditiona Tail Expectation. Assumes more positive values are higher risk measures, so a higher p will return a more positive number, but this can be reversed if `rev` is `true`.
+
+# May also be called with `ConditionalTailExpectation(...)`.
+
+# Also known as Tail Value at Risk (TVaR), or Tail Conditional Expectation (TCE)
+# """
+# function CTE(v::T, p; sorted=false) where {T<:AbstractArray}
+#     if sorted
+#         _CTE_sorted(v, p)
+#     else
+#         _CTE_sorted(sort(v), p)
+#     end
+# end
+
+# # Core CTE assumes v is sorted
+# function _CTE_sorted(v, p)
+#     i = 1
+#     n = length(v)
+#     q_prior = 0.0
+#     x_prior = first(v)
+#     sub_total = zero(eltype(v))
+#     in_range = false
+#     for (i, x) in enumerate(v)
+#         q = i / n
+#         if in_range || q >= p
+#             # return weighted between two points
+#             # return x * (p - q_prior) / (q - q_prior) + x_prior * (q - p) / (q - q_prior)
+#             in_range = true
+
+#         end
+#         x_prior = x
+#         q_prior = q
+#     end
+
+#     return last(v)
+# end
 
-"""
-[`CTE`](@ref)
-"""
-ConditionalTailExpectation = CTE
+# """
+# [`CTE`](@ref)
+# """
+# ConditionalTailExpectation = CTE
diff --git a/test/risk_measures.jl b/test/risk_measures.jl
index afb84be..6d3279f 100644
--- a/test/risk_measures.jl
+++ b/test/risk_measures.jl
@@ -2,18 +2,32 @@
     sample = [0:100...]
 
     @testset "VaR" begin
-        @test VaR(sample,.9) ≈ 90
-        @test VaR(sample,.9,rev=true) ≈ 10
-        @test VaR(sample,.1) ≈ VaR(sample,1-.9)
-        @test VaR(sample,.1) == ValueAtRisk(sample,.1) 
+        @test VaR(sample, 0.9) ≈ 90
+        @test VaR(sample, 0.9, rev=true) ≈ 10
+        @test VaR(sample, 0.1) ≈ VaR(sample, 1 - 0.9)
+        @test VaR(sample, 0.1) == ValueAtRisk(sample, 0.1)
     end
-    
+
     @testset "CTE" begin
-        @test CTE(sample,.9) ≈ sum(90:100) / length(90:100)
-        @test CTE(sample,.1) ≈ sum(10:100) / length(10:100)
-        @test CTE(sample,.15) >= CTE(sample,.1) # monotonicity  
-        @test CTE(sample,.9,rev=true) ≈ sum(0:10) / length(0:10)
-        @test CTE(sample,.9) == ConditionalTailExpectation(sample,.9)
+        @test CTE(sample, 0.9) ≈ sum(90:100) / length(90:100)
+        @test CTE(sample, 0.1) ≈ sum(10:100) / length(10:100)
+        @test CTE(sample, 0.15) >= CTE(sample, 0.1) # monotonicity  
+        @test CTE(sample, 0.9, rev=true) ≈ sum(0:10) / length(0:10)
+        @test CTE(sample, 0.9) == ConditionalTailExpectation(sample, 0.9)
+    end
+
+
+    @testset "duplicated values" begin
+        sample = zeros(100)
+        sample[end] = 100
+
+        @test CTE(sample, 0) ≈ 100 / 100
+        @test CTE(sample, 0.5) ≈ 100 / 50
+
+        @test VaR(sample, 0) ≈ 0
+        @test VaR(sample, 0.5) ≈ 0
+        @test VaR(sample, 0.99) ≈ 0
+
     end
 
 end