diff --git a/PRAS.jl/examples/pras_walkthrough.jl b/PRAS.jl/examples/pras_walkthrough.jl
index 081f982a..40b5fc53 100644
--- a/PRAS.jl/examples/pras_walkthrough.jl
+++ b/PRAS.jl/examples/pras_walkthrough.jl
@@ -176,3 +176,30 @@ println("Surplus in")
 # performed on the subset of samples in which the event was observed, using the
 # `ShortfallSamples`, `UtilizationSamples`, and
 # `StorageEnergySamples` result specifications instead.
+
+## [Assess tail-risk metrics](@id assess_tail_risks)
+
+# PRAS can calculate the risk-tail adequacy metrics, which are useful when analyzing extreme
+# occurences of shortfall events. The Conditional Value at Risk (CVAR) measures the weighted 
+# average outcomes beyond a pre-defined percentile $\alpha$. Therefore, the CVAR calculates
+# the expected shortfalls in the worst $(1-\alpha)$ outcomes. 
+# In general, CVAR can be used for calculating tail-risk outcomes across any of the resource
+# adequacy planning dimensions (energy shortfall, duration, and magnitude).
+# Currently, PRAS implements only unserved energy CVAR. 
+
+# Start by defining the confidence level $\alpha$, which defines the tail threshold beyond which we
+# want measure the expected shortfall.
+
+alpha = 0.95 # nth percentile for tail events
+
+# We can check the overall system unserved energy CVAR
+
+ue_cvar = CVAR(:energy, shortfall, alpha)
+println("System CVAR: $(ue_cvar)")
+
+# And we can calculate the noramzlied CVAR (NCVAR)
+ue_ncvar = NCVAR(shortfall, ue_cvar)
+println("System NCVAR: $(ue_ncvar)")
+
+# We can also find the unserved energy CVAR by region
+regional_ue_cvar = CVAR.(:energy, shortfall, alpha, sys.regions.names)
\ No newline at end of file
diff --git a/PRAS.jl/test/runtests.jl b/PRAS.jl/test/runtests.jl
index 7189f62e..9e9d0bb9 100644
--- a/PRAS.jl/test/runtests.jl
+++ b/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(:energy, 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(:energy, 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
\ No newline at end of file
diff --git a/PRASCore.jl/src/Results/Results.jl b/PRASCore.jl/src/Results/Results.jl
index 3e4e3300..b50731ff 100644
--- a/PRASCore.jl/src/Results/Results.jl
+++ b/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,
@@ -30,6 +30,25 @@ export
 include("utils.jl")
 include("metrics.jl")
 
+function _cvar(estimate::AbstractVector{<:Real}, alpha::Float64)
+    var = quantile(estimate, alpha)
+    tail = estimate[estimate .> var]   
+    cvar = isempty(tail) ? MeanEstimate(0.) : MeanEstimate(tail)
+    return cvar, var
+end
+
+function _ncvar(cvar::CVAR, demand::Real)
+    if demand > 0
+        scale = demand / 1e6
+        ncvar = div(cvar.cvar, scale)
+        var = cvar.var / scale
+    else
+        ncvar = MeanEstimate(0.)
+        var = 0.0
+    end
+    return ncvar, var
+end
+
 abstract type ResultSpec end
 
 abstract type ResultAccumulator{R<:ResultSpec} end
@@ -80,6 +99,16 @@ NEUE(x::AbstractShortfallResult, r::AbstractString, ::Colon) =
 NEUE(x::AbstractShortfallResult, ::Colon, ::Colon) =
     NEUE.(x, x.regions.names, permutedims(x.timestamps))
 
+CVAR(dim::Symbol, x::AbstractShortfallResult, alpha::Float64, args...) =
+    CVAR(Val(dim), x, alpha, args...)
+
+
+function CVAR(::Val, ::AbstractShortfallResult, ::Float64, args...)
+    throw(ArgumentError("Invalid dim, use one of $CVAR_METRICS"))
+end
+
+NCVAR(x::AbstractShortfallResult, cvar::CVAR, ::Colon) =
+    NCVAR.(x, cvar, x.regions.names)
 include("Shortfall.jl")
 include("ShortfallSamples.jl")
 
diff --git a/PRASCore.jl/src/Results/Shortfall.jl b/PRASCore.jl/src/Results/Shortfall.jl
index 0b00d7d2..71d8fe38 100644
--- a/PRASCore.jl/src/Results/Shortfall.jl
+++ b/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,6 +94,8 @@ 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,
@@ -97,7 +103,8 @@ function accumulator(
         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,7 +151,6 @@ 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
 
     shortfall_std::Float64
@@ -149,6 +158,9 @@ struct ShortfallResult{N, L, T <: Period, E <: EnergyUnit, S} <:
     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,
@@ -165,8 +177,10 @@ struct ShortfallResult{N, L, T <: Period, E <: EnergyUnit, S} <:
         shortfall_std::Float64,
         shortfall_region_std::Vector{Float64},
         shortfall_period_std::Vector{Float64},
-        shortfall_regionperiod_std::Matrix{Float64}
-    ) where {N,L,T<:Period,E<:EnergyUnit,S <: Union{Shortfall,DemandResponseShortfall}}
+        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 ||
             throw(DomainError("Sample count must be positive or `nothing`."))
@@ -185,7 +199,9 @@ 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,
@@ -195,7 +211,8 @@ struct ShortfallResult{N, L, T <: Period, E <: EnergyUnit, S} <:
             eventperiod_regionperiod_mean, eventperiod_regionperiod_std,
             shortfall_mean, shortfall_std,
             shortfall_region_std, shortfall_period_std,
-            shortfall_regionperiod_std)
+            shortfall_regionperiod_std, shortfall_samples,
+            shortfall_region_samples)
 
     end
 
