Stochastic Walkthrough

Usage and Examples

This page demonstrates how to use the Stochastic simulation mode. Unlike EAD mode (which gives you a single expected value), Stochastic mode reveals the distribution of outcomes – including tail risk that expected-value analysis hides. This is critical when evaluating strategies that concentrate risk, like dike-heavy portfolios.

Prerequisites

If you are new to iCOW, start with the Quickstart for a gentler introduction. This page assumes familiarity with the basic concepts (Config, Scenario, Policy, Outcome).

Setup

Code

using ICOW.Stochastic
using Distributions
using Random
using CairoMakie
CairoMakie.activate!(; type="svg")
using Statistics

Configuring the City

StochasticConfig holds the same 28 city parameters as EADConfig. The defaults describe a Manhattan-like coastal city.

config = StochasticConfig()

StochasticConfig{Float64}(1.5e12, 30.0, 17.0, 2000.0, 43000.0, 1.75, 2.0, 3.0, 0.5, 1000.0, 1.1, 0.95, 1.0, 0.01, 1.25, 0.35, 0.115, 0.4, 3.0, 0.39, 0.03, 1.5, 0.95, 0.05, 1.1, 4.0e9, 1.0, 1.01)

Key parameters:

City value: $1.5T
City elevation: 17.0 m
Building height: 30.0 m
Seawall height: 1.75 m

Defining the Surge Scenario

Unlike EAD mode (which takes GEV distribution parameters), Stochastic mode takes a realized time series of annual maximum surge heights – one value per year. This represents a single “what if” trajectory of storms. Here we sample from the GEV distribution calibrated to The Battery, NYC per Ceres, Forest, and Keller (2019).

rng = MersenneTwister(123)
n_years = 100
surge_dist = GeneralizedExtremeValue(0.935, 0.260, 0.030)
surges = rand(rng, surge_dist, n_years)

mean_sea_level = collect(range(0.0; step=0.01, length=n_years))
scenario = StochasticScenario(; surges=surges, discount_rate=0.02, mean_sea_level=mean_sea_level)

StochasticScenario{Float64}(SimOptDecisions.TimeSeriesParameter{Float64, Int64}([1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  91, 92, 93, 94, 95, 96, 97, 98, 99, 100], [1.1339625979696317, 0.7481591972960504, 1.217981334268734, 1.1694619304987077, 0.9136931218354459, 0.830627039252128, 1.0282403363617671, 0.8503580515651172, 0.7091096434751376, 1.0585959659901716  …  1.2747293261432038, 0.6506622574547744, 1.1290517052957716, 1.030525101992608, 1.2353120253992311, 1.2212410088779981, 0.9132689259158091, 1.3664177784150835, 0.7058810531547546, 1.1045918098604222]), SimOptDecisions.TimeSeriesParameter{Float64, Int64}([1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  91, 92, 93, 94, 95, 96, 97, 98, 99, 100], [0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09  …  0.9, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99]), SimOptDecisions.ContinuousParameter{Float64}(0.02, (0.0, 1.0)))

fig = Figure(; size=(600, 300))
ax = Axis(fig[1, 1]; xlabel="Year", ylabel="Surge height (m)", title="Surge time series")
barplot!(ax, 1:n_years, surges; color=:steelblue, gap=0.1)
hlines!(ax, [config.H_seawall]; color=:red, linestyle=:dash, linewidth=2, label="Seawall ($(config.H_seawall)m)")
axislegend(ax; position=:rt)
fig

Surges above the seawall line cause flooding. The dashed red line shows the existing seawall protection – everything above it contributes to flood damage (after wave runup adjustment).

Tip

In practice, you generate many surge sequences (e.g., from a GEV distribution) and simulate each one to build a distribution of outcomes.

Choosing a Protection Strategy

StaticPolicy uses the same reparameterized fractions as EAD mode. Here are several archetypal strategies:

archetypes = [
    ("No additional protection", StaticPolicy(; a_frac=0.0, w_frac=0.0, b_frac=0.0, r_frac=0.0, P=0.0)),
    ("Dike only",                StaticPolicy(; a_frac=0.5, w_frac=0.0, b_frac=0.0, r_frac=0.0, P=0.0)),
    ("Withdraw + resist",        StaticPolicy(; a_frac=0.2, w_frac=0.5, b_frac=1.0, r_frac=0.3, P=0.7)),
    ("Mixed",                    StaticPolicy(; a_frac=0.3, w_frac=0.1, b_frac=0.2, r_frac=0.1, P=0.5)),
    ("High dike + resist",       StaticPolicy(; a_frac=0.5, w_frac=0.0, b_frac=0.1, r_frac=0.1, P=0.3)),
]
policies = last.(archetypes)
names = first.(archetypes)

5-element Vector{String}:
 "No additional protection"
 "Dike only"
 "Withdraw + resist"
 "Mixed"
 "High dike + resist"

Running a Simulation

policy = policies[4]  # Mixed strategy
rng = MersenneTwister(42)
outcome = simulate(config, scenario, policy, rng)
outcome

StochasticOutcome{Float64}(SimOptDecisions.ContinuousParameter{Float64}(5.455527272132463e10, (-Inf, Inf)), SimOptDecisions.ContinuousParameter{Float64}(6.572674410956048e8, (-Inf, Inf)))

Why Stochastic? The Role of Storm Variability

EAD mode integrates over the surge distribution analytically, giving a single expected damage number. But real cities experience a specific sequence of storms – some lucky (mostly small surges), some unlucky (a rare extreme event early on). Stochastic mode captures this by simulating many independent surge sequences, revealing the full distribution of outcomes.

n_scenarios = 200
rng_scenarios = MersenneTwister(42)
scenarios = [
    StochasticScenario(;
        surges=rand(rng_scenarios, surge_dist, n_years),
        discount_rate=0.02,
        mean_sea_level=mean_sea_level,
    )
    for _ in 1:n_scenarios
]

damages = [
    value(simulate(config, s, policy, MersenneTwister(i)).damage)
    for (i, s) in enumerate(scenarios)
]

200-element Vector{Float64}:
 8.081095184741836e8
 1.72222658843448e9
 1.518042868083947e9
 1.9485620011771386e9
 2.0866185831997738e9
 1.3076489119572062e9
 1.8572629766290944e9
 1.1507323415666895e9
 3.9580771288925967e9
 1.2799604351120772e9
 ⋮
 2.0962817172517223e9
 1.0918618951900034e9
 1.7501628562891257e9
 1.6228433183008301e9
 4.754463359129002e8
 1.1112133285527663e9
 3.6385450122275e9
 2.0017881937773347e9
 1.4953139733404613e9

fig = Figure(; size=(600, 350))
ax = Axis(fig[1, 1]; xlabel="Damage (\$B)", ylabel="Count", title="Damage distribution across $n_scenarios surge sequences")
hist!(ax, damages ./ 1e9; bins=30, color=(:steelblue, 0.7))
vlines!(ax, [mean(damages) / 1e9]; color=:red, linewidth=2, label="Mean")
vlines!(ax, [quantile(damages, 0.05) / 1e9]; color=:orange, linewidth=2, linestyle=:dash, label="5th pctile")
vlines!(ax, [quantile(damages, 0.95) / 1e9]; color=:purple, linewidth=2, linestyle=:dash, label="95th pctile")
axislegend(ax; position=:rt)
fig

Mean damage: $2.02B
5th percentile: $0.81B
95th percentile: $3.48B

Warning

The gap between the 5th and 95th percentiles reveals hidden risk that expected-value analysis misses. A policy that looks optimal under EAD may have unacceptable tail losses.

Comparing Policies

We use explore() to evaluate all five policies across many independent surge sequences. Each scenario is a different realization of storm history over the planning horizon:

result = explore(config, scenarios, policies; progress=false)
result

