Chapter 10: Interest Rate Term Structure Modeling Practice

Knowledge Summary

1. Fundamentals of Interest Rate Term Structure

The interest rate term structure describes the relationship between risk-free interest rates of different maturities. Within the Markov framework, we model interest rates as a continuous-time Markov process.

Spot Rate Curve: $r(T) = \text{spot rate from time 0 to maturity T}$

Forward Rate: $f(t,T) = \text{instantaneous forward rate at time T as seen from time t}$

Relationship: $f(0,T) = r(T) + T \frac{dr(T)}{dT}$

2. Hull-White Model

The Hull-White model is the most commonly used single-factor interest rate model, describing the dynamics of short-term interest rates:

Stochastic Differential Equation: $dr_t = (\theta(t) - a r_t)dt + \sigma dW_t$

Where:

• $r_t$ is the short-term interest rate
• $\theta(t)$ is the time-varying mean reversion level
• $a$ is the mean reversion speed
• $\sigma$ is the volatility
• $W_t$ is Brownian motion

graph TD
    A[Short rate r₀] --> B[Time t₁]
    B --> C[Time t₂]
    C --> D[Time t₃]
    D --> E[Time T]

    A -.->|θ(t)| F[Mean reversion level]
    B -.->|θ(t)| F
    C -.->|θ(t)| F
    D -.->|θ(t)| F

    style F fill:#e8f5e8
    style A fill:#e1f5fe
    style E fill:#fff3e0

3. Bond Pricing Formula

Under the Hull-White model, zero-coupon bond prices have analytical solutions:

$P(t,T) = A(t,T) e^{-B(t,T) r_t}$

Where: $B(t,T) = \frac{1-e^{-a(T-t)}}{a}$

$A(t,T) = \exp\left(\int_t^T \theta(s) B(s,T) ds - \frac{\sigma^2}{2a^2}(B(t,T) - T + t)\left(1 - e^{-2a(T-t)}\right)\right)$

4. Interest Rate Tree Construction

To price complex derivatives, we construct a discrete interest rate tree:

Trinomial Tree Method:

• Each node has three branches: up, flat, down
• Branch probabilities are calculated based on model parameters
• Ensure tree convergence and numerical stability

Example Code

Example 1: Hull-White Model Implementation

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from scipy.interpolate import interp1d
import pandas as pd

class HullWhiteModel:
    """
    Hull-White single-factor interest rate model
    """

    def __init__(self, a=0.1, sigma=0.01):
        """
        Initialize Hull-White model

        Parameters:
        a: Mean reversion speed
        sigma: Volatility
        """
        self.a = a
        self.sigma = sigma
        self.theta = None  # Time-varying function, determined through calibration

    def set_theta_function(self, theta_func):
        """Set theta function"""
        self.theta = theta_func

    def bond_price(self, r, t, T):
        """
        Calculate zero-coupon bond price

        Parameters:
        r: Current short-term interest rate
        t: Current time
        T: Maturity time

        Returns:
        Bond price
        """
        if t >= T:
            return 1.0

        tau = T - t
        B = (1 - np.exp(-self.a * tau)) / self.a if self.a != 0 else tau

        # Simplified A(t,T) calculation (assuming constant theta)
        if self.theta is None:
            theta_avg = 0.03  # Assume average theta is 3%
        else:
            theta_avg = self.theta(t) if callable(self.theta) else self.theta

        A = np.exp(theta_avg * B - self.sigma**2 / (4 * self.a) * B**2 * (1 - np.exp(-2 * self.a * tau)))

        return A * np.exp(-B * r)

    def simulate_paths(self, r0, T, n_steps, n_paths, theta_func=None):
        """
        Monte Carlo simulation of interest rate paths

        Parameters:
        r0: Initial interest rate
        T: Total time
        n_steps: Number of time steps
        n_paths: Number of paths
        theta_func: theta function

        Returns:
        Time grid and interest rate paths
        """
        dt = T / n_steps
        times = np.linspace(0, T, n_steps + 1)
        rates = np.zeros((n_paths, n_steps + 1))
        rates[:, 0] = r0

        for i in range(n_steps):
            t = times[i]
            theta_t = theta_func(t) if theta_func else 0.03

            # Euler-Markov method
            dW = np.random.normal(0, np.sqrt(dt), n_paths)
            dr = (theta_t - self.a * rates[:, i]) * dt + self.sigma * dW
            rates[:, i + 1] = rates[:, i] + dr

        return times, rates

    def calibrate_to_yield_curve(self, market_yields, maturities):
        """
        Calibrate model to market yield curve

        Parameters:
        market_yields: Market yields
        maturities: Corresponding maturities

        Returns:
        Calibrated parameters
        """
        def objective(params):
            a_new, sigma_new = params
            model_yields = []

            for T in maturities:
                # Calculate model yields using current parameters
                # This is simplified; actual implementation should use more complex calibration
                r0 = market_yields[0]  # Assume r0 is the shortest-term rate
                price = self.bond_price_analytical(r0, 0, T, a_new, sigma_new)
                yield_model = -np.log(price) / T if T > 0 else r0
                model_yields.append(yield_model)

            # Calculate error
            model_yields = np.array(model_yields)
            return np.sum((market_yields - model_yields)**2)

        # Optimize
        initial_guess = [self.a, self.sigma]
        bounds = [(0.01, 1.0), (0.001, 0.1)]
        result = minimize(objective, initial_guess, bounds=bounds, method='L-BFGS-B')

        if result.success:
            self.a, self.sigma = result.x
            return result.x
        else:
            print("Calibration failed")
            return None

    def bond_price_analytical(self, r, t, T, a=None, sigma=None):
        """Calculate bond price using analytical solution (simplified version)"""
        a = a or self.a
        sigma = sigma or self.sigma

        tau = T - t
        if tau <= 0:
            return 1.0

        B = (1 - np.exp(-a * tau)) / a if a != 0 else tau
        A = np.exp(-sigma**2 / (4 * a) * B**2 * (1 - np.exp(-2 * a * tau)))

        return A * np.exp(-B * r)

