From bd3c07af947df87802e9c3601484e3c0e3a7019b Mon Sep 17 00:00:00 2001
From: Abdelrahman Ayad <abdelrahman.ayad@gmail.com>
Date: Fri, 30 Jan 2026 16:58:22 -0700
Subject: [PATCH 01/23] initial CVAR functions commit

---
 PRASCore.jl/src/Results/Results.jl          |  4 +-
 PRASCore.jl/src/Results/ShortfallSamples.jl | 55 +++++++++++++++++++++
 PRASCore.jl/src/Results/metrics.jl          | 32 ++++++++++++
 3 files changed, 89 insertions(+), 2 deletions(-)

diff --git a/PRASCore.jl/src/Results/Results.jl b/PRASCore.jl/src/Results/Results.jl
index 3e4e3300..c3356b56 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,
 
     # Result specifications
     Shortfall, ShortfallSamples,
diff --git a/PRASCore.jl/src/Results/ShortfallSamples.jl b/PRASCore.jl/src/Results/ShortfallSamples.jl
index 99106433..61d7b1cd 100644
--- a/PRASCore.jl/src/Results/ShortfallSamples.jl
+++ b/PRASCore.jl/src/Results/ShortfallSamples.jl
@@ -185,6 +185,61 @@ function NEUE(x::ShortfallSamplesResult, r::AbstractString)
 
 end
 
+function CVAR(x::ShortfallSamplesResult{N,L,T,P,E}, alpha::Float64) where {N,L,T,P,E}
+    estimate = x[]
+    tail_losses = estimate[estimate .>= quantile(estimate, alpha)]
+
+    cvar = if !isempty(tail_losses)
+        MeanEstimate(tail_losses)
+    else
+        MeanEstimate(0.)
+    end
+    
+    return CVAR{N,L,T,E}(cvar, alpha)
+
+end
+
+function CVAR(x::ShortfallSamplesResult{N,L,T,P,E}, alpha::Float64, r::AbstractString) where {N,L,T,P,E}
+    estimate = x[r]
+    tail_losses = estimate[estimate .>= quantile(estimate, alpha)]
+
+    cvar = if !isempty(tail_losses)
+        MeanEstimate(tail_losses)
+    else
+        MeanEstimate(0.)
+    end
+    
+    return CVAR{N,L,T,E}(cvar, alpha)
+
+end
+
+function CVAR(x::ShortfallSamplesResult{N,L,T,P,E}, alpha::Float64, t::ZonedDateTime) where {N,L,T,P,E}
+    estimate = x[t]
+    tail_losses = estimate[estimate .>= quantile(estimate, alpha)]
+
+    cvar = if !isempty(tail_losses)
+        MeanEstimate(tail_losses)
+    else
+        MeanEstimate(0.)
+    end
+    
+    return CVAR{N,L,T,E}(cvar, alpha)
+
+end
+
+function CVAR(x::ShortfallSamplesResult{N,L,T,P,E}, alpha::Float64, r::AbstractString, t::ZonedDateTime) where {N,L,T,P,E}
+    estimate = x[r, t]
+    tail_losses = estimate[estimate .>= quantile(estimate, alpha)]
+
+    cvar = if !isempty(tail_losses)
+        MeanEstimate(tail_losses)
+    else
+        MeanEstimate(0.)
+    end
+    
+    return CVAR{N,L,T,E}(cvar, alpha)
+
+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..0eea834c 100644
--- a/PRASCore.jl/src/Results/metrics.jl
+++ b/PRASCore.jl/src/Results/metrics.jl
@@ -164,3 +164,35 @@ function Base.show(io::IO, x::NEUE)
     print(io, "NEUE = ", x.neue, " ppm")
 
 end
