API Reference

This section provides auto-generated documentation directly extracted from the Aethel codebase.

---

Simulator Configuration

aethel.engine.parameters.SimulatorConfig dataclass

Configuration parameters for the Economic Scenario Generator. Encapsulates all mathematical and structural constants.

Source code in src/aethel/engine/parameters.py

@dataclass
class SimulatorConfig:
    """
    Configuration parameters for the Economic Scenario Generator.
    Encapsulates all mathematical and structural constants.
    """
    duration_years: int = 60
    num_scenarios: int = 1000
    seed: int = 42

    # Chunking and thread worker optional overrides
    chunk_size: Optional[int] = None
    max_workers: Optional[int] = None

    # Structural yield/tenor constants
    tenors: np.ndarray = field(
        default_factory=lambda: np.array([0.25, 0.5, 1.0, 2.0, 5.0, 10.0, 20.0, 30.0])
    )

    # Policy and inflation boundary thresholds
    mu_min: float = 0.010
    pi_min: float = -0.02

    # Model dynamics parameters
    alpha_smooth: float = 0.20
    beta_drag: float = 0.15
    eta_erp: float = 0.10
    lambda_irp: float = 0.05
    kappa_irp: float = 0.15

    # Merton Jump-Diffusion parameters
    lambda_J: float = 0.20
    mu_J: float = -0.15
    sigma_J: float = 0.10

    # Initial state values (Standard international economic terms)
    initial_rate: float = 0.10
    initial_inflation: float = 0.045

    # Calibrated parameter overrides (filled dynamically by the Calibrator)
    ou_mu: Optional[float] = None
    cir_theta_val: Optional[float] = None
    cir_sigma_val: Optional[float] = None
    cir_mu_val: Optional[float] = None
    ou_theta_val: Optional[float] = None
    ou_sigma_val: Optional[float] = None
    gbm_sigma_val: Optional[float] = None

    @property
    def steps(self) -> int:
        return self.duration_years * 12

    @property
    def dt(self) -> float:
        return 1.0 / 12.0

    @property
    def k_jump(self) -> float:
        return np.exp(self.mu_J + 0.5 * (self.sigma_J ** 2)) - 1.0

    def to_dict(self) -> dict:
        """
        Converts the configuration instance into a dictionary.
        Automatically converts NumPy arrays to standard lists for JSON serializability.
        """
        import dataclasses
        d = dataclasses.asdict(self)
        for key, value in d.items():
            if isinstance(value, np.ndarray):
                d[key] = value.tolist()
        return d

    @classmethod
    def from_dict(cls, d: dict) -> "SimulatorConfig":
        """
        Constructs a SimulatorConfig instance from a dictionary.
        Handles both static constructor fields and dynamic attributes added by calibration.
        """
        valid_fields = {f.name for f in fields(cls)}

        init_kwargs = {}
        dynamic_attrs = {}
        for k, v in d.items():
            if k in valid_fields:
                init_kwargs[k] = v
            else:
                dynamic_attrs[k] = v

        if "tenors" in init_kwargs and isinstance(init_kwargs["tenors"], list):
            init_kwargs["tenors"] = np.array(init_kwargs["tenors"])

        config = cls(**init_kwargs)

        for k, v in dynamic_attrs.items():
            setattr(config, k, v)

        return config

from_dict(d) classmethod

Constructs a SimulatorConfig instance from a dictionary. Handles both static constructor fields and dynamic attributes added by calibration.

Source code in src/aethel/engine/parameters.py

@classmethod
def from_dict(cls, d: dict) -> "SimulatorConfig":
    """
    Constructs a SimulatorConfig instance from a dictionary.
    Handles both static constructor fields and dynamic attributes added by calibration.
    """
    valid_fields = {f.name for f in fields(cls)}

    init_kwargs = {}
    dynamic_attrs = {}
    for k, v in d.items():
        if k in valid_fields:
            init_kwargs[k] = v
        else:
            dynamic_attrs[k] = v

    if "tenors" in init_kwargs and isinstance(init_kwargs["tenors"], list):
        init_kwargs["tenors"] = np.array(init_kwargs["tenors"])

    config = cls(**init_kwargs)

    for k, v in dynamic_attrs.items():
        setattr(config, k, v)

    return config

to_dict()

Converts the configuration instance into a dictionary. Automatically converts NumPy arrays to standard lists for JSON serializability.

Source code in src/aethel/engine/parameters.py

def to_dict(self) -> dict:
    """
    Converts the configuration instance into a dictionary.
    Automatically converts NumPy arrays to standard lists for JSON serializability.
    """
    import dataclasses
    d = dataclasses.asdict(self)
    for key, value in d.items():
        if isinstance(value, np.ndarray):
            d[key] = value.tolist()
    return d

---

Simulation Engine

aethel.engine.simulator.MarketSimulator

An Economic Scenario Generator (ESG) that orchestrates parameters, concurrent random path generation, and analytical yield derivation.

Source code in src/aethel/engine/simulator.py