def generate_sample_yield_curve():
    """Generate sample yield curve"""
    maturities = np.array([0.25, 0.5, 1, 2, 3, 5, 7, 10, 15, 20, 30])

    # Nelson-Siegel model parameters
    beta0, beta1, beta2, tau = 0.04, -0.01, 0.02, 2.0

    yields = beta0 + beta1 * (1 - np.exp(-maturities/tau)) / (maturities/tau) + \
             beta2 * ((1 - np.exp(-maturities/tau)) / (maturities/tau) - np.exp(-maturities/tau))

    return maturities, yields

# Generate sample yield curve
maturities, market_yields = generate_sample_yield_curve()

print("Sample Yield Curve:")
yield_df = pd.DataFrame({
    'Maturity (years)': maturities,
    'Yield (%)': market_yields * 100
})
print(yield_df)

# Create Hull-White model
hw_model = HullWhiteModel(a=0.1, sigma=0.01)

print(f"\nInitial parameters:")
print(f"Mean reversion speed a: {hw_model.a}")
print(f"Volatility σ: {hw_model.sigma}")

# Calibrate model
print(f"\nCalibrating model to market yield curve...")
calibrated_params = hw_model.calibrate_to_yield_curve(market_yields, maturities)

if calibrated_params is not None:
    print(f"Calibrated parameters:")
    print(f"Mean reversion speed a: {calibrated_params[0]:.4f}")
    print(f"Volatility σ: {calibrated_params[1]:.4f}")

Example 2: Interest Rate Path Simulation and Visualization

# Set theta function (simplified as constant)
def theta_function(t):
    return 0.03 + 0.01 * np.sin(2 * np.pi * t)  # Periodically varying theta

hw_model.set_theta_function(theta_function)

# Simulate interest rate paths
r0 = market_yields[0]  # Initial rate
T = 10  # Simulate 10 years
n_steps = 1000
n_paths = 1000

print(f"\nSimulation parameters:")
print(f"Initial rate: {r0:.2%}")
print(f"Simulation period: {T} years")
print(f"Time steps: {n_steps}")
print(f"Number of paths: {n_paths}")

times, rate_paths = hw_model.simulate_paths(r0, T, n_steps, n_paths, theta_function)

# Visualize interest rate paths
fig, axes = plt.subplots(2, 2, figsize=(15, 12))