+
+"""
+    CVAR
+
+`CVAR` reports conditional value at risk of 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
+
+    cvar::MeanEstimate
+    alpha::Float64
+    
+    function CVAR{N,L,T,E}(cvar::MeanEstimate, alpha::Float64) where {N,L,T<:Period,E<:EnergyUnit}
+        val(cvar) >= 0 || throw(DomainError(
+            "$val is not a valid CVAR"))
+        new{N,L,T,E}(cvar, alpha)
+    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
\ No newline at end of file

From b1f105ac1ba63fc4f003ad1ce4a805a3a1197ca6 Mon Sep 17 00:00:00 2001
From: Abdelrahman Ayad <abdelrahman.ayad@gmail.com>
Date: Mon, 2 Feb 2026 00:32:28 -0700
Subject: [PATCH 02/23] added normalized CVAR (NCVAR)

---
 PRASCore.jl/src/Results/Results.jl          |  2 +-
 PRASCore.jl/src/Results/ShortfallSamples.jl | 46 +++++++++++++++++++++
 PRASCore.jl/src/Results/metrics.jl          | 30 ++++++++++++++
 3 files changed, 77 insertions(+), 1 deletion(-)

diff --git a/PRASCore.jl/src/Results/Results.jl b/PRASCore.jl/src/Results/Results.jl
index c3356b56..93f99f27 100644
--- a/PRASCore.jl/src/Results/Results.jl
+++ b/PRASCore.jl/src/Results/Results.jl
@@ -13,7 +13,7 @@ export
 
     # Metrics
     ReliabilityMetric, LOLE, EUE, NEUE,
-    val, stderror, CVAR,
+    val, stderror, CVAR, NCVAR,
 
     # Result specifications
     Shortfall, ShortfallSamples,
diff --git a/PRASCore.jl/src/Results/ShortfallSamples.jl b/PRASCore.jl/src/Results/ShortfallSamples.jl
index 61d7b1cd..fdefc495 100644
--- a/PRASCore.jl/src/Results/ShortfallSamples.jl
+++ b/PRASCore.jl/src/Results/ShortfallSamples.jl
@@ -240,6 +240,52 @@ function CVAR(x::ShortfallSamplesResult{N,L,T,P,E}, alpha::Float64, r::AbstractS
     return CVAR{N,L,T,E}(cvar, alpha)
 
 end
+
+function NCVAR(x::ShortfallSamplesResult{N,L,T,P}, alpha::Float64) where {N,L,T,P}
+    demand = sum(x.regions.load)
+
+    estimate = x[]
+    tail_losses = estimate[estimate .>= quantile(estimate, alpha)]
+
+    cvar = if !isempty(tail_losses)
+        MeanEstimate(tail_losses)
+    else
+        MeanEstimate(0.)
+    end
+
+    ncvar = if demand > 0
+        div(cvar, demand/1e6)
+    else
+        MeanEstimate(0.)
+    end
+
+    return NCVAR{N,L,T}(ncvar, alpha)
+
+end
+
+function NCVAR(x::ShortfallSamplesResult{N,L,T,P}, alpha::Float64, r::AbstractString) where {N,L,T,P}
+    i_r = findfirstunique(x.regions.names, r)
+    demand = sum(x.regions.load[i_r, :])
+
+    estimate = x[r]
+    tail_losses = estimate[estimate .>= quantile(estimate, alpha)]
+
+    cvar = if !isempty(tail_losses)
+        MeanEstimate(tail_losses)
+    else
+        MeanEstimate(0.)
+    end
+
+    ncvar = if demand > 0
+        div(cvar, demand/1e6)
+    else
+        MeanEstimate(0.)
+    end
+    
+    return NCVAR{N,L,T}(ncvar, alpha)
+
+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 0eea834c..b856e733 100644
--- a/PRASCore.jl/src/Results/metrics.jl
+++ b/PRASCore.jl/src/Results/metrics.jl
@@ -195,4 +195,34 @@ 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.
+
+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{N,L,T<:Period} <: ReliabilityMetric
+
+    ncvar::MeanEstimate
+    alpha::Float64
+    
+    function NCVAR{N,L,T}(ncvar::MeanEstimate, alpha::Float64) where {N,L,T<:Period}
+        val(ncvar) >= 0 || throw(DomainError(
+            "$val is not a valid NCVAR"))
+        new{N,L,T}(ncvar, alpha)
+    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

From e633074925b5699eba95b8f6175cc3e7e5a34043 Mon Sep 17 00:00:00 2001
From: Abdelrahman Ayad <abdelrahman.ayad@gmail.com>
Date: Mon, 2 Feb 2026 11:01:15 -0700
Subject: [PATCH 03/23] formulate normalized CVAR based on CVAR struct

---
 PRASCore.jl/src/Results/ShortfallSamples.jl | 30 +++++----------------
 PRASCore.jl/src/Results/metrics.jl          |  6 ++---
 2 files changed, 9 insertions(+), 27 deletions(-)

diff --git a/PRASCore.jl/src/Results/ShortfallSamples.jl b/PRASCore.jl/src/Results/ShortfallSamples.jl
index fdefc495..3e78e087 100644
--- a/PRASCore.jl/src/Results/ShortfallSamples.jl
+++ b/PRASCore.jl/src/Results/ShortfallSamples.jl
@@ -241,48 +241,30 @@ function CVAR(x::ShortfallSamplesResult{N,L,T,P,E}, alpha::Float64, r::AbstractS
 
 end
 
-function NCVAR(x::ShortfallSamplesResult{N,L,T,P}, alpha::Float64) where {N,L,T,P}
+function NCVAR(x::ShortfallSamplesResult{N,L,T,P}, cvar::CVAR) where {N,L,T,P}
     demand = sum(x.regions.load)
 
-    estimate = x[]
-    tail_losses = estimate[estimate .>= quantile(estimate, alpha)]
-
-    cvar = if !isempty(tail_losses)
-        MeanEstimate(tail_losses)
-    else
-        MeanEstimate(0.)
-    end
-
     ncvar = if demand > 0
-        div(cvar, demand/1e6)
+        div(cvar.cvar, demand/1e6)
     else
         MeanEstimate(0.)
     end
 
-    return NCVAR{N,L,T}(ncvar, alpha)
+    return NCVAR(ncvar, cvar.alpha)
 
 end
 
-function NCVAR(x::ShortfallSamplesResult{N,L,T,P}, alpha::Float64, r::AbstractString) where {N,L,T,P}
+function NCVAR(x::ShortfallSamplesResult{N,L,T,P}, cvar::CVAR, r::AbstractString) where {N,L,T,P}
     i_r = findfirstunique(x.regions.names, r)
     demand = sum(x.regions.load[i_r, :])
 
-    estimate = x[r]
-    tail_losses = estimate[estimate .>= quantile(estimate, alpha)]
-
-    cvar = if !isempty(tail_losses)
-        MeanEstimate(tail_losses)
-    else
-        MeanEstimate(0.)
-    end
-
     ncvar = if demand > 0
-        div(cvar, demand/1e6)
+        div(cvar.cvar, demand/1e6)
     else
         MeanEstimate(0.)
     end
     
-    return NCVAR{N,L,T}(ncvar, alpha)
+    return NCVAR(ncvar, cvar.alpha)
 
 end
 
diff --git a/PRASCore.jl/src/Results/metrics.jl b/PRASCore.jl/src/Results/metrics.jl
index b856e733..3acf14ac 100644
--- a/PRASCore.jl/src/Results/metrics.jl
+++ b/PRASCore.jl/src/Results/metrics.jl
@@ -206,15 +206,15 @@ 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{N,L,T<:Period} <: ReliabilityMetric
+struct NCVAR <: ReliabilityMetric
 
     ncvar::MeanEstimate
     alpha::Float64
     
-    function NCVAR{N,L,T}(ncvar::MeanEstimate, alpha::Float64) where {N,L,T<:Period}
+    function NCVAR(ncvar::MeanEstimate, alpha::Float64)
         val(ncvar) >= 0 || throw(DomainError(
             "$val is not a valid NCVAR"))
-        new{N,L,T}(ncvar, alpha)
+        new(ncvar, alpha)
     end
 
 end

From fac883c4a5d5d05d4b46ba77de9f77c01a4e7b57 Mon Sep 17 00:00:00 2001
From: Abdelrahman Ayad <abdelrahman.ayad@gmail.com>
Date: Mon, 2 Feb 2026 11:01:40 -0700
Subject: [PATCH 04/23] broadcast CVAR and NCVAR - to check

---
 PRASCore.jl/src/Results/Results.jl | 18 +++++++++++++++++-
 1 file changed, 17 insertions(+), 1 deletion(-)

diff --git a/PRASCore.jl/src/Results/Results.jl b/PRASCore.jl/src/Results/Results.jl
index 93f99f27..1025bd4c 100644
--- a/PRASCore.jl/src/Results/Results.jl
+++ b/PRASCore.jl/src/Results/Results.jl
@@ -13,7 +13,7 @@ export
 
     # Metrics
     ReliabilityMetric, LOLE, EUE, NEUE,
-    val, stderror, CVAR, NCVAR,
+    val, stderror, CVAR, NCVAR, 
 
     # Result specifications
     Shortfall, ShortfallSamples,
@@ -80,6 +80,22 @@ NEUE(x::AbstractShortfallResult, r::AbstractString, ::Colon) =
 NEUE(x::AbstractShortfallResult, ::Colon, ::Colon) =
     NEUE.(x, x.regions.names, permutedims(x.timestamps))
 
+##TODO: Need to check these CVAR and NCVAR implementations
+CVAR(x::AbstractShortfallResult, ::Colon, t::ZonedDateTime) =
+    CVAR.(x, x.regions.names, t)
+
+CVAR(x::AbstractShortfallResult, r::AbstractString, ::Colon) =
+    CVAR.(x, r, x.timestamps)
+
+CVAR(x::AbstractShortfallResult, ::Colon, ::Colon) =
+    CVAR.(x, x.regions.names, permutedims(x.timestamps))    
+
+NCVAR(x::AbstractShortfallResult, r::AbstractString, ::Colon) =
+    NCVAR.(x, r, x.timestamps)
+
+NCVAR(x::AbstractShortfallResult, ::Colon, ::Colon) =
+    NCVAR.(x, x.regions.names, permutedims(x.timestamps))
+
 include("Shortfall.jl")
 include("ShortfallSamples.jl")
 

From 37ab22c7ee056977128e57f2336701eb7cc52c13 Mon Sep 17 00:00:00 2001
From: Abdelrahman Ayad <abdelrahman.ayad@gmail.com>
Date: Mon, 2 Feb 2026 11:02:01 -0700
Subject: [PATCH 05/23] broadcast CVAR and NCVAR - to check

---
 PRASCore.jl/src/Results/Results.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/PRASCore.jl/src/Results/Results.jl b/PRASCore.jl/src/Results/Results.jl
index 1025bd4c..d78c46ab 100644
--- a/PRASCore.jl/src/Results/Results.jl
+++ b/PRASCore.jl/src/Results/Results.jl
@@ -13,7 +13,7 @@ export
 
     # Metrics
     ReliabilityMetric, LOLE, EUE, NEUE,
-    val, stderror, CVAR, NCVAR, 
+    val, stderror, CVAR, NCVAR,
 
     # Result specifications
     Shortfall, ShortfallSamples,

From 2091f496b69542f25c2c9370665820b2d23e0d55 Mon Sep 17 00:00:00 2001
From: Abdelrahman Ayad <abdelrahman.ayad@gmail.com>
Date: Mon, 2 Feb 2026 11:10:47 -0700
Subject: [PATCH 06/23] add tests for CVAR and NCVAR

---
 PRAS.jl/test/runtests.jl | 43 +++++++++++++++++++++++++++++++---------
 1 file changed, 34 insertions(+), 9 deletions(-)

diff --git a/PRAS.jl/test/runtests.jl b/PRAS.jl/test/runtests.jl
index 7189f62e..d6fd641c 100644
--- a/PRAS.jl/test/runtests.jl
+++ b/PRAS.jl/test/runtests.jl
@@ -1,15 +1,40 @@
 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
+    @test val(eue) isa Float64
+    @test stderror(eue) isa Float64
+    @test val(neue) isa Float64
+    @test stderror(neue) 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
\ No newline at end of file

From 04e6f25684ec6e2cabc288dd4470fd6679591e74 Mon Sep 17 00:00:00 2001
From: Abdelrahman Ayad <abdelrahman.ayad@gmail.com>
Date: Mon, 2 Feb 2026 11:49:27 -0700
Subject: [PATCH 07/23] CVAR and NCVAR shortfall and metrics test

---
 PRASCore.jl/test/Results/metrics.jl   | 24 ++++++++++++++++
 PRASCore.jl/test/Results/shortfall.jl | 40 +++++++++++++++++++++++++++
 2 files changed, 64 insertions(+)

diff --git a/PRASCore.jl/test/Results/metrics.jl b/PRASCore.jl/test/Results/metrics.jl
index fb7bf796..174ed0c9 100644
--- a/PRASCore.jl/test/Results/metrics.jl
+++ b/PRASCore.jl/test/Results/metrics.jl
@@ -72,4 +72,28 @@
 
     end
 
+    @testset "CVAR" begin
+
+        cvar1 = CVAR{2,1,Hour,MWh}(MeanEstimate(1.2), 0.95)
+        @test string(cvar1) == "CVAR@0.95 = 1.20000 MWh/2h"
+
+        cvar2 = CVAR{1,2,Year,GWh}(MeanEstimate(17.2, 1.3), 0.95)
+        @test string(cvar2) == "CVAR@0.95 = 17±1 GWh/2y"
+
+        @test_throws DomainError CVAR{1,1,Hour,MWh}(MeanEstimate(-1.2), 0.95)
+
+    end    
+
+    @testset "NCVAR" begin
+
+        ncvar1 = NCVAR(MeanEstimate(1.2), 0.95)
+        @test string(ncvar1) == "NCVAR@0.95 = 1.20000 ppm"
+
+        ncvar2 = NCVAR(MeanEstimate(17.2, 1.3), 0.95)
+        @test string(ncvar2) == "NCVAR@0.95 = 17±1 ppm"
+
+        @test_throws DomainError NCVAR(MeanEstimate(-1.2), 0.95)
+
+    end
+
 end
diff --git a/PRASCore.jl/test/Results/shortfall.jl b/PRASCore.jl/test/Results/shortfall.jl
index 9f5d4a01..8e39ce6b 100644
--- a/PRASCore.jl/test/Results/shortfall.jl
+++ b/PRASCore.jl/test/Results/shortfall.jl
@@ -100,6 +100,7 @@ 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
+    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 +124,16 @@ end
     @test val(neue) ≈ mean(result[]) / load*1e6
     @test stderror(neue) ≈ std(result[]) / sqrt(DD.nsamples) / load*1e6
 
+    cvar = CVAR(result, alpha)
+    estimate = result[];
+    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 length(result[r]) == DD.nsamples
@@ -142,10 +153,22 @@ end
     @test val(region_neue) ≈ mean(result[r]) / load*1e6
     @test stderror(region_neue) ≈ std(result[r]) / sqrt(DD.nsamples) / load*1e6
 
+    region_cvar = CVAR(result, alpha, r)
+    region_estimate = result[r];
+    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(result, alpha, r_bad)
+    @test_throws BoundsError NCVAR(result, region_cvar, r_bad)
 
     # Period-specific
 
@@ -161,9 +184,16 @@ end
     @test val(period_eue) ≈ mean(result[t])
     @test stderror(period_eue) ≈ std(result[t]) / sqrt(DD.nsamples)
 
+    period_cvar = CVAR(result, alpha, t)
+    period_estimate = result[t];
+    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(result, alpha, t_bad)
 
     # Region + period-specific
 
@@ -180,6 +210,12 @@ end
     @test val(regionperiod_eue) ≈ mean(result[r, t])
     @test stderror(regionperiod_eue) ≈ std(result[r, t]) / sqrt(DD.nsamples)
 
+    regionperiod_cvar = CVAR(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 +228,8 @@ end
     @test_throws BoundsError EUE(result, r_bad, t)
     @test_throws BoundsError EUE(result, r_bad, t_bad)
 
+    @test_throws BoundsError CVAR(result, alpha, r, t_bad)
+    @test_throws BoundsError CVAR(result, alpha, r_bad, t)
+    @test_throws BoundsError CVAR(result, alpha, r_bad, t_bad)
+
 end

From 3547a87b7f776e8101594df07e7617e6083717ef Mon Sep 17 00:00:00 2001
From: Abdelrahman Ayad <abdelrahman.ayad@gmail.com>
Date: Mon, 2 Feb 2026 12:17:33 -0700
Subject: [PATCH 08/23] broadcast CVAR and NCVAR functions

---
 PRASCore.jl/src/Results/Results.jl | 22 ++++++++++------------
 1 file changed, 10 insertions(+), 12 deletions(-)

diff --git a/PRASCore.jl/src/Results/Results.jl b/PRASCore.jl/src/Results/Results.jl
index d78c46ab..881bc872 100644
--- a/PRASCore.jl/src/Results/Results.jl
+++ b/PRASCore.jl/src/Results/Results.jl
@@ -80,22 +80,20 @@ NEUE(x::AbstractShortfallResult, r::AbstractString, ::Colon) =
 NEUE(x::AbstractShortfallResult, ::Colon, ::Colon) =
     NEUE.(x, x.regions.names, permutedims(x.timestamps))
 
-##TODO: Need to check these CVAR and NCVAR implementations
-CVAR(x::AbstractShortfallResult, ::Colon, t::ZonedDateTime) =
-    CVAR.(x, x.regions.names, t)
+CVAR(x::AbstractShortfallResult, alpha::Float64, ::Colon, t::ZonedDateTime) =
+    CVAR.(x, alpha, x.regions.names, t)
 
-CVAR(x::AbstractShortfallResult, r::AbstractString, ::Colon) =
-    CVAR.(x, r, x.timestamps)
+CVAR(x::AbstractShortfallResult, alpha::Float64, r::AbstractString, ::Colon) =
+    CVAR.(x, alpha, r, x.timestamps)
 
-CVAR(x::AbstractShortfallResult, ::Colon, ::Colon) =
-    CVAR.(x, x.regions.names, permutedims(x.timestamps))    
+CVAR(x::AbstractShortfallResult, alpha::Float64, ::Colon, ::Colon) =
+    CVAR.(x, alpha, x.regions.names, permutedims(x.timestamps))    
 
-NCVAR(x::AbstractShortfallResult, r::AbstractString, ::Colon) =
-    NCVAR.(x, r, x.timestamps)
-
-NCVAR(x::AbstractShortfallResult, ::Colon, ::Colon) =
-    NCVAR.(x, 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")
 

From 1140338d809100c4fcb2c35c986378efddebed6d Mon Sep 17 00:00:00 2001
From: Abdelrahman Ayad <abdelrahman.ayad@gmail.com>
Date: Wed, 4 Feb 2026 17:51:55 -0700
Subject: [PATCH 09/23] report var value and add field to distinguish between
 period and sample CVAR

---
 PRASCore.jl/src/Results/metrics.jl | 21 +++++++++++++++------
 1 file changed, 15 insertions(+), 6 deletions(-)

diff --git a/PRASCore.jl/src/Results/metrics.jl b/PRASCore.jl/src/Results/metrics.jl
index 3acf14ac..e8a42147 100644
--- a/PRASCore.jl/src/Results/metrics.jl
+++ b/PRASCore.jl/src/Results/metrics.jl
@@ -178,11 +178,13 @@ struct CVAR{N,L,T<:Period,E<:EnergyUnit} <: ReliabilityMetric
 
     cvar::MeanEstimate
     alpha::Float64
+    var::Float64
+    cvar_type::String
     
-    function CVAR{N,L,T,E}(cvar::MeanEstimate, alpha::Float64) where {N,L,T<:Period,E<:EnergyUnit}
+    function CVAR{N,L,T,E}(cvar::MeanEstimate, alpha::Float64, var::Float64, cvar_type::String) where {N,L,T<:Period,E<:EnergyUnit}
         val(cvar) >= 0 || throw(DomainError(
             "$val is not a valid CVAR"))
-        new{N,L,T,E}(cvar, alpha)
+        new{N,L,T,E}(cvar, alpha, var, cvar_type)
     end
 
 end
@@ -191,9 +193,15 @@ 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, " ",
+    if x.cvar_type == "period"
+        print(io, "CVAR@$(x.alpha) = ", x.cvar, " ",
+          unitsymbol(E))
+    else
+        print(io, "CVAR@$(x.alpha) = ", x.cvar, " ",
           unitsymbol(E), "/", N*L == 1 ? "" : N*L, unitsymbol(T))
+    end
+    # print(io, "CVAR@$(x.alpha) = ", x.cvar, " ",
+    #       unitsymbol(E), "/", N*L == 1 ? "" : N*L, unitsymbol(T))
 
 end
 
@@ -210,11 +218,12 @@ struct NCVAR <: ReliabilityMetric
 
     ncvar::MeanEstimate
     alpha::Float64
+    var::Float64
     
-    function NCVAR(ncvar::MeanEstimate, alpha::Float64)
+    function NCVAR(ncvar::MeanEstimate, alpha::Float64, var::Float64)
         val(ncvar) >= 0 || throw(DomainError(
             "$val is not a valid NCVAR"))
-        new(ncvar, alpha)
+        new(ncvar, alpha, var)
     end
 
 end

From dac226257f45012d2c325f9bcfa4f9f247cc0968 Mon Sep 17 00:00:00 2001
From: Abdelrahman Ayad <abdelrahman.ayad@gmail.com>
Date: Wed, 4 Feb 2026 17:52:22 -0700
Subject: [PATCH 10/23] Shortfall based CVAR - worst periods

---
 PRASCore.jl/src/Results/Shortfall.jl | 101 +++++++++++++++++++++++++++
 1 file changed, 101 insertions(+)

diff --git a/PRASCore.jl/src/Results/Shortfall.jl b/PRASCore.jl/src/Results/Shortfall.jl
index 0b00d7d2..d960623e 100644
--- a/PRASCore.jl/src/Results/Shortfall.jl
+++ b/PRASCore.jl/src/Results/Shortfall.jl
@@ -292,6 +292,107 @@ 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_mean[:]
+    var = quantile(estimate, alpha)
+    tail_losses = estimate[estimate .>= var]
+
+    cvar = if !isempty(tail_losses)
+        MeanEstimate(tail_losses)
+    else
+        MeanEstimate(0.)
+    end
+
+    cvar_type = "period"
+
+    return CVAR{N,L,T,E}(cvar, alpha, var, cvar_type)
+  
+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_mean[i_r, :]
+    var = quantile(estimate, alpha)
+    tail_losses = estimate[estimate .>= var]
+
+    cvar = if !isempty(tail_losses)
+        MeanEstimate(tail_losses)
+    else
+        MeanEstimate(0.)
+    end
+
+    cvar_type = "period"
+
+    return CVAR{N,L,T,E}(cvar, alpha, var, cvar_type)
+  
+end
+
+function CVAR(x::ShortfallResult{N,L,T,E}, alpha::Float64, t::ZonedDateTime) where {N,L,T,E}
+    estimate = x.shortfall_mean[t, :]
+    var = quantile(estimate, alpha)
+    tail_losses = estimate[estimate .>= var]
+
+    cvar = if !isempty(tail_losses)
+        MeanEstimate(tail_losses)
+    else
+        MeanEstimate(0.)
+    end
+
+    cvar_type = "period"
+
+    return CVAR{N,L,T,E}(cvar, alpha, var, cvar_type)
+  
+end
+
+function CVAR(x::ShortfallResult{N,L,T,E}, alpha::Float64, r::AbstractString, t::ZonedDateTime) where {N,L,T,E}
+    estimate = x.shortfall_mean[r, t, :]
+    var = quantile(estimate, alpha)
+    tail_losses = estimate[estimate .>= var]
+
+    cvar = if !isempty(tail_losses)
+        MeanEstimate(tail_losses)
+    else
+        MeanEstimate(0.)
+    end
+
+    cvar_type = "period"
+
+    return CVAR{N,L,T,E}(cvar, alpha, var, cvar_type)
+  
+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},

From a8dd8dd421dde3662103b4edf57913c6fab59def Mon Sep 17 00:00:00 2001
From: Abdelrahman Ayad <abdelrahman.ayad@gmail.com>
Date: Wed, 4 Feb 2026 17:52:59 -0700
Subject: [PATCH 11/23] add type field for CVAR and store var

---
 PRASCore.jl/src/Results/ShortfallSamples.jl | 47 ++++++++++++++-------
 1 file changed, 31 insertions(+), 16 deletions(-)

diff --git a/PRASCore.jl/src/Results/ShortfallSamples.jl b/PRASCore.jl/src/Results/ShortfallSamples.jl
index 3e78e087..79b1627e 100644
--- a/PRASCore.jl/src/Results/ShortfallSamples.jl
+++ b/PRASCore.jl/src/Results/ShortfallSamples.jl
@@ -187,21 +187,25 @@ end
 
 function CVAR(x::ShortfallSamplesResult{N,L,T,P,E}, alpha::Float64) where {N,L,T,P,E}
     estimate = x[]
-    tail_losses = estimate[estimate .>= quantile(estimate, alpha)]
+    var = quantile(estimate, alpha)
+    tail_losses = estimate[estimate .>= var]
 
     cvar = if !isempty(tail_losses)
         MeanEstimate(tail_losses)
     else
         MeanEstimate(0.)
     end
+
+    cvar_type = "sample"
     
-    return CVAR{N,L,T,E}(cvar, alpha)
+    return CVAR{N,L,T,E}(cvar, alpha, var, cvar_type)
 
 end
 
 function CVAR(x::ShortfallSamplesResult{N,L,T,P,E}, alpha::Float64, r::AbstractString) where {N,L,T,P,E}
     estimate = x[r]
-    tail_losses = estimate[estimate .>= quantile(estimate, alpha)]
+    var = quantile(estimate, alpha)
+    tail_losses = estimate[estimate .>= var]
 
     cvar = if !isempty(tail_losses)
         MeanEstimate(tail_losses)
@@ -209,13 +213,16 @@ function CVAR(x::ShortfallSamplesResult{N,L,T,P,E}, alpha::Float64, r::AbstractS
         MeanEstimate(0.)
     end
     
-    return CVAR{N,L,T,E}(cvar, alpha)
+    cvar_type = "sample"
+
+    return CVAR{N,L,T,E}(cvar, alpha, var, cvar_type)
 
 end
 
 function CVAR(x::ShortfallSamplesResult{N,L,T,P,E}, alpha::Float64, t::ZonedDateTime) where {N,L,T,P,E}
     estimate = x[t]
-    tail_losses = estimate[estimate .>= quantile(estimate, alpha)]
+    var = quantile(estimate, alpha)
+    tail_losses = estimate[estimate .>= var]
 
     cvar = if !isempty(tail_losses)
         MeanEstimate(tail_losses)
@@ -223,13 +230,16 @@ function CVAR(x::ShortfallSamplesResult{N,L,T,P,E}, alpha::Float64, t::ZonedDate
         MeanEstimate(0.)
     end
     
-    return CVAR{N,L,T,E}(cvar, alpha)
+    cvar_type = "sample"
+    
+    return CVAR{N,L,T,E}(cvar, alpha, var, cvar_type)
 
 end
 
 function CVAR(x::ShortfallSamplesResult{N,L,T,P,E}, alpha::Float64, r::AbstractString, t::ZonedDateTime) where {N,L,T,P,E}
     estimate = x[r, t]
-    tail_losses = estimate[estimate .>= quantile(estimate, alpha)]
+    var = quantile(estimate, alpha)
+    tail_losses = estimate[estimate .>= var]
 
     cvar = if !isempty(tail_losses)
         MeanEstimate(tail_losses)
@@ -237,20 +247,23 @@ function CVAR(x::ShortfallSamplesResult{N,L,T,P,E}, alpha::Float64, r::AbstractS
         MeanEstimate(0.)
     end
     
-    return CVAR{N,L,T,E}(cvar, alpha)
+    cvar_type = "sample"
 
+    return CVAR{N,L,T,E}(cvar, alpha, var, cvar_type)
 end
 
 function NCVAR(x::ShortfallSamplesResult{N,L,T,P}, cvar::CVAR) where {N,L,T,P}
     demand = sum(x.regions.load)
 
-    ncvar = if demand > 0
-        div(cvar.cvar, demand/1e6)
+    if demand > 0
+        ncvar = div(cvar.cvar, demand/1e6)
+        var = div(cvar.var, demand/1e6)
     else
-        MeanEstimate(0.)
+        ncvar = MeanEstimate(0.)
+        var = MeanEstimate(0.)
     end
 
-    return NCVAR(ncvar, cvar.alpha)
+    return NCVAR(ncvar, cvar.alpha, var)
 
 end
 
@@ -258,13 +271,15 @@ function NCVAR(x::ShortfallSamplesResult{N,L,T,P}, cvar::CVAR, r::AbstractString
     i_r = findfirstunique(x.regions.names, r)
     demand = sum(x.regions.load[i_r, :])
 
-    ncvar = if demand > 0
-        div(cvar.cvar, demand/1e6)
+    if demand > 0
+        ncvar = div(cvar.cvar, demand/1e6)
+        var = div(cvar.var, demand/1e6)
     else
-        MeanEstimate(0.)
+        ncvar = MeanEstimate(0.)
+        var = MeanEstimate(0.)
     end
     
-    return NCVAR(ncvar, cvar.alpha)
+    return NCVAR(ncvar, cvar.alpha, var)
 
 end
 

From a8c3e7fea6d0b348910af932f21446feeb725255 Mon Sep 17 00:00:00 2001
From: Abdelrahman Ayad <abdelrahman.ayad@gmail.com>
Date: Fri, 20 Feb 2026 01:05:10 -0700
Subject: [PATCH 12/23] add ShortfallResult tests

---
 PRAS.jl/test/runtests.jl | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/PRAS.jl/test/runtests.jl b/PRAS.jl/test/runtests.jl
index d6fd641c..db44b014 100644
--- a/PRAS.jl/test/runtests.jl
+++ b/PRAS.jl/test/runtests.jl
@@ -10,10 +10,19 @@ using Test
     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
 
 @testset "ShortfallSamplesResult" begin

From 04395999529b6f1d20aa2fffa4ab3fb0111d2b48 Mon Sep 17 00:00:00 2001
From: Abdelrahman Ayad <abdelrahman.ayad@gmail.com>
Date: Fri, 20 Feb 2026 01:05:57 -0700
Subject: [PATCH 13/23] metrics.jl update

---
 PRASCore.jl/src/Results/metrics.jl | 27 +++++++++++++++------------
 1 file changed, 15 insertions(+), 12 deletions(-)

diff --git a/PRASCore.jl/src/Results/metrics.jl b/PRASCore.jl/src/Results/metrics.jl
index e8a42147..7bcc544d 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
@@ -179,12 +183,17 @@ struct CVAR{N,L,T<:Period,E<:EnergyUnit} <: ReliabilityMetric
     cvar::MeanEstimate
     alpha::Float64
     var::Float64
-    cvar_type::String
+    period_cvar::MeanEstimate
+    period_var::Float64
     
-    function CVAR{N,L,T,E}(cvar::MeanEstimate, alpha::Float64, var::Float64, cvar_type::String) where {N,L,T<:Period,E<:EnergyUnit}
+    function CVAR{N,L,T,E}(cvar::MeanEstimate,
+                           alpha::Float64,
+                           var::Float64,
+                           period_cvar::MeanEstimate,
+                           period_var::Float64) where {N,L,T<:Period,E<:EnergyUnit}
         val(cvar) >= 0 || throw(DomainError(
             "$val is not a valid CVAR"))
-        new{N,L,T,E}(cvar, alpha, var, cvar_type)
+        new{N,L,T,E}(cvar, alpha, var, period_cvar, period_var)
     end
 
 end
@@ -193,15 +202,9 @@ 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}
-    if x.cvar_type == "period"
-        print(io, "CVAR@$(x.alpha) = ", x.cvar, " ",
-          unitsymbol(E))
-    else
-        print(io, "CVAR@$(x.alpha) = ", x.cvar, " ",
+    
+    print(io, "CVAR@$(x.alpha) = ", x.cvar, " ",
           unitsymbol(E), "/", N*L == 1 ? "" : N*L, unitsymbol(T))
-    end
-    # print(io, "CVAR@$(x.alpha) = ", x.cvar, " ",
-    #       unitsymbol(E), "/", N*L == 1 ? "" : N*L, unitsymbol(T))
 
 end
 

From 6dded9d7cb3e1ec24ac242e755e80b11de22da3c Mon Sep 17 00:00:00 2001
From: Abdelrahman Ayad <abdelrahman.ayad@gmail.com>
Date: Fri, 20 Feb 2026 01:06:37 -0700
Subject: [PATCH 14/23] add unservedload_sample and unservedload_region_sample
 to ShortfallAccumulator

---
 PRASCore.jl/src/Results/Shortfall.jl | 90 +++++++++++++++-------------
 1 file changed, 47 insertions(+), 43 deletions(-)

diff --git a/PRASCore.jl/src/Results/Shortfall.jl b/PRASCore.jl/src/Results/Shortfall.jl
index d960623e..f0c32c26 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::Union{Vector{Float64}, Vector{Int}}
+    unservedload_region_sample::Union{Matrix{Float64}, Matrix{Int}}
+
 end
 
 function accumulator(
@@ -90,6 +94,8 @@ function accumulator(
 
     unservedload_total_currentsim = 0
     unservedload_region_currentsim = zeros(Int, nregions)
+    unservedload_sample = zeros(Float64, nsamples)
+    unservedload_region_sample = zeros(Float64, 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
@@ -149,6 +159,9 @@ struct ShortfallResult{N, L, T <: Period, E <: EnergyUnit, S} <:
     shortfall_period_std::Vector{Float64}
     shortfall_regionperiod_std::Matrix{Float64}
 
+    shortfall_samples::Union{Vector{Float64}, Vector{Int}}
+    shortfall_region_samples::Union{Matrix{Float64}, Matrix{Int}}
+
     function ShortfallResult{N,L,T,E,S}(
         nsamples::Union{Int,Nothing},
         regions::Regions,
@@ -165,7 +178,9 @@ 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}
+        shortfall_regionperiod_std::Matrix{Float64},
+        shortfall_samples::Union{Vector{Float64}, Vector{Int}},
+        shortfall_region_samples::Union{Matrix{Float64}, Matrix{Int}}
     ) where {N,L,T<:Period,E<:EnergyUnit,S <: Union{Shortfall,DemandResponseShortfall}}
 
         isnothing(nsamples) || nsamples > 0 ||
@@ -185,7 +200,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 +212,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
 
@@ -293,7 +311,8 @@ 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_mean[:]
+
+    estimate = x.shortfall_samples
     var = quantile(estimate, alpha)
     tail_losses = estimate[estimate .>= var]
 
@@ -301,35 +320,26 @@ function CVAR(x::ShortfallResult{N,L,T,E}, alpha::Float64) where {N,L,T,E}
         MeanEstimate(tail_losses)
     else
         MeanEstimate(0.)
-    end
-
-    cvar_type = "period"
-
-    return CVAR{N,L,T,E}(cvar, alpha, var, cvar_type)
-  
-end
-
-function CVAR(x::ShortfallResult{N,L,T,E}, alpha::Float64, r::AbstractString) where {N,L,T,E}
+    end                                                           
 
-    i_r = findfirstunique(x.regions.names, r)
-    estimate = x.shortfall_mean[i_r, :]
-    var = quantile(estimate, alpha)
-    tail_losses = estimate[estimate .>= var]
+    period_estimate = x.shortfall_mean[:]
+    period_var = quantile(period_estimate, alpha)
+    period_tail_losses = period_estimate[period_estimate .>= period_var]
 
-    cvar = if !isempty(tail_losses)
-        MeanEstimate(tail_losses)
+    period_cvar = if !isempty(period_tail_losses)
+        MeanEstimate(period_tail_losses)
     else
         MeanEstimate(0.)
-    end
-
-    cvar_type = "period"
+    end                                                               
 
-    return CVAR{N,L,T,E}(cvar, alpha, var, cvar_type)
+    return CVAR{N,L,T,E}(cvar, alpha, var, period_cvar, period_var)
   
 end
 
-function CVAR(x::ShortfallResult{N,L,T,E}, alpha::Float64, t::ZonedDateTime) where {N,L,T,E}
-    estimate = x.shortfall_mean[t, :]
+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]
 
@@ -339,26 +349,17 @@ function CVAR(x::ShortfallResult{N,L,T,E}, alpha::Float64, t::ZonedDateTime) whe
         MeanEstimate(0.)
     end
 
-    cvar_type = "period"
-
-    return CVAR{N,L,T,E}(cvar, alpha, var, cvar_type)
-  
-end
+    period_estimate = x.shortfall_mean[i_r, :]
+    period_var = quantile(period_estimate, alpha)
+    period_tail_losses = period_estimate[period_estimate .>= period_var]
 
-function CVAR(x::ShortfallResult{N,L,T,E}, alpha::Float64, r::AbstractString, t::ZonedDateTime) where {N,L,T,E}
-    estimate = x.shortfall_mean[r, t, :]
-    var = quantile(estimate, alpha)
-    tail_losses = estimate[estimate .>= var]
-
-    cvar = if !isempty(tail_losses)
-        MeanEstimate(tail_losses)
+    period_cvar = if !isempty(period_tail_losses)
+        MeanEstimate(period_tail_losses)
     else
         MeanEstimate(0.)
-    end
-
-    cvar_type = "period"
+    end   
 
-    return CVAR{N,L,T,E}(cvar, alpha, var, cvar_type)
+    return CVAR{N,L,T,E}(cvar, alpha, var, period_cvar, period_var)
   
 end
 
@@ -418,6 +419,8 @@ function finalize(
     ue_region_std .*= p2e
     ue_period_std .*= p2e
     ue_regionperiod_std .*= p2e
+    ue_sample = acc.unservedload_sample .* p2e
+    ue_region_sample = acc.unservedload_region_sample .* p2e
 
     return ShortfallResult{N,L,T,E,S}(
         nsamples, system.regions, system.timestamps,
@@ -425,6 +428,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

From 5613f0b284faebbb1a04ff5ee9c6edfddc8ea33d Mon Sep 17 00:00:00 2001
From: Abdelrahman Ayad <abdelrahman.ayad@gmail.com>
Date: Fri, 20 Feb 2026 01:07:07 -0700
Subject: [PATCH 15/23] update CVAR implementation in ShortfallSamples to
 include energy and capacity calculations

---
 PRASCore.jl/src/Results/ShortfallSamples.jl | 39 ++++++++++++++-------
 1 file changed, 27 insertions(+), 12 deletions(-)

diff --git a/PRASCore.jl/src/Results/ShortfallSamples.jl b/PRASCore.jl/src/Results/ShortfallSamples.jl
index 79b1627e..eda2134f 100644
--- a/PRASCore.jl/src/Results/ShortfallSamples.jl
+++ b/PRASCore.jl/src/Results/ShortfallSamples.jl
@@ -186,9 +186,10 @@ function NEUE(x::ShortfallSamplesResult, r::AbstractString)
 end
 
 function CVAR(x::ShortfallSamplesResult{N,L,T,P,E}, alpha::Float64) where {N,L,T,P,E}
+    
     estimate = x[]
     var = quantile(estimate, alpha)
-    tail_losses = estimate[estimate .>= var]
+    tail_losses = estimate[estimate .> var]
 
     cvar = if !isempty(tail_losses)
         MeanEstimate(tail_losses)
@@ -196,9 +197,18 @@ function CVAR(x::ShortfallSamplesResult{N,L,T,P,E}, alpha::Float64) where {N,L,T
         MeanEstimate(0.)
     end
 
-    cvar_type = "sample"
+    period_estimate = x.shortfall[:]
+    period_var = quantile(period_estimate, alpha)
+    period_tail_losses = period_estimate[period_estimate .>= period_var]
+
+    period_cvar = if !isempty(period_tail_losses)
+        MeanEstimate(period_tail_losses)
+    else
+        MeanEstimate(0.)
+    end     
+
     
-    return CVAR{N,L,T,E}(cvar, alpha, var, cvar_type)
+    return CVAR{N,L,T,E}(cvar, alpha, var, period_cvar, period_var)
 
 end
 
@@ -213,9 +223,18 @@ function CVAR(x::ShortfallSamplesResult{N,L,T,P,E}, alpha::Float64, r::AbstractS
         MeanEstimate(0.)
     end
     
-    cvar_type = "sample"
+    i_r = findfirstunique(x.regions.names, r)
+    period_estimate = x.shortfall[i_r, :, :][:]
+    period_var = quantile(period_estimate, alpha)
+    period_tail_losses = period_estimate[period_estimate .>= period_var]
 
-    return CVAR{N,L,T,E}(cvar, alpha, var, cvar_type)
+    period_cvar = if !isempty(period_tail_losses)
+        MeanEstimate(period_tail_losses)
+    else
+        MeanEstimate(0.)
+    end  
+
+    return CVAR{N,L,T,E}(cvar, alpha, var, period_cvar, period_var)
 
 end
 
@@ -229,10 +248,8 @@ function CVAR(x::ShortfallSamplesResult{N,L,T,P,E}, alpha::Float64, t::ZonedDate
     else
         MeanEstimate(0.)
     end
-    
-    cvar_type = "sample"
-    
-    return CVAR{N,L,T,E}(cvar, alpha, var, cvar_type)
+
+    return CVAR{N,L,T,E}(cvar, alpha, var, cvar, var)
 
 end
 
@@ -247,9 +264,7 @@ function CVAR(x::ShortfallSamplesResult{N,L,T,P,E}, alpha::Float64, r::AbstractS
         MeanEstimate(0.)
     end
     
-    cvar_type = "sample"
-
-    return CVAR{N,L,T,E}(cvar, alpha, var, cvar_type)
+    return CVAR{N,L,T,E}(cvar, alpha, var, cvar, var)
 end
 
 function NCVAR(x::ShortfallSamplesResult{N,L,T,P}, cvar::CVAR) where {N,L,T,P}

From 5befde49102772029e2518ee1a9d68bcde47686b Mon Sep 17 00:00:00 2001
From: Abdelrahman Ayad <abdelrahman.ayad@gmail.com>
Date: Fri, 20 Feb 2026 01:07:27 -0700
Subject: [PATCH 16/23] update to the recording and reset! functions

---
 PRASCore.jl/src/Simulations/recording.jl | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

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

From 914e54124832c86877892cfe470dcd97b4e7e3d5 Mon Sep 17 00:00:00 2001
From: Abdelrahman Ayad <abdelrahman.ayad@gmail.com>
Date: Fri, 20 Feb 2026 01:07:43 -0700
Subject: [PATCH 17/23] udpate TestData with the CVAR values

---
 PRASCore.jl/src/Systems/TestData.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/PRASCore.jl/src/Systems/TestData.jl b/PRASCore.jl/src/Systems/TestData.jl
index 6b5ba264..07c2cfc2 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 = 11.84
 
 ## 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 = 30.44
 
 
 # Single-Region System B, with storage

From 5bf96cf53ea2ce6e5dc452036382c650a7e51de2 Mon Sep 17 00:00:00 2001
From: Abdelrahman Ayad <abdelrahman.ayad@gmail.com>
Date: Fri, 20 Feb 2026 01:08:01 -0700
Subject: [PATCH 18/23] update DummyData with CVAR values

---
 PRASCore.jl/test/dummydata.jl | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/PRASCore.jl/test/dummydata.jl b/PRASCore.jl/test/dummydata.jl
index 06f7d139..49bbcd00 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)
@@ -32,6 +33,8 @@ d1 = rand()
 d1_resource = rand(nresources)
 d1_period = rand(nperiods)
 d1_resourceperiod = rand(nresources, nperiods)
+d1_sample = rand(0:999, nsamples)
+d1_resourcesample = rand(0:999, nresources, nsamples)
 
 d2 = rand()
 d2_resource = rand(nresources)

From a97980126e861771f4cd3a204e94eb22dbc33903 Mon Sep 17 00:00:00 2001
From: Abdelrahman Ayad <abdelrahman.ayad@gmail.com>
Date: Fri, 20 Feb 2026 01:08:21 -0700
Subject: [PATCH 19/23] include CVAR tests

---
 PRASCore.jl/test/Results/metrics.jl      | 12 +++++------
 PRASCore.jl/test/Results/shortfall.jl    | 27 +++++++++++++++++++++++-
 PRASCore.jl/test/Simulations/runtests.jl |  7 ++++++
 3 files changed, 39 insertions(+), 7 deletions(-)

diff --git a/PRASCore.jl/test/Results/metrics.jl b/PRASCore.jl/test/Results/metrics.jl
index 174ed0c9..08232c04 100644
--- a/PRASCore.jl/test/Results/metrics.jl
+++ b/PRASCore.jl/test/Results/metrics.jl
@@ -74,25 +74,25 @@
 
     @testset "CVAR" begin
 
-        cvar1 = CVAR{2,1,Hour,MWh}(MeanEstimate(1.2), 0.95)
+        cvar1 = CVAR{2,1,Hour,MWh}(MeanEstimate(1.2), 0.95, 1.2, MeanEstimate(1.0), 1.2)
         @test string(cvar1) == "CVAR@0.95 = 1.20000 MWh/2h"
 
-        cvar2 = CVAR{1,2,Year,GWh}(MeanEstimate(17.2, 1.3), 0.95)
+        cvar2 = CVAR{1,2,Year,GWh}(MeanEstimate(17.2, 1.3), 0.95, 1.2, MeanEstimate(17.2, 1.3), 1.2)
         @test string(cvar2) == "CVAR@0.95 = 17±1 GWh/2y"
 
-        @test_throws DomainError CVAR{1,1,Hour,MWh}(MeanEstimate(-1.2), 0.95)
+        @test_throws DomainError CVAR{1,1,Hour,MWh}(MeanEstimate(-1.2), 0.95, 1.2, MeanEstimate(-1.2), 1.2)
 
     end    
 
     @testset "NCVAR" begin
 
-        ncvar1 = NCVAR(MeanEstimate(1.2), 0.95)
+        ncvar1 = NCVAR(MeanEstimate(1.2), 0.95, 1.2)
         @test string(ncvar1) == "NCVAR@0.95 = 1.20000 ppm"
 
-        ncvar2 = NCVAR(MeanEstimate(17.2, 1.3), 0.95)
+        ncvar2 = NCVAR(MeanEstimate(17.2, 1.3), 0.95, 1.3)
         @test string(ncvar2) == "NCVAR@0.95 = 17±1 ppm"
 
-        @test_throws DomainError NCVAR(MeanEstimate(-1.2), 0.95)
+        @test_throws DomainError NCVAR(MeanEstimate(-1.2), 0.95, -1.2)
 
     end
 
diff --git a/PRASCore.jl/test/Results/shortfall.jl b/PRASCore.jl/test/Results/shortfall.jl
index 8e39ce6b..296e8118 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.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(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,22 @@
     @test val(region_neue) ≈ first(result[r]) / load*1e6
     @test stderror(region_neue) ≈ last(result[r]) / sqrt(DD.nsamples) / load*1e6
 
+    region_cvar = CVAR(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(result, alpha, r_bad)
+    @test_throws BoundsError NCVAR(result, region_cvar, r_bad)
 
     # Period-specific
 
diff --git a/PRASCore.jl/test/Simulations/runtests.jl b/PRASCore.jl/test/Simulations/runtests.jl
index 3c518634..fe1ebf8f 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(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(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(shortfall_1b, alpha, "Region"),
+                          TestData.singlenode_b_cvar, nstderr_tol)                                                    
 
         @test all(LOLE.(shortfall_1b, timestamps_b) .≈
                   LOLE.(shortfall2_1b, timestamps_b))

From 35579e23b001df1295708207a7c67372bd5b0cf8 Mon Sep 17 00:00:00 2001
From: Abdelrahman Ayad <abdelrahman.ayad@gmail.com>
Date: Mon, 23 Feb 2026 12:19:55 -0700
Subject: [PATCH 20/23] rename capacity_cvar

---
 PRASCore.jl/src/Results/metrics.jl | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/PRASCore.jl/src/Results/metrics.jl b/PRASCore.jl/src/Results/metrics.jl
index 7bcc544d..91a5dbc1 100644
--- a/PRASCore.jl/src/Results/metrics.jl
+++ b/PRASCore.jl/src/Results/metrics.jl
@@ -172,7 +172,7 @@ end
 """
     CVAR
 
