Navigation :

How to aggregate ranvars across echelons with a Gaussian copula

This guide shows how to aggregate store-level ranvars into a warehouse-level ranvar while preserving shared demand shocks. It compares three cases: independence, full dependence as a stress test, and a Gaussian copula in between, then extends the calculation to many items, warehouses, and periods.

Use this pattern when the downstream ranvars already include local starting conditions such as store stock, and the upstream echelon must inherit that information.

Table of contents

Step 1: Build the downstream marginals

Model one ranvar per downstream location after local stock has already been consumed. In the example below, Stores.UncoveredDemand is the quantity that must still be supplied by the upstream echelon.

Step 2: Aggregate with three dependence assumptions

Paste the script below into a new Envision project and run it.

preview = show region { 1..8, 1..24 }

table Stores[store] = with
  [| as store, as StockOnHand, as Mu, as Dispersion |]
  [| "Paris" , 3, 6, 2.0 |]
  [| "Lyon"  , 1, 4, 1.5 |]
  [| "Lille" , 0, 2, 1.2 |]

Stores.Demand = negativeBinomial(Stores.Mu, Stores.Dispersion)
Stores.UncoveredDemand = max(Stores.Demand - Stores.StockOnHand, 0)

// Baseline: independent downstream locations.
warehouseCov0 = sum(Stores.UncoveredDemand)

// Stress test: fully dependent locations, using a common quantile draw.
nbSamples = 1000
montecarlo nbSamples with
  u = random.uniform(0, 1)
  total = sum(quantile(Stores.UncoveredDemand, u))
  sample warehouseCov1 = ranvar(total)

// Core method: Gaussian copula between the store marginals.
rho = 0.55
epsilon = 0.0001 // `normal.quantile` expects p in ]0, 1[

montecarlo nbSamples with
  commonU = max(epsilon, min(1 - epsilon, random.uniform(0, 1)))
  globalShock = normal.quantile(0, 1, commonU)

  Stores.LocalNoise = random.normal(0 into Stores, 1)
  Stores.Latent = rho * globalShock + sqrt(1 - rho^2) * Stores.LocalNoise
  Stores.U = normal.cdf(0, 1, Stores.Latent)

  total = sum(quantile(Stores.UncoveredDemand, Stores.U))
  sample warehouseCopula = ranvar(total)

show table "Store marginals" { .., 1..8 in preview } with
  Stores.store
  Stores.StockOnHand
  mean(Stores.UncoveredDemand) as "Mean uncovered"
  variance(Stores.UncoveredDemand) as "Var uncovered"
  quantile(Stores.UncoveredDemand, 0.9) as "P90 uncovered"

show summary "Warehouse comparison" { .., 10..21 in preview } with
  mean(warehouseCov0) as "Mean cov=0" // ~8.26
  variance(warehouseCov0) as "Var cov=0" // ~18.33
  quantile(warehouseCov0, 0.9) as "P90 cov=0" // 14
  mean(warehouseCov1) as "Mean cov=1" // ~7.82
  variance(warehouseCov1) as "Var cov=1" // ~44.40
  quantile(warehouseCov1, 0.9) as "P90 cov=1" // 17
  mean(warehouseCopula) as "Mean copula" // ~8.19
  variance(warehouseCopula) as "Var copula" // ~27.64
  quantile(warehouseCopula, 0.9) as "P90 copula" // 15

Run the script and verify:

Subtracting a store’s own StockOnHand cannot increase its demand quantiles.
The three warehouse ranvars should have similar means, while their spread differs markedly.
In the displayed example, the aggregate variance and P90 follow cov=0 < copula < cov=1. Check the quantiles for your own marginals; this P90 ordering is specific to the example.

This uses normal.quantile and normal.cdf to generate correlated uniforms without going through a discretized zedfunc representation of the standard normal CDF.

Step 3: Aggregate many item-warehouse-period groups efficiently

Use one table of downstream marginals and a grouping table for the upstream totals. Each input row should represent one item’s uncovered demand at one store for one period, together with the warehouse supplying it. Select the contributing locations and required periods before creating the groups.

Create Groups with by over the actual (Item, Warehouse, Period) keys. This gives each input row one upstream group and avoids building combinations that have no contributors. Use the same keys for any independence or full dependence comparisons.

Run the following standalone example. Its marginals already represent demand after local stock has been consumed; periods 1 and 2 represent two weeks. The dirac marginals are fixed quantities, so their contributions can be added outside the simulation.

