25  Scenario Generation

“Scenarios are the rehearsal of possible futures. They are not predictions, but pathways that help us prepare for the unknown.” — Peter Schwartz

25.1 Chapter Overview

Scenario generation is the discipline of simulating coherent paths for the economic and market variables that feed valuation, risk, and planning models. You will often combine several stochastic sub-models (interest rates, equity returns, credit spreads, inflation, demographics) to create realistic joint scenarios. This chapter reviews commonly used building blocks and illustrates their Julia implementations so you can stitch together scenario engines for portfolio analytics, risk transfer pricing, or asset-liability management (ALM).

25.2 Common Use Cases of Scenario Generators

Scenario generators are widely used in risk management, investment analysis, and regulatory compliance to model potential future outcomes. For forecasting under the real-world (physical) measure P, use real-world (RW) scenarios. For valuation and hedging, use risk‑neutral (RN) scenarios under the martingale measure Q (e.g., with drift adjusted to the short rate minus dividends for equities and to fit the initial curve for interest rates).

25.2.1 Risk management, especially Value at Risk (VaR) & Expected Shortfall (ES).

RW scenario generators are used to simulate market movements to estimate potential portfolio losses. Basel III regulatory capital requirements have adopted these approaches.

25.2.2 Stress Testing & Regulatory Compliance

RW scenario generators can also be used to generate extreme but plausible market conditions to assess resilience, which is required by central banks and financial regulators (e.g., Federal Reserve and ECB).

25.2.3 Portfolio Optimization & Asset Allocation

RW scenario generators are used to simulate thousands of market conditions to determine optimal portfolio allocations which is commonly used in modern portfolio theory (MPT) and Black-Litterman models.

25.2.4 Pension & Insurance Risk Modeling

RW scenario generators can be used to simulate longevity risk, policyholder behavior, and interest rate movements. They are also used for economic capital estimation under uncertain economic scenarios.

25.2.5 Economic & Macro-Financial Forecasting

Central banks and institutions (e.g., IMF, World Bank) use RW scenario generators to predict macroeconomic trends.

25.2.6 Asset Pricing & Hedging

RN scenario generators help value options using stochastic models (e.g., Black-Scholes model). They can help simulate future stock price movements under different volatility conditions. They can also be used for hedging purposes to test how a portfolio performs under different inflation, interest rate, or commodity price scenarios.

25.2.7 Fixed Income & Interest Rate Modeling

Yield curve modeling uses RN scenarios to value bonds and interest rate derivatives. Swaps, swaptions, and credit default swaps (CDS) also rely on RN pricing. RN scenario generators can also simulate yield curves for bond and fixed-income pricing. Models like Cox-Ingersoll-Ross (CIR) or Hull-White generate future interest rate paths.

25.2.8 Regulatory & Accounting Valuations

IFRS 13 & fair value accounting uses RN models determine the market-consistent value of liabilities. Solvency II for insurers requires valuation of policyholder guarantees using RN scenarios.

25.3 Common Economic Scenario Generation Approaches

Economic scenario generation involves the development of plausible future economic scenarios to assess the potential impact on financial portfolios, investments, or decision-making processes. Various approaches are used to generate economic scenarios, such as adapting underlying stochastic differential equations (SDEs) for Monte Carlo scenario generation techniques.

25.3.1 Interest Rate Models

Short-rate models produce positive, mean-reverting simulated rates that can be bootstrapped into full yield curves. They remain popular in ALM where projecting discount curves, reinvestment rates, and collateral flows is critical.

25.3.1.1 Vasicek and Cox Ingersoll Ross (CIR)

The Vasicek model is a one-factor model commonly used for simulating interest rate scenarios. It describes the dynamics of short-term interest rates using a stochastic differential equation (SDE). In a Monte Carlo simulation, we can use the Vasicek model to generate multiple interest rate paths. The CIR model is an extension of the Vasicek model with non-constant volatility. It is a square‑root diffusion with state‑dependent volatility that helps keep rates positive. It addresses the issue of negative interest rates by ensuring that interest rates remain positive.