-`CVAR` reports conditional value at risk of shortfalls.
+`CVAR` reports conditional value at risk of shortfalls, for total unserved energy and capacity shortfalls.
 
 Contains both the estimated value itself as well as the standard error
 of that estimate, which can be extracted with `val` and `stderror`,
@@ -183,17 +183,17 @@ struct CVAR{N,L,T<:Period,E<:EnergyUnit} <: ReliabilityMetric
     cvar::MeanEstimate
     alpha::Float64
     var::Float64
-    period_cvar::MeanEstimate
-    period_var::Float64
+    capacity_cvar::MeanEstimate
+    capacity_var::Float64
     
     function CVAR{N,L,T,E}(cvar::MeanEstimate,
                            alpha::Float64,
                            var::Float64,
-                           period_cvar::MeanEstimate,
-                           period_var::Float64) where {N,L,T<:Period,E<:EnergyUnit}
+                           capacity_cvar::MeanEstimate,
+                           capacity_var::Float64) where {N,L,T<:Period,E<:EnergyUnit}
         val(cvar) >= 0 || throw(DomainError(
             "$val is not a valid CVAR"))
-        new{N,L,T,E}(cvar, alpha, var, period_cvar, period_var)
+        new{N,L,T,E}(cvar, alpha, var, capacity_cvar, capacity_var)
     end
 
 end