table Rows small 100 = with
  [| as Item, as Warehouse, as Period, as Store, as UncoveredDemand |]
  [| "A", "North", 1, "N1", negativeBinomial(4, 2) |]
  [| "A", "North", 1, "N2", negativeBinomial(2, 2) |]
  [| "A", "North", 1, "N3", dirac(3) |]
  [| "A", "North", 2, "N1", negativeBinomial(8, 2) |]
  [| "A", "North", 2, "N2", negativeBinomial(4, 2) |]
  [| "A", "South", 1, "S1", negativeBinomial(3, 2) |]
  [| "A", "South", 1, "S2", dirac(2) |]
  [| "B", "North", 1, "N1", dirac(2) |]
  [| "B", "North", 1, "N2", dirac(0) |]

// Define the output groups before filtering out constant contributions.
table Groups = by [Rows.Item, Rows.Warehouse, Rows.Period]
Groups.Item = same(Rows.Item)
Groups.Warehouse = same(Rows.Warehouse)
Groups.Period = same(Rows.Period)
Groups.ExpectedMean = sum(mean(Rows.UncoveredDemand))

Rows.IsStochastic = variance(Rows.UncoveredDemand) > 0
Groups.Constant = sum(mean(Rows.UncoveredDemand)) when not Rows.IsStochastic

rho = 0.55
nbSamples = 1000
where Rows.IsStochastic
  montecarlo nbSamples with
    // One common shock per group, and independent noise per contributing row.
    Groups.Shock = random.normal(0 into Groups, 1)
    Rows.LocalNoise = random.normal(0 into Rows, 1)
    Rows.U = normal.cdf(0, 1, rho * Groups.Shock + sqrt(1 - rho^2) * Rows.LocalNoise)

    // Sum across locations, then collect each group's totals across samples.
    Groups.Total = sum(quantile(Rows.UncoveredDemand, Rows.U))
    sample Groups.SimulatedDemand = ranvar(Groups.Total)

Groups.Demand = Groups.SimulatedDemand + Groups.Constant
show table "Warehouse demand by item and period" with
  Groups.Item
  Groups.Warehouse
  Groups.Period
  Groups.ExpectedMean
  mean(Groups.Demand) as "Simulated mean"
  quantile(Groups.Demand, 0.9) as "P90"
  order by Groups.Item, Groups.Warehouse, Groups.Period

One run gives the following results, with means rounded to two decimals. Sampled values can vary between runs.

Item	Warehouse	Period	Expected mean	Simulated mean	P90
A	North	1	9.00	8.88	15
A	North	2	12.00	12.11	19
A	South	1	5.00	5.02	8
B	North	1	2.00	2.00	2

Check that the result has these four groups and that the simulated means are close to the sums of the input means. Item B’s demand is exactly 2 units, including at P90, because both of its contributors are deterministic.

Use the table relationship for both shocks and totals

The by declaration makes Groups upstream of Rows. Inside each sample, Groups.Shock broadcasts the appropriate common shock to each row, and Groups.Total = sum(...) aggregates sampled quantities back to their group. The sample statement then collects one empirical ranvar per group over all iterations. See data model and relationships for the underlying table rules.

Keep 0 into Groups and 0 into Rows in the random calls. The arguments determine where the draws occur: random.normal(0, 1) produces one scalar draw, even when assigned to a column or written inside a sum. Each sample needs one shared shock per group and an independent local draw per row. Here, rho is the loading on the common shock; the latent correlation between two different rows in the same group is rho^2. Step 4 addresses calibration.

Sum the sampled quantities inside the simulation. Summing the input ranvars directly with sum(Rows.UncoveredDemand) would apply the independence assumption used in Step 2.

Remove repeated scans and constant work

Prepare group membership once, before montecarlo, and sample all group totals in that block. An outer each over output groups, containing a Monte Carlo block that scans all input rows with matching-key predicates, repeats that scan for every group and every sample. Grouping first removes this group-count multiplier.

For S samples, N contributing rows, and G groups, the traversal work therefore changes roughly from S * G * N to S * (N + G), in addition to preparing the groups and processing the distributions. Actual runtime also depends on the marginals and other work in the script. Compare formulations with the same sample count and probability model.

When many marginals have zero variance, precompute their group totals as in the example and filter them out before entering montecarlo. Their quantiles always return the same quantity. Add those constants back after sampling. Keep marginals with positive variance even when most of their mass is at zero.

