Expectile (#9)

* implemented expectile need tests

* woops left the test cases commented in the last commit

* fix dumb mistakes

* typos

* tests for expectation

* ignore Manifest.toml files

* ignore Manifest.toml since it is machine specific

* fix docs

* Test that the expectile is coherent for \alpha <= 0.5

* remove _e and change the direction of alpha

* change wording on return values

* update index and remove mention of _e methods

* use sortperm instead of sort eachindex

* remove mention of _e methods in risk measures

* fix tests :)

* fix documentation

* updated var risk level

* arg min instead of min

---------

Co-authored-by: Marek Petrik <[email protected]>

.gitignore

@@ -0,0 +1 @@
+Manifest.toml

README.md

@@ -10,14 +10,15 @@ RiskMeasures
 
 Julia library for computing risk measures for random variables. The random variable represents *profits* or *rewards* that are to be maximized. Also the computed risk value is preferable when it is greater.
 
-All risk measures get more conservative with an *increasing* risk level alpha.
+All risk measures, except ERM, are non-decreasing in risk level alpha. The ERM is non-increasing in level beta.
 
 The following risk measures are currently supported
 
 - VaR: Value at risk
 - CVaR: Conditional value at risk
 - ERM: Entropic risk measure
 - EVaR: Entropic value at risk
+- expectile: Expectile
 
 The focus is currently on random variables with categorical (discrete) probability distributions, but continuous probabilty distributions may be supported in the future too. 
 
@@ -36,6 +37,7 @@ VaR(x̃, 0.1)   # value at risk
 CVaR(x̃, 0.1)  # conditional value at risk
 EVaR(x̃, 0.1)  # entropic value at risk
 ERM(x̃, 0.1)   # entropic risk measure
+expectile(x̃, 0.1)   # entropic risk measure
 ```
 
 We can also compute risk measures of transformed random variables
@@ -45,14 +47,15 @@ VaR(5*x̃ + 10, 0.1)   # value at risk
 CVaR(x̃ - 10, 0.1)  # conditional value at risk
 ```
 
-Extended methods `VaR_e`, `CVaR_e`, and `EVaR_e` also return additional statistics and values, such as the distribution that attains the risk value and the optimal `β` in EVaR.
-
 Please see the unit tests for examples of how this package can be used to compute the risk. 
 
 ## Future development plans:
 
 - Analytical computation for special distributions, like Normal and others
 - Add an optional intergration with Mosek's exponential cones to support computation of EVaR. 
+- Coquet capacity risk measures
+- General risk measure construction from utility functions, such as CE, OCE, utility shortfall risk measures. 
+- Phi-divergence risk mesures for any phi-divergence function
 
 ## See Also
 

docs/src/index.md

@@ -35,6 +35,7 @@ VaR(x̃, 0.1)   # value at risk
 CVaR(x̃, 0.1)  # conditional value at risk
 EVaR(x̃, 0.1)  # entropic value at risk
 ERM(x̃, 0.1)   # entropic risk measure
+expectile(x̃, 0.1)  # expectile
 ```
 
 We can also compute risk measures of transformed random variables
@@ -44,7 +45,7 @@ VaR(5*x̃ + 10, 0.1)   # value at risk
 CVaR(x̃ - 10, 0.1)  # conditional value at risk
 ```
 
-Extended methods `VaR_e`, `CVaR_e`, and `EVaR_e` also return additional statistics and values, such as the distribution that attains the risk value and the optimal `β` in EVaR.
+Many methods methods `VaR`, `CVaR`, and `EVaR` also return additional statistics and values, such as the distribution that attains the risk value and the optimal `β` in EVaR. These are returned as named tuples.
 
 ## See Also
 
@@ -59,32 +60,19 @@ Extended methods `VaR_e`, `CVaR_e`, and `EVaR_e` also return additional statisti
 VaR
 ```
 
-```@docs
-VaR_e
-```
-
-
 ## Conditional Value at Risk
 
 ```@docs
 CVaR
 ```
 
-```@docs
-CVaR_e
-```
-
 ## Entropic Value at Risk
 
 
 ```@docs
 EVaR
 ```
 
-```@docs
-EVaR_e
-```
-
 ## Entropic Risk Measure
 
 ```@docs
@@ -95,12 +83,14 @@ ERM
 softmin
 ```
 
-## Essential Infimum
+## Expectile
 
 ```@docs
-essinf
+expectile
 ```
 
+## Essential Infimum
+
 ```@docs