class MarketSimulator:
    """
    An Economic Scenario Generator (ESG) that orchestrates parameters,
    concurrent random path generation, and analytical yield derivation.
    """
    def __init__(self, config: Optional[SimulatorConfig] = None):
        self.config = config if config is not None else SimulatorConfig()

    def run(self) -> LazyScenarioList:
        """
        Executes the simulation engine, dynamically choosing between single-block
        or memory-safe chunked execution based on available host RAM.
        """
        cfg = self.config

        # 1. Deterministic generation of scenario-level parameters
        master_rng = np.random.default_rng(cfg.seed)

        ou_mu_mean = cfg.ou_mu if cfg.ou_mu is not None else 0.055
        ou_mu = np.clip(master_rng.normal(ou_mu_mean, 0.008, cfg.num_scenarios), ou_mu_mean - 0.02, ou_mu_mean + 0.02)
        pi_target = ou_mu

        r_real_mean = cfg.cir_mu_val - ou_mu_mean if (cfg.cir_mu_val is not None) else 0.050
        r_real_mean = np.clip(r_real_mean, 0.02, 0.08)
        r_real_target = np.clip(master_rng.normal(r_real_mean, 0.008, cfg.num_scenarios), r_real_mean - 0.02, r_real_mean + 0.02)

        gamma = master_rng.uniform(0.3, 0.7, cfg.num_scenarios)
        base_erp = master_rng.uniform(0.005, 0.025, cfg.num_scenarios)

        cir_sigma_mean = cfg.cir_sigma_val if cfg.cir_sigma_val is not None else 0.08
        cir_sigma = np.clip(master_rng.normal(cir_sigma_mean, 0.01, cfg.num_scenarios), cir_sigma_mean - 0.02, cir_sigma_mean + 0.02)

        ou_sigma_mean = cfg.ou_sigma_val if cfg.ou_sigma_val is not None else 0.015
        ou_sigma = np.clip(master_rng.normal(ou_sigma_mean, 0.003, cfg.num_scenarios), ou_sigma_mean - 0.005, ou_sigma_mean + 0.007)

        gbm_sigma_mean = cfg.gbm_sigma_val if cfg.gbm_sigma_val is not None else 0.18
        gbm_sigma = np.clip(master_rng.normal(gbm_sigma_mean, 0.02, cfg.num_scenarios), gbm_sigma_mean - 0.04, gbm_sigma_mean + 0.04)

        cir_theta_val = cfg.cir_theta_val if cfg.cir_theta_val is not None else 0.25
        cir_theta = np.full(cfg.num_scenarios, cir_theta_val)

        ou_theta_val = cfg.ou_theta_val if cfg.ou_theta_val is not None else 0.35
        ou_theta = np.full(cfg.num_scenarios, ou_theta_val)

        # 2. Derive global random variables correlation matrix
        rho_12 = np.clip(master_rng.normal(0.40, 0.05), 0.25, 0.55)
        rho_13 = np.clip(master_rng.normal(-0.15, 0.05), -0.25, -0.05)
        rho_23 = np.clip(master_rng.normal(-0.10, 0.05), -0.20, 0.00)

        correlation_matrix = np.array([
            [1.00, rho_12, rho_13],
            [rho_12, 1.00, rho_23],
            [rho_13, rho_23, 1.00]
        ])
        L = np.linalg.cholesky(correlation_matrix)

        # 3. Dynamic Hardware-Aware Concurrency & Safety Buffer
        avail_ram = get_available_system_ram()

        os_safety_buffer = 4 * 1024 * 1024 * 1024
        usable_ram = max(1 * 1024 * 1024 * 1024, avail_ram - os_safety_buffer)

        peak_mem_needed = 136 * cfg.steps * cfg.num_scenarios

        if cfg.max_workers is not None:
            max_workers = cfg.max_workers
        else:
            system_cpus = os.cpu_count() or 1
            if usable_ram < 8 * 1024 * 1024 * 1024:
                max_workers = min(2, system_cpus)
            else:
                max_workers = min(4, system_cpus)

        # 4. Defensive Triggering Rules
        if cfg.chunk_size is not None:
            chunk_size = cfg.chunk_size
            use_chunking = True
            print(f"[ESG] Manual override: Chunked execution activated (chunk size: {chunk_size}).")
        else:
            use_chunking = (cfg.num_scenarios > 20000) or (peak_mem_needed > (0.25 * usable_ram))
            if use_chunking:
                temp_mem_per_scenario = 64 * cfg.steps
                total_temp_limit = 256 * 1024 * 1024
                temp_limit_per_worker = total_temp_limit / max_workers

                chunk_size = int(temp_limit_per_worker / temp_mem_per_scenario)
                chunk_size = min(5000, max(100, chunk_size))

                print(f"[ESG] High memory projection detected ({peak_mem_needed / (1024**3):.2f} GB estimated, {avail_ram / (1024**3):.2f} GB system available).")
                print(f"[ESG] Defensive dynamic chunking active (chunk size: {chunk_size}, active workers: {max_workers}).")
            else:
                chunk_size = cfg.num_scenarios

        num_chunks = int(np.ceil(cfg.num_scenarios / chunk_size))

        # 5. Symmetrical Master Outputs Pre-Allocation (In-Memory)
        rate_paths = np.zeros((cfg.steps + 1, cfg.num_scenarios), dtype=np.float64)
        inflation_paths = np.zeros((cfg.steps + 1, cfg.num_scenarios), dtype=np.float64)
        y_target_paths = np.zeros((cfg.steps + 1, cfg.num_scenarios), dtype=np.float64)
        mu_rate_paths = np.zeros((cfg.steps + 1, cfg.num_scenarios), dtype=np.float64)

        equity_returns = np.zeros((cfg.steps, cfg.num_scenarios), dtype=np.float64)
        cpis = np.zeros((cfg.steps, cfg.num_scenarios), dtype=np.float64)
        deposit_rates = np.zeros((cfg.steps, cfg.num_scenarios), dtype=np.float64)

        # 6. Dispatch threads filling master array views directly
        futures = []
        with ThreadPoolExecutor(max_workers=max_workers) as executor:
            for j in range(num_chunks):
                start_idx = j * chunk_size
                end_idx = min(start_idx + chunk_size, cfg.num_scenarios)

                futures.append(
                    executor.submit(
                        self._run_chunk,
                        start_idx=start_idx,
                        end_idx=end_idx,
                        seed_offset=cfg.seed + start_idx,
                        ou_mu_chunk=ou_mu[start_idx:end_idx],
                        r_real_target_chunk=r_real_target[start_idx:end_idx],
                        gamma_chunk=gamma[start_idx:end_idx],
                        base_erp_chunk=base_erp[start_idx:end_idx],
                        cir_sigma_chunk=cir_sigma[start_idx:end_idx],
                        ou_sigma_chunk=ou_sigma[start_idx:end_idx],
                        gbm_sigma_chunk=gbm_sigma[start_idx:end_idx],
                        cir_theta_chunk=cir_theta[start_idx:end_idx],
                        ou_theta_chunk=ou_theta[start_idx:end_idx],
                        L=L,
                        equity_returns_view=equity_returns[:, start_idx:end_idx],
                        cpis_view=cpis[:, start_idx:end_idx],
                        deposit_rates_view=deposit_rates[:, start_idx:end_idx],
                        rate_paths_view=rate_paths[:, start_idx:end_idx],
                        inflation_paths_view=inflation_paths[:, start_idx:end_idx],
                        y_target_paths_view=y_target_paths[:, start_idx:end_idx],
                        mu_rate_paths_view=mu_rate_paths[:, start_idx:end_idx]
                    )
                )

            for fut in futures:
                fut.result()

        print("[ESG] Core processing phase complete.")
        return LazyScenarioList(
            equity_returns=np.ascontiguousarray(equity_returns.T),
            cpis=np.ascontiguousarray(cpis.T),
            deposit_rates=np.ascontiguousarray(deposit_rates.T),
            rate_paths=rate_paths,
            mu_rate_paths=mu_rate_paths,
            inflation_paths=inflation_paths,
            y_target_paths=y_target_paths,
            cir_theta=cir_theta,
            cir_sigma=cir_sigma,
            ou_theta=ou_theta,
            ou_sigma=ou_sigma,
            tenors=cfg.tenors,
            pi_min=cfg.pi_min,
            lambda_irp=cfg.lambda_irp,
            kappa_irp=cfg.kappa_irp
        )

    def _run_chunk(
        self,
        start_idx: int,
        end_idx: int,
        seed_offset: int,
        ou_mu_chunk: np.ndarray,
        r_real_target_chunk: np.ndarray,
        gamma_chunk: np.ndarray,
        base_erp_chunk: np.ndarray,
        cir_sigma_chunk: np.ndarray,
        ou_sigma_chunk: np.ndarray,
        gbm_sigma_chunk: np.ndarray,
        cir_theta_chunk: np.ndarray,
        ou_theta_chunk: np.ndarray,
        L: np.ndarray,
        equity_returns_view: np.ndarray,
        cpis_view: np.ndarray,
        deposit_rates_view: np.ndarray,
        rate_paths_view: np.ndarray,
        inflation_paths_view: np.ndarray,
        y_target_paths_view: np.ndarray,
        mu_rate_paths_view: np.ndarray
    ):
        """Worker thread compiling and running an isolated block of scenarios."""
        cfg = self.config
        chunk_size = end_idx - start_idx
        sqrt_dt = np.sqrt(cfg.dt)

        Z_raw_flat = np.empty((3, cfg.steps * chunk_size))
        jump_shocks_all = np.zeros((cfg.steps, chunk_size), dtype=np.float64)

        for i in range(chunk_size):
            scenario_seed = seed_offset + i
            local_rng = np.random.default_rng(scenario_seed)

            Z_raw_flat[:, i * cfg.steps : (i + 1) * cfg.steps] = local_rng.normal(0.0, 1.0, (3, cfg.steps))

            jumps = local_rng.poisson(cfg.lambda_J * cfg.dt, cfg.steps)
            active_idx = np.where(jumps > 0)[0]
            for t_step in active_idx:
                num_jumps = jumps[t_step]
                jump_shocks_all[t_step, i] = np.sum(local_rng.normal(cfg.mu_J, cfg.sigma_J, num_jumps))

        Z_corr = (L @ Z_raw_flat).reshape(3, chunk_size, cfg.steps).transpose(2, 0, 1)
        Z_corr = np.ascontiguousarray(Z_corr)

        y_paths = np.empty((cfg.steps + 1, chunk_size))
        smoothed_inflation_paths = np.empty((cfg.steps + 1, chunk_size))

        # Direct writes initialization into master view pointers
        rate_paths_view[0] = cfg.initial_rate
        inflation_paths_view[0] = cfg.initial_inflation
        y_paths[0] = inflation_paths_view[0] - cfg.pi_min
        smoothed_inflation_paths[0] = inflation_paths_view[0]

        cir_theta_dt = cir_theta_chunk * cfg.dt
        cir_sigma_sqrt_dt = cir_sigma_chunk * sqrt_dt
        ou_theta_dt = ou_theta_chunk * cfg.dt
        ou_sigma_sqrt_dt = ou_sigma_chunk * sqrt_dt

        gbm_sigma_sq_dt = 0.5 * (gbm_sigma_chunk ** 2) * cfg.dt
        gbm_sigma_sqrt_dt = gbm_sigma_chunk * sqrt_dt
        drift_adjustment = - (cfg.lambda_J * cfg.k_jump) * cfg.dt - gbm_sigma_sq_dt

        if HAS_NUMBA:
            run_simulation_loop_numba(
                cfg.steps, chunk_size, cfg.dt,
                Z_corr, jump_shocks_all,
                ou_mu_chunk, ou_mu_chunk, r_real_target_chunk, gamma_chunk, base_erp_chunk,
                cir_theta_dt, cir_sigma_sqrt_dt, ou_theta_dt, ou_sigma_sqrt_dt,
                gbm_sigma_sqrt_dt, drift_adjustment,
                rate_paths_view, inflation_paths_view, y_paths, smoothed_inflation_paths,
                mu_rate_paths_view, y_target_paths_view, equity_returns_view,
                cfg.mu_min, cfg.pi_min, cfg.beta_drag, cfg.alpha_smooth, cfg.eta_erp
            )
        else:
            run_simulation_loop_numpy(
                cfg.steps, chunk_size, cfg.dt,
                Z_corr, jump_shocks_all,
                ou_mu_chunk, ou_mu_chunk, r_real_target_chunk, gamma_chunk, base_erp_chunk,
                cir_theta_dt, cir_sigma_sqrt_dt, ou_theta_dt, ou_sigma_sqrt_dt,
                gbm_sigma_sqrt_dt, drift_adjustment,
                rate_paths_view, inflation_paths_view, y_paths, smoothed_inflation_paths,
                mu_rate_paths_view, y_target_paths_view, equity_returns_view,
                cfg.mu_min, cfg.pi_min, cfg.beta_drag, cfg.alpha_smooth, cfg.eta_erp
            )

        mu_rate_paths_view[cfg.steps] = mu_rate_paths_view[cfg.steps - 1]
        y_target_paths_view[cfg.steps] = y_target_paths_view[cfg.steps - 1]

        # Direct in-place writes for monthly outputs (Zero memory copies)
        deposit_rates_view[:] = (1.0 + np.maximum(-0.99, rate_paths_view[:-1, :])) ** (1.0 / 12.0) - 1.0
        inflation_monthly_all = (1.0 + np.maximum(-0.99, inflation_paths_view[:-1, :])) ** (1.0 / 12.0) - 1.0
        cpis_view[:] = np.cumprod(1.0 + inflation_monthly_all, axis=0)

run()

Executes the simulation engine, dynamically choosing between single-block or memory-safe chunked execution based on available host RAM.

Source code in src/aethel/engine/simulator.py