# Subplot 1: Sample interest rate paths
axes[0, 0].plot(times, rate_paths[:50].T, alpha=0.3, linewidth=0.5)
axes[0, 0].plot(times, np.mean(rate_paths, axis=0), 'r-', linewidth=2, label='Average path')
axes[0, 0].set_title('Hull-White Model Interest Rate Path Samples')
axes[0, 0].set_xlabel('Time (years)')
axes[0, 0].set_ylabel('Interest Rate')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# Subplot 2: Evolution of interest rate distribution
time_points = [1, 3, 5, 10]
colors = ['blue', 'green', 'orange', 'red']

for i, (t_point, color) in enumerate(zip(time_points, colors)):
    t_idx = int(t_point * n_steps / T)
    rates_at_t = rate_paths[:, t_idx]
    axes[0, 1].hist(rates_at_t, bins=50, alpha=0.6, density=True,
                   color=color, label=f't={t_point} years')

axes[0, 1].set_title('Interest Rate Distribution at Different Times')
axes[0, 1].set_xlabel('Interest Rate')
axes[0, 1].set_ylabel('Density')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# Subplot 3: Interest rate statistics over time
mean_rates = np.mean(rate_paths, axis=0)
std_rates = np.std(rate_paths, axis=0)
percentile_5 = np.percentile(rate_paths, 5, axis=0)
percentile_95 = np.percentile(rate_paths, 95, axis=0)

axes[1, 0].plot(times, mean_rates, 'b-', linewidth=2, label='Mean')
axes[1, 0].fill_between(times, percentile_5, percentile_95, alpha=0.3,
                       color='blue', label='90% Confidence Interval')
axes[1, 0].plot(times, mean_rates + 2*std_rates, 'r--', alpha=0.7, label='Mean±2σ')
axes[1, 0].plot(times, mean_rates - 2*std_rates, 'r--', alpha=0.7)
axes[1, 0].set_title('Evolution of Interest Rate Statistics')
axes[1, 0].set_xlabel('Time (years)')
axes[1, 0].set_ylabel('Interest Rate')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)

# Subplot 4: Yield curve comparison
model_maturities = np.linspace(0.25, 10, 20)
model_yields = []

for T in model_maturities:
    # Calculate model-predicted yields
    price = hw_model.bond_price(r0, 0, T)
    model_yield = -np.log(price) / T if T > 0 else r0
    model_yields.append(model_yield)

axes[1, 1].plot(maturities, market_yields * 100, 'bo-', label='Market Yields', markersize=6)
axes[1, 1].plot(model_maturities, np.array(model_yields) * 100, 'r-',
               label='Hull-White Model', linewidth=2)
axes[1, 1].set_title('Yield Curve Comparison')
axes[1, 1].set_xlabel('Maturity (years)')
axes[1, 1].set_ylabel('Yield (%)')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Calculate interest rate statistics
print(f"\nInterest Rate Path Statistics:")
print(f"Initial rate: {r0:.2%}")
print(f"Average rate after 10 years: {mean_rates[-1]:.2%}")
print(f"Rate standard deviation after 10 years: {std_rates[-1]:.2%}")
print(f"Rate range after 10 years: [{percentile_5[-1]:.2%}, {percentile_95[-1]:.2%}]")

Example 3: Trinomial Tree Interest Rate Model