Vasicek is defined as

\[ dr(t) = \kappa (\theta - r(t)) \, dt + \sigma \, dW(t) \]

where

• \(r(t)\) is the short-term interest rate at time \(t\).
• \(κ\) is the speed of mean reversion, representing how quickly the interest rate reverts to its long-term mean.
• \(θ\) is the long-term mean or equilibrium level of the interest rate.
• \(σ\) is the volatility of the interest rate.
• \(dW(t)\) is a Wiener process or Brownian motion, representing a random shock.

And CIR is defined as

\[ dr(t) = \kappa (\theta - r(t)) \, dt + \sigma \sqrt{r(t)} \, dW(t) \]

where

• \(r(t)\) is the short-term interest rate at time \(t\).
• \(κ\) is the speed of mean reversion, representing how quickly the interest rate reverts to its long-term mean.
• \(θ\) is the long-term mean or equilibrium level of the interest rate.
• \(σ\) is the volatility of the interest rate.
• \(dW(t)\) is a Wiener process or Brownian motion, representing a random shock.

The following code shows a simplified implementation of a CIR model. The specification of \(dr\) can be changed to make it a Vasicek model.

using Random, CairoMakie

# Function to simulate CIR process
function cir_simulation(rng::AbstractRNG, κ, θ, σ, r₀, Δt, num_steps, num_sims)
    interest_rate_paths = zeros(num_steps, num_sims)
    st = sqrt(Δt)
    for j in 1:num_sims
        shocks = randn(rng, num_steps - 1)
        interest_rate_paths[1, j] = r₀
        for i in 2:num_steps
            # dr = κ * (θ - interest_rate_paths[i-1, j]) * Δt + σ * dW # for Vasicek
            dr = κ * (θ - interest_rate_paths[i-1, j]) * Δt + σ * sqrt(max(interest_rate_paths[i-1, j], 0.0)) * st * shocks[i - 1]
            interest_rate_paths[i, j] = max(interest_rate_paths[i-1, j] + dr, 0.0)  # Ensure non-negativity
        end
    end
    return interest_rate_paths
end

let
    # Set seed for reproducibility
    rng = MersenneTwister(1234)

    # CIR model parameters
    κ = 0.2 # Speed of mean reversion
    θ = 0.05 # Long-term mean
    σ = 0.1 # Volatility

    # Initial short-term interest rate
    r₀ = 0.03

    # Number of time steps and simulations
    num_steps = 252
    num_simulations = 1_000

    # Time increment
    Δt = 1 / 252
    # Run CIR simulation
    cir_paths = cir_simulation(rng, κ, θ, σ, r₀, Δt, num_steps, num_simulations)

    # Plot the simulated interest rate paths
    f = Figure()
    ax = Axis(f[1, 1], xlabel="Trading days", ylabel="Short rate", title="CIR paths")
    ti = 0:(num_steps-1)
    for i in 1:num_simulations
        lines!(ax, ti, cir_paths[:, i], color=:steelblue, linewidth=0.8)
    end
    f
end

The exact CIR process has a non-central chi-squared distribution for \(r(t)\) in closed form.

