Bond Risk Measures and Immunization › Immunization

def immunize_portfolio(liability_pv, target_duration, available_bonds, yield_rate):
    """Match a positive liability PV and its Macaulay horizon using two bonds."""
    from math import isclose, isfinite

    if not isfinite(liability_pv) or liability_pv <= 0:
        raise ValueError("liability_pv must be positive and finite")
    if not isfinite(target_duration) or target_duration <= 0:
        raise ValueError("target_duration must be positive and finite")
    if len(available_bonds) != 2:
        raise ValueError("This simple immunization requires exactly 2 bonds")
    durations = [
        macaulay_duration(b["face_value"], b["coupon_rate"], b["maturity"], yield_rate)
        for b in available_bonds
    ]
    prices = [
        bond_price(b["face_value"], b["coupon_rate"], b["maturity"], yield_rate)
        for b in available_bonds
    ]
    d1, d2 = durations
    if isclose(d1, d2, rel_tol=1e-12, abs_tol=1e-12):
        raise ValueError("Two distinct durations are required")
    w1 = (target_duration - d2) / (d1 - d2)
    w2 = 1 - w1
    if min(w1, w2) < 0:
        raise ValueError("The target horizon requires a short position")
    return [w1, w2], durations, prices

# Immunization example: $1,000,000 liability due in 7 years
liability_amount = 1_000_000
liability_years = 7
yield_rate = 0.05

# Present value of liability
frequency = 2
liability_pv = liability_amount / (1 + yield_rate / frequency) ** (
    frequency * liability_years
)

# Available bonds for immunization
immunization_bonds = [
    {
        "face_value": 1000,
        "coupon_rate": 0.04,
        "maturity": 5,
        "name": "5-year, 4%",
    },
    {
        "face_value": 1000,
        "coupon_rate": 0.06,
        "maturity": 15,
        "name": "15-year, 6%",
    },
]

weights, durations, prices = immunize_portfolio(
    liability_pv, liability_years, immunization_bonds, yield_rate
)

print(f"Liability: ${liability_amount:,.0f} due in {liability_years} years")
print(f"Present Value of Liability: ${liability_pv:,.2f}")
print(f"Target Macaulay Duration: {liability_years} years")
print()
print("Immunization Portfolio:")
print("-" * 60)

for i, (bond, w, d, p) in enumerate(
    zip(immunization_bonds, weights, durations, prices)
):
    investment = liability_pv * w
    num_bonds = investment / p
    print(f"{bond['name']}:")
    print(f"  Duration: {d:.2f} years")
    print(f"  Weight: {w:.2%}")
    print(f"  Investment: ${investment:,.2f}")
    print(f"  Number of bonds: {num_bonds:.1f}")
    print()

# Verify portfolio duration
port_duration = weights[0] * durations[0] + weights[1] * durations[1]
print(f"Portfolio Macaulay Duration: {port_duration:.2f} years")

def simulate_immunization(
    liability_amount,
    liability_years,
    bonds,
    initial_weights,
    initial_yield,
    yield_shock,
):
    """Horizon value of fixed semiannual cash flows after a permanent flat shock.

    Payments before the coupon-date horizon compound forward. Payments after
    the horizon are discounted back, giving the residual bond's sale value.
    Assume no default, frictionless reinvestment and fractional holdings.
    """
    from math import isclose, isfinite

    frequency = 2
    if not isfinite(liability_amount) or liability_amount <= 0:
        raise ValueError("liability_amount must be positive and finite")
    if not isfinite(liability_years) or liability_years <= 0:
        raise ValueError("liability_years must be positive and finite")
    horizon_periods = round(liability_years * frequency)
    if not isclose(
        horizon_periods, liability_years * frequency, rel_tol=0, abs_tol=1e-10
    ):
        raise ValueError("The horizon must lie on the semiannual payment grid")
    if len(bonds) != len(initial_weights) or not bonds:
        raise ValueError("Provide one weight per bond")
    if any(not isfinite(w) or w < 0 for w in initial_weights):
        raise ValueError("Weights must be finite and nonnegative")
    if not isclose(sum(initial_weights), 1.0, rel_tol=0, abs_tol=1e-12):
        raise ValueError("Weights must sum to one")
    if not isfinite(initial_yield) or not isfinite(yield_shock):
        raise ValueError("Yields must be finite")
    initial_base = 1 + initial_yield / frequency
    shocked_base = 1 + (initial_yield + yield_shock) / frequency
    if min(initial_base, shocked_base) <= 0:
        raise ValueError("Periodic discount bases must be positive")
    initial_pv = liability_amount / initial_base**horizon_periods
    terminal_value = 0.0
    for bond, weight in zip(bonds, initial_weights):
        old_price = bond_price(
            bond["face_value"],
            bond["coupon_rate"],
            bond["maturity"],
            initial_yield,
            frequency,
        )
        quantity = initial_pv * weight / old_price
        periods = round(bond["maturity"] * frequency)
        if not isclose(periods, bond["maturity"] * frequency, rel_tol=0, abs_tol=1e-10):
            raise ValueError("Bond maturity must lie on the payment grid")
        coupon = bond["face_value"] * bond["coupon_rate"] / frequency
        for period in range(1, periods + 1):
            cash_flow = coupon + (bond["face_value"] if period == periods else 0)
            terminal_value += (
                quantity * cash_flow * shocked_base ** (horizon_periods - period)
            )
    return terminal_value

# Test immunization under different yield scenarios
yield_shocks = [-0.02, -0.01, 0, 0.01, 0.02]

print("Immunization Test: Terminal Values Under Different Yield Scenarios")
print("=" * 65)
print(f"{'Yield Shock':>12} {'New Yield':>12} {'Terminal Value':>16} {'vs Target':>12}")
print("-" * 65)

for shock in yield_shocks:
    terminal = simulate_immunization(
        liability_amount,
        liability_years,
        immunization_bonds,
        weights,
        yield_rate,
        shock,
    )
    diff_pct = (terminal - liability_amount) / liability_amount * 100
    print(
        f"{shock * 100:>+12.1f}% {(yield_rate + shock) * 100:>12.1f}% ${terminal:>15,.0f} {diff_pct:>+11.2f}%"
    )