YAXArray Dataset
Shared Axes: 
None
Variables with additional axes:
  Additional Axes: 
  (↓ policy Sampled{Int64} 1:5 ForwardOrdered Regular Points,
  → scenario Sampled{Int64} 1:200 ForwardOrdered Regular Points)
  Variables: 
  damage, investment

  Additional Axes: 
  (↓ policy Sampled{Int64} 1:5 ForwardOrdered Regular Points)
  Variables: 
  policy_P, policy_a_frac, policy_b_frac, policy_r_frac, policy_w_frac

  Additional Axes: 
  (↓ scenario Sampled{Int64} 1:200 ForwardOrdered Regular Points)
  Variables: 
  scenario_discount_rate, scenario_mean_sea_level[100], scenario_mean_sea_level[10], scenario_mean_sea_level[11], scenario_mean_sea_level[12], scenario_mean_sea_level[13], scenario_mean_sea_level[14], scenario_mean_sea_level[15], scenario_mean_sea_level[16], scenario_mean_sea_level[17], scenario_mean_sea_level[18], scenario_mean_sea_level[19], scenario_mean_sea_level[1], scenario_mean_sea_level[20], scenario_mean_sea_level[21], scenario_mean_sea_level[22], scenario_mean_sea_level[23], scenario_mean_sea_level[24], scenario_mean_sea_level[25], scenario_mean_sea_level[26], scenario_mean_sea_level[27], scenario_mean_sea_level[28], scenario_mean_sea_level[29], scenario_mean_sea_level[2], scenario_mean_sea_level[30], scenario_mean_sea_level[31], scenario_mean_sea_level[32], scenario_mean_sea_level[33], scenario_mean_sea_level[34], scenario_mean_sea_level[35], scenario_mean_sea_level[36], scenario_mean_sea_level[37], scenario_mean_sea_level[38], scenario_mean_sea_level[39], scenario_mean_sea_level[3], scenario_mean_sea_level[40], scenario_mean_sea_level[41], scenario_mean_sea_level[42], scenario_mean_sea_level[43], scenario_mean_sea_level[44], scenario_mean_sea_level[45], scenario_mean_sea_level[46], scenario_mean_sea_level[47], scenario_mean_sea_level[48], scenario_mean_sea_level[49], scenario_mean_sea_level[4], scenario_mean_sea_level[50], scenario_mean_sea_level[51], scenario_mean_sea_level[52], scenario_mean_sea_level[53], scenario_mean_sea_level[54], scenario_mean_sea_level[55], scenario_mean_sea_level[56], scenario_mean_sea_level[57], scenario_mean_sea_level[58], scenario_mean_sea_level[59], scenario_mean_sea_level[5], scenario_mean_sea_level[60], scenario_mean_sea_level[61], scenario_mean_sea_level[62], scenario_mean_sea_level[63], scenario_mean_sea_level[64], scenario_mean_sea_level[65], scenario_mean_sea_level[66], scenario_mean_sea_level[67], scenario_mean_sea_level[68], scenario_mean_sea_level[69], scenario_mean_sea_level[6], scenario_mean_sea_level[70], scenario_mean_sea_level[71], scenario_mean_sea_level[72], scenario_mean_sea_level[73], scenario_mean_sea_level[74], scenario_mean_sea_level[75], scenario_mean_sea_level[76], scenario_mean_sea_level[77], scenario_mean_sea_level[78], scenario_mean_sea_level[79], scenario_mean_sea_level[7], scenario_mean_sea_level[80], scenario_mean_sea_level[81], scenario_mean_sea_level[82], scenario_mean_sea_level[83], scenario_mean_sea_level[84], scenario_mean_sea_level[85], scenario_mean_sea_level[86], scenario_mean_sea_level[87], scenario_mean_sea_level[88], scenario_mean_sea_level[89], scenario_mean_sea_level[8], scenario_mean_sea_level[90], scenario_mean_sea_level[91], scenario_mean_sea_level[92], scenario_mean_sea_level[93], scenario_mean_sea_level[94], scenario_mean_sea_level[95], scenario_mean_sea_level[96], scenario_mean_sea_level[97], scenario_mean_sea_level[98], scenario_mean_sea_level[99], scenario_mean_sea_level[9], scenario_surges[100], scenario_surges[10], scenario_surges[11], scenario_surges[12], scenario_surges[13], scenario_surges[14], scenario_surges[15], scenario_surges[16], scenario_surges[17], scenario_surges[18], scenario_surges[19], scenario_surges[1], scenario_surges[20], scenario_surges[21], scenario_surges[22], scenario_surges[23], scenario_surges[24], scenario_surges[25], scenario_surges[26], scenario_surges[27], scenario_surges[28], scenario_surges[29], scenario_surges[2], scenario_surges[30], scenario_surges[31], scenario_surges[32], scenario_surges[33], scenario_surges[34], scenario_surges[35], scenario_surges[36], scenario_surges[37], scenario_surges[38], scenario_surges[39], scenario_surges[3], scenario_surges[40], scenario_surges[41], scenario_surges[42], scenario_surges[43], scenario_surges[44], scenario_surges[45], scenario_surges[46], scenario_surges[47], scenario_surges[48], scenario_surges[49], scenario_surges[4], scenario_surges[50], scenario_surges[51], scenario_surges[52], scenario_surges[53], scenario_surges[54], scenario_surges[55], scenario_surges[56], scenario_surges[57], scenario_surges[58], scenario_surges[59], scenario_surges[5], scenario_surges[60], scenario_surges[61], scenario_surges[62], scenario_surges[63], scenario_surges[64], scenario_surges[65], scenario_surges[66], scenario_surges[67], scenario_surges[68], scenario_surges[69], scenario_surges[6], scenario_surges[70], scenario_surges[71], scenario_surges[72], scenario_surges[73], scenario_surges[74], scenario_surges[75], scenario_surges[76], scenario_surges[77], scenario_surges[78], scenario_surges[79], scenario_surges[7], scenario_surges[80], scenario_surges[81], scenario_surges[82], scenario_surges[83], scenario_surges[84], scenario_surges[85], scenario_surges[86], scenario_surges[87], scenario_surges[88], scenario_surges[89], scenario_surges[8], scenario_surges[90], scenario_surges[91], scenario_surges[92], scenario_surges[93], scenario_surges[94], scenario_surges[95], scenario_surges[96], scenario_surges[97], scenario_surges[98], scenario_surges[99], scenario_surges[9]