def run(self) -> LazyScenarioList:
    """
    Executes the simulation engine, dynamically choosing between single-block
    or memory-safe chunked execution based on available host RAM.
    """
    cfg = self.config

    # 1. Deterministic generation of scenario-level parameters
    master_rng = np.random.default_rng(cfg.seed)

    ou_mu_mean = cfg.ou_mu if cfg.ou_mu is not None else 0.055
    ou_mu = np.clip(master_rng.normal(ou_mu_mean, 0.008, cfg.num_scenarios), ou_mu_mean - 0.02, ou_mu_mean + 0.02)
    pi_target = ou_mu

    r_real_mean = cfg.cir_mu_val - ou_mu_mean if (cfg.cir_mu_val is not None) else 0.050
    r_real_mean = np.clip(r_real_mean, 0.02, 0.08)
    r_real_target = np.clip(master_rng.normal(r_real_mean, 0.008, cfg.num_scenarios), r_real_mean - 0.02, r_real_mean + 0.02)

    gamma = master_rng.uniform(0.3, 0.7, cfg.num_scenarios)
    base_erp = master_rng.uniform(0.005, 0.025, cfg.num_scenarios)

    cir_sigma_mean = cfg.cir_sigma_val if cfg.cir_sigma_val is not None else 0.08
    cir_sigma = np.clip(master_rng.normal(cir_sigma_mean, 0.01, cfg.num_scenarios), cir_sigma_mean - 0.02, cir_sigma_mean + 0.02)

    ou_sigma_mean = cfg.ou_sigma_val if cfg.ou_sigma_val is not None else 0.015
    ou_sigma = np.clip(master_rng.normal(ou_sigma_mean, 0.003, cfg.num_scenarios), ou_sigma_mean - 0.005, ou_sigma_mean + 0.007)

    gbm_sigma_mean = cfg.gbm_sigma_val if cfg.gbm_sigma_val is not None else 0.18
    gbm_sigma = np.clip(master_rng.normal(gbm_sigma_mean, 0.02, cfg.num_scenarios), gbm_sigma_mean - 0.04, gbm_sigma_mean + 0.04)

    cir_theta_val = cfg.cir_theta_val if cfg.cir_theta_val is not None else 0.25
    cir_theta = np.full(cfg.num_scenarios, cir_theta_val)

    ou_theta_val = cfg.ou_theta_val if cfg.ou_theta_val is not None else 0.35
    ou_theta = np.full(cfg.num_scenarios, ou_theta_val)

    # 2. Derive global random variables correlation matrix
    rho_12 = np.clip(master_rng.normal(0.40, 0.05), 0.25, 0.55)
    rho_13 = np.clip(master_rng.normal(-0.15, 0.05), -0.25, -0.05)
    rho_23 = np.clip(master_rng.normal(-0.10, 0.05), -0.20, 0.00)

    correlation_matrix = np.array([
        [1.00, rho_12, rho_13],
        [rho_12, 1.00, rho_23],
        [rho_13, rho_23, 1.00]
    ])
    L = np.linalg.cholesky(correlation_matrix)

    # 3. Dynamic Hardware-Aware Concurrency & Safety Buffer
    avail_ram = get_available_system_ram()

    os_safety_buffer = 4 * 1024 * 1024 * 1024
    usable_ram = max(1 * 1024 * 1024 * 1024, avail_ram - os_safety_buffer)

    peak_mem_needed = 136 * cfg.steps * cfg.num_scenarios

    if cfg.max_workers is not None:
        max_workers = cfg.max_workers
    else:
        system_cpus = os.cpu_count() or 1
        if usable_ram < 8 * 1024 * 1024 * 1024:
            max_workers = min(2, system_cpus)
        else:
            max_workers = min(4, system_cpus)

    # 4. Defensive Triggering Rules
    if cfg.chunk_size is not None:
        chunk_size = cfg.chunk_size
        use_chunking = True
        print(f"[ESG] Manual override: Chunked execution activated (chunk size: {chunk_size}).")
    else:
        use_chunking = (cfg.num_scenarios > 20000) or (peak_mem_needed > (0.25 * usable_ram))
        if use_chunking:
            temp_mem_per_scenario = 64 * cfg.steps
            total_temp_limit = 256 * 1024 * 1024
            temp_limit_per_worker = total_temp_limit / max_workers

            chunk_size = int(temp_limit_per_worker / temp_mem_per_scenario)
            chunk_size = min(5000, max(100, chunk_size))

            print(f"[ESG] High memory projection detected ({peak_mem_needed / (1024**3):.2f} GB estimated, {avail_ram / (1024**3):.2f} GB system available).")
            print(f"[ESG] Defensive dynamic chunking active (chunk size: {chunk_size}, active workers: {max_workers}).")
        else:
            chunk_size = cfg.num_scenarios

    num_chunks = int(np.ceil(cfg.num_scenarios / chunk_size))

    # 5. Symmetrical Master Outputs Pre-Allocation (In-Memory)
    rate_paths = np.zeros((cfg.steps + 1, cfg.num_scenarios), dtype=np.float64)
    inflation_paths = np.zeros((cfg.steps + 1, cfg.num_scenarios), dtype=np.float64)
    y_target_paths = np.zeros((cfg.steps + 1, cfg.num_scenarios), dtype=np.float64)
    mu_rate_paths = np.zeros((cfg.steps + 1, cfg.num_scenarios), dtype=np.float64)

    equity_returns = np.zeros((cfg.steps, cfg.num_scenarios), dtype=np.float64)
    cpis = np.zeros((cfg.steps, cfg.num_scenarios), dtype=np.float64)
    deposit_rates = np.zeros((cfg.steps, cfg.num_scenarios), dtype=np.float64)

    # 6. Dispatch threads filling master array views directly
    futures = []
    with ThreadPoolExecutor(max_workers=max_workers) as executor:
        for j in range(num_chunks):
            start_idx = j * chunk_size
            end_idx = min(start_idx + chunk_size, cfg.num_scenarios)

            futures.append(
                executor.submit(
                    self._run_chunk,
                    start_idx=start_idx,
                    end_idx=end_idx,
                    seed_offset=cfg.seed + start_idx,
                    ou_mu_chunk=ou_mu[start_idx:end_idx],
                    r_real_target_chunk=r_real_target[start_idx:end_idx],
                    gamma_chunk=gamma[start_idx:end_idx],
                    base_erp_chunk=base_erp[start_idx:end_idx],
                    cir_sigma_chunk=cir_sigma[start_idx:end_idx],
                    ou_sigma_chunk=ou_sigma[start_idx:end_idx],
                    gbm_sigma_chunk=gbm_sigma[start_idx:end_idx],
                    cir_theta_chunk=cir_theta[start_idx:end_idx],
                    ou_theta_chunk=ou_theta[start_idx:end_idx],
                    L=L,
                    equity_returns_view=equity_returns[:, start_idx:end_idx],
                    cpis_view=cpis[:, start_idx:end_idx],
                    deposit_rates_view=deposit_rates[:, start_idx:end_idx],
                    rate_paths_view=rate_paths[:, start_idx:end_idx],
                    inflation_paths_view=inflation_paths[:, start_idx:end_idx],
                    y_target_paths_view=y_target_paths[:, start_idx:end_idx],
                    mu_rate_paths_view=mu_rate_paths[:, start_idx:end_idx]
                )
            )

        for fut in futures:
            fut.result()

    print("[ESG] Core processing phase complete.")
    return LazyScenarioList(
        equity_returns=np.ascontiguousarray(equity_returns.T),
        cpis=np.ascontiguousarray(cpis.T),
        deposit_rates=np.ascontiguousarray(deposit_rates.T),
        rate_paths=rate_paths,
        mu_rate_paths=mu_rate_paths,
        inflation_paths=inflation_paths,
        y_target_paths=y_target_paths,
        cir_theta=cir_theta,
        cir_sigma=cir_sigma,
        ou_theta=ou_theta,
        ou_sigma=ou_sigma,
        tenors=cfg.tenors,
        pi_min=cfg.pi_min,
        lambda_irp=cfg.lambda_irp,
        kappa_irp=cfg.kappa_irp
    )

---

Performance Optimizations

aethel.engine.simulator.LazyScenarioList

A list-like adapter wrapping pre-allocated contiguous arrays. Saves memory by loading or slicing array regions on demand, fully matching standard Python list expectations.

Source code in src/aethel/engine/simulator.py