\[ \begin{aligned} d &= \frac{4 \kappa \theta}{\sigma^2} \\ c &= \frac{\sigma^2 \left( 1 - e^{-\kappa \Delta t} \right)}{4 \kappa} \\ \lambda &= \frac{4 \kappa e^{-\kappa \Delta t} r_t}{\sigma^2 \left( 1 - e^{-\kappa \Delta t} \right)} \\ r_{t+\Delta t} &\sim c \cdot \chi'^2_{d}(\lambda) \\ \end{aligned} \]

where \(\chi'^2_{d}(\lambda)\) means the non-central chi-square distribution with \(d\) degrees of freedom and a non-centrality parameter \(λ\). The following code shows an exact implementation of a CIR model.

using Random, CairoMakie, Distributions

function cir_exact_paths(rng::AbstractRNG, κ, θ, σ, r₀, T, num_steps, num_sims)
    Δt = T / num_steps
    d = 4 * κ * θ / (σ * σ)
    c = (σ * σ * (1 - exp(-κ * Δt))) / (4 * κ)
    rates = zeros(Float64, num_steps + 1, num_sims)
    time = range(0, T; length=num_steps + 1)
    rates[1, :] .= r₀
    for j in 1:num_sims
        r = r₀
        for i in 2:length(time)
            λ = (4 * κ * exp(-κ * Δt) * r) / (σ * σ * (1 - exp(-κ * Δt)))
            r = c * rand(rng, NoncentralChisq(d, λ))
            rates[i, j] = r
        end
    end
    return time, rates
end

let
    rng = MersenneTwister(1234)

    # CIR parameters
    κ = 0.5 # mean reversion speed
    θ = 0.04 # long-term mean
    σ = 0.1 # volatility coefficient
    r₀ = 0.03 # initial rate
    T = 5.0 # years
    nsteps = 500 # time steps
    nsims = 1000 # number of simulated paths

    tgrid, rpaths = cir_exact_paths(rng, κ, θ, σ, r₀, T, nsteps, nsims)
    
    # Plot
    fig = Figure(size=(800, 500))
    ax = Axis(fig[1, 1], xlabel="Time (years)", ylabel="Short rate", title="Exact CIR Process Simulation (Non-central Chi-square)")
    for path in 1:nsims
        lines!(ax, tgrid, rpaths[:, path], color=:steelblue, linewidth=0.8)
    end
    fig
end

25.3.1.2 Hull-White

The Hull-White model is a one-factor model that extends the Vasicek model by allowing the mean reversion and volatility parameters to be time-dependent. It is commonly used for pricing interest rate derivatives. Brace-Gatarek-Musiela (BGM) Model extends the Hull-White model to incorporate more factors. It is one of the Libor Market Model (LMM) that describes the evolution of forward rates. It allows for the modeling of both the short-rate and the entire yield curve. It is defined as

\[ dr(t) = (\theta(t) - a r(t)) \, dt + \sigma(t) \, dW(t) \]

where

• \(r(t)\) is the short-term interest rate at time \(t\).
• \(θ(t)\) is the long-term mean or equilibrium level of the interest rate.
• \(a\) is the speed of mean reversion.
• \(σ(t)\) is the time-dependent volatility of the interest rate.
• \(dW(t)\) is a Wiener process or Brownian motion, representing a random shock.

using DataInterpolations, Random, CairoMakie

function hull_white_theta_timevar_sigma(years, prices, σt; a::Float64)
    @assert length(years) == length(prices) == length(σt) "times and prices and sigmas must have same length"

    # Interpolate logP for smooth derivatives (cubic spline)
    logP = log.(prices)
    interp_logP = CubicSpline(logP, years)

    # instantaneous forward f(0,t) = -∂_t log P(0,t)
    f0(t) = -DataInterpolations.derivative(interp_logP, t)
    df0(t) = -DataInterpolations.derivative(interp_logP, t, 2)

    # Precompute Δs_j (left Riemann widths). assumes years[1] == 0 or handles first width as years[1].
    Δs = Vector{Float64}(undef, length(years))
    Δs[1] = years[1] # if years[1]==0 this is zero; else left interval from 0
    for j in 2:length(years)
        Δs[j] = years[j] - years[j-1]
    end

    # Compute integral term I(t_i) = 1/2 * ∫_0^{t_i} e^{-2a(t_i - s)} σ(s)^2 ds
    I = zeros(length(years))
    for i in 1:length(years)
        ti = years[i]
        s = 0.0
        acc = 0.0
        # left Riemann: use sigma at s = years[j] for interval [years[j], years[j+1])
        for j in 1:i
            sj = years[j]
            # width Δs[j] approximates length of [sj, sj+1) (or [0, years[1]] for j=1)
            w = Δs[j]
            if w <= 0
                # if grid is pure points (years[1]==0 and Δs[1]==0), skip zero-width
                continue
            end
            # integrand evaluated at left point sj
            integrand = exp(-2a * (ti - sj)) * (σt[j]^2)
            acc += integrand * w
        end
        I[i] = 0.5 * acc
    end

    # theta vector: df0 + a*f0 + I
    θt = [df0(years[i]) + a * f0(years[i]) + I[i] for i in 1:length(years)]

    return θt, f0, df0
end

let
    # Seed for reproducibility
    rng = MersenneTwister(1234)

    # Market curve setup (synthetic)
    years = collect(0.0:0.1:10.0)
    yields = 0.02 .+ 0.002 .* years
    prices = exp.(-yields .* years)

    # Parameters
    a = 0.1 # mean reversion speed
    σt = 0.008 .+ 0.004 .* sin.(0.5 .* years) # sinusoidal time-varying sigma
    θt, f0, df0 = hull_white_theta_timevar_sigma(years, prices, σt; a=a)
    r0 = 0.03 # initial short rate
    nsim = 1000 # number of simulated paths

    # Preallocate
    rpaths = zeros(length(years), nsim)

    # Euler–Maruyama simulation
    for path in 1:nsim
        r = r0
        rpaths[1, path] = r
        for i in 2:length(years)
            dt = years[i] - years[i - 1]
            dW = sqrt(dt) * randn(rng)
            drift = (θt[i - 1] - a * r) * dt
            diffusion = σt[i - 1] * dW
            r += drift + diffusion
            rpaths[i, path] = r
        end
    end

    # Plot
    fig = Figure(size=(800, 500))
    ax = Axis(fig[1, 1], xlabel="Time (years)", ylabel="Short rate", title="Hull–White Simulation with time-dependent θ(t) and σ(t)")
    for path in 1:nsim
        lines!(ax, years, rpaths[:, path])
    end
    fig
end

25.3.2 Equity Models

Equity scenario blocks usually mix a diffusion for price levels with a separate volatility process (to match option smiles) or with fat-tailed jump structures. We start with GBM for its closed-form intuition and then add a simple GARCH(1,1) volatility process.

25.3.2.1 Geometric Brownian Motion (GBM)

GBM is a stochastic process commonly used to model the price movement of financial instruments, including stocks. It assumes constant volatility and is characterized by a log-normal distribution. It is defined as

\[ dS(t) = \mu S(t) \, dt + \sigma S(t) \, dW(t) \]

where

• \(S(t)\) is the stock price at time \(t\).
• \(μ\) is the drift coefficient (expected return).
• \(σ\) is the volatility coefficient.
• \(dW(t)\) is a Wiener process or Brownian motion, representing a random shock.

using Random, CairoMakie

# Function to simulate GBM
function gbm_simulation(rng::AbstractRNG, μ, σ, S₀, Δt, num_steps, num_simulations)
    stock_price_paths = zeros(num_steps, num_simulations)
    drift = (μ - 0.5 * σ * σ) * Δt
    vol = σ * sqrt(Δt)
    for j in 1:num_simulations
        stock_price_paths[1, j] = S₀
        shocks = randn(rng, num_steps - 1)
        for i in 2:num_steps
            stock_price_paths[i, j] = stock_price_paths[i-1, j] * exp(drift + vol * shocks[i - 1])
        end
    end
    return stock_price_paths
end

let
    # Set seed for reproducibility
    rng = MersenneTwister(1234)

    # GBM parameters
    μ = 0.05 # Drift (expected return)
    σ = 0.2 # Volatility

    # Initial stock price
    S₀ = 100.0

    # Number of time steps and simulations
    num_steps = 252
    num_simulations = 1_000

    # Time increment
    Δt = 1 / 252

    # Run GBM simulation
    gbm_paths = gbm_simulation(rng, μ, σ, S₀, Δt, num_steps, num_simulations)

    # Plot the simulated stock price paths
    f = Figure()
    ax = Axis(f[1, 1], xlabel="Year", ylabel="Price", title="GBM price scenarios")
    for i in 1:num_simulations
        lines!(ax, 1:num_steps, gbm_paths[:, i])
    end
    f
end

25.3.2.2 Generalized Autoregressive Conditional Heteroskedasticity (GARCH)

GARCH models capture time-varying volatility. They are often used in conjunction with other models to forecast volatility. It is defined as

\[ \sigma^2_t = \omega + \alpha_1 r^2_{t-1} + \beta_1 \sigma^2_{t-1} \]

\[ r_t = \varepsilon_t \sqrt{\sigma^2_t} \]

• \(σ^2_t\) is the conditional variance at time \(t\)
• \(r_t\) is the return at time \(t\)
• \(\varepsilon_t\) is a white noise or innovation process
• \(\omega\), \(\alpha_1\), \(\beta_1\) are model parameters

using Random, CairoMakie

# Function to simulate GARCH(1,1) volatility
function garch_simulation(rng::AbstractRNG, ω, α₁, β₁, num_steps, num_sims)
    volatility_paths = zeros(num_steps, num_sims)
    σ²_ss = ω / (1 - α₁ - β₁)
    returns = zeros(Float64, num_steps, num_sims)
    for j in 1:num_sims
        returns[1, j] = volatility_paths[1, j] * randn(rng)
        for i in 2:num_steps
            σ² = ω + α₁ * returns[i - 1, j] * returns[i - 1, j] + β₁ * volatility_paths[i - 1, j] * volatility_paths[i - 1, j]
            volatility_paths[i, j] = sqrt(σ²)
            returns[i, j] = volatility_paths[i, j] * randn(rng)
        end
    end
    return volatility_paths
end

let
    # Set seed for reproducibility
    rng = MersenneTwister(1234)

    # GARCH(1,1) parameters
    ω = 0.01 # Constant term
    α₁ = 0.1 # Coefficient for lagged squared returns
    β₁ = 0.8 # Coefficient for lagged conditional volatility
    @assert α₁ + β₁ < 1 "Stationarity requires α₁ + β₁ < 1"

    # Number of time steps and simulations
    num_steps = 252
    num_simulations = 1_000

    # Time increment
    Δt = 1 / 252

    # Run GARCH simulation
    garch_paths = garch_simulation(rng, ω, α₁, β₁, num_steps, num_simulations)

    # Plot the simulated volatility paths
    f = Figure()
    ax = Axis(f[1, 1], xlabel="Trading days", ylabel="Volatility", title="GARCH(1, 1) conditional volatility")
    for i in 1:num_simulations
        lines!(ax, 1:num_steps, garch_paths[:, i])
    end
    f
end

25.3.3 Copulas

Copulas give us a flexible way to model dependence structures between random variables, independently from their marginal distributions. With copulas, we can model each variable’s distribution (marginal) separately, and then “glue” them together with a copula to describe the dependence structure. This is useful when variables have very different distributions (e.g. stock returns vs. interest rates). Besides, copulas can capture nonlinear dependence and tail dependence (how extreme events co-occur), v.s. traditional correlation (e.g. Pearson) only captures linear dependence.

Simulating data using copulas involves generating multivariate samples with specified marginal distributions and a copula structure.

using Random, Copulas, CairoMakie

let
    Random.seed!(1234)

    # 2D Gaussian copula with correlation ρ = 0.8
    ρ = 0.8
    Σ = [1.0 ρ;
        ρ 1.0]
    Cg = GaussianCopula(Σ)
    U = rand(Cg, 10_000)  # n × d matrix with U(0,1) margins
    f = Figure()
    ax = Axis(f[1, 1], xlabel="U₁", ylabel="U₂")
    scatter!(ax, U[1, :], U[2, :], markersize=2, color=:steelblue, alpha=0.4)
    f
end

Copulas can also be used to infer combined distributions from data samples.

using Random, Copulas, Distributions, CairoMakie

let
    Random.seed!(1234)
    X₁ = Gamma(2.0, 3.0)
    X₂ = Pareto(3.0, 1.0)
    X₃ = LogNormal(0.0, 1.0)
    C = ClaytonCopula(3, 0.7) # A 3-variate Clayton Copula with θ = 0.7
    D = SklarDist(C, (X₁, X₂, X₃)) # The final distribution

    # Generate a dataset
    simu = rand(D, 1000)
    # Estimate marginals by MLE
    d̂1 = fit(Gamma, simu[:, 1])
    d̂2 = fit(Pareto, simu[:, 2])
    d̂3 = fit(LogNormal, simu[:, 3])
    # Probability integral transform to uniforms
    Û = hcat(cdf.(d̂1, simu[:, 1]), cdf.(d̂2, simu[:, 2]), cdf.(d̂3, simu[:, 3]))
    # Fit copula on uniforms
    Ĉ = fit(ClaytonCopula, Û)
    # Compose the estimated Sklar distribution
    D̂ = SklarDist(Ĉ, (d̂1, d̂2, d̂3))
    
    # Generate new samples from the fitted distribution
    simu_fitted = rand(D̂, 1000)
    
    # Create visualizations
    # fig = Figure(size=(1000, 800))
    
    # # Plot pairwise scatter plots
    # ax1 = Axis(fig[1, 1], xlabel="X₁ (Gamma)", ylabel="X₂ (Pareto)", title="Original Data")
    # scatter!(ax1, simu[:, 1], simu[:, 2], markersize=3, color=:steelblue, alpha=0.5)
    
    # ax2 = Axis(fig[1, 2], xlabel="X₁ (Gamma)", ylabel="X₂ (Pareto)", title="Fitted Distribution")
    # scatter!(ax2, simu_fitted[:, 1], simu_fitted[:, 2], markersize=3, color=:coral, alpha=0.5)
    
    # ax3 = Axis(fig[2, 1], xlabel="X₁ (Gamma)", ylabel="X₃ (LogNormal)")
    # scatter!(ax3, simu[:, 1], simu[:, 3], markersize=3, color=:steelblue, alpha=0.5)
    
    # ax4 = Axis(fig[2, 2], xlabel="X₁ (Gamma)", ylabel="X₃ (LogNormal)")
    # scatter!(ax4, simu_fitted[:, 1], simu_fitted[:, 3], markersize=3, color=:coral, alpha=0.5)
    
    # ax5 = Axis(fig[3, 1], xlabel="X₂ (Pareto)", ylabel="X₃ (LogNormal)")
    # scatter!(ax5, simu[:, 2], simu[:, 3], markersize=3, color=:steelblue, alpha=0.5)
    
    # ax6 = Axis(fig[3, 2], xlabel="X₂ (Pareto)", ylabel="X₃ (LogNormal)")
    # scatter!(ax6, simu_fitted[:, 2], simu_fitted[:, 3], markersize=3, color=:coral, alpha=0.5)
    
    # fig
end

3×1000 Matrix{Float64}:
 8.0567   8.46045   0.926878  0.660617  …  2.78892   4.28651   2.45952
 1.84058  1.7161    2.17833   1.63618      2.51906   2.02895   3.71619
 0.39327  0.217874  5.94546   4.49684      0.858615  0.768066  0.00817772

25.4 Where to Go Next

Realistic enterprise scenario engines incorporate macro regime switches, stochastic volatility surfaces, credit migration, and balance-sheet feedback loops. Julia’s multiple dispatch and composability make it straightforward to encapsulate each component (e.g., AbstractInterestRateProcess, EquityBlock, CopulaLink) and then assemble them into reproducible scenario trees or Monte Carlo streams. In subsequent chapters we will look at calibrating these models and feeding the resulting scenarios into portfolio analytics and ALM dashboards.