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 use RN models to 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 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 (for the continuous-time process, strictly so when the Feller condition \(2\kappa\theta \geq \sigma^2\) holds and \(r(0) > 0\)).

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, discretized with a truncated Euler–Maruyama scheme. The specification of \(dr\) can be changed to make it a Vasicek model.

using Random, Statistics, 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 = κ*(θ - r)*Δt + σ*dW for Vasicek
            r_prev = interest_rate_paths[i-1, j]
            drift = κ * (θ - r_prev) * Δt
            diffusion = σ * sqrt(max(r_prev, 0.0)) * st * shocks[i-1]
            dr = drift + diffusion
            # Ensure non-negativity
            interest_rate_paths[i, j] = max(r_prev + dr, 0.0)
        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)

    # Compute percentiles across simulations at each time step
    ti = 0:(num_steps-1)
    p05 = [quantile(cir_paths[i, :], 0.05) for i in 1:num_steps]
    p25 = [quantile(cir_paths[i, :], 0.25) for i in 1:num_steps]
    p50 = [quantile(cir_paths[i, :], 0.50) for i in 1:num_steps]
    p75 = [quantile(cir_paths[i, :], 0.75) for i in 1:num_steps]
    p95 = [quantile(cir_paths[i, :], 0.95) for i in 1:num_steps]

    # Plot shaded bands
    f = Figure()
    ax = Axis(f[1, 1], xlabel="Trading days",
        ylabel="Short rate", title="CIR paths")
    band!(ax, ti, p05, p95, color=(:steelblue, 0.2), label="5th–95th percentile")
    band!(ax, ti, p25, p75, color=(:steelblue, 0.4), label="25th–75th percentile")
    lines!(ax, ti, p50, color=:steelblue, linewidth=2, label="Median")
    axislegend(ax, position=:rb)
    f
end

Precompiling packages...
   3682.6 ms  ✓ QuartoNotebookWorkerMakieExt (serial)
  1 dependency successfully precompiled in 4 seconds
Precompiling packages...
   3199.9 ms  ✓ QuartoNotebookWorkerCairoMakieExt (serial)
  1 dependency successfully precompiled in 3 seconds

Note

The positivity guarantee above belongs to the continuous-time process; the discretized simulation is another matter. Over a finite time step the Gaussian shock can push the rate negative even when the Feller condition holds, which is why the code clamps at zero (max(r_prev, 0.0) inside the square root, and truncation of the updated rate). This truncation is convenient and common in practice, but it means the code is an approximation of CIR: the clamping slightly biases the simulated distribution near zero, and the discretization error shrinks only as the time step is refined. The exact simulation scheme below avoids both issues.

The CIR process can also be simulated exactly, because conditional on \(r_t\) the distribution of \(r_{t+\Delta t}\) is a scaled non-central chi-squared distribution.

\[ \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, Statistics, 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)

    # Compute percentiles across simulations at each time step
    nsteps_total = length(tgrid)
    p05 = [quantile(rpaths[i, :], 0.05) for i in 1:nsteps_total]
    p25 = [quantile(rpaths[i, :], 0.25) for i in 1:nsteps_total]
    p50 = [quantile(rpaths[i, :], 0.50) for i in 1:nsteps_total]
    p75 = [quantile(rpaths[i, :], 0.75) for i in 1:nsteps_total]
    p95 = [quantile(rpaths[i, :], 0.95) for i in 1:nsteps_total]

    # Plot shaded bands
    fig = Figure(size=(800, 500))
    ax = Axis(fig[1, 1], xlabel="Time (years)",
        ylabel="Short rate", title="Exact CIR Simulation")
    band!(ax, collect(tgrid), p05, p95, color=(:steelblue, 0.2), label="5th–95th percentile")
    band!(ax, collect(tgrid), p25, p75, color=(:steelblue, 0.4), label="25th–75th percentile")
    lines!(ax, collect(tgrid), p50, color=:steelblue, linewidth=2, label="Median")
    axislegend(ax, position=:rb)
    fig
end

25.3.1.2 Hull-White