The result is a YAXArrays Dataset with dimensions (policy, scenario). Each scenario is an independent surge sequence, so the scenario dimension captures storm variability:

# Mean and 95th percentile damage across scenarios for each policy
mean_dmg = [mean(result[:damage].data[p, :]) for p in 1:length(policies)] ./ 1e9
p95_dmg = [quantile(result[:damage].data[p, :], 0.95) for p in 1:length(policies)] ./ 1e9

fig = Figure(; size=(600, 400))
ax = Axis(fig[1, 1]; xlabel="Mean Damage (\$B)", ylabel="95th Percentile Damage (\$B)",
    title="Mean vs Tail Risk by Policy")
scatter!(ax, mean_dmg, p95_dmg; markersize=14, color=:steelblue)
for (i, name) in enumerate(names)
    text!(ax, mean_dmg[i], p95_dmg[i]; text=name, fontsize=10, offset=(8, 4))
end
fig

Mean vs. Tail Risk

Policies that reduce mean damage do not always reduce tail risk proportionally. A dike-heavy strategy may have low average damage but catastrophic outcomes when the dike fails. The mixed strategy balances mean and tail performance by diversifying across defense types.

This distinction matters for decision making: a risk-neutral decision-maker cares only about expected cost and would use EAD mode. A risk-averse decision-maker weights bad outcomes more heavily and needs the distributional information that Stochastic mode provides. Which perspective is appropriate depends on the decision context – see Keller, Helgeson, and Srikrishnan (2021) for a broader discussion.

References

Ceres, Robert L., Chris E. Forest, and Klaus Keller. 2019. “Optimization of Multiple Storm Surge Risk Mitigation Strategies for an Island City On a Wedge.” Environmental Modelling & Software 119 (September): 341–53. https://doi.org/10.1016/j.envsoft.2019.06.011.

Keller, Klaus, Casey Helgeson, and Vivek Srikrishnan. 2021. “Climate Risk Management.” Annual Review of Earth and Planetary Sciences 49 (1): 95–116. https://doi.org/10.1146/annurev-earth-080320-055847.