Add Conditional Value at Risk (CVAR) metric

PRAS.jl/test/runtests.jl

@@ -1,15 +1,49 @@
 using PRAS
 using Test
 
-sys = PRAS.rts_gmlc()
+@testset "ShortfallResult" begin
+    sys = PRAS.rts_gmlc()
 
-sf, = assess(sys, SequentialMonteCarlo(samples=100), Shortfall())
+    sf, = assess(sys, SequentialMonteCarlo(samples=100), Shortfall())
 
-eue = EUE(sf)
-lole = LOLE(sf)
-neue = NEUE(sf)
+    eue = EUE(sf)
+    lole = LOLE(sf)
+    neue = NEUE(sf)
 
-@test val(eue) isa Float64
-@test stderror(eue) isa Float64
-@test val(neue) isa Float64
-@test stderror(neue) isa Float64
+    alpha = 0.95
+    cvar = CVAR(sf, alpha)
+    ncvar = NCVAR(sf, cvar)
+
+    @test val(eue) isa Float64
+    @test stderror(eue) isa Float64
+    @test val(neue) isa Float64
+    @test stderror(neue) isa Float64
+    @test val(cvar) isa Float64
+    @test stderror(cvar) isa Float64
+    @test val(ncvar) isa Float64
+    @test stderror(ncvar) isa Float64
+
+end
+
+@testset "ShortfallSamplesResult" begin
+    sys = PRAS.rts_gmlc()
+
+    sf, = assess(sys, SequentialMonteCarlo(samples=100), ShortfallSamples())
+
+    eue = EUE(sf)
+    lole = LOLE(sf)
+    neue = NEUE(sf)
+
+    alpha = 0.95
+    cvar = CVAR(sf, alpha)
+    ncvar = NCVAR(sf, cvar)
+
+    @test val(eue) isa Float64
+    @test stderror(eue) isa Float64
+    @test val(neue) isa Float64
+    @test stderror(neue) isa Float64
+    @test val(cvar) isa Float64
+    @test stderror(cvar) isa Float64
+    @test val(ncvar) isa Float64
+    @test stderror(ncvar) isa Float64
+end

PRASCore.jl/src/Results/Results.jl

@@ -4,7 +4,7 @@ import Base: broadcastable, getindex, merge!
 import OnlineStats: Series
 import OnlineStatsBase: EqualWeight, Mean, Variance, value
 import Printf: @sprintf
-import StatsBase: mean, std, stderror
+import StatsBase: mean, std, stderror, quantile
 
 import ..Systems: SystemModel, ZonedDateTime, Period,
                   PowerUnit, EnergyUnit, conversionfactor,
@@ -13,7 +13,7 @@ export
 
     # Metrics
     ReliabilityMetric, LOLE, EUE, NEUE,
-    val, stderror,
+    val, stderror, CVAR, NCVAR,
 
     # Result specifications
     Shortfall, ShortfallSamples,
@@ -80,6 +80,20 @@ NEUE(x::AbstractShortfallResult, r::AbstractString, ::Colon) =
 NEUE(x::AbstractShortfallResult, ::Colon, ::Colon) =
     NEUE.(x, x.regions.names, permutedims(x.timestamps))
 
+CVAR(x::AbstractShortfallResult, alpha::Float64, ::Colon, t::ZonedDateTime) =
+    CVAR.(x, alpha, x.regions.names, t)
+
+CVAR(x::AbstractShortfallResult, alpha::Float64, r::AbstractString, ::Colon) =
+    CVAR.(x, alpha, r, x.timestamps)
+
+CVAR(x::AbstractShortfallResult, alpha::Float64, ::Colon, ::Colon) =
+    CVAR.(x, alpha, x.regions.names, permutedims(x.timestamps))    
+
+NCVAR(x::AbstractShortfallResult, alpha::Float64, r::AbstractString, ::Colon) =
+    NCVAR.(x, alpha, r, x.timestamps)
+
+NCVAR(x::AbstractShortfallResult, alpha::Float64, ::Colon, ::Colon) =
+    NCVAR.(x, alpha, x.regions.names, permutedims(x.timestamps))
 include("Shortfall.jl")
 include("ShortfallSamples.jl")
 

PRASCore.jl/src/Results/Shortfall.jl

@@ -67,6 +67,10 @@ mutable struct ShortfallAccumulator{S} <: ResultAccumulator{Shortfall}
     unservedload_total_currentsim::Int
     unservedload_region_currentsim::Vector{Int}
 