The Hull-White model is a one-factor short-rate model that extends the Vasicek model by allowing the mean reversion level and volatility to be time-dependent. It is commonly used for pricing interest rate derivatives. A separate class of models, the Libor Market Models (LMMs) such as the Brace-Gatarek-Musiela (BGM) model, instead model the evolution of forward rates directly rather than the short rate. The Hull-White model 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 calibrated market-implied 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, Statistics, CairoMakie

function hull_white_theta_timevar_sigma(years, prices, σt; a::Float64)
    n = length(years)
    @assert n == length(prices) == length(σt) "mismatched lengths"

    # 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) = ∫_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] = 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
    # sinusoidal time-varying sigma
    σt = 0.008 .+ 0.004 .* sin.(0.5 .* years)
    θ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

    # Compute percentiles across simulations at each time step
    nsteps_total = length(years)
    p05 = [quantile(rpaths[i, :], 0.05) for i in 1:nsteps_total]
    p25 = [quantile(rpaths[i, :], 0.25) for i in 1:nsteps_total]
    p50 = [quantile(rpaths[i, :], 0.50) for i in 1:nsteps_total]
    p75 = [quantile(rpaths[i, :], 0.75) for i in 1:nsteps_total]
    p95 = [quantile(rpaths[i, :], 0.95) for i in 1:nsteps_total]

    # Plot shaded bands
    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)")
    band!(ax, years, p05, p95, color=(:steelblue, 0.2), label="5th–95th percentile")
    band!(ax, years, p25, p75, color=(:steelblue, 0.4), label="25th–75th percentile")
    lines!(ax, years, p50, color=:steelblue, linewidth=2, label="Median")
    axislegend(ax, position=:rb)
    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)

Geometric Brownian Motion (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, Statistics, 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
            S_prev = stock_price_paths[i-1, j]
            stock_price_paths[i, j] = S_prev * 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)

    # Compute percentiles across simulations at each time step
    ti = 1:num_steps
    p05 = [quantile(gbm_paths[i, :], 0.05) for i in 1:num_steps]
    p25 = [quantile(gbm_paths[i, :], 0.25) for i in 1:num_steps]
    p50 = [quantile(gbm_paths[i, :], 0.50) for i in 1:num_steps]
    p75 = [quantile(gbm_paths[i, :], 0.75) for i in 1:num_steps]
    p95 = [quantile(gbm_paths[i, :], 0.95) for i in 1:num_steps]

    # Plot shaded bands
    f = Figure()
    ax = Axis(f[1, 1], xlabel="Trading days", ylabel="Price", title="GBM price scenarios")
    band!(ax, collect(ti), p05, p95, color=(:steelblue, 0.2), label="5th–95th percentile")
    band!(ax, collect(ti), p25, p75, color=(:steelblue, 0.4), label="25th–75th percentile")
    lines!(ax, collect(ti), p50, color=:steelblue, linewidth=2, label="Median")
    axislegend(ax, position=:lt)
    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. The GARCH(1,1) model 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, Statistics, 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
        volatility_paths[1, j] = sqrt(σ²_ss)
        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)

    # Compute percentiles across simulations at each time step
    ti = 1:num_steps
    p05 = [quantile(garch_paths[i, :], 0.05) for i in 1:num_steps]
    p25 = [quantile(garch_paths[i, :], 0.25) for i in 1:num_steps]
    p50 = [quantile(garch_paths[i, :], 0.50) for i in 1:num_steps]
    p75 = [quantile(garch_paths[i, :], 0.75) for i in 1:num_steps]
    p95 = [quantile(garch_paths[i, :], 0.95) for i in 1:num_steps]

    # Plot shaded bands
    f = Figure()
    ax = Axis(
        f[1, 1],
        xlabel="Trading days", ylabel="Volatility",
        title="GARCH(1, 1) conditional volatility"
    )
    band!(ax, collect(ti), p05, p95, color=(:steelblue, 0.2), label="5th–95th percentile")
    band!(ax, collect(ti), p25, p75, color=(:steelblue, 0.4), label="25th–75th percentile")
    lines!(ax, collect(ti), p50, color=:steelblue, linewidth=2, label="Median")
    axislegend(ax, position=:rt)
    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). Moreover, copulas can capture nonlinear dependence and tail dependence (how extreme events co-occur), whereas 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)  # d × n 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
    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 (returns dimensions × observations)
    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
    Û = permutedims(
        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

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.