class LazyScenarioList:
    """
    A list-like adapter wrapping pre-allocated contiguous arrays.
    Saves memory by loading or slicing array regions on demand,
    fully matching standard Python list expectations.
    """
    def __init__(
        self,
        equity_returns: np.ndarray,
        cpis: np.ndarray,
        deposit_rates: np.ndarray,
        rate_paths: np.ndarray,
        mu_rate_paths: np.ndarray,
        inflation_paths: np.ndarray,
        y_target_paths: np.ndarray,
        cir_theta: np.ndarray,
        cir_sigma: np.ndarray,
        ou_theta: np.ndarray,
        ou_sigma: np.ndarray,
        tenors: np.ndarray,
        pi_min: float,
        lambda_irp: float,
        kappa_irp: float
    ):
        self.equity_returns = equity_returns
        self.cpis = cpis
        self.deposit_rates = deposit_rates
        self.rate_paths = rate_paths
        self.mu_rate_paths = mu_rate_paths
        self.inflation_paths = inflation_paths
        self.y_target_paths = y_target_paths
        self.cir_theta = cir_theta
        self.cir_sigma = cir_sigma
        self.ou_theta = ou_theta
        self.ou_sigma = ou_sigma
        self.tenors = tenors
        self.pi_min = pi_min
        self.lambda_irp = lambda_irp
        self.kappa_irp = kappa_irp
        self.num_scenarios = len(equity_returns)

    def __len__(self) -> int:
        return self.num_scenarios

    def __getitem__(self, idx: Union[int, slice]) -> Any:
        if isinstance(idx, slice):
            start, stop, step = idx.indices(self.num_scenarios)
            return [self[i] for i in range(start, stop, step)]

        if idx < 0:
            idx += self.num_scenarios
        if idx < 0 or idx >= self.num_scenarios:
            raise IndexError("Scenario index out of range.")

        # Dynamically generate yield curves for the requested scenario index on-the-fly
        nominal_yields = self._generate_scenario_nominal_yields(idx)
        real_yields = self._generate_scenario_real_yields(idx, nominal_yields)

        return {
            "stock_returns": self.equity_returns[idx],
            "cpis": self.cpis[idx],
            "deposit_rates": self.deposit_rates[idx],
            "nominal_yield_curves": nominal_yields,
            "real_yield_curves": real_yields,
            "tenors": self.tenors
        }

    def __iter__(self):
        for i in range(self.num_scenarios):
            yield self[i]

    def cleanup(self) -> None:
        """Symmetrical clean-up method."""
        pass

    def _generate_scenario_nominal_yields(self, idx: int) -> np.ndarray:
        r_path = self.rate_paths[:, idx]
        mu_path = self.mu_rate_paths[:, idx]
        theta = self.cir_theta[idx]
        sigma = self.cir_sigma[idx]
        tenors = self.tenors

        h = np.sqrt(theta ** 2 + 2.0 * (sigma ** 2))
        denominator = (theta + h) * (np.exp(h * tenors) - 1.0) + 2.0 * h
        base_A = (2.0 * h * np.exp((theta + h) * tenors / 2.0)) / denominator
        B_tau = (2.0 * (np.exp(h * tenors) - 1.0)) / denominator

        log_base_A_div_tau = np.log(base_A) / tenors
        B_tau_div_tau = B_tau / tenors
        safe_sigma_sq = max(1e-6, sigma ** 2)
        power_factor = (2.0 * theta * mu_path) / safe_sigma_sq

        yields = r_path[:, np.newaxis] * B_tau_div_tau[np.newaxis, :]
        yields -= power_factor[:, np.newaxis] * log_base_A_div_tau[np.newaxis, :]
        return yields

    def _generate_scenario_real_yields(self, idx: int, nominal_yields: np.ndarray) -> np.ndarray:
        inflation_path = self.inflation_paths[:, idx]
        mu_local_path = self.y_target_paths[:, idx] + self.pi_min
        theta = self.ou_theta[idx]
        sigma = self.ou_sigma[idx]
        tenors = self.tenors

        theta_tau = theta * tenors
        factor = np.where(
            theta_tau > 1e-4,
            (1.0 - np.exp(-theta_tau)) / theta_tau,
            1.0 - 0.5 * theta_tau + (theta_tau ** 2) / 6.0
        )

        diff = inflation_path - mu_local_path
        irp = (self.lambda_irp * sigma) * (1.0 - np.exp(-self.kappa_irp * tenors))

        yields_real = nominal_yields - mu_local_path[:, np.newaxis]
        yields_real -= diff[:, np.newaxis] * factor[np.newaxis, :]
        yields_real -= irp[np.newaxis, :]
        return yields_real

    def __repr__(self) -> str:
        return f"LazyScenarioList(scenarios={self.num_scenarios}, steps={self.equity_returns.shape[1]})"

cleanup()

Symmetrical clean-up method.

Source code in src/aethel/engine/simulator.py

def cleanup(self) -> None:
    """Symmetrical clean-up method."""
    pass

---

Parameter Calibration

Unified Calibrator

aethel.calibration.calibrator.MarketCalibrator

Unified API orchestrating model-parameter calibration from raw input historical series.

Source code in src/aethel/calibration/calibrator.py

class MarketCalibrator:
    """
    Unified API orchestrating model-parameter calibration
    from raw input historical series.
    """

    def __init__(self, base_config: Optional[SimulatorConfig] = None):
        self.config = base_config if base_config is not None else SimulatorConfig()

    def fit(
        self,
        historical_inflation: np.ndarray,
        historical_rates: np.ndarray,
        historical_equity_returns: np.ndarray,
        historical_yield_curve: Optional[np.ndarray] = None,
        tenors: Optional[np.ndarray] = None
    ) -> SimulatorConfig:
        """
        Calibrates model configurations based on historical data.
        Returns a new SimulatorConfig containing the calibrated parameters.
        """
        # 1. Calibrate Ornstein-Uhlenbeck (OU) Inflation parameters
        inf_params = InflationCalibrator.calibrate(historical_inflation, dt=self.config.dt)

        # 2. Calibrate CIR Short Rate Parameters
        if historical_yield_curve is not None and tenors is not None:
            # Match yield curve structure
            initial_rate = historical_rates[-1]
            rate_params = RatesCalibrator.fit_yield_curve_to_target(
                target_yields=historical_yield_curve,
                tenors=tenors,
                initial_rate=initial_rate
            )
        else:
            # Fallback to time-series analysis
            rate_params = RatesCalibrator.calibrate_short_rate_series(historical_rates, dt=self.config.dt)

        # 3. Calibrate Merton Jump Diffusion parameters
        eq_params = EquityCalibrator.calibrate(historical_equity_returns, dt=self.config.dt)

        # 4. Construct a new configuration inheriting non-calibrated variables
        new_config = SimulatorConfig(
            duration_years=self.config.duration_years,
            num_scenarios=self.config.num_scenarios,
            seed=self.config.seed,
            tenors=self.config.tenors,
            mu_min=self.config.mu_min,
            pi_min=self.config.pi_min,
            alpha_smooth=self.config.alpha_smooth,
            beta_drag=self.config.beta_drag,
            eta_erp=self.config.eta_erp,
            lambda_irp=self.config.lambda_irp,
            kappa_irp=self.config.kappa_irp,
            initial_rate=float(historical_rates[-1]),
            initial_inflation=float(historical_inflation[-1]),

            # Calibrated variables
            ou_mu=inf_params["ou_mu"],
            lambda_J=eq_params["lambda_J"],
            mu_J=eq_params["mu_J"],
            sigma_J=eq_params["sigma_J"]
        )

        # Storing directly as attributes to allow the Simulator to read them
        new_config.cir_theta_val = rate_params["cir_theta"]
        new_config.cir_sigma_val = rate_params["cir_sigma"]
        new_config.cir_mu_val = rate_params["cir_mu"]

        new_config.ou_theta_val = inf_params["ou_theta"]
        new_config.ou_sigma_val = inf_params["ou_sigma"]

        new_config.gbm_sigma_val = eq_params["gbm_sigma"]

        return new_config

fit(historical_inflation, historical_rates, historical_equity_returns, historical_yield_curve=None, tenors=None)

Calibrates model configurations based on historical data. Returns a new SimulatorConfig containing the calibrated parameters.

Source code in src/aethel/calibration/calibrator.py

def fit(
    self,
    historical_inflation: np.ndarray,
    historical_rates: np.ndarray,
    historical_equity_returns: np.ndarray,
    historical_yield_curve: Optional[np.ndarray] = None,
    tenors: Optional[np.ndarray] = None
) -> SimulatorConfig:
    """
    Calibrates model configurations based on historical data.
    Returns a new SimulatorConfig containing the calibrated parameters.
    """
    # 1. Calibrate Ornstein-Uhlenbeck (OU) Inflation parameters
    inf_params = InflationCalibrator.calibrate(historical_inflation, dt=self.config.dt)

    # 2. Calibrate CIR Short Rate Parameters
    if historical_yield_curve is not None and tenors is not None:
        # Match yield curve structure
        initial_rate = historical_rates[-1]
        rate_params = RatesCalibrator.fit_yield_curve_to_target(
            target_yields=historical_yield_curve,
            tenors=tenors,
            initial_rate=initial_rate
        )
    else:
        # Fallback to time-series analysis
        rate_params = RatesCalibrator.calibrate_short_rate_series(historical_rates, dt=self.config.dt)

    # 3. Calibrate Merton Jump Diffusion parameters
    eq_params = EquityCalibrator.calibrate(historical_equity_returns, dt=self.config.dt)

    # 4. Construct a new configuration inheriting non-calibrated variables
    new_config = SimulatorConfig(
        duration_years=self.config.duration_years,
        num_scenarios=self.config.num_scenarios,
        seed=self.config.seed,
        tenors=self.config.tenors,
        mu_min=self.config.mu_min,
        pi_min=self.config.pi_min,
        alpha_smooth=self.config.alpha_smooth,
        beta_drag=self.config.beta_drag,
        eta_erp=self.config.eta_erp,
        lambda_irp=self.config.lambda_irp,
        kappa_irp=self.config.kappa_irp,
        initial_rate=float(historical_rates[-1]),
        initial_inflation=float(historical_inflation[-1]),

        # Calibrated variables
        ou_mu=inf_params["ou_mu"],
        lambda_J=eq_params["lambda_J"],
        mu_J=eq_params["mu_J"],
        sigma_J=eq_params["sigma_J"]
    )

    # Storing directly as attributes to allow the Simulator to read them
    new_config.cir_theta_val = rate_params["cir_theta"]
    new_config.cir_sigma_val = rate_params["cir_sigma"]
    new_config.cir_mu_val = rate_params["cir_mu"]

    new_config.ou_theta_val = inf_params["ou_theta"]
    new_config.ou_sigma_val = inf_params["ou_sigma"]

    new_config.gbm_sigma_val = eq_params["gbm_sigma"]

    return new_config