-essinf_e
+essinf
 ```

src/RiskMeasures.jl

@@ -3,19 +3,21 @@ module RiskMeasures
 include("general.jl")
 
 include("essinf.jl")
-export essinf, essinf_e
+export essinf, essinf
 
 include("erm.jl")
 export ERM, softmin
 
 include("var.jl")
-export VaR, VaR_e
+export VaR, VaR
 
 include("cvar.jl")
-export CVaR, CVaR_e
-    
+export CVaR, CVaR
+
 include("evar.jl")
-export EVaR, EVaR_e
+export EVaR, EVaR
 
+include("expectile.jl")
+export expectile
 
 end

src/cvar.jl

@@ -5,15 +5,16 @@ using Distributions
 
 Compute the conditional value at risk at level `α` for the random variable `x̃`.
 
-The risk level `α` must satisfy the ``α ∈ [0,1]``. Risk aversion increases with an increasing `α` and, `α = 0` represents the expectation,  `α = 1` computes the essential infimum (smallest value with positive probability).
+The risk level `α` must satisfy the ``α ∈ [0,1]``. Risk aversion descreses with an increasing `α` and, `α = 1` represents the expectation,
+`α = 0` computes the essential infimum (smallest value with positive probability).
 
 Assumes a reward maximization setting and solves the dual form
 ```math
 \\min_{q ∈ \\mathcal{Q}} q^T x̃
 ```
 where
 ```math