+    # Sample-level UE for current simulation
+    unservedload_sample::Vector{Int}
+    unservedload_region_sample::Matrix{Int}
+
 end
 
 function accumulator(
@@ -90,14 +94,17 @@ function accumulator(
 
     unservedload_total_currentsim = 0
     unservedload_region_currentsim = zeros(Int, nregions)
+    unservedload_sample = zeros(Int, nsamples)
+    unservedload_region_sample = zeros(Int, nregions, nsamples)
 
     return ShortfallAccumulator{S}(
         periodsdropped_total, periodsdropped_region,
         periodsdropped_period, periodsdropped_regionperiod,
         periodsdropped_total_currentsim, periodsdropped_region_currentsim,
         unservedload_total, unservedload_region,
         unservedload_period, unservedload_regionperiod,
-        unservedload_total_currentsim, unservedload_region_currentsim)
+        unservedload_total_currentsim, unservedload_region_currentsim,
+        unservedload_sample, unservedload_region_sample)
 
 end
 
@@ -115,6 +122,9 @@ function merge!(
     foreach(merge!, x.unservedload_period, y.unservedload_period)
     foreach(merge!, x.unservedload_regionperiod, y.unservedload_regionperiod)
 
+    x.unservedload_sample .+= y.unservedload_sample
+    x.unservedload_region_sample .+= y.unservedload_region_sample
+
     return
 
 end
@@ -141,14 +151,17 @@ struct ShortfallResult{N, L, T <: Period, E <: EnergyUnit, S} <:
     eventperiod_regionperiod_mean::Matrix{Float64}
     eventperiod_regionperiod_std::Matrix{Float64}
 
-
     shortfall_mean::Matrix{Float64} # r x t
+    capacity_shortfall_mean::Matrix{Float64} # r x t
 
     shortfall_std::Float64
     shortfall_region_std::Vector{Float64}
     shortfall_period_std::Vector{Float64}
     shortfall_regionperiod_std::Matrix{Float64}
 
+    shortfall_samples::Vector{Float64}
+    shortfall_region_samples::Matrix{Float64}
+
     function ShortfallResult{N,L,T,E,S}(
         nsamples::Union{Int,Nothing},
         regions::Regions,
@@ -162,10 +175,13 @@ struct ShortfallResult{N, L, T <: Period, E <: EnergyUnit, S} <:
         eventperiod_regionperiod_mean::Matrix{Float64},
         eventperiod_regionperiod_std::Matrix{Float64},
         shortfall_mean::Matrix{Float64},
+        capacity_shortfall_mean::Matrix{Float64},
         shortfall_std::Float64,
         shortfall_region_std::Vector{Float64},
         shortfall_period_std::Vector{Float64},
-        shortfall_regionperiod_std::Matrix{Float64}
+        shortfall_regionperiod_std::Matrix{Float64},
+        shortfall_samples::Vector{Float64},
+        shortfall_region_samples::Matrix{Float64},
     ) where {N,L,T<:Period,E<:EnergyUnit,S <: Union{Shortfall,DemandResponseShortfall}}
 
         isnothing(nsamples) || nsamples > 0 ||
@@ -185,17 +201,20 @@ struct ShortfallResult{N, L, T <: Period, E <: EnergyUnit, S} <:
         size(eventperiod_regionperiod_std) == (nregions, N) &&
         length(shortfall_region_std) == nregions &&
         length(shortfall_period_std) == N &&
-        size(shortfall_regionperiod_std) == (nregions, N) ||
+        size(shortfall_regionperiod_std) == (nregions, N) &&
+        size(shortfall_samples) == (nsamples,) &&
+        size(shortfall_region_samples) == (nregions, nsamples) ||
             error("Inconsistent input data sizes")
 
         new{N,L,T,E,S}(nsamples, regions, timestamps,
             eventperiod_mean, eventperiod_std,
             eventperiod_region_mean, eventperiod_region_std,
             eventperiod_period_mean, eventperiod_period_std,
             eventperiod_regionperiod_mean, eventperiod_regionperiod_std,
-            shortfall_mean, shortfall_std,
+            shortfall_mean, capacity_shortfall_mean, shortfall_std,
             shortfall_region_std, shortfall_period_std,
-            shortfall_regionperiod_std)
+            shortfall_regionperiod_std, shortfall_samples,
+            shortfall_region_samples)
 
     end
 
@@ -292,6 +311,90 @@ function NEUE(x::ShortfallResult, r::AbstractString)
 
 end
 
+function CVAR(x::ShortfallResult{N,L,T,E}, alpha::Float64) where {N,L,T,E}
+
+    estimate = x.shortfall_samples
+    var = quantile(estimate, alpha)
+    tail_losses = estimate[estimate .>= var]
+
+    cvar = if !isempty(tail_losses)
+        MeanEstimate(tail_losses)
+    else
+        MeanEstimate(0.)
+    end                                                           
+
+    capacity_estimate = x.capacity_shortfall_mean[:]
+    capacity_var = quantile(capacity_estimate, alpha)
+    capacity_tail_losses = capacity_estimate[capacity_estimate .>= capacity_var]
+
+    capacity_cvar = if !isempty(capacity_tail_losses)
+        MeanEstimate(capacity_tail_losses)
+    else
+        MeanEstimate(0.)
+    end                                                               
+
+    return CVAR{N,L,T,E}(cvar, alpha, var, capacity_cvar, capacity_var)
+
+end
+
+function CVAR(x::ShortfallResult{N,L,T,E}, alpha::Float64, r::AbstractString) where {N,L,T,E}
+
+    i_r = findfirstunique(x.regions.names, r)
+    estimate = x.shortfall_region_samples[i_r, :]
+    var = quantile(estimate, alpha)
+    tail_losses = estimate[estimate .>= var]
+
+    cvar = if !isempty(tail_losses)
+        MeanEstimate(tail_losses)
+    else
+        MeanEstimate(0.)
+    end
+
+    capacity_estimate = x.capacity_shortfall_mean[i_r, :]
+    capacity_var = quantile(capacity_estimate, alpha)
+    capacity_tail_losses = capacity_estimate[capacity_estimate .>= capacity_var]
+
+    capacity_cvar = if !isempty(capacity_tail_losses)
+        MeanEstimate(capacity_tail_losses)
+    else
+        MeanEstimate(0.)
+    end 
+
+    return CVAR{N,L,T,E}(cvar, alpha, var, capacity_cvar, capacity_var)
+
+end
+
+function NCVAR(x::ShortfallResult{N,L,T,E}, cvar::CVAR) where {N,L,T,E}
+    demand = sum(x.regions.load)
+
+    if demand > 0
+        ncvar = div(cvar.cvar, demand/1e6)
+        var = div(cvar.var, demand/1e6)
+    else
+        ncvar = MeanEstimate(0.)
+        var = MeanEstimate(0.)
+    end
+
+    return NCVAR(ncvar, cvar.alpha, var)
+
+end
+
+function NCVAR(x::ShortfallResult{N,L,T,E}, cvar::CVAR, r::AbstractString) where {N,L,T,E}
+    i_r = findfirstunique(x.regions.names, r)
+    demand = sum(x.regions.load[i_r, :])
+
+    if demand > 0
+        ncvar = div(cvar.cvar, demand/1e6)
+        var = div(cvar.var, demand/1e6)
+    else
+        ncvar = MeanEstimate(0.)
+        var = MeanEstimate(0.)
+    end
+
+    return NCVAR(ncvar, cvar.alpha, var)
+
+end
+
 function finalize(
     acc::ShortfallAccumulator{S},
     system::SystemModel{N,L,T,P,E},
@@ -312,18 +415,22 @@ function finalize(
     nsamples = first(acc.unservedload_total.stats).n
 
     p2e = conversionfactor(L,T,P,E)
+    capacity_shortfall_mean = ue_regionperiod_mean
     ue_regionperiod_mean .*= p2e
     ue_total_std *= p2e
     ue_region_std .*= p2e
     ue_period_std .*= p2e
     ue_regionperiod_std .*= p2e
+    ue_sample = float(acc.unservedload_sample .* p2e)
+    ue_region_sample = float(acc.unservedload_region_sample .* p2e)
 
     return ShortfallResult{N,L,T,E,S}(
         nsamples, system.regions, system.timestamps,
         ep_total_mean, ep_total_std, ep_region_mean, ep_region_std,
         ep_period_mean, ep_period_std,
         ep_regionperiod_mean, ep_regionperiod_std,
-        ue_regionperiod_mean, ue_total_std,
-        ue_region_std, ue_period_std, ue_regionperiod_std)
+        ue_regionperiod_mean, capacity_shortfall_mean, ue_total_std,
+        ue_region_std, ue_period_std, ue_regionperiod_std,
+        ue_sample, ue_region_sample, )
 
 end