Core Model Sub-Calibrators

aethel.calibration.inflation.InflationCalibrator

Calibrates Shifted-CIR / Ornstein-Uhlenbeck (OU) inflation parameters using historical monthly inflation series via OLS AR(1) mapping.

Source code in src/aethel/calibration/inflation.py

class InflationCalibrator:
    """
    Calibrates Shifted-CIR / Ornstein-Uhlenbeck (OU) inflation parameters
    using historical monthly inflation series via OLS AR(1) mapping.
    """

    @staticmethod
    def calibrate(historical_inflation: np.ndarray, dt: float = 1.0 / 12.0) -> dict:
        """
        Fits an AR(1) process to the inflation series and maps the parameters
        to continuous-time Ornstein-Uhlenbeck parameters.
        """
        inflation = np.asarray(historical_inflation, dtype=np.float64)
        if len(inflation) < 3:
            raise ValueError("Historical inflation series must contain at least 3 historical points.")

        # Prepare lagged series
        x_t = inflation[1:]
        x_lag = inflation[:-1]

        # Fit OLS: x_t = a + b * x_lag + e
        poly = np.polyfit(x_lag, x_t, deg=1)
        b = poly[0]
        a = poly[1]

        residuals = x_t - (a + b * x_lag)
        residual_var = np.var(residuals, ddof=2)

        # Enforce stability constraints
        if b <= 0.0 or b >= 1.0:
            # Fallback to realistic bounds if non-stationary
            b = np.clip(b, 0.01, 0.99)

        # Map AR(1) to continuous-time OU
        theta_ou = -np.log(b) / dt
        mu_ou = a / (1.0 - b)

        # Continuous variance mapping
        sigma_ou_sq = residual_var * (2.0 * theta_ou) / (1.0 - b**2)
        sigma_ou = np.sqrt(np.maximum(1e-6, sigma_ou_sq))

        return {
            "ou_theta": float(theta_ou),
            "ou_mu": float(mu_ou),
            "ou_sigma": float(sigma_ou)
        }

calibrate(historical_inflation, dt=1.0 / 12.0) staticmethod

Fits an AR(1) process to the inflation series and maps the parameters to continuous-time Ornstein-Uhlenbeck parameters.

Source code in src/aethel/calibration/inflation.py

@staticmethod
def calibrate(historical_inflation: np.ndarray, dt: float = 1.0 / 12.0) -> dict:
    """
    Fits an AR(1) process to the inflation series and maps the parameters
    to continuous-time Ornstein-Uhlenbeck parameters.
    """
    inflation = np.asarray(historical_inflation, dtype=np.float64)
    if len(inflation) < 3:
        raise ValueError("Historical inflation series must contain at least 3 historical points.")

    # Prepare lagged series
    x_t = inflation[1:]
    x_lag = inflation[:-1]

    # Fit OLS: x_t = a + b * x_lag + e
    poly = np.polyfit(x_lag, x_t, deg=1)
    b = poly[0]
    a = poly[1]

    residuals = x_t - (a + b * x_lag)
    residual_var = np.var(residuals, ddof=2)

    # Enforce stability constraints
    if b <= 0.0 or b >= 1.0:
        # Fallback to realistic bounds if non-stationary
        b = np.clip(b, 0.01, 0.99)

    # Map AR(1) to continuous-time OU
    theta_ou = -np.log(b) / dt
    mu_ou = a / (1.0 - b)

    # Continuous variance mapping
    sigma_ou_sq = residual_var * (2.0 * theta_ou) / (1.0 - b**2)
    sigma_ou = np.sqrt(np.maximum(1e-6, sigma_ou_sq))

    return {
        "ou_theta": float(theta_ou),
        "ou_mu": float(mu_ou),
        "ou_sigma": float(sigma_ou)
    }

aethel.calibration.rates.RatesCalibrator

Calibrates CIR short-rate parameter matrices and expectations using historical short-rates and yield curve term structures.

Source code in src/aethel/calibration/rates.py

class RatesCalibrator:
    """
    Calibrates CIR short-rate parameter matrices and expectations
    using historical short-rates and yield curve term structures.
    """

    @staticmethod
    def calibrate_short_rate_series(historical_rates: np.ndarray, dt: float = 1.0 / 12.0) -> dict:
        """
        Performs discrete-time estimation of the CIR drift and diffusion
        parameters based on historical short rate transitions.
        """
        r = np.asarray(historical_rates, dtype=np.float64)
        if len(r) < 3:
            raise ValueError("Historical short rate series must have at least 3 observations.")

        # CIR: dx_t = theta * (mu - x_t) * dt + sigma * sqrt(x_t) * dW_t
        r_t = r[:-1]
        r_next = r[1:]

        y = (r_next - r_t) / np.sqrt(np.maximum(1e-5, r_t))
        x1 = dt / np.sqrt(np.maximum(1e-5, r_t))
        x2 = -dt * np.sqrt(np.maximum(1e-5, r_t))

        X = np.column_stack((x1, x2))

        beta, _, _, _ = np.linalg.lstsq(X, y, rcond=None)

        theta_mu = beta[0]
        theta = beta[1]

        if theta <= 0.0:
            theta = 0.25
        mu = np.maximum(0.01, theta_mu / theta)

        predicted = theta * (mu - r_t) * dt
        actual_diff = r_next - r_t
        residuals = (actual_diff - predicted) / np.sqrt(np.maximum(1e-5, r_t))
        sigma = np.sqrt(np.var(residuals, ddof=1) / dt)

        return {
            "cir_theta": float(theta),
            "cir_mu": float(mu),
            "cir_sigma": float(np.clip(sigma, 0.01, 0.20))
        }

    @staticmethod
    def fit_yield_curve_to_target(
        target_yields: np.ndarray,
        tenors: np.ndarray,
        initial_rate: float
    ) -> dict:
        """
        Calibrates CIR parameters by matching analytical yield formulas
        to an observed target yield curve using optimization.
        """
        target = np.asarray(target_yields, dtype=np.float64)
        tau = np.asarray(tenors, dtype=np.float64)

        def objective(params):
            theta, mu, sigma = params
            if theta <= 0.01 or mu <= 0.001 or sigma <= 0.001:
                return 1e10

            h = np.sqrt(theta**2 + 2.0 * sigma**2)
            denominator = (theta + h) * (np.exp(h * tau) - 1.0) + 2.0 * h
            base_A = (2.0 * h * np.exp((theta + h) * tau / 2.0)) / denominator
            B_tau = (2.0 * (np.exp(h * tau) - 1.0)) / denominator

            log_base_A_div_tau = np.log(base_A) / tau
            B_tau_div_tau = B_tau / tau

            power_factor = (2.0 * theta * mu) / (sigma**2)
            model_yields = initial_rate * B_tau_div_tau - power_factor * log_base_A_div_tau

            return np.mean((target - model_yields) ** 2)

        init_guess = [0.25, np.mean(target), 0.08]
        bounds = [(0.01, 2.0), (0.005, 0.25), (0.005, 0.25)]

        res = minimize(objective, init_guess, method="L-BFGS-B", bounds=bounds)
        if not res.success:
            return {"cir_theta": 0.25, "cir_mu": float(np.mean(target)), "cir_sigma": 0.08}

        return {
            "cir_theta": float(res.x[0]),
            "cir_mu": float(res.x[1]),
            "cir_sigma": float(res.x[2])
        }

calibrate_short_rate_series(historical_rates, dt=1.0 / 12.0) staticmethod

Performs discrete-time estimation of the CIR drift and diffusion parameters based on historical short rate transitions.

Source code in src/aethel/calibration/rates.py

@staticmethod
def calibrate_short_rate_series(historical_rates: np.ndarray, dt: float = 1.0 / 12.0) -> dict:
    """
    Performs discrete-time estimation of the CIR drift and diffusion
    parameters based on historical short rate transitions.
    """
    r = np.asarray(historical_rates, dtype=np.float64)
    if len(r) < 3:
        raise ValueError("Historical short rate series must have at least 3 observations.")

    # CIR: dx_t = theta * (mu - x_t) * dt + sigma * sqrt(x_t) * dW_t
    r_t = r[:-1]
    r_next = r[1:]

    y = (r_next - r_t) / np.sqrt(np.maximum(1e-5, r_t))
    x1 = dt / np.sqrt(np.maximum(1e-5, r_t))
    x2 = -dt * np.sqrt(np.maximum(1e-5, r_t))

    X = np.column_stack((x1, x2))

    beta, _, _, _ = np.linalg.lstsq(X, y, rcond=None)

    theta_mu = beta[0]
    theta = beta[1]

    if theta <= 0.0:
        theta = 0.25
    mu = np.maximum(0.01, theta_mu / theta)

    predicted = theta * (mu - r_t) * dt
    actual_diff = r_next - r_t
    residuals = (actual_diff - predicted) / np.sqrt(np.maximum(1e-5, r_t))
    sigma = np.sqrt(np.var(residuals, ddof=1) / dt)

    return {
        "cir_theta": float(theta),
        "cir_mu": float(mu),
        "cir_sigma": float(np.clip(sigma, 0.01, 0.20))
    }

fit_yield_curve_to_target(target_yields, tenors, initial_rate) staticmethod

Calibrates CIR parameters by matching analytical yield formulas to an observed target yield curve using optimization.

Source code in src/aethel/calibration/rates.py