@@ -211,7 +211,7 @@ end
 """
     NCVAR
 
-`NCVAR` reports normalized conditional value at risk of shortfalls.
+`NCVAR` reports normalized conditional value at risk of shortfalls, for total unserved energy and capacity shortfalls.
 
 Contains both the estimated value itself as well as the standard error
 of that estimate, which can be extracted with `val` and `stderror`,

From 0959c2339cfd38a590019d8ba08dcafbd63d522d Mon Sep 17 00:00:00 2001
From: Abdelrahman Ayad <abdelrahman.ayad@gmail.com>
Date: Mon, 23 Feb 2026 12:21:20 -0700
Subject: [PATCH 21/23] add capacity_cvar to shortfall and shortfallsamples

---
 PRASCore.jl/src/Results/Shortfall.jl        | 56 +++++++++++----------
 PRASCore.jl/src/Results/ShortfallSamples.jl | 33 ++++++------
 2 files changed, 47 insertions(+), 42 deletions(-)

diff --git a/PRASCore.jl/src/Results/Shortfall.jl b/PRASCore.jl/src/Results/Shortfall.jl
index f0c32c26..1d95fd1a 100644
--- a/PRASCore.jl/src/Results/Shortfall.jl
+++ b/PRASCore.jl/src/Results/Shortfall.jl
@@ -68,8 +68,8 @@ mutable struct ShortfallAccumulator{S} <: ResultAccumulator{Shortfall}
     unservedload_region_currentsim::Vector{Int}
 
     # Sample-level UE for current simulation
-    unservedload_sample::Union{Vector{Float64}, Vector{Int}}
-    unservedload_region_sample::Union{Matrix{Float64}, Matrix{Int}}
+    unservedload_sample::Vector{Int}
+    unservedload_region_sample::Matrix{Int}
 
 end
 
@@ -94,8 +94,8 @@ function accumulator(
 
     unservedload_total_currentsim = 0
     unservedload_region_currentsim = zeros(Int, nregions)
-    unservedload_sample = zeros(Float64, nsamples)
-    unservedload_region_sample = zeros(Float64, nregions, nsamples)
+    unservedload_sample = zeros(Int, nsamples)
+    unservedload_region_sample = zeros(Int, nregions, nsamples)
 
     return ShortfallAccumulator{S}(
         periodsdropped_total, periodsdropped_region,
@@ -151,16 +151,16 @@ 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::Union{Vector{Float64}, Vector{Int}}
-    shortfall_region_samples::Union{Matrix{Float64}, Matrix{Int}}
+    shortfall_samples::Vector{Float64}
+    shortfall_region_samples::Matrix{Float64}
 
     function ShortfallResult{N,L,T,E,S}(
         nsamples::Union{Int,Nothing},
@@ -175,12 +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_samples::Union{Vector{Float64}, Vector{Int}},
-        shortfall_region_samples::Union{Matrix{Float64}, Matrix{Int}}
+        shortfall_samples::Vector{Float64},
+        shortfall_region_samples::Matrix{Float64},
     ) where {N,L,T<:Period,E<:EnergyUnit,S <: Union{Shortfall,DemandResponseShortfall}}
 
         isnothing(nsamples) || nsamples > 0 ||
@@ -210,7 +211,7 @@ struct ShortfallResult{N, L, T <: Period, E <: EnergyUnit, S} <:
             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_samples,
             shortfall_region_samples)
@@ -322,17 +323,17 @@ function CVAR(x::ShortfallResult{N,L,T,E}, alpha::Float64) where {N,L,T,E}
         MeanEstimate(0.)
     end                                                           
 
-    period_estimate = x.shortfall_mean[:]
-    period_var = quantile(period_estimate, alpha)
-    period_tail_losses = period_estimate[period_estimate .>= period_var]
+    capacity_estimate = x.capacity_shortfall_mean[:]
+    capacity_var = quantile(capacity_estimate, alpha)
+    capacity_tail_losses = capacity_estimate[capacity_estimate .>= capacity_var]
 
-    period_cvar = if !isempty(period_tail_losses)
-        MeanEstimate(period_tail_losses)
+    capacity_cvar = if !isempty(capacity_tail_losses)
+        MeanEstimate(capacity_tail_losses)
     else
         MeanEstimate(0.)
     end                                                               
 
-    return CVAR{N,L,T,E}(cvar, alpha, var, period_cvar, period_var)
+    return CVAR{N,L,T,E}(cvar, alpha, var, capacity_cvar, capacity_var)
   
 end
 
@@ -349,17 +350,17 @@ function CVAR(x::ShortfallResult{N,L,T,E}, alpha::Float64, r::AbstractString) wh
         MeanEstimate(0.)
     end
 
-    period_estimate = x.shortfall_mean[i_r, :]
-    period_var = quantile(period_estimate, alpha)
-    period_tail_losses = period_estimate[period_estimate .>= period_var]
+    capacity_estimate = x.capacity_shortfall_mean[i_r, :]
+    capacity_var = quantile(capacity_estimate, alpha)
+    capacity_tail_losses = capacity_estimate[capacity_estimate .>= capacity_var]
 
-    period_cvar = if !isempty(period_tail_losses)
-        MeanEstimate(period_tail_losses)
+    capacity_cvar = if !isempty(capacity_tail_losses)
+        MeanEstimate(capacity_tail_losses)
     else
         MeanEstimate(0.)
-    end   
+    end 
 
-    return CVAR{N,L,T,E}(cvar, alpha, var, period_cvar, period_var)
+    return CVAR{N,L,T,E}(cvar, alpha, var, capacity_cvar, capacity_var)
   
 end
 
@@ -414,21 +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 = acc.unservedload_sample .* p2e
-    ue_region_sample = acc.unservedload_region_sample .* 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_regionperiod_mean, capacity_shortfall_mean, ue_total_std,
         ue_region_std, ue_period_std, ue_regionperiod_std,
-        ue_sample, ue_region_sample)
+        ue_sample, ue_region_sample, )
 
 end
diff --git a/PRASCore.jl/src/Results/ShortfallSamples.jl b/PRASCore.jl/src/Results/ShortfallSamples.jl
index eda2134f..e0ab7342 100644
--- a/PRASCore.jl/src/Results/ShortfallSamples.jl
+++ b/PRASCore.jl/src/Results/ShortfallSamples.jl
@@ -82,6 +82,7 @@ struct ShortfallSamplesResult{N,L,T<:Period,P<:PowerUnit,E<:EnergyUnit, S} <: Ab
     timestamps::StepRange{ZonedDateTime,T}
 
     shortfall::Array{Int,3} # r x t x s
+    capacity_shortfall::Array{Int,2} # r x s
 
 end
 
@@ -197,18 +198,17 @@ function CVAR(x::ShortfallSamplesResult{N,L,T,P,E}, alpha::Float64) where {N,L,T
         MeanEstimate(0.)
     end
 
-    period_estimate = x.shortfall[:]
-    period_var = quantile(period_estimate, alpha)
-    period_tail_losses = period_estimate[period_estimate .>= period_var]
+    capacity_estimate = x.capacity_shortfall[:]
+    capacity_var = quantile(capacity_estimate, alpha)
+    capacity_tail_losses = capacity_estimate[capacity_estimate .>= capacity_var]
 
-    period_cvar = if !isempty(period_tail_losses)
-        MeanEstimate(period_tail_losses)
+    capacity_cvar = if !isempty(capacity_tail_losses)
+        MeanEstimate(capacity_tail_losses)
     else
         MeanEstimate(0.)
     end     
 
-    
-    return CVAR{N,L,T,E}(cvar, alpha, var, period_cvar, period_var)
+    return CVAR{N,L,T,E}(cvar, alpha, var, capacity_cvar, capacity_var)
 
 end
 
@@ -224,17 +224,17 @@ function CVAR(x::ShortfallSamplesResult{N,L,T,P,E}, alpha::Float64, r::AbstractS
     end
     
     i_r = findfirstunique(x.regions.names, r)
-    period_estimate = x.shortfall[i_r, :, :][:]
-    period_var = quantile(period_estimate, alpha)
-    period_tail_losses = period_estimate[period_estimate .>= period_var]
+    capacity_estimate = x.capacity_shortfall[i_r, :]
+    capacity_var = quantile(capacity_estimate, alpha)
+    capacity_tail_losses = capacity_estimate[capacity_estimate .>= capacity_var]
 
-    period_cvar = if !isempty(period_tail_losses)
-        MeanEstimate(period_tail_losses)
+    capacity_cvar = if !isempty(capacity_tail_losses)
+        MeanEstimate(capacity_tail_losses)
     else
         MeanEstimate(0.)
-    end  
+    end     
 
-    return CVAR{N,L,T,E}(cvar, alpha, var, period_cvar, period_var)
+    return CVAR{N,L,T,E}(cvar, alpha, var, capacity_cvar, capacity_var)
 
 end
 
@@ -303,7 +303,10 @@ function finalize(
     system::SystemModel{N,L,T,P,E},
 ) where {N,L,T,P,E,S<:Union{ShortfallSamples,DemandResponseShortfallSamples}}
 
+    n_regions = length(system.regions)
+    p2e = conversionfactor(L,T,P,E)
+    max_capacity_shortfall = reshape(maximum(acc.shortfall, dims=[2]) ./ p2e, n_regions, :)
     return ShortfallSamplesResult{N,L,T,P,E,S}(
-        system.regions, system.timestamps, acc.shortfall)
+        system.regions, system.timestamps, acc.shortfall, max_capacity_shortfall)
 
 end

From 59e7983afdae55a49a5672cfa1dd989763e6f510 Mon Sep 17 00:00:00 2001
From: Abdelrahman Ayad <abdelrahman.ayad@gmail.com>
Date: Mon, 23 Feb 2026 12:21:47 -0700
Subject: [PATCH 22/23] update Dummydata

---
 PRASCore.jl/test/dummydata.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/PRASCore.jl/test/dummydata.jl b/PRASCore.jl/test/dummydata.jl
index 49bbcd00..12f96184 100644
--- a/PRASCore.jl/test/dummydata.jl
+++ b/PRASCore.jl/test/dummydata.jl
@@ -28,13 +28,14 @@ testperiod = periods[testperiod_idx]
 notaperiod = ZonedDateTime(2010,1,1,0,tz)
 
 d = rand(0:999, nresources, nperiods, nsamples)
+cap_d = rand(0:999, nresources, nsamples)
 
 d1 = rand()
 d1_resource = rand(nresources)
 d1_period = rand(nperiods)
 d1_resourceperiod = rand(nresources, nperiods)
-d1_sample = rand(0:999, nsamples)
-d1_resourcesample = rand(0:999, nresources, nsamples)
+d1_sample = rand(nsamples) * 999
+d1_resourcesample = rand(nresources, nsamples) * 999
 
 d2 = rand()
 d2_resource = rand(nresources)

From 7a610f66c615b777e5f2e7e4fd74d02b5dcee1d0 Mon Sep 17 00:00:00 2001
From: Abdelrahman Ayad <abdelrahman.ayad@gmail.com>
Date: Mon, 23 Feb 2026 12:34:24 -0700
Subject: [PATCH 23/23] CVAR capacity tests

---
 PRASCore.jl/test/Results/shortfall.jl | 36 +++++++++++++++++++++++++--
 1 file changed, 34 insertions(+), 2 deletions(-)

diff --git a/PRASCore.jl/test/Results/shortfall.jl b/PRASCore.jl/test/Results/shortfall.jl
index 296e8118..c9df817e 100644
--- a/PRASCore.jl/test/Results/shortfall.jl
+++ b/PRASCore.jl/test/Results/shortfall.jl
@@ -10,7 +10,7 @@
         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.d3_resourceperiod, DD.d3_resourceperiod, 
         DD.d4, DD.d4_resource, DD.d4_period, DD.d4_resourceperiod,
         DD.d1_sample, DD.d1_resourcesample)
 
@@ -37,6 +37,14 @@
     @test val(cvar) ≈ mean(tail_losses)
     @test stderror(cvar) ≈ std(tail_losses) / sqrt(length(tail_losses))
 
+    capacity_cvar = cvar.capacity_cvar
+    cap_shortfal = vec(reshape(result.capacity_shortfall_mean, 1, :))
+    var = quantile(cap_shortfal, alpha)
+    tail_losses = cap_shortfal[cap_shortfal .>= var]
+    capacity_cvar_check = mean(tail_losses)
+    @test val(capacity_cvar) ≈ capacity_cvar_check
+    @test stderror(capacity_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
@@ -64,6 +72,14 @@
     @test val(region_cvar) ≈ mean(region_tail_losses)
     @test stderror(region_cvar) ≈ std(region_tail_losses) / sqrt(length(region_tail_losses))
 
+    region_capacity_cvar = region_cvar.capacity_cvar
+    region_cap_shortfal = result.capacity_shortfall_mean[r_idx, :]
+    var = quantile(region_cap_shortfal, alpha)
+    tail_losses = region_cap_shortfal[region_cap_shortfal .>= var]
+    region_capacity_cvar_check = mean(tail_losses)
+    @test val(region_capacity_cvar) ≈ region_capacity_cvar_check
+    @test stderror(region_capacity_cvar) ≈ std(tail_losses) / sqrt(length(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
@@ -128,7 +144,7 @@ end
     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)
+        Regions{N,MW}(DD.resourcenames, DD.resource_vals), DD.periods, DD.d, DD.cap_d)
 
     # Overall
 
@@ -155,6 +171,14 @@ end
     @test val(cvar) ≈ mean(tail_losses)
     @test stderror(cvar) ≈ std(tail_losses) / sqrt(length(tail_losses))
 
+    capacity_cvar = cvar.capacity_cvar
+    cap_shortfal = vec(reshape(result.capacity_shortfall, 1, :))
+    var = quantile(cap_shortfal, alpha)
+    tail_losses = cap_shortfal[cap_shortfal .>= var]
+    capacity_cvar_check = mean(tail_losses)
+    @test val(capacity_cvar) ≈ capacity_cvar_check
+    @test stderror(capacity_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
@@ -184,6 +208,14 @@ end
     @test val(region_cvar) ≈ mean(region_tail_losses)
     @test stderror(region_cvar) ≈ std(region_tail_losses) / sqrt(length(region_tail_losses))
 
+    region_capacity_cvar = region_cvar.capacity_cvar
+    region_cap_shortfal = result.capacity_shortfall[r_idx, :]
+    var = quantile(region_cap_shortfal, alpha)
+    tail_losses = region_cap_shortfal[region_cap_shortfal .>= var]
+    region_capacity_cvar_check = mean(tail_losses)
+    @test val(region_capacity_cvar) ≈ region_capacity_cvar_check
+    @test stderror(region_capacity_cvar) ≈ std(tail_losses) / sqrt(length(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