Create the groups before this stochastic filter. Filtering Rows leaves its upstream Groups table intact, including groups with no stochastic rows; their sampled sum is zero and their final ranvar contains only the constant total. If your output also needs groups with no input rows at all, retain a separate output table and fill unmatched demand with dirac(0).

Batch complete groups and check the quantities used downstream

Set bounds on the tables used by the simulation; see the small table limits. Size batches using the expanded item-location-period row count. Keep every contributor to an (Item, Warehouse, Period) group in the same batch. If earlier calculations split those contributors across files, bring their marginals together before the copula aggregation. Independently aggregating fragments and summing their ranvars would lose the intended dependence across fragments.

Check that output keys and constant totals are preserved. With rho = 0, the simulated group means and variances should approach those obtained by summing independent marginals as the sample count increases. Different execution layouts can use different draws, so compare their distributions with sampling tolerance.

Choose the sample count from the quantiles or decisions consumed downstream; follow How to choose a Monte Carlo sample count. Increasing this count reduces aggregation sampling noise, but cannot recover tail information absent from previously simulated input marginals. This recipe produces a separate marginal for each period; preserve or model joint scenarios separately if a decision requires dependence across periods.

Step 4: Recover `rho` from censored downstream time series

Use a separate script when you want to estimate the dependence parameter from a panel of synchronized downstream histories. The script below generates 100 store time series with negative-binomial marginals and a Gaussian copula, then applies stock caps to create days where demand is only partially observable.

After generation, the calibration operates on the observed sales only. The estimated parameter is the average pair correlation of the latent Gaussian layer. The rho used in the copula construction is recovered afterward as its square root.

trueRho = 0.55
truePairCorr = trueRho^2

numStores = 100
epsilon = 0.0001

startDate = date(2024, 1, 1)
endDate = date(2026, 9, 26)
keep span date = [startDate .. endDate]

table Stores max 100 = extend.range(numStores)
Stores.Mu = 2 + Stores.N / 20
Stores.Dispersion = max(1, Stores.Mu * 2)
Stores.DemandRanvar = negativeBinomial(Stores.Mu, Stores.Dispersion)
Stores.Cdf = cdf(Stores.DemandRanvar)

Day.CommonU = random.uniform(epsilon into Day, 1 - epsilon)
Day.Global = normal.quantile(0, 1, Day.CommonU)

table StoreDays max 200k = cross(Stores, Day)
StoreDays.Local = random.normal(0 into StoreDays, 1)
StoreDays.Latent = trueRho * Day.Global + sqrt(1 - truePairCorr) * StoreDays.Local
StoreDays.U = normal.cdf(0, 1, StoreDays.Latent)
StoreDays.Demand = quantile(Stores.DemandRanvar, StoreDays.U)

// On capped days, only the sales are observed, not the full demand.
StoreDays.StockCap = quantile(Stores.DemandRanvar, random.uniform(0.3 into StoreDays, 0.8))
StoreDays.Sales = min(StoreDays.Demand, StoreDays.StockCap)
StoreDays.IsCensored = StoreDays.Sales >= StoreDays.StockCap
StoreDays.IsObserved = not StoreDays.IsCensored

StoreDays.F = valueAt(Stores.Cdf, StoreDays.Sales)
StoreDays.P = int(Stores.DemandRanvar, StoreDays.Sales, StoreDays.Sales)
StoreDays.Uhat = max(epsilon, min(1 - epsilon, StoreDays.F - 0.5 * StoreDays.P))
StoreDays.Zhat = if StoreDays.IsObserved then normal.quantile(0, 1, StoreDays.Uhat) else 0

Day.ObservedStoreCount = count(StoreDays.IsObserved)
Day.SumZ = sum(StoreDays.Zhat) when StoreDays.IsObserved
Day.SumZ2 = sum(StoreDays.Zhat^2) when StoreDays.IsObserved
Day.PairSignal = if Day.ObservedStoreCount < 2 then 0
                 else (Day.SumZ^2 - Day.SumZ2) /
                      (Day.ObservedStoreCount * (Day.ObservedStoreCount - 1))

pairCorr = avg(Day.PairSignal) when (Day.ObservedStoreCount >= 2)

rho = sqrt(pairCorr)

show summary "Recovered dependence" with
  trueRho as "True rho" // 0.55
  rho as "Learned rho" // ~0.60
  truePairCorr as "True pair corr" // 0.3025
  pairCorr as "Estimated pair corr" // ~0.3599
  avg(Day.ObservedStoreCount) as "Avg observed stores" // ~49.7