@staticmethod
def fit_yield_curve_to_target(
    target_yields: np.ndarray,
    tenors: np.ndarray,
    initial_rate: float
) -> dict:
    """
    Calibrates CIR parameters by matching analytical yield formulas
    to an observed target yield curve using optimization.
    """
    target = np.asarray(target_yields, dtype=np.float64)
    tau = np.asarray(tenors, dtype=np.float64)

    def objective(params):
        theta, mu, sigma = params
        if theta <= 0.01 or mu <= 0.001 or sigma <= 0.001:
            return 1e10

        h = np.sqrt(theta**2 + 2.0 * sigma**2)
        denominator = (theta + h) * (np.exp(h * tau) - 1.0) + 2.0 * h
        base_A = (2.0 * h * np.exp((theta + h) * tau / 2.0)) / denominator
        B_tau = (2.0 * (np.exp(h * tau) - 1.0)) / denominator

        log_base_A_div_tau = np.log(base_A) / tau
        B_tau_div_tau = B_tau / tau

        power_factor = (2.0 * theta * mu) / (sigma**2)
        model_yields = initial_rate * B_tau_div_tau - power_factor * log_base_A_div_tau

        return np.mean((target - model_yields) ** 2)

    init_guess = [0.25, np.mean(target), 0.08]
    bounds = [(0.01, 2.0), (0.005, 0.25), (0.005, 0.25)]

    res = minimize(objective, init_guess, method="L-BFGS-B", bounds=bounds)
    if not res.success:
        return {"cir_theta": 0.25, "cir_mu": float(np.mean(target)), "cir_sigma": 0.08}

    return {
        "cir_theta": float(res.x[0]),
        "cir_mu": float(res.x[1]),
        "cir_sigma": float(res.x[2])
    }

aethel.calibration.equity.EquityCalibrator

Calibrates Merton Jump Diffusion parameters using Maximum Likelihood Estimation (MLE) from historical asset return paths.

Source code in src/aethel/calibration/equity.py

class EquityCalibrator:
    """
    Calibrates Merton Jump Diffusion parameters using Maximum Likelihood Estimation (MLE) 
    from historical asset return paths.
    """

    @staticmethod
    def calibrate(historical_returns: np.ndarray, dt: float = 1.0 / 12.0) -> dict:
        """
        Finds drift, diffusion, and jump parameters using historical returns.
        Uses a truncated infinite series representation of the Merton density function.
        """
        returns = np.asarray(historical_returns, dtype=np.float64)
        log_returns = np.log(1.0 + returns)

        if len(log_returns) < 10:
            raise ValueError("Merton MLE requires a minimum of 10 return observations.")

        # Numerical approximation truncation level for Poisson jump events
        K_max = 5 

        # Parameters to fit: [mu, sigma (diffusion), lambda_J, mu_J, sigma_J]
        def negative_log_likelihood(params):
            mu, sigma, lmbda, mu_J, sigma_J = params

            # Boundary constraints for physical realism
            if sigma <= 0.001 or lmbda < 0.0 or sigma_J <= 0.001:
                return 1e10
            if lmbda > 5.0:  # Restrict extreme unphysical jump rates
                return 1e10

            # Calculate Merton density for each observation
            likelihood_sum = np.zeros_like(log_returns)
            for k in range(K_max):
                poisson_prob = np.exp(-lmbda * dt) * ((lmbda * dt) ** k) / math.factorial(k)

                # Distribution parameters conditional on k jumps
                mean_k = (mu - 0.5 * sigma**2 - lmbda * (np.exp(mu_J + 0.5 * sigma_J**2) - 1.0)) * dt + k * mu_J
                variance_k = (sigma**2) * dt + k * (sigma_J**2)

                likelihood_sum += poisson_prob * norm.pdf(log_returns, loc=mean_k, scale=np.sqrt(variance_k))

            # Numerical stability adjustment
            likelihood_sum = np.maximum(likelihood_sum, 1e-12)
            return -np.sum(np.log(likelihood_sum))

        # Initial assumptions via standard sample moments
        sample_mean = np.mean(log_returns) / dt
        sample_std = np.std(log_returns) / np.sqrt(dt)

        # Assumptions: [mu, sigma, lambda_J, mu_J, sigma_J]
        init_guess = [sample_mean, sample_std * 0.80, 0.20, -0.10, 0.10]
        bounds = [
            (-0.5, 0.5),          # Drift range
            (0.01, 0.5),          # Continuous volatility
            (0.0, 2.0),           # Annual jump rate (0 to 2 occurrences per year)
            (-0.4, 0.1),          # Average jump return impact
            (0.01, 0.3)           # Jump impact uncertainty
        ]

        res = minimize(negative_log_likelihood, init_guess, method="L-BFGS-B", bounds=bounds)

        if not res.success:
            print("Warning: Fallback applied, solver failed to converge")
            # FALLBACK
            # Identify "outlier" months that are likely jump events (e.g., movements > 1.5 standard deviations)
            threshold = 1.5 * sample_std
            jumps = log_returns[np.abs(log_returns - np.mean(log_returns)) > threshold]

            if len(jumps) > 0:
                # Estimate jump characteristics directly from these outlier events
                est_lambda = len(jumps) / (len(log_returns) * dt)
                est_mu_J = float(np.mean(jumps))
                est_sigma_J = float(np.clip(np.std(jumps), 0.01, 0.20))
            else:
                # Defaults if no outliers are found
                est_lambda = 0.20
                est_mu_J = -0.10
                est_sigma_J = 0.05

            return {
                "gbm_sigma": float(sample_std * 0.90), # Leave room for the jump variance
                "lambda_J": float(np.clip(est_lambda, 0.05, 1.0)),
                "mu_J": float(np.clip(est_mu_J, -0.30, 0.0)),
                "sigma_J": est_sigma_J
            }

        return {
            "gbm_sigma": float(res.x[1]),
            "lambda_J": float(res.x[2]),
            "mu_J": float(res.x[3]),
            "sigma_J": float(res.x[4])
        }

calibrate(historical_returns, dt=1.0 / 12.0) staticmethod

Finds drift, diffusion, and jump parameters using historical returns. Uses a truncated infinite series representation of the Merton density function.

Source code in src/aethel/calibration/equity.py

@staticmethod
def calibrate(historical_returns: np.ndarray, dt: float = 1.0 / 12.0) -> dict:
    """
    Finds drift, diffusion, and jump parameters using historical returns.
    Uses a truncated infinite series representation of the Merton density function.
    """
    returns = np.asarray(historical_returns, dtype=np.float64)
    log_returns = np.log(1.0 + returns)

    if len(log_returns) < 10:
        raise ValueError("Merton MLE requires a minimum of 10 return observations.")

    # Numerical approximation truncation level for Poisson jump events
    K_max = 5 

    # Parameters to fit: [mu, sigma (diffusion), lambda_J, mu_J, sigma_J]
    def negative_log_likelihood(params):
        mu, sigma, lmbda, mu_J, sigma_J = params

        # Boundary constraints for physical realism
        if sigma <= 0.001 or lmbda < 0.0 or sigma_J <= 0.001:
            return 1e10
        if lmbda > 5.0:  # Restrict extreme unphysical jump rates
            return 1e10

        # Calculate Merton density for each observation
        likelihood_sum = np.zeros_like(log_returns)
        for k in range(K_max):
            poisson_prob = np.exp(-lmbda * dt) * ((lmbda * dt) ** k) / math.factorial(k)

            # Distribution parameters conditional on k jumps
            mean_k = (mu - 0.5 * sigma**2 - lmbda * (np.exp(mu_J + 0.5 * sigma_J**2) - 1.0)) * dt + k * mu_J
            variance_k = (sigma**2) * dt + k * (sigma_J**2)

            likelihood_sum += poisson_prob * norm.pdf(log_returns, loc=mean_k, scale=np.sqrt(variance_k))

        # Numerical stability adjustment
        likelihood_sum = np.maximum(likelihood_sum, 1e-12)
        return -np.sum(np.log(likelihood_sum))

    # Initial assumptions via standard sample moments
    sample_mean = np.mean(log_returns) / dt
    sample_std = np.std(log_returns) / np.sqrt(dt)

    # Assumptions: [mu, sigma, lambda_J, mu_J, sigma_J]
    init_guess = [sample_mean, sample_std * 0.80, 0.20, -0.10, 0.10]
    bounds = [
        (-0.5, 0.5),          # Drift range
        (0.01, 0.5),          # Continuous volatility
        (0.0, 2.0),           # Annual jump rate (0 to 2 occurrences per year)
        (-0.4, 0.1),          # Average jump return impact
        (0.01, 0.3)           # Jump impact uncertainty
    ]

    res = minimize(negative_log_likelihood, init_guess, method="L-BFGS-B", bounds=bounds)

    if not res.success:
        print("Warning: Fallback applied, solver failed to converge")
        # FALLBACK
        # Identify "outlier" months that are likely jump events (e.g., movements > 1.5 standard deviations)
        threshold = 1.5 * sample_std
        jumps = log_returns[np.abs(log_returns - np.mean(log_returns)) > threshold]

        if len(jumps) > 0:
            # Estimate jump characteristics directly from these outlier events
            est_lambda = len(jumps) / (len(log_returns) * dt)
            est_mu_J = float(np.mean(jumps))
            est_sigma_J = float(np.clip(np.std(jumps), 0.01, 0.20))
        else:
            # Defaults if no outliers are found
            est_lambda = 0.20
            est_mu_J = -0.10
            est_sigma_J = 0.05

        return {
            "gbm_sigma": float(sample_std * 0.90), # Leave room for the jump variance
            "lambda_J": float(np.clip(est_lambda, 0.05, 1.0)),
            "mu_J": float(np.clip(est_mu_J, -0.30, 0.0)),
            "sigma_J": est_sigma_J
        }

    return {
        "gbm_sigma": float(res.x[1]),
        "lambda_J": float(res.x[2]),
        "mu_J": float(res.x[3]),
        "sigma_J": float(res.x[4])
    }