@@ -292,6 +309,51 @@ function NEUE(x::ShortfallResult, r::AbstractString)
 
 end
 
+function CVAR(::Val{:energy}, x::ShortfallResult{N,L,T,E},alpha::Float64) where {N,L,T,E}
+    cvar, var = _cvar(x.shortfall_samples, alpha)
+    return CVAR{N,L,T,E}(:energy, cvar, alpha, var)
+end
+
+function CVAR(::Val{:energy}, x::ShortfallResult{N,L,T,E},alpha::Float64, r::AbstractString) where {N,L,T,E}
+    i_r = findfirstunique(x.regions.names, r)
+    cvar, var = _cvar(x.shortfall_region_samples[i_r, :], alpha)
+    return CVAR{N,L,T,E}(:energy, cvar, alpha, var)
+end
+
+function CVAR(::Val, ::ShortfallResult, ::Float64, ::StepRange{ZonedDateTime})
+    # To support this, add unservedload_period_sample::Matrix{Int} (N × nsamples)
+    # to ShortfallAccumulator, record it in recording.jl, and slice [i_t0:i_tf, :] here.
+    throw(ArgumentError(
+        "Time-sliced CVAR is not supported for ShortfallResult. " *
+        "Use ShortfallSamplesResult instead."))
+end
+
+function CVAR(::Val, ::ShortfallResult, ::Float64, ::AbstractString, ::ZonedDateTime)
+    # To support this, add unservedload_regionperiod_sample::Array{Int,3} (nregions × N × nsamples)
+    # to ShortfallAccumulator, record it in recording.jl, and slice [i_r, i_t, :] here.
+    throw(ArgumentError(
+        "Region+timestep CVAR is not supported for ShortfallResult. " *
+        "Use ShortfallSamplesResult instead."))
+end
+
+function NCVAR(x::ShortfallResult, cvar::CVAR)
+    demand = sum(x.regions.load)
+
+    ncvar, var = _ncvar(cvar, demand)
+
+    return NCVAR(cvar.dim, ncvar, cvar.alpha, var)
+end
+
+function NCVAR(x::ShortfallResult, cvar::CVAR, r::AbstractString)
+    i_r = findfirstunique(x.regions.names, r)
+    demand = sum(x.regions.load[i_r, :])
+
+    ncvar, var = _ncvar(cvar, demand)
+
+    return NCVAR(cvar.dim, ncvar, cvar.alpha, var)
+  
+end
+
 function finalize(
     acc::ShortfallAccumulator{S},
     system::SystemModel{N,L,T,P,E},
@@ -317,6 +379,8 @@ function finalize(
     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,
@@ -324,6 +388,7 @@ function finalize(
         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_region_std, ue_period_std, ue_regionperiod_std,
+        ue_sample, ue_region_sample)
 
 end
diff --git a/PRASCore.jl/src/Results/ShortfallSamples.jl b/PRASCore.jl/src/Results/ShortfallSamples.jl
index 99106433..c8d982ef 100644
--- a/PRASCore.jl/src/Results/ShortfallSamples.jl
+++ b/PRASCore.jl/src/Results/ShortfallSamples.jl
@@ -117,6 +117,25 @@ function getindex(
     return vec(p2e * x.shortfall[i_r, i_t, :])
 end
 
+function getindex(
+    x::ShortfallSamplesResult{N,L,T,P,E}, t::StepRange{ZonedDateTime}
+) where {N,L,T,P,E}
+    i_t0 = findfirstunique(x.timestamps, first(t))
+    i_tf = findlastunique(x.timestamps, last(t))
+    p2e = conversionfactor(L, T, P, E)
+    return vec(p2e * sum(x.shortfall[:, i_t0:i_tf, :], dims=1:2))
+end
+
+function getindex(
+    x::ShortfallSamplesResult{N,L,T,P,E}, r::AbstractString, t::StepRange{ZonedDateTime}
+) where {N,L,T,P,E}
+    i_r = findfirstunique(x.regions.names, r)
+    i_t0 = findfirstunique(x.timestamps, first(t))
+    i_tf = findlastunique(x.timestamps, last(t))
+    p2e = conversionfactor(L, T, P, E)
+    return vec(p2e * sum(x.shortfall[i_r, i_t0:i_tf, :], dims=1))
+end
+
 
 function LOLE(x::ShortfallSamplesResult{N,L,T}) where {N,L,T}
     eventperiods = sum(sum(x.shortfall, dims=1) .> 0, dims=2)
@@ -185,6 +204,52 @@ function NEUE(x::ShortfallSamplesResult, r::AbstractString)
 
 end
 
+function CVAR(::Val{:energy}, x::ShortfallSamplesResult{N,L,T,P,E}, alpha::Float64) where {N,L,T,P,E}
+    cvar, var = _cvar(x[], alpha)
+    return CVAR{N,L,T,E}(:energy, cvar, alpha, var)
+end
+
+function CVAR(::Val{:energy}, x::ShortfallSamplesResult{N,L,T,P,E}, alpha::Float64, r::AbstractString) where {N,L,T,P,E}
+    cvar, var = _cvar(x[r], alpha)
+    return CVAR{N,L,T,E}(:energy, cvar, alpha, var)
+end
+
+function CVAR(::Val{:energy}, x::ShortfallSamplesResult{N,L,T,P,E}, alpha::Float64, t::StepRange{ZonedDateTime}) where {N,L,T,P,E}
+    cvar, var = _cvar(x[t], alpha)
+    return CVAR{N,L,T,E}(:energy, cvar, alpha, var)
+end
+
+function CVAR(::Val{:energy}, x::ShortfallSamplesResult{N,L,T,P,E}, alpha::Float64, r::AbstractString, t::ZonedDateTime) where {N,L,T,P,E}
+    cvar, var = _cvar(x[r, t], alpha)
+    return CVAR{N,L,T,E}(:energy, cvar, alpha, var)
+end
+
+function CVAR(::Val{:energy}, x::ShortfallSamplesResult{N,L,T,P,E}, alpha::Float64, r::AbstractString, t::StepRange{ZonedDateTime}) where {N,L,T,P,E}
+    cvar, var = _cvar(x[r, t], alpha)
+    return CVAR{N,L,T,E}(:energy, cvar, alpha, var)
+end
+
+CVAR(dim::Symbol, x::ShortfallSamplesResult, alpha::Float64, ::Colon, t::StepRange{ZonedDateTime}) =
+    CVAR.(dim, x, alpha, x.regions.names, Ref(t))
+
+function NCVAR(x::ShortfallSamplesResult, cvar::CVAR)
+    demand = sum(x.regions.load)
+
+    ncvar, var = _ncvar(cvar, demand)
+
+    return NCVAR(cvar.dim, ncvar, cvar.alpha, var)
+
+end
+
+function NCVAR(x::ShortfallSamplesResult, cvar::CVAR, r::AbstractString)
+    i_r = findfirstunique(x.regions.names, r)
+    demand = sum(x.regions.load[i_r, :])
+
+    ncvar, var = _ncvar(cvar, demand)
+    return NCVAR(cvar.dim, ncvar, cvar.alpha, var)
+
+end
+
 function finalize(
     acc::ShortfallSamplesAccumulator{S},
     system::SystemModel{N,L,T,P,E},
diff --git a/PRASCore.jl/src/Results/metrics.jl b/PRASCore.jl/src/Results/metrics.jl
index f44c7096..bc7f6d8d 100644
--- a/PRASCore.jl/src/Results/metrics.jl
+++ b/PRASCore.jl/src/Results/metrics.jl
@@ -22,7 +22,11 @@ MeanEstimate(mu::Real, sigma::Real, n::Int) = MeanEstimate(mu, sigma / sqrt(n))
 
 function MeanEstimate(xs::AbstractArray{<:Real})
     est = mean(xs)
-    return MeanEstimate(est, std(xs, mean=est), length(xs))
+    if length(xs) > 1
+        MeanEstimate(est, std(xs, mean=est), length(xs))
+    else
+        MeanEstimate(est)
+    end
 end
 
 val(est::MeanEstimate) = est.estimate
@@ -164,3 +168,82 @@ function Base.show(io::IO, x::NEUE)
     print(io, "NEUE = ", x.neue, " ppm")
 
 end
+
+const CVAR_METRICS = (:energy,)
+
+function _check_dim(dim::Symbol, valid::Tuple{Vararg{Symbol}}=CVAR_METRICS)
+    dim in valid ||
+        throw(ArgumentError("$dim is not a valid dim, use one of $valid"))
+end
+
+"""
+    CVAR
+
+`CVAR` reports conditional value at risk of shortfalls, for total unserved energy shortfalls.
+
+Contains both the estimated value itself as well as the standard error
+of that estimate, which can be extracted with `val` and `stderror`,
+respectively.
+"""
+struct CVAR{N,L,T<:Period,E<:EnergyUnit} <: ReliabilityMetric
+
+    dim::Symbol
+    cvar::MeanEstimate
+    alpha::Float64
+    var::Float64
+
+    function CVAR{N,L,T,E}(dim::Symbol,
+                           cvar::MeanEstimate,
+                           alpha::Float64,
+                           var::Float64) where {N,L,T<:Period,E<:EnergyUnit}
+        _check_dim(dim)
+        val(cvar) >= 0 || throw(DomainError(val(cvar),
+            "$(val(cvar)) is not a valid CVAR"))
+        0 <= alpha < 1 || throw(DomainError(alpha,
+            "$alpha is not a valid confidence level"))
+        new{N,L,T,E}(dim, cvar, alpha, var)
+    end
+
+end
+
+val(x::CVAR) = val(x.cvar)
+stderror(x::CVAR) = stderror(x.cvar)
+
+function Base.show(io::IO, x::CVAR{N,L,T,E}) where {N,L,T,E}
+    print(io, "CVAR@$(x.alpha) = ", x.cvar, " ",
+          unitsymbol(E), "/", N*L == 1 ? "" : N*L, unitsymbol(T))
+end
+
+"""
+    NCVAR
+
+`NCVAR` reports normalized conditional value at risk of shortfalls, for total unserved energy shortfalls.
+
+Contains both the estimated value itself as well as the standard error
+of that estimate, which can be extracted with `val` and `stderror`,
+respectively.
+"""
+struct NCVAR <: ReliabilityMetric
+
+    dim::Symbol
+    ncvar::MeanEstimate
+    alpha::Float64
+    var::Float64
+    
+    function NCVAR(dim::Symbol, ncvar::MeanEstimate, alpha::Float64, var::Float64)
+
+        _check_dim(dim)
+        val(ncvar) >= 0 || throw(DomainError(val(ncvar),
+            "$(val(ncvar)) is not a valid NCVAR"))
+        new(dim, ncvar, alpha, var)
+    end
+
+end
+
+val(x::NCVAR) = val(x.ncvar)
+stderror(x::NCVAR) = stderror(x.ncvar)
+
+function Base.show(io::IO, x::NCVAR)
+    print(io, "NCVAR@$(x.alpha) = ", x.ncvar, " ppm")
+
+end
\ No newline at end of file
diff --git a/PRASCore.jl/src/Results/utils.jl b/PRASCore.jl/src/Results/utils.jl
index 41782995..ba202c35 100644
--- a/PRASCore.jl/src/Results/utils.jl
+++ b/PRASCore.jl/src/Results/utils.jl
@@ -40,3 +40,10 @@ function findfirstunique(a::AbstractVector{T}, i::T) where T
     i_idx === nothing && throw(BoundsError(a))
     return i_idx
 end
+
+function findlastunique(a::AbstractVector{T}, i::T) where T
+    i_idx = findlast(isequal(i), a)
+    i_idx === nothing && throw(BoundsError(a))
+    return i_idx
+end
+
diff --git a/PRASCore.jl/src/Simulations/recording.jl b/PRASCore.jl/src/Simulations/recording.jl
index a156eb32..44f6ff7d 100644
--- a/PRASCore.jl/src/Simulations/recording.jl
+++ b/PRASCore.jl/src/Simulations/recording.jl
@@ -55,7 +55,9 @@ function record!(
 end
 
 function reset!(acc::Results.ShortfallAccumulator, sampleid::Int)
-
+    # Store total UE for each sample
+    acc.unservedload_sample[sampleid] = acc.unservedload_total_currentsim 
+    
     # Store regional / total sums for current simulation
     fit!(acc.periodsdropped_total, acc.periodsdropped_total_currentsim)
     fit!(acc.unservedload_total, acc.unservedload_total_currentsim)
@@ -63,6 +65,7 @@ function reset!(acc::Results.ShortfallAccumulator, sampleid::Int)
     for r in eachindex(acc.periodsdropped_region)
         fit!(acc.periodsdropped_region[r], acc.periodsdropped_region_currentsim[r])
         fit!(acc.unservedload_region[r], acc.unservedload_region_currentsim[r])
+        acc.unservedload_region_sample[r, sampleid] = acc.unservedload_region_currentsim[r]
     end
 
     # Reset for new simulation
diff --git a/PRASCore.jl/src/Systems/TestData.jl b/PRASCore.jl/src/Systems/TestData.jl
index 6b5ba264..f77dba5a 100644
--- a/PRASCore.jl/src/Systems/TestData.jl
+++ b/PRASCore.jl/src/Systems/TestData.jl
@@ -36,6 +36,7 @@ singlenode_a_lole = 0.355
 singlenode_a_lolps = [0.028, 0.271, 0.028, 0.028]
 singlenode_a_eue = 1.59
 singlenode_a_eues = [0.29, 0.832, 0.29, 0.178]
+singlenode_a_cvar = 13.34
 
 ## Single-Region System A - 5 minute version
 
@@ -99,6 +100,7 @@ singlenode_b_lole = 0.96
 singlenode_b_lolps = [0.19, 0.19, 0.19, 0.1, 0.1, 0.19]
 singlenode_b_eue = 7.11
 singlenode_b_eues = [1.29, 1.29, 1.29, 0.85, 1.05, 1.34]
+singlenode_b_cvar = 31.71
 
 
 # Single-Region System B, with storage
diff --git a/PRASCore.jl/test/Results/metrics.jl b/PRASCore.jl/test/Results/metrics.jl
index fb7bf796..0172c1f7 100644
--- a/PRASCore.jl/test/Results/metrics.jl
+++ b/PRASCore.jl/test/Results/metrics.jl
@@ -71,5 +71,29 @@
         @test_throws DomainError NEUE(MeanEstimate(-1.2))
 
     end
+  
+    @testset "CVAR" begin
+
+        cvar1 = CVAR{2,1,Hour,MWh}(:energy,MeanEstimate(1.2), 0.95, 2.0)
+        @test string(cvar1) == "CVAR@0.95 = 1.20000 MWh/2h"
+
+        cvar2 = CVAR{1,2,Year,GWh}(:energy,MeanEstimate(17.2, 1.3), 0.95, 1.2)
+        @test string(cvar2) == "CVAR@0.95 = 17±1 GWh/2y"
+
+        @test_throws DomainError CVAR{1,1,Hour,MWh}(:energy,MeanEstimate(-1.2), 0.95, 1.2)
+
+    end    
+
+    @testset "NCVAR" begin
+
+        ncvar1 = NCVAR(:energy, MeanEstimate(1.2), 0.95, 1.2)
+        @test string(ncvar1) == "NCVAR@0.95 = 1.20000 ppm"
+
+        ncvar2 = NCVAR(:energy, MeanEstimate(17.2, 1.3), 0.95, 1.3)
+        @test string(ncvar2) == "NCVAR@0.95 = 17±1 ppm"
+
+        @test_throws DomainError NCVAR(:energy, MeanEstimate(-1.2), 0.95, -1.2)
+
+    end
 
 end
diff --git a/PRASCore.jl/test/Results/shortfall.jl b/PRASCore.jl/test/Results/shortfall.jl
index 9f5d4a01..503d6b72 100644
--- a/PRASCore.jl/test/Results/shortfall.jl
+++ b/PRASCore.jl/test/Results/shortfall.jl
@@ -4,13 +4,15 @@
     N = DD.nperiods
     r, r_idx, r_bad = DD.testresource, DD.testresource_idx, DD.notaresource
     t, t_idx, t_bad = DD.testperiod, DD.testperiod_idx, DD.notaperiod
+    alpha = DD.alpha
 
     result = PRASCore.Results.ShortfallResult{N,1,Hour,MWh,Shortfall}(
-        DD.nsamples, Regions{N,MW}(DD.resourcenames, DD.resource_vals), DD.periods,
-        DD.d1, DD.d2, DD.d1_resource, DD.d2_resource,
-        DD.d1_period, DD.d2_period, DD.d1_resourceperiod, DD.d2_resourceperiod,
-        DD.d3_resourceperiod,
-        DD.d4, DD.d4_resource, DD.d4_period, DD.d4_resourceperiod)
+        DD.nsamples, Regions{N,MW}(DD.resourcenames, DD.resource_vals),
+        DD.periods, DD.d1, DD.d2, DD.d1_resource, DD.d2_resource,
+        DD.d1_period, DD.d2_period, DD.d1_resourceperiod,
+        DD.d2_resourceperiod, DD.d3_resourceperiod, DD.d4,
+        DD.d4_resource, DD.d4_period, DD.d4_resourceperiod,
+        DD.d1_sample, DD.d1_resourcesample)
 
     # Overall
 
@@ -28,6 +30,17 @@
     load = sum(DD.resource_vals)
     @test val(neue) ≈ first(result[]) / load*1e6
     @test stderror(neue) ≈ last(result[]) / sqrt(DD.nsamples) / load*1e6
+
+    cvar = CVAR(:energy, result, alpha)
+    estimate = result.shortfall_samples;
+    tail_losses = estimate[estimate .> quantile(estimate, alpha)];
+    @test val(cvar) ≈ mean(tail_losses)
+    @test stderror(cvar) ≈ std(tail_losses) / sqrt(length(tail_losses))
+
+    ncvar = NCVAR(result, cvar)
+    @test val(ncvar) ≈ val(cvar) / load*1e6
+    @test stderror(ncvar) ≈ stderror(cvar) / load*1e6
+
     # Region-specific
 
     @test result[r] ≈ (sum(DD.d3_resourceperiod[r_idx,:]), DD.d4_resource[r_idx])
@@ -45,10 +58,23 @@
     @test val(region_neue) ≈ first(result[r]) / load*1e6
     @test stderror(region_neue) ≈ last(result[r]) / sqrt(DD.nsamples) / load*1e6
 
+    region_cvar = CVAR(:energy, result, alpha, r)
+    region_estimate = result.shortfall_region_samples[r_idx, :];
+    region_tail_losses = region_estimate[region_estimate .>= quantile(region_estimate, alpha)];
+    @test val(region_cvar) ≈ mean(region_tail_losses)
+    @test stderror(region_cvar) ≈ std(region_tail_losses) / sqrt(length(region_tail_losses))
+
+    region_ncvar = NCVAR(result, region_cvar, r)
+    @test val(region_ncvar) ≈ val(region_cvar) / load*1e6
+    @test stderror(region_ncvar) ≈ stderror(region_cvar) / load*1e6
+
     @test_throws BoundsError result[r_bad]
     @test_throws BoundsError LOLE(result, r_bad)
     @test_throws BoundsError EUE(result, r_bad)
     @test_throws BoundsError NEUE(result, r_bad)
+    @test_throws BoundsError CVAR(:energy,result, alpha, r_bad)
+    @test_throws BoundsError NCVAR(result, region_cvar, r_bad)
+    @test_throws ArgumentError CVAR(:power, result, alpha)
 
     # Period-specific
 
@@ -62,6 +88,8 @@
     @test val(period_eue) ≈ first(result[t])
     @test stderror(period_eue) ≈ last(result[t]) / sqrt(DD.nsamples)
 
+    @test_throws ArgumentError CVAR(:energy, result, alpha, DD.periods)
+
     @test_throws BoundsError result[t_bad]
     @test_throws BoundsError LOLE(result, t_bad)
     @test_throws BoundsError EUE(result, t_bad)
@@ -99,7 +127,8 @@ end
 
     N = DD.nperiods
     r, r_idx, r_bad = DD.testresource, DD.testresource_idx, DD.notaresource
-    t, t_idx, t_bad = DD.testperiod, DD.testperiod_idx, DD.notaperiod
+    t, t_idx, t_bad, badperiods = DD.testperiod, DD.testperiod_idx, DD.notaperiod, DD.badperiods
+    alpha = 0.95
 
     result = PRASCore.Results.ShortfallSamplesResult{N,1,Hour,MW,MWh,ShortfallSamples}(
         Regions{N,MW}(DD.resourcenames, DD.resource_vals), DD.periods, DD.d)
@@ -123,6 +152,16 @@ end
     @test val(neue) ≈ mean(result[]) / load*1e6
     @test stderror(neue) ≈ std(result[]) / sqrt(DD.nsamples) / load*1e6
 
+    ue_cvar = CVAR(:energy, result, alpha)
+    estimate = result[];
+    tail_losses = estimate[estimate .>= quantile(estimate, alpha)];
+    @test val(ue_cvar) ≈ mean(tail_losses)
+    @test stderror(ue_cvar) ≈ std(tail_losses) / sqrt(length(tail_losses))
+
+    ncvar = NCVAR(result, ue_cvar)
+    @test val(ncvar) ≈ val(ue_cvar) / load*1e6
+    @test stderror(ncvar) ≈ stderror(ue_cvar) / load*1e6
+
     # Region-specific
 
     @test length(result[r]) == DD.nsamples
@@ -142,10 +181,22 @@ end
     @test val(region_neue) ≈ mean(result[r]) / load*1e6
     @test stderror(region_neue) ≈ std(result[r]) / sqrt(DD.nsamples) / load*1e6
 
+    region_ue_cvar = CVAR(:energy, result, alpha, r)
+    region_estimate = result[r];
+    region_tail_losses = region_estimate[region_estimate .>= quantile(region_estimate, alpha)];
+    @test val(region_ue_cvar) ≈ mean(region_tail_losses)
+    @test stderror(region_ue_cvar) ≈ std(region_tail_losses) / sqrt(length(region_tail_losses))
+
+    region_ncvar = NCVAR(result, region_ue_cvar, r)
+    @test val(region_ncvar) ≈ val(region_ue_cvar) / load*1e6
+    @test stderror(region_ncvar) ≈ stderror(region_ue_cvar) / load*1e6
+
     @test_throws BoundsError result[r_bad]
     @test_throws BoundsError LOLE(result, r_bad)
     @test_throws BoundsError EUE(result, r_bad)
     @test_throws BoundsError NEUE(result, r_bad)
+    @test_throws BoundsError CVAR(:energy, result, alpha, r_bad)
+    @test_throws BoundsError NCVAR(result, region_ue_cvar, r_bad)
 
     # Period-specific
 
@@ -161,9 +212,16 @@ end
     @test val(period_eue) ≈ mean(result[t])
     @test stderror(period_eue) ≈ std(result[t]) / sqrt(DD.nsamples)
 
+    period_cvar = CVAR(:energy, result, alpha, DD.periods)
+    period_estimate = result[DD.periods];
+    period_tail_losses = period_estimate[period_estimate .> quantile(period_estimate, alpha)];
+    @test val(period_cvar) ≈ mean(period_tail_losses)
+    @test stderror(period_cvar) ≈ std(period_tail_losses) / sqrt(length(period_tail_losses))
+
     @test_throws BoundsError result[t_bad]
     @test_throws BoundsError LOLE(result, t_bad)
     @test_throws BoundsError EUE(result, t_bad)
+    @test_throws BoundsError CVAR(:energy, result, alpha, badperiods)
 
     # Region + period-specific
 
@@ -180,6 +238,12 @@ end
     @test val(regionperiod_eue) ≈ mean(result[r, t])
     @test stderror(regionperiod_eue) ≈ std(result[r, t]) / sqrt(DD.nsamples)
 
+    regionperiod_cvar = CVAR(:energy, result, alpha, r, t)
+    regionperiod_estimate = result[r, t];
+    regionperiod_tail_losses = regionperiod_estimate[regionperiod_estimate .>= quantile(regionperiod_estimate, alpha)];
+    @test val(regionperiod_cvar) ≈ mean(regionperiod_tail_losses)
+    @test stderror(regionperiod_cvar) ≈ std(regionperiod_tail_losses) / sqrt(length(regionperiod_tail_losses))
+
     @test_throws BoundsError result[r, t_bad]
     @test_throws BoundsError result[r_bad, t]
     @test_throws BoundsError result[r_bad, t_bad]
@@ -192,4 +256,8 @@ end
     @test_throws BoundsError EUE(result, r_bad, t)
     @test_throws BoundsError EUE(result, r_bad, t_bad)
 
+    @test_throws BoundsError CVAR(:energy, result, alpha, r, t_bad)
+    @test_throws BoundsError CVAR(:energy, result, alpha, r_bad, t)
+    @test_throws BoundsError CVAR(:energy, result, alpha, r_bad, t_bad)
+
 end
diff --git a/PRASCore.jl/test/Simulations/runtests.jl b/PRASCore.jl/test/Simulations/runtests.jl
index 3c518634..a904e522 100644
--- a/PRASCore.jl/test/Simulations/runtests.jl
+++ b/PRASCore.jl/test/Simulations/runtests.jl
@@ -5,6 +5,7 @@
     end
 
     nstderr_tol = 3
+    alpha = 0.95
 
     simspec = SequentialMonteCarlo(samples=100_000, seed=1, threaded=false)
     smallsample = SequentialMonteCarlo(samples=10, seed=123, threaded=false)
@@ -66,6 +67,8 @@
                           TestData.singlenode_a_lole, nstderr_tol)
         @test withinrange(EUE(shortfall_1a),
                           TestData.singlenode_a_eue, nstderr_tol)
+        @test withinrange(CVAR(:energy, shortfall_1a, alpha),
+                          TestData.singlenode_a_cvar, nstderr_tol)                          
         @test withinrange(LOLE(shortfall_1a, "Region"),
                           TestData.singlenode_a_lole, nstderr_tol)
         @test withinrange(EUE(shortfall_1a, "Region"),
@@ -134,10 +137,14 @@
                           TestData.singlenode_b_lole, nstderr_tol)
         @test withinrange(EUE(shortfall_1b),
                           TestData.singlenode_b_eue, nstderr_tol)
+        @test withinrange(CVAR(:energy, shortfall_1b, alpha),
+                          TestData.singlenode_b_cvar, nstderr_tol)                          
         @test withinrange(LOLE(shortfall_1b, "Region"),
                           TestData.singlenode_b_lole, nstderr_tol)
         @test withinrange(EUE(shortfall_1b, "Region"),
                           TestData.singlenode_b_eue, nstderr_tol)
+        @test withinrange(CVAR(:energy, shortfall_1b, alpha, "Region"),
+                          TestData.singlenode_b_cvar, nstderr_tol)                                                    
 
         @test all(LOLE.(shortfall_1b, timestamps_b) .≈
                   LOLE.(shortfall2_1b, timestamps_b))
diff --git a/PRASCore.jl/test/dummydata.jl b/PRASCore.jl/test/dummydata.jl
index 06f7d139..f28606c4 100644
--- a/PRASCore.jl/test/dummydata.jl
+++ b/PRASCore.jl/test/dummydata.jl
@@ -6,6 +6,7 @@ using TimeZones
 const tz = tz"UTC"
 
 nsamples = 100
+alpha = 0.95
 
 resourcenames = ["A", "B", "C"]
 nresources = length(resourcenames)
@@ -20,6 +21,7 @@ testinterface = interfacenames[testinterface_idx]
 notaninterface = "X"=>"Y"
 
 periods = ZonedDateTime(2012,4,1,0,tz):Hour(1):ZonedDateTime(2012,4,7,23,tz)
+badperiods = ZonedDateTime(2013,4,1,0,tz):Hour(1):ZonedDateTime(2013,4,7,23,tz)
 nperiods = length(periods)
 resource_vals = rand(0:999, nresources, nperiods)
 testperiod_idx = 29
@@ -32,6 +34,8 @@ d1 = rand()
 d1_resource = rand(nresources)
 d1_period = rand(nperiods)
 d1_resourceperiod = rand(nresources, nperiods)
+d1_sample = rand(nsamples) * 999
+d1_resourcesample = rand(nresources, nsamples) * 999
 
 d2 = rand()
 d2_resource = rand(nresources)
diff --git a/docs/src/PRAS/results.md b/docs/src/PRAS/results.md
index 893af5d7..2dd88f92 100644
--- a/docs/src/PRAS/results.md
+++ b/docs/src/PRAS/results.md
@@ -100,7 +100,7 @@ the total shortfall across all regions and time periods.
 
 Shortfall results are unique among the built-in result types in that the raw
 results can also be converted to specific probabilistic risk metrics
-(**EUE** and **LOLE**). For sampling-based methods, both metric
+(**EUE**,  **LOLE**, or **CVAR**). For sampling-based methods, both metric
 estimates and the standard error of those estimates are provided. For example,
 after assessing the system, metrics across all regions and the full simulation
 horizon can be extracted as:
@@ -109,6 +109,7 @@ horizon can be extracted as:
 shortfall, = assess(sys, SequentialMonteCarlo(), Shortfall())
 eue_overall = EUE(shortfall)
 lole_overall = LOLE(shortfall)
+ue_cvar_overall = CVAR(:energy, shortfall, 0.95)
 ```
 
 More specific metrics can be obtained as well:
diff --git a/docs/src/PRASCore/api.md b/docs/src/PRASCore/api.md
index a555af07..41865b63 100644
--- a/docs/src/PRASCore/api.md
+++ b/docs/src/PRASCore/api.md
@@ -42,4 +42,6 @@ PRASCore.Results.DemandResponseAvailability
 PRASCore.Results.DemandResponseEnergy
 PRASCore.Results.DemandResponseEnergySamples
 PRASCore.Results.LineAvailability
+PRASCore.Results.CVAR
+PRASCore.Results.NCVAR
 ```