class TrinomialTree:
    """
    Trinomial tree implementation for Hull-White model
    """

    def __init__(self, hw_model, T, n_steps):
        """
        Initialize trinomial tree

        Parameters:
        hw_model: Hull-White model
        T: Total time
        n_steps: Number of time steps
        """
        self.hw_model = hw_model
        self.T = T
        self.n_steps = n_steps
        self.dt = T / n_steps
        self.tree = {}
        self.probabilities = {}

    def build_tree(self, r0):
        """Build interest rate tree"""
        # Time step
        dt = self.dt
        a = self.hw_model.a
        sigma = self.hw_model.sigma

        # Space step
        dr = sigma * np.sqrt(3 * dt)

        # Initialize tree
        self.tree[0] = {0: r0}
        self.probabilities[0] = {0: 1.0}

        # Build each level of the tree
        for i in range(self.n_steps):
            self.tree[i + 1] = {}
            self.probabilities[i + 1] = {}

            for j in self.tree[i]:
                r = self.tree[i][j]
                prob = self.probabilities[i][j]

                # Calculate mean reversion drift
                theta_t = 0.03  # Simplified as constant
                drift = (theta_t - a * r) * dt

                # Three branch rate values
                r_up = r + drift + dr
                r_mid = r + drift
                r_down = r + drift - dr

                # Calculate branch probabilities
                # This is a simplified version, actual should match moments more precisely
                p_up = 1/6 + (r * a * dt - drift)**2 / (2 * dr**2) + (r * a * dt - drift) / (2 * dr)
                p_down = 1/6 + (r * a * dt - drift)**2 / (2 * dr**2) - (r * a * dt - drift) / (2 * dr)
                p_mid = 2/3 - (r * a * dt - drift)**2 / (dr**2)

                # Ensure probabilities are positive and sum to 1
                p_up = max(0, min(1, p_up))
                p_down = max(0, min(1, p_down))
                p_mid = max(0, min(1, 1 - p_up - p_down))

                # Normalize
                total_prob = p_up + p_mid + p_down
                if total_prob > 0:
                    p_up /= total_prob
                    p_mid /= total_prob
                    p_down /= total_prob

                # Add to tree
                for k, (r_new, p_new) in enumerate([(r_up, p_up), (r_mid, p_mid), (r_down, p_down)]):
                    node_id = 3 * j + k

                    if node_id in self.tree[i + 1]:
                        self.probabilities[i + 1][node_id] += prob * p_new
                    else:
                        self.tree[i + 1][node_id] = r_new
                        self.probabilities[i + 1][node_id] = prob * p_new

    def price_zero_coupon_bond(self, face_value=1.0):
        """
        Price zero-coupon bond using tree

        Parameters:
        face_value: Face value

        Returns:
        Bond price
        """
        # Backward induction from maturity
        bond_values = {}

        # Maturity value
        for j in self.tree[self.n_steps]:
            bond_values[(self.n_steps, j)] = face_value

        # Backward recursion
        for i in range(self.n_steps - 1, -1, -1):
            for j in self.tree[i]:
                r = self.tree[i][j]

                # Calculate expected value
                expected_value = 0
                for k in range(3):  # Three branches
                    next_node = 3 * j + k
                    if next_node in self.tree[i + 1]:
                        # Branch probability (simplified calculation)
                        p = 1/3  # Simplified as equal probability
                        expected_value += p * bond_values.get((i + 1, next_node), 0)

                # Discount
                bond_values[(i, j)] = expected_value * np.exp(-r * self.dt)

        return bond_values[(0, 0)]

    def price_bond_option(self, strike, option_type='call', bond_maturity=None):
        """
        Price bond option

        Parameters:
        strike: Strike price
        option_type: Option type ('call' or 'put')
        bond_maturity: Underlying bond maturity

        Returns:
        Option price
        """
        if bond_maturity is None:
            bond_maturity = self.T

        # Option values
        option_values = {}

        # Maturity value
        for j in self.tree[self.n_steps]:
            r = self.tree[self.n_steps][j]

            # Calculate underlying bond value
            bond_price = self.hw_model.bond_price(r, self.T, bond_maturity)

            # Option intrinsic value
            if option_type == 'call':
                payoff = max(0, bond_price - strike)
            else:
                payoff = max(0, strike - bond_price)

            option_values[(self.n_steps, j)] = payoff

        # Backward recursion
        for i in range(self.n_steps - 1, -1, -1):
            for j in self.tree[i]:
                r = self.tree[i][j]

                # Calculate expected value
                expected_value = 0
                for k in range(3):
                    next_node = 3 * j + k
                    if next_node in self.tree[i + 1]:
                        p = 1/3  # Simplified as equal probability
                        expected_value += p * option_values.get((i + 1, next_node), 0)

                # Discount
                option_values[(i, j)] = expected_value * np.exp(-r * self.dt)

        return option_values[(0, 0)]

# Create trinomial tree
tree = TrinomialTree(hw_model, T=2, n_steps=20)
tree.build_tree(r0)

# Price zero-coupon bond
zcb_price = tree.price_zero_coupon_bond()
analytical_price = hw_model.bond_price(r0, 0, 2)

print(f"\nZero-Coupon Bond Pricing Comparison (2-year):")
print(f"Trinomial tree price: {zcb_price:.6f}")
print(f"Analytical solution price: {analytical_price:.6f}")
print(f"Price difference: {abs(zcb_price - analytical_price):.6f}")