---

Simulation Results & Analytics

aethel.output.results.SimulationResults

A unified analysis, post-processing, and visualization engine for ESG outputs. Maintains a robust state-caching framework and delegates operations to submodules.

Source code in src/aethel/output/results.py

class SimulationResults:
    """
    A unified analysis, post-processing, and visualization engine for ESG outputs.
    Maintains a robust state-caching framework and delegates operations to submodules.
    """

    def __init__(self, raw_scenarios: List[Dict[str, np.ndarray]]):
        if not raw_scenarios:
            raise ValueError("Scenario list cannot be empty.")

        self.scenarios = raw_scenarios
        self.num_scenarios = len(raw_scenarios)
        self.steps = len(raw_scenarios[0]["stock_returns"])
        self.tenors = list(raw_scenarios[0]["tenors"])
        self.steps_per_year = 12

        # Performance Optimization: Cache base matrices to prevent redundant calculations
        self._cache: Dict[str, np.ndarray] = {}

    def query(
        self,
        metric: str,
        stat: str = "raw",
        time: Optional[Union[str, float, int, List[float]]] = None,
        year: Optional[Union[str, float, int, List[float]]] = None,
        month: Optional[Union[str, int, List[int]]] = None,
        step: Optional[Union[str, int, List[int]]] = None,
        tenor: Optional[float] = None,
        annualized: bool = False,
        portfolio_weights: Optional[Dict[str, float]] = None
    ) -> Union[np.ndarray, float]:
        """
        Extracts, transforms, and calculates statistical parameters from the simulation.
        """
        # 1. Fetch base 2D matrix (steps, scenarios)
        matrix = self._extract_base_matrix(metric, tenor, portfolio_weights)

        # 2. Apply annualization transformations if requested
        if annualized:
            matrix = time_utils.apply_annualization(matrix, metric, self.steps_per_year)

        # 3. Resolve time inputs to indices
        step_indices = time_utils.resolve_time_to_indices(
            time_legacy=time,
            year=year,
            month=month,
            step=step,
            max_idx=len(matrix) - 1,
            steps_per_year=self.steps_per_year
        )

        # 4. Handle time horizon slicing
        if step_indices is not None:
            matrix = matrix[step_indices, :]

        # 5. Collapse dimension 1 (scenarios) according to the requested statistic
        return self._apply_statistics(matrix, stat)

    # --- Portfolio Blending Logic (Delegated) ---

    def calculate_portfolio_returns(self, weights: Dict[str, float]) -> np.ndarray:
        """
        Calculates the blended monthly return series for a given portfolio allocation.
        """
        return portfolio.calculate_portfolio_returns(self, weights)

    def calculate_portfolio_growth(self, weights: Dict[str, float]) -> np.ndarray:
        """
        Calculates the blended cumulative growth series of $1.00 for a given portfolio allocation.
        """
        return portfolio.calculate_portfolio_growth(self, weights)

    # --- Decumulation (Withdrawal) Simulator (Delegated) ---

    def simulate_decumulation(
        self,
        initial_balance: float,
        initial_monthly_withdrawal: float,
        portfolio_weights: Dict[str, float],
        withdrawal_timing: str = "beginning",
        inflate_withdrawals: bool = True,
        frictional_drag_annual: float = 0.0,
        tax_on_gains_rate: float = 0.0,
        liquidation_strategy: str = "constant_mix",
        withdrawal_policy: Optional[Callable[[np.ndarray, np.ndarray, int, np.ndarray], np.ndarray]] = None
    ) -> Dict[str, np.ndarray]:
        """
        Simulates the monthly decumulation of a starting balance.
        """
        return decumulation.simulate_decumulation(
            self,
            initial_balance,
            initial_monthly_withdrawal,
            portfolio_weights,
            withdrawal_timing,
            inflate_withdrawals,
            frictional_drag_annual,
            tax_on_gains_rate,
            liquidation_strategy,
            withdrawal_policy
        )

    # --- Optimized Base Matrix Retrieval (With Memoization) ---

    def _extract_base_matrix(
        self,
        metric: str,
        tenor: Optional[float] = None,
        portfolio_weights: Optional[Dict[str, float]] = None
    ) -> np.ndarray:
        m = metric.lower().strip()

        # Build stable cache key for portfolio blends
        if portfolio_weights is not None:
            sorted_weights = sorted(portfolio_weights.items())
            weights_str = "_".join([f"{k}:{v:.4f}" for k, v in sorted_weights])
            cache_key = f"{m}_{weights_str}"
        else:
            cache_key = f"{m}_{tenor}" if tenor is not None else m

        # Return cached array if already computed to avoid redundant loops
        if cache_key in self._cache:
            return self._cache[cache_key]

        # Check for optimized LazyScenarioList processing
        from aethel.engine.simulator import LazyScenarioList
        is_lazy = isinstance(self.scenarios, LazyScenarioList)

        if m == "portfolio_returns":
            res = self.calculate_portfolio_returns(portfolio_weights)

        elif m == "portfolio_growth":
            res = self.calculate_portfolio_growth(portfolio_weights)

        elif m == "decumulation_balance":
            if "decumulation_balance" not in self._cache:
                raise ValueError(
                    "Decumulation balance not simulated yet. Please run "
                    "simulate_decumulation(...) first to populate the cache."
                )
            return self._cache["decumulation_balance"]

        elif m == "decumulation_withdrawal":
            if "decumulation_withdrawal" not in self._cache:
                raise ValueError(
                    "Decumulation withdrawal not simulated yet. Please run "
                    "simulate_decumulation(...) first to populate the cache."
                )
            return self._cache["decumulation_withdrawal"]

        elif m in {"equity_returns", "returns"}:
            if is_lazy:
                res = self.scenarios.equity_returns.T
            else:
                res = np.column_stack([s["stock_returns"] for s in self.scenarios])

        elif m in {"equity_growth", "growth"}:
            returns = self._extract_base_matrix("returns")
            growth = np.vstack([np.zeros(self.num_scenarios), returns])
            res = np.cumprod(1.0 + growth, axis=0)

        elif m in {"cpi", "inflation_index"}:
            if is_lazy:
                res = self.scenarios.cpis.T
            else:
                res = np.column_stack([s["cpis"] for s in self.scenarios])

        elif m in {"inflation_rate", "inflation", "ipca"}:
            cpis_matrix = self._extract_base_matrix("cpi")
            cpis_padded = np.vstack([np.ones(self.num_scenarios), cpis_matrix])
            with np.errstate(divide='ignore', invalid='ignore'):
                monthly_rates = (cpis_padded[1:, :] / cpis_padded[:-1, :]) - 1.0
            res = (1.0 + monthly_rates) ** 12 - 1.0

        elif m in {"short_rate", "rate", "cdi", "deposit_rates"}:
            if is_lazy:
                res = (1.0 + self.scenarios.deposit_rates.T) ** 12 - 1.0
            else:
                monthly_rates = np.column_stack([s["deposit_rates"] for s in self.scenarios])
                res = (1.0 + monthly_rates) ** 12 - 1.0

        elif m in {"nominal_yield", "real_yield"}:
            if tenor is None:
                raise ValueError(f"A 'tenor' must be specified when querying '{metric}'. Available: {self.tenors}")
            if tenor not in self.tenors:
                raise ValueError(f"Tenor {tenor} not available. Choose from: {self.tenors}")

            if is_lazy:
                if m == "real_yield":
                    res = self._generate_real_yields_for_tenor(tenor)[:-1, :]
                else:
                    res = self._generate_cir_yields_for_tenor(tenor)[:-1, :]
            else:
                tenor_idx = self.tenors.index(tenor)
                key = "real_yield_curves" if m == "real_yield" else "nominal_yield_curves"
                res = np.column_stack([s[key][:-1, tenor_idx] for s in self.scenarios])

        else:
            raise ValueError(f"Unknown metric '{metric}'. Choose from correct keywords.")

        self._cache[cache_key] = res
        return res

    # --- On-The-Fly Vectorized Tenor Yield Derivations ---

    def _generate_cir_yields_for_tenor(self, tenor: float) -> np.ndarray:
        """Derives CIR nominal yields for a single tenor across all scenarios instantly."""
        r = self.scenarios.rate_paths.T
        mu = self.scenarios.mu_rate_paths.T
        theta = self.scenarios.cir_theta
        sigma = self.scenarios.cir_sigma

        h = np.sqrt(theta ** 2 + 2.0 * (sigma ** 2))
        denominator = (theta + h) * (np.exp(h * tenor) - 1.0) + 2.0 * h
        base_A = (2.0 * h * np.exp((theta + h) * tenor / 2.0)) / denominator
        B_tau = (2.0 * (np.exp(h * tenor) - 1.0)) / denominator

        log_base_A_div_tenor = np.log(base_A) / tenor
        B_tau_div_tenor = B_tau / tenor
        safe_sigma_sq = np.maximum(1e-6, sigma ** 2)
        power_factor = (2.0 * theta[:, np.newaxis] * mu) / safe_sigma_sq[:, np.newaxis]

        yields = r * B_tau_div_tenor[:, np.newaxis]
        yields -= power_factor * log_base_A_div_tenor[:, np.newaxis]
        return yields.T

    def _generate_real_yields_for_tenor(self, tenor: float) -> np.ndarray:
        """Derives Fisher real yields for a single tenor across all scenarios on-the-fly."""
        yields_nominal = self._generate_cir_yields_for_tenor(tenor)

        inflation_rates = self.scenarios.inflation_paths
        mu_local = self.scenarios.y_target_paths + self.scenarios.pi_min

        theta = self.scenarios.ou_theta
        sigma = self.scenarios.ou_sigma

        theta_tau = theta * tenor
        factor = np.where(
            theta_tau > 1e-4,
            (1.0 - np.exp(-theta_tau)) / theta_tau,
            1.0 - 0.5 * theta_tau + (theta_tau ** 2) / 6.0
        )

        diff = inflation_rates - mu_local
        irp = (self.scenarios.lambda_irp * sigma) * (1.0 - np.exp(-self.scenarios.kappa_irp * tenor))

        yields_real = yields_nominal - mu_local
        yields_real -= diff * factor[np.newaxis, :]
        yields_real -= irp[np.newaxis, :]
        return yields_real

    def _apply_statistics(self, matrix: np.ndarray, stat: str) -> Union[np.ndarray, float]:
        s = stat.lower().strip()

        if s == "raw":
            return matrix[0] if matrix.shape[0] == 1 else matrix

        if s == "mean":
            aggregated = np.mean(matrix, axis=1)
        elif s == "median":
            aggregated = np.percentile(matrix, 50.0, axis=1)
        elif s in {"std", "vol", "volatility"}:
            aggregated = np.std(matrix, axis=1)
        elif s.startswith("p"):
            try:
                percentile_val = float(s[1:])
                if not (0.0 <= percentile_val <= 100.0):
                    raise ValueError
                aggregated = np.percentile(matrix, percentile_val, axis=1)
            except ValueError:
                raise ValueError(f"Could not parse valid percentile from '{stat}'.")
        else:
            raise ValueError(f"Unknown statistic '{stat}'.")

        return float(aggregated[0]) if len(aggregated) == 1 else aggregated

    # --- Backward Compatibility Methods ---

    def get_matrix(self, metric: str) -> np.ndarray:
        valid_map = {
            "stock_returns": "returns",
            "cpis": "cpi",
            "deposit_rates": "deposit_rates"
        }
        mapped = valid_map.get(metric)
        if not mapped:
            raise ValueError(f"Metric '{metric}' must be one of {list(valid_map.keys())}")
        return self._extract_base_matrix(mapped, tenor=None)

    def get_yield_curve_matrix(self, real: bool = False, tenor_idx: int = 0) -> np.ndarray:
        tenor = self.tenors[tenor_idx]
        metric = "real_yield" if real else "nominal_yield"
        return self._extract_base_matrix(metric, tenor=tenor)

    def get_summary_statistics(self, metric: str, percentiles: List[float] = [5.0, 50.0, 95.0]) -> Dict[str, np.ndarray]:
        valid_map = {"cpis": "cpi", "deposit_rates": "deposit_rates", "stock_returns": "returns"}
        mapped = valid_map.get(metric)
        if not mapped:
            raise ValueError(f"Metric must be one of {list(valid_map.keys())}")

        stats = {}
        for p in percentiles:
            stats[f"p{int(p)}"] = self.query(mapped, stat=f"p{p}", step="all")
        stats["mean"] = self.query(mapped, stat="mean", step="all")
        return stats

    def to_pandas(self, scenario_idx: int) -> pd.DataFrame:
        if scenario_idx < 0 or scenario_idx >= self.num_scenarios:
            raise IndexError(f"Scenario index must be between 0 and {self.num_scenarios - 1}")

        scenario = self.scenarios[scenario_idx]
        df = pd.DataFrame({
            "step": np.arange(self.steps),
            "stock_returns": scenario["stock_returns"],
            "cpis": scenario["cpis"],
            "deposit_rates": scenario["deposit_rates"]
        })

        for idx, tenor in enumerate(self.tenors):
            df[f"nominal_yield_{tenor}y"] = scenario["nominal_yield_curves"][:-1, idx]
            df[f"real_yield_{tenor}y"] = scenario["real_yield_curves"][:-1, idx]

        return df.set_index("step")

    def cleanup(self) -> None:
        """Clears cached matrices."""
        self._cache.clear()
        if hasattr(self.scenarios, "cleanup"):
            self.scenarios.cleanup()

    # --- Plotly Interactive Plotting Helpers (Delegated) ---

    def plot_fan_chart(
        self,
        metric: str,
        tenor: Optional[float] = None,
        annualized: bool = False,
        title: Optional[str] = None
    ) -> Any:
        return visualizer.plot_fan_chart(self, metric, tenor, annualized, title)

    def plot_scenario_paths(
        self,
        metric: str,
        num_paths: int = 15,
        tenor: Optional[float] = None,
        annualized: bool = False,
        title: Optional[str] = None
    ) -> Any:
        return visualizer.plot_scenario_paths(self, metric, num_paths, tenor, annualized, title)

    def plot_yield_curve_evolution(
        self,
        years_milestones: List[float] = [0.0, 1.0, 5.0, 15.0, 30.0],
        real: bool = False,
        title: Optional[str] = None
    ) -> Any:
        return visualizer.plot_yield_curve_evolution(self, years_milestones, real, title)

    def plot_horizon_distribution(
        self,
        metric: str,
        target_year: float,
        tenor: Optional[float] = None,
        annualized: bool = False,
        bins: int = 40,
        title: Optional[str] = None
    ) -> Any:
        return visualizer.plot_horizon_distribution(self, metric, target_year, tenor, annualized, bins, title)