-\\mathcal{Q} = \\left\\{q ∈ Δ^n : q_i ≤ \\frac{p_i}{1-α}\\right\\}
+\\mathcal{Q} = \\left\\{q ∈ Δ^n : q_i ≤ \\frac{p_i}{α}\\right\\}
 ```
 and ``Δ^n`` is the probability simplex, and ``p`` is the distribution of ``x̃``. 
 
@@ -23,7 +24,7 @@ More details: https://en.wikipedia.org/wiki/Expected_shortfall
 function CVaR end
 
 """
-    CVaR_e(x̃, α)
+    CVaR(x̃, α)
 
 Compute the conditional value at risk at level `α` for the random variable `x̃`. Also
 compute the equivalent random variable with the same support but a different distribution.
@@ -32,61 +33,56 @@ Returns a named tuple with `value` and the `solution` random variable that solve
 CVaR formulation. That is if `ỹ` is the solution, then the support of `x̃` and `ỹ` are the same
 and  ``\\mathbb{E}[ỹ] = \\operatorname{CVaR}_{α}[x̃]``.
 """
-function CVaR_e end
+function CVaR end
 
 
 
 """
-    CVaR_e(values, pmf, α; ...) 
+    CVaR(values, pmf, α; ...) 
 
 Compute CVaR for a discrete random variable with `values` and the probability mass
 function `pmf`. See `CVaR(x̃, α)` for more details.
 """
-function CVaR_e(values::AbstractVector{<:Real}, pmf::AbstractVector{<:Real}, α::Real;
-                check_inputs = true) 
-    _check_α(α)
+function CVaR(values::AbstractVector{<:Real}, pmf::AbstractVector{<:Real}, α::Real; check_inputs=true)
+    check_inputs && _check_α(α)
     check_inputs && _check_pmf(values, pmf)
-    
+
     T = eltype(pmf)
 
     # handle special cases
-    if iszero(α)
-        return (value = values'* pmf, pmf = Vector(pmf))
-    elseif isone(α)
-        minval = essinf_e(values, pmf; check_inputs = false)
+    if isone(α)
+        return (value=values' * pmf, pmf=Vector(pmf))
+    elseif iszero(α)
+        minval = essinf(values, pmf; check_inputs=false)
         minpmf = zeros(T, length(pmf))
         minpmf[minval.index] = one(T)
-        return (value = minval.value, pmf = minpmf)
+        return (value=minval.value, pmf=minpmf)
     end
-    
+
     # Here on: α ∈ (0,1)
     pc = zeros(T, length(pmf))  # this is the new distribution
     value = zero(T)                  # CVaR value
     p_left = one(T)           # probabilities left for allocation
-    α̂ = one(α) - α                      # probabilities to allocate
+    α̂ = α                      # probabilities to allocate
 
     # Efficiency note: sorting by values is O(n*log n);
     # quickselect is O(n) and would suffice but would need be based on quantile
-    sortedi = sort(eachindex(values, pmf); by=(i->@inbounds values[i]))
+    sortedi = sortperm(values)
     @inbounds for i ∈ sortedi
         # update index's probability and probability left to sum to 1.0
         increment = min(pmf[i] / α̂, p_left)
         # update  return values
         pc[i] = increment
         value += increment * values[i]
         p_left -= increment
-        p_left ≤ zero(p_left) && break 
+        p_left ≤ zero(p_left) && break
     end
-    return (value = value, pmf = pc)
+    return (value=value, pmf=pc)
 end
 
-# Definition for DiscreteNonparametric
-CVaR(x̃, α::Real; kwargs...) =
-    CVaR_e(rv2pmf(x̃)..., α; kwargs...).value
-
-function CVaR_e(x̃, α::Real; kwargs...)
+function CVaR(x̃, α::Real; kwargs...)
     supp, pmf = rv2pmf(x̃)
-    v1 = CVaR_e(supp, pmf, α; kwargs...)
+    v1 = CVaR(supp, pmf, α; kwargs...)
     ỹ = DiscreteNonParametric(supp, v1.pmf)
-    (value = v1.value, solution = ỹ)
+    (value=v1.value, pmf=ỹ)
 end

src/erm.jl

@@ -17,23 +17,23 @@ More details: https://en.wikipedia.org/wiki/Entropic_risk_measure
 function ERM end
 
 function ERM(values::AbstractVector{<:Real}, pmf::AbstractVector{<:Real}, β::Real;
-             x̃min::Real = -Inf, check_inputs = true)
-    
+    x̃min::Real=-Inf, check_inputs=true)
+
     check_inputs && _check_pmf(values, pmf)
     β < zero(β) && _bad_risk("Risk level β must be non-negative.")
 
     if iszero(β)
-        return values' * pmf 
+        return values' * pmf
     elseif isinf(β) && β > zero(β)
-        return essinf_e(values, pmf; check_inputs = false).value
-    end 
+        return essinf(values, pmf; check_inputs=false).value
+    end
     # because entropic risk measure is translation equivariant, we can change values
     # so that it is positive. That change makes it less likely that it overflows
     x̃min = isfinite(x̃min) ? x̃min : minimum(values)
-    
+
     #@fastmath x̃min-one(β) / β*log(sum(pmf .* exp.(-β .* (values .- x̃min) ))) |> float
-    a = Iterators.map((x,p) -> (@fastmath p * exp(-β * (x - x̃min))), values, pmf) |> sum
-    @fastmath x̃min - one(β) / β * log(a) 
+    a = Iterators.map((x, p) -> (@fastmath p * exp(-β * (x - x̃min))), values, pmf) |> sum
+    @fastmath x̃min - one(β) / β * log(a)
 end
 
 ERM(x̃, β::Real; kwargs...) = ERM(rv2pmf(x̃)..., β; kwargs...)
@@ -58,17 +58,17 @@ The value β must be positive
 function softmin end
 
 function softmin(values::AbstractVector{<:Real}, pmf::AbstractVector{<:Real},
-                 β::Real; x̃min::Real = -Inf, check_inputs = true) 
-    
+    β::Real; x̃min::Real=-Inf, check_inputs=true)
+
     check_inputs && _check_pmf(values, pmf)
 
     if β > zero(β)
         # TODO: only needs to look at values with pos prob.
         # because softmin is translation invariant add the subtract the smallest value
         # ensures that values ≥ 0
-        if isfinite(x̃min) 
+        if isfinite(x̃min)
             np = similar(pmf)
-            np .= @fastmath pmf .* exp.(-β .* (values .- x̃min) )
+            np .= @fastmath pmf .* exp.(-β .* (values .- x̃min))
             np .= @fastmath inv(sum(np)) .* np
             mapreduce(isfinite, &, np) || error("Overflow, reduce β.")
             np

src/essinf.jl

@@ -7,14 +7,14 @@ value with positive probability
 function essinf end
 
 """
-    essinf_e(values, pmf; ...)
+    essinf(values, pmf; ...)
 
 Compute the value for a random variable with `values` and the probability mass
 function `pmf`.
 
-See `essinf_e(x̃)` for more details.
+See `essinf(x̃)` for more details.
 """
-function essinf_e(values::AbstractVector{Tval}, pmf::AbstractVector{<:Real};
+function essinf(values::AbstractVector{Tval}, pmf::AbstractVector{<:Real};
                   check_inputs = true) :: @NamedTuple{value::Tval, index::Int} where
     {Tval <: Real}