# Price bond option
call_price = tree.price_bond_option(strike=0.9, option_type='call', bond_maturity=3)
put_price = tree.price_bond_option(strike=0.9, option_type='put', bond_maturity=3)

print(f"\nBond Option Pricing (strike=0.9, underlying 3-year bond):")
print(f"Call option price: {call_price:.6f}")
print(f"Put option price: {put_price:.6f}")

# Verify put-call parity
current_bond_price = hw_model.bond_price(r0, 0, 3)
strike = 0.9
discount_factor = hw_model.bond_price(r0, 0, 2)

put_call_parity_diff = call_price - put_price - (current_bond_price - strike * discount_factor)
print(f"Put-call parity verification error: {put_call_parity_diff:.8f}")

Example 4: Interest Rate Risk Management

class InterestRateRiskManager:
    """
    Interest rate risk management tools
    """

    def __init__(self, hw_model):
        self.hw_model = hw_model

    def calculate_duration(self, bond_cashflows, bond_times, current_rate):
        """
        Calculate modified duration

        Parameters:
        bond_cashflows: Bond cash flows
        bond_times: Cash flow times
        current_rate: Current interest rate

        Returns:
        Modified duration
        """
        # Current bond price
        current_price = sum(cf * self.hw_model.bond_price(current_rate, 0, t)
                           for cf, t in zip(bond_cashflows, bond_times))

        # Calculate price change with small rate change
        dr = 0.0001  # 1 basis point
        price_up = sum(cf * self.hw_model.bond_price(current_rate + dr, 0, t)
                      for cf, t in zip(bond_cashflows, bond_times))
        price_down = sum(cf * self.hw_model.bond_price(current_rate - dr, 0, t)
                        for cf, t in zip(bond_cashflows, bond_times))

        # Modified duration
        modified_duration = -(price_up - price_down) / (2 * dr * current_price)
        return modified_duration, current_price

    def calculate_convexity(self, bond_cashflows, bond_times, current_rate):
        """Calculate convexity"""
        dr = 0.0001
        current_price = sum(cf * self.hw_model.bond_price(current_rate, 0, t)
                           for cf, t in zip(bond_cashflows, bond_times))

        price_up = sum(cf * self.hw_model.bond_price(current_rate + dr, 0, t)
                      for cf, t in zip(bond_cashflows, bond_times))
        price_down = sum(cf * self.hw_model.bond_price(current_rate - dr, 0, t)
                        for cf, t in zip(bond_cashflows, bond_times))

        convexity = (price_up + price_down - 2 * current_price) / (dr**2 * current_price)
        return convexity

    def portfolio_risk_analysis(self, portfolio):
        """
        Portfolio interest rate risk analysis

        Parameters:
        portfolio: Portfolio information [{'weight': w, 'cashflows': cf, 'times': t}, ...]

        Returns:
        Risk metrics
        """
        current_rate = 0.03  # Assume current rate

        portfolio_duration = 0
        portfolio_convexity = 0
        portfolio_value = 0

        bond_details = []

        for bond in portfolio:
            weight = bond['weight']
            cashflows = bond['cashflows']
            times = bond['times']

            # Calculate individual bond risk metrics
            duration, price = self.calculate_duration(cashflows, times, current_rate)
            convexity = self.calculate_convexity(cashflows, times, current_rate)

            # Calculate weight
            bond_value = weight * price
            portfolio_value += bond_value

            bond_details.append({
                'weight': weight,
                'price': price,
                'duration': duration,
                'convexity': convexity,
                'value': bond_value
            })

        # Calculate portfolio weighted average
        for bond in bond_details:
            value_weight = bond['value'] / portfolio_value
            portfolio_duration += value_weight * bond['duration']
            portfolio_convexity += value_weight * bond['convexity']

        return {
            'portfolio_duration': portfolio_duration,
            'portfolio_convexity': portfolio_convexity,
            'portfolio_value': portfolio_value,
            'bond_details': bond_details
        }

    def hedge_ratio(self, hedged_instrument, hedging_instrument, current_rate):
        """Calculate hedge ratio"""
        # Calculate duration of hedged instrument
        duration_hedged, price_hedged = self.calculate_duration(
            hedged_instrument['cashflows'],
            hedged_instrument['times'],
            current_rate
        )

        # Calculate duration of hedging instrument
        duration_hedging, price_hedging = self.calculate_duration(
            hedging_instrument['cashflows'],
            hedging_instrument['times'],
            current_rate
        )

        # Hedge ratio
        hedge_ratio = (duration_hedged * price_hedged) / (duration_hedging * price_hedging)

        return hedge_ratio