calculate_portfolio_growth(weights)

Calculates the blended cumulative growth series of $1.00 for a given portfolio allocation.

Source code in src/aethel/output/results.py

def calculate_portfolio_growth(self, weights: Dict[str, float]) -> np.ndarray:
    """
    Calculates the blended cumulative growth series of $1.00 for a given portfolio allocation.
    """
    return portfolio.calculate_portfolio_growth(self, weights)

calculate_portfolio_returns(weights)

Calculates the blended monthly return series for a given portfolio allocation.

Source code in src/aethel/output/results.py

def calculate_portfolio_returns(self, weights: Dict[str, float]) -> np.ndarray:
    """
    Calculates the blended monthly return series for a given portfolio allocation.
    """
    return portfolio.calculate_portfolio_returns(self, weights)

cleanup()

Clears cached matrices.

Source code in src/aethel/output/results.py

def cleanup(self) -> None:
    """Clears cached matrices."""
    self._cache.clear()
    if hasattr(self.scenarios, "cleanup"):
        self.scenarios.cleanup()

query(metric, stat='raw', time=None, year=None, month=None, step=None, tenor=None, annualized=False, portfolio_weights=None)

Extracts, transforms, and calculates statistical parameters from the simulation.

Source code in src/aethel/output/results.py

def query(
    self,
    metric: str,
    stat: str = "raw",
    time: Optional[Union[str, float, int, List[float]]] = None,
    year: Optional[Union[str, float, int, List[float]]] = None,
    month: Optional[Union[str, int, List[int]]] = None,
    step: Optional[Union[str, int, List[int]]] = None,
    tenor: Optional[float] = None,
    annualized: bool = False,
    portfolio_weights: Optional[Dict[str, float]] = None
) -> Union[np.ndarray, float]:
    """
    Extracts, transforms, and calculates statistical parameters from the simulation.
    """
    # 1. Fetch base 2D matrix (steps, scenarios)
    matrix = self._extract_base_matrix(metric, tenor, portfolio_weights)

    # 2. Apply annualization transformations if requested
    if annualized:
        matrix = time_utils.apply_annualization(matrix, metric, self.steps_per_year)

    # 3. Resolve time inputs to indices
    step_indices = time_utils.resolve_time_to_indices(
        time_legacy=time,
        year=year,
        month=month,
        step=step,
        max_idx=len(matrix) - 1,
        steps_per_year=self.steps_per_year
    )

    # 4. Handle time horizon slicing
    if step_indices is not None:
        matrix = matrix[step_indices, :]

    # 5. Collapse dimension 1 (scenarios) according to the requested statistic
    return self._apply_statistics(matrix, stat)

simulate_decumulation(initial_balance, initial_monthly_withdrawal, portfolio_weights, withdrawal_timing='beginning', inflate_withdrawals=True, frictional_drag_annual=0.0, tax_on_gains_rate=0.0, liquidation_strategy='constant_mix', withdrawal_policy=None)

Simulates the monthly decumulation of a starting balance.

Source code in src/aethel/output/results.py

def simulate_decumulation(
    self,
    initial_balance: float,
    initial_monthly_withdrawal: float,
    portfolio_weights: Dict[str, float],
    withdrawal_timing: str = "beginning",
    inflate_withdrawals: bool = True,
    frictional_drag_annual: float = 0.0,
    tax_on_gains_rate: float = 0.0,
    liquidation_strategy: str = "constant_mix",
    withdrawal_policy: Optional[Callable[[np.ndarray, np.ndarray, int, np.ndarray], np.ndarray]] = None
) -> Dict[str, np.ndarray]:
    """
    Simulates the monthly decumulation of a starting balance.
    """
    return decumulation.simulate_decumulation(
        self,
        initial_balance,
        initial_monthly_withdrawal,
        portfolio_weights,
        withdrawal_timing,
        inflate_withdrawals,
        frictional_drag_annual,
        tax_on_gains_rate,
        liquidation_strategy,
        withdrawal_policy
    )