|
| 1 | +import numpy as np |
| 2 | +import cvxpy as cp |
| 3 | +import matplotlib.pyplot as plt |
| 4 | + |
| 5 | +def nelson_siegel(tau, beta0, beta1, beta2, lambda_param): |
| 6 | + """ |
| 7 | + Calculate Nelson-Siegel yield curve values. |
| 8 | + |
| 9 | + Parameters: |
| 10 | + tau: maturity times |
| 11 | + beta0, beta1, beta2: model parameters |
| 12 | + lambda_param: decay parameter |
| 13 | + """ |
| 14 | + exp_term = cp.exp(-tau / lambda_param) |
| 15 | + factor1 = (1 - exp_term) / (tau / lambda_param) |
| 16 | + factor2 = factor1 - exp_term |
| 17 | + |
| 18 | + return beta0 + beta1 * factor1 + beta2 * factor2 |
| 19 | + |
| 20 | +def fit_nelson_siegel_cvxpy(maturities, yields, lambda_init=1.0, |
| 21 | + beta_bounds=(-10, 10), lambda_bounds=(0.1, 10)): |
| 22 | + """ |
| 23 | + Fit Nelson-Siegel model using CVXPY with NLP solver. |
| 24 | + |
| 25 | + Parameters: |
| 26 | + maturities: array of maturity times |
| 27 | + yields: observed yields |
| 28 | + lambda_init: initial value for lambda parameter |
| 29 | + beta_bounds: bounds for beta parameters |
| 30 | + lambda_bounds: bounds for lambda parameter |
| 31 | + """ |
| 32 | + n = len(maturities) |
| 33 | + |
| 34 | + # Define variables |
| 35 | + beta0 = cp.Variable() |
| 36 | + beta1 = cp.Variable() |
| 37 | + beta2 = cp.Variable() |
| 38 | + lambda_param = cp.Variable(pos=True) |
| 39 | + |
| 40 | + # Calculate model predictions using Nelson-Siegel formula |
| 41 | + # Note: We need to handle the division by zero case when tau approaches 0 |
| 42 | + predictions = [] |
| 43 | + for tau in maturities: |
| 44 | + if tau < 1e-6: # Handle near-zero maturity |
| 45 | + pred = beta0 + beta1 |
| 46 | + else: |
| 47 | + exp_term = cp.exp(-tau / lambda_param) |
| 48 | + factor1 = (1 - exp_term) / (tau / lambda_param) |
| 49 | + factor2 = factor1 - exp_term |
| 50 | + pred = beta0 + beta1 * factor1 + beta2 * factor2 |
| 51 | + predictions.append(pred) |
| 52 | + |
| 53 | + predictions = cp.vstack(predictions) |
| 54 | + |
| 55 | + # Define objective: minimize sum of squared errors |
| 56 | + objective = cp.Minimize(cp.sum_squares(predictions - yields.reshape(-1, 1))) |
| 57 | + |
| 58 | + # Define constraints |
| 59 | + constraints = [ |
| 60 | + beta0 >= beta_bounds[0], beta0 <= beta_bounds[1], |
| 61 | + beta1 >= beta_bounds[0], beta1 <= beta_bounds[1], |
| 62 | + beta2 >= beta_bounds[0], beta2 <= beta_bounds[1], |
| 63 | + lambda_param >= lambda_bounds[0], |
| 64 | + lambda_param <= lambda_bounds[1] |
| 65 | + ] |
| 66 | + |
| 67 | + # Set initial values (important for NLP solvers) |
| 68 | + beta0.value = np.mean(yields) |
| 69 | + lambda_param.value = lambda_init |
| 70 | + |
| 71 | + # Create and solve problem |
| 72 | + problem = cp.Problem(objective, constraints) |
| 73 | + |
| 74 | + # Solve using NLP solver (e.g., IPOPT through CVXPY) |
| 75 | + # Note: You need to have an NLP solver installed |
| 76 | + problem.solve(solver=cp.IPOPT, verbose=True, nlp=True) |
| 77 | + |
| 78 | + if problem.status not in ["infeasible", "unbounded"]: |
| 79 | + return { |
| 80 | + 'beta0': beta0.value, |
| 81 | + 'beta1': beta1.value, |
| 82 | + 'beta2': beta2.value, |
| 83 | + 'lambda': lambda_param.value, |
| 84 | + 'objective': problem.value, |
| 85 | + 'status': problem.status |
| 86 | + } |
| 87 | + else: |
| 88 | + raise ValueError(f"Optimization failed with status: {problem.status}") |
| 89 | + |
| 90 | +def plot_results(maturities, yields, fitted_params, title="Nelson-Siegel Fit"): |
| 91 | + """Plot observed vs fitted yield curve.""" |
| 92 | + |
| 93 | + # Generate smooth curve for plotting |
| 94 | + tau_smooth = np.linspace(min(maturities), max(maturities), 100) |
| 95 | + |
| 96 | + # Calculate fitted values |
| 97 | + beta0 = fitted_params['beta0'] |
| 98 | + beta1 = fitted_params['beta1'] |
| 99 | + beta2 = fitted_params['beta2'] |
| 100 | + lambda_val = fitted_params['lambda'] |
| 101 | + |
| 102 | + # Nelson-Siegel formula for smooth curve |
| 103 | + exp_term = np.exp(-tau_smooth / lambda_val) |
| 104 | + factor1 = np.where(tau_smooth < 1e-6, 1, (1 - exp_term) / (tau_smooth / lambda_val)) |
| 105 | + factor2 = factor1 - exp_term |
| 106 | + y_fitted_smooth = beta0 + beta1 * factor1 + beta2 * factor2 |
| 107 | + |
| 108 | + # Calculate fitted values at observed points |
| 109 | + exp_term_obs = np.exp(-maturities / lambda_val) |
| 110 | + factor1_obs = np.where(maturities < 1e-6, 1, |
| 111 | + (1 - exp_term_obs) / (maturities / lambda_val)) |
| 112 | + factor2_obs = factor1_obs - exp_term_obs |
| 113 | + y_fitted_obs = beta0 + beta1 * factor1_obs + beta2 * factor2_obs |
| 114 | + |
| 115 | + plt.figure(figsize=(10, 6)) |
| 116 | + plt.scatter(maturities, yields, color='blue', label='Observed', s=50) |
| 117 | + plt.plot(tau_smooth, y_fitted_smooth, 'r-', label='Fitted NS Curve', linewidth=2) |
| 118 | + plt.scatter(maturities, y_fitted_obs, color='red', alpha=0.5, s=30) |
| 119 | + |
| 120 | + plt.xlabel('Maturity (years)') |
| 121 | + plt.ylabel('Yield (%)') |
| 122 | + plt.title(title) |
| 123 | + plt.legend() |
| 124 | + plt.grid(True, alpha=0.3) |
| 125 | + |
| 126 | + # Add parameter text |
| 127 | + param_text = f'β₀={beta0:.4f}, β₁={beta1:.4f}, β₂={beta2:.4f}, λ={lambda_val:.4f}' |
| 128 | + plt.text(0.02, 0.98, param_text, transform=plt.gca().transAxes, |
| 129 | + verticalalignment='top', fontsize=10, |
| 130 | + bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5)) |
| 131 | + |
| 132 | + plt.show() |
| 133 | + |
| 134 | +# Example usage |
| 135 | +if __name__ == "__main__": |
| 136 | + # Generate sample yield curve data |
| 137 | + np.random.seed(42) |
| 138 | + |
| 139 | + # Maturities in years |
| 140 | + maturities = np.array([0.25, 0.5, 1, 2, 3, 5, 7, 10, 15, 20, 30]) |
| 141 | + |
| 142 | + # True parameters for generating synthetic data |
| 143 | + true_beta0 = 5.0 |
| 144 | + true_beta1 = -2.0 |
| 145 | + true_beta2 = 1.5 |
| 146 | + true_lambda = 2.0 |
| 147 | + |
| 148 | + # Generate synthetic yields with some noise |
| 149 | + exp_term = np.exp(-maturities / true_lambda) |
| 150 | + factor1 = (1 - exp_term) / (maturities / true_lambda) |
| 151 | + factor2 = factor1 - exp_term |
| 152 | + true_yields = true_beta0 + true_beta1 * factor1 + true_beta2 * factor2 |
| 153 | + |
| 154 | + # Add noise |
| 155 | + noise = np.random.normal(0, 0.1, len(maturities)) |
| 156 | + observed_yields = true_yields + noise |
| 157 | + |
| 158 | + print("Fitting Nelson-Siegel model using CVXPY...") |
| 159 | + print("-" * 50) |
| 160 | + |
| 161 | + # Fit the model |
| 162 | + fitted_params = fit_nelson_siegel_cvxpy( |
| 163 | + maturities, |
| 164 | + observed_yields, |
| 165 | + lambda_init=2.0, |
| 166 | + beta_bounds=(-10, 10), |
| 167 | + lambda_bounds=(0.5, 5.0) |
| 168 | + ) |
| 169 | + |
| 170 | + print(f"Optimization Status: {fitted_params['status']}") |
| 171 | + print(f"Objective Value (SSE): {fitted_params['objective']:.6f}") |
| 172 | + print("\nFitted Parameters:") |
| 173 | + print(f" β₀ (level): {fitted_params['beta0']:.4f} (true: {true_beta0:.4f})") |
| 174 | + print(f" β₁ (slope): {fitted_params['beta1']:.4f} (true: {true_beta1:.4f})") |
| 175 | + print(f" β₂ (curvature): {fitted_params['beta2']:.4f} (true: {true_beta2:.4f})") |
| 176 | + print(f" λ (decay): {fitted_params['lambda']:.4f} (true: {true_lambda:.4f})") |
| 177 | + |
| 178 | + # Calculate R-squared |
| 179 | + exp_term_fit = np.exp(-maturities / fitted_params['lambda']) |
| 180 | + factor1_fit = (1 - exp_term_fit) / (maturities / fitted_params['lambda']) |
| 181 | + factor2_fit = factor1_fit - exp_term_fit |
| 182 | + y_fitted = (fitted_params['beta0'] + |
| 183 | + fitted_params['beta1'] * factor1_fit + |
| 184 | + fitted_params['beta2'] * factor2_fit) |
| 185 | + |
| 186 | + ss_res = np.sum((observed_yields - y_fitted) ** 2) |
| 187 | + ss_tot = np.sum((observed_yields - np.mean(observed_yields)) ** 2) |
| 188 | + r_squared = 1 - (ss_res / ss_tot) |
| 189 | + print(f"\nR-squared: {r_squared:.4f}") |
| 190 | + |
| 191 | + # Plot results |
| 192 | + plot_results(maturities, observed_yields, fitted_params, |
| 193 | + "Nelson-Siegel Model Fit with CVXPY") |
| 194 | + |
| 195 | + # Alternative: Using different lambda initialization |
| 196 | + print("\n" + "=" * 50) |
| 197 | + print("Testing sensitivity to lambda initialization...") |
| 198 | + |
| 199 | + for lambda_init in [0.5, 1.0, 3.0, 4.5]: |
| 200 | + try: |
| 201 | + params = fit_nelson_siegel_cvxpy( |
| 202 | + maturities, observed_yields, |
| 203 | + lambda_init=lambda_init |
| 204 | + ) |
| 205 | + print(f"λ_init={lambda_init:.1f}: λ_final={params['lambda']:.4f}, " |
| 206 | + f"SSE={params['objective']:.4f}") |
| 207 | + except Exception as e: |
| 208 | + print(f"λ_init={lambda_init:.1f}: Failed - {e}") |
0 commit comments