# Create risk manager
risk_manager = InterestRateRiskManager(hw_model)

# Define sample portfolio
portfolio = [
    {
        'weight': 1000000,  # 1 million face value
        'cashflows': [50000, 50000, 50000, 1050000],  # 5% coupon, 3-year
        'times': [1, 2, 3, 3]
    },
    {
        'weight': 2000000,  # 2 million face value
        'cashflows': [40000, 40000, 40000, 40000, 40000, 1040000],  # 4% coupon, 5-year
        'times': [1, 2, 3, 4, 5, 5]
    },
    {
        'weight': 1500000,  # 1.5 million face value zero-coupon bond
        'cashflows': [1500000],  # 10-year zero-coupon bond
        'times': [10]
    }
]

# Perform risk analysis
risk_analysis = risk_manager.portfolio_risk_analysis(portfolio)

print(f"\nPortfolio Interest Rate Risk Analysis:")
print("=" * 50)
print(f"Portfolio total value: {risk_analysis['portfolio_value']:,.0f}")
print(f"Portfolio modified duration: {risk_analysis['portfolio_duration']:.2f}")
print(f"Portfolio convexity: {risk_analysis['portfolio_convexity']:.2f}")

print(f"\nDetailed bond information:")
for i, bond in enumerate(risk_analysis['bond_details']):
    print(f"Bond {i+1}:")
    print(f"  Face value: {bond['weight']:,.0f}")
    print(f"  Current value: {bond['value']:,.0f}")
    print(f"  Modified duration: {bond['duration']:.2f}")
    print(f"  Convexity: {bond['convexity']:.2f}")

# Interest rate shock analysis
rate_shocks = [-0.02, -0.01, -0.005, 0, 0.005, 0.01, 0.02]  # Rate changes
portfolio_values = []

current_rate = 0.03
base_value = risk_analysis['portfolio_value']

for shock in rate_shocks:
    new_rate = current_rate + shock
    shocked_value = 0

    for bond in portfolio:
        bond_value = sum(cf * hw_model.bond_price(new_rate, 0, t)
                        for cf, t in zip(bond['cashflows'], bond['times']))
        shocked_value += bond['weight'] / bond['cashflows'][-1] * bond_value  # Adjust weight

    portfolio_values.append(shocked_value)

# Visualize interest rate risk
plt.figure(figsize=(12, 8))

# Subplot 1: Interest rate shock analysis
plt.subplot(2, 2, 1)
pnl_values = [(val - base_value) / base_value * 100 for val in portfolio_values]
plt.plot(np.array(rate_shocks) * 100, pnl_values, 'bo-', linewidth=2, markersize=6)
plt.title('Portfolio Interest Rate Sensitivity Analysis')
plt.xlabel('Rate Change (basis points)')
plt.ylabel('Portfolio Value Change (%)')
plt.grid(True, alpha=0.3)

# Duration and convexity approximation formula
duration_approx = []
for shock in rate_shocks:
    approx_change = (-risk_analysis['portfolio_duration'] * shock +
                    0.5 * risk_analysis['portfolio_convexity'] * shock**2) * 100
    duration_approx.append(approx_change)

plt.plot(np.array(rate_shocks) * 100, duration_approx, 'r--',
         linewidth=2, label='Duration+Convexity Approximation')
plt.legend()

# Subplot 2: Bond duration distribution
plt.subplot(2, 2, 2)
durations = [bond['duration'] for bond in risk_analysis['bond_details']]
weights = [bond['value']/risk_analysis['portfolio_value'] for bond in risk_analysis['bond_details']]
labels = [f'Bond {i+1}' for i in range(len(durations))]

plt.bar(labels, durations, color=['blue', 'green', 'orange'], alpha=0.7)
plt.title('Modified Duration by Bond')
plt.ylabel('Modified Duration')
plt.grid(True, alpha=0.3, axis='y')

# Add weight information to bar chart
for i, (dur, weight) in enumerate(zip(durations, weights)):
    plt.text(i, dur + 0.1, f'Weight: {weight:.1%}', ha='center', va='bottom')

# Subplot 3: Yield curve shift scenarios
plt.subplot(2, 2, 3)
scenarios = {
    'Parallel shift +100bp': 0.01,
    'Parallel shift -100bp': -0.01,
    'Steepening': 0,  # Simplified
    'Flattening': 0   # Simplified
}

scenario_impacts = []
for scenario, shock in scenarios.items():
    if scenario in ['Parallel shift +100bp', 'Parallel shift -100bp']:
        impact = (-risk_analysis['portfolio_duration'] * shock +
                 0.5 * risk_analysis['portfolio_convexity'] * shock**2) * 100
    else:
        impact = 0  # Simplified, actual requires more complex calculation
    scenario_impacts.append(impact)

colors = ['red' if x < 0 else 'green' for x in scenario_impacts]
plt.bar(range(len(scenarios)), scenario_impacts, color=colors, alpha=0.7)
plt.title('Portfolio Impact Under Different Rate Scenarios')
plt.xticks(range(len(scenarios)), scenarios.keys(), rotation=45)
plt.ylabel('Value Change (%)')
plt.grid(True, alpha=0.3, axis='y')

# Subplot 4: Hedging strategy analysis
plt.subplot(2, 2, 4)
hedge_instruments = ['2-year Treasury', '5-year Treasury', '10-year Treasury', '30-year Treasury']
hedge_durations = [1.9, 4.5, 8.5, 18.0]  # Assumed durations
hedge_effectiveness = []

for dur in hedge_durations:
    # Simplified hedge effectiveness calculation
    effectiveness = 1 - abs(risk_analysis['portfolio_duration'] - dur) / risk_analysis['portfolio_duration']
    effectiveness = max(0, effectiveness)
    hedge_effectiveness.append(effectiveness)

plt.bar(hedge_instruments, hedge_effectiveness, color='purple', alpha=0.7)
plt.title('Effectiveness of Different Hedging Instruments')
plt.ylabel('Hedge Effectiveness')
plt.xticks(rotation=45)
plt.grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.show()

print(f"\nInterest Rate Risk Management Recommendations:")
print(f"1. Portfolio duration is {risk_analysis['portfolio_duration']:.2f}, relatively sensitive to rate changes")
print(f"2. A 100bp rate increase will lead to approximately {risk_analysis['portfolio_duration']:.1f}% decline in portfolio value")
print(f"3. Recommend hedging with instruments of similar duration")
print(f"4. Positive convexity is more favorable for rate declines")

Theoretical Analysis

Mathematical Properties of Hull-White Model

Mean Reversion Property: The solution of the Hull-White model is: $r_t = r_0 e^{-at} + \int_0^t \theta(s) e^{-a(t-s)} ds + \sigma \int_0^t e^{-a(t-s)} dW_s$

Variance Structure: $\text{Var}(r_t) = \frac{\sigma^2}{2a}(1-e^{-2at})$

Long-term Variance: $\lim_{t \to \infty} \text{Var}(r_t) = \frac{\sigma^2}{2a}$

Affine Property of Bond Pricing

The Hull-White model belongs to the affine term structure model class, with bond prices having exponential-affine form: $P(t,T) = \exp(A(t,T) - B(t,T) r_t)$

This enables:

1. Analytical solutions for bond prices
2. Relatively simple interest rate derivative pricing
3. More efficient calibration process

Convergence of Interest Rate Trees

Convergence conditions for trinomial tree method:

1. Stability condition: $\Delta t \leq \frac{1}{a}$
2. Non-negative probability: All branch probabilities $\geq 0$
3. Moment matching: Correctly match first two moments

Mathematical Formula Summary

1. Hull-White SDE: $dr_t = (\theta(t) - ar_t)dt + \sigma dW_t$

2. Bond Price Formula: $P(t,T) = A(t,T)e^{-B(t,T)r_t}$

3. B Function: $B(t,T) = \frac{1-e^{-a(T-t)}}{a}$

4. Modified Duration: $D_{mod} = -\frac{1}{P}\frac{\partial P}{\partial r}$

5. Convexity: $C = \frac{1}{P}\frac{\partial^2 P}{\partial r^2}$

6. Hedge Ratio: $h = \frac{D_H \cdot P_H}{D_I \cdot P_I}$