Calibrating an Arbitrage-Free Nelson–Siegel Yield-Curve Model

# Arbitrage Free Nelson Siegel Model Python


# Arbitrage Free Nelson Siegel Model Python












I am currently trying to implement the Arbitrage Free Nelson Siegel (AFNS) model in Python. However, I am encountering the problem that my results do not match the current yield curve at all. Is there anyone here who has experience with the AFNS and/or can help with the Kalman filter for determining the variabes?

I follow the work of Christensen, Diebold and Rudebusch (https://www.frbsf.org/wp-content/uploads/wp07-20bk.pdf)

Many thanks in advance! Enclosed is my code

```
import numpy as np
from scipy.linalg import expm
from scipy.optimize import minimize
from scipy.special import expit  # Sigmoid function for the transformation
from scipy.stats import norm
import pandas as pd

# DNS factor loading matrix
def NS_B(lambda_, tau):
    col1 = np.ones_like(tau)
    col2 = (1 - np.exp(-lambda_ * tau)) / (lambda_ * tau)
    col3 = col2 - np.exp(-lambda_ * tau)
    return np.column_stack((col1, col2, col3))

# yield adjustment term in AFNS
def AFNS_C(sigma, lambda_, tau):
    s = sigma
    A = s[0, 0] ** 2 + s[0, 1] ** 2 + s[0, 2] ** 2
    B = s[1, 0] ** 2 + s[1, 1] ** 2 + s[1, 2] ** 2
    C = s[2, 0] ** 2 + s[2, 1] ** 2 + s[2, 2] ** 2
    D = s[0, 0] * s[1, 0] + s[0, 1] * s[1, 1] + s[0, 2] * s[1, 2]
    E = s[0, 0] * s[2, 0] + s[0, 1] * s[2, 1] + s[0, 2] * s[2, 2]
    F = s[1, 0] * s[2, 0] + s[1, 1] * s[2, 1] + s[1, 2] * s[2, 2]

    t = tau
    la = lambda_
    r1 = -A * t ** 2 / 6
    r2 = -B * (1 / (2 * la ** 2) - (1 - np.exp(-la * t)) / (la ** 3 * t) + (1 - np.exp(-2 * la * t)) / (
                4 * la ** 3 * t))
    r3 = -C * (1 / (2 * la ** 2) + np.exp(-la * t) / (la ** 2) - t * np.exp(-2 * la * t) / (4 * la) -
               3 * np.exp(-2 * la * t) / (4 * la ** 2) - 2 * (1 - np.exp(-la * t)) / (la ** 3 * t) +
               5 * (1 - np.exp(-2 * la * t)) / (8 * la ** 3 * t))
    r4 = -D * (t / (2 * la) + np.exp(-la * t) / (la ** 2) - (1 - np.exp(-la * t)) / (la ** 3 * t))
    r5 = -E * (3 * np.exp(-la * t) / (la ** 2) + t / (2 * la) + t * np.exp(-la * t) / la -
               3 * (1 - np.exp(-la * t)) / (la ** 3 * t))
    r6 = -F * (1 / (la ** 2) + np.exp(-la * t) / (la ** 2) - np.exp(-2 * la * t) / (2 * la ** 2) -
               3 * (1 - np.exp(-la * t)) / (la ** 3 * t) + 3 * (1 - np.exp(-2 * la * t)) / (4 * la ** 3 * t))
    return r1 + r2 + r3 + r4 + r5 + r6

# Parameter restrictions
def transform_parameters(b):
    bb = b.copy()
    bb[0] = 1 / (1 + np.exp(b[0]))  # kappa11
    bb[12] = b[12] ** 2  # lambda
    bb[13:] = b[13:] ** 2  # measurement error
    return bb

# Log-likelihood function
def log_likelihood(para_un, m_spot, v_mat, dt, nf, nobs):
    para = transform_parameters(para_un)

    kappa = np.diag(para[:3])
    sigma = np.array([[para[3], 0, 0], [para[5], para[6], 0], [para[7], para[8], para[9]]])
    theta = para[10:13]
    lambda_ = para[13]
    H = np.diag(para[14:14 + len(v_mat)]) 
    B = NS_B(lambda_, v_mat)
    C = AFNS_C(sigma, lambda_, v_mat)

    # Eigen decomposition for kappa
    evals, evecs = np.linalg.eig(kappa)
    Smat = np.linalg.inv(evecs) @ sigma @ sigma.T @ np.linalg.inv(evecs.T)

    Vdt = np.zeros((nf, nf))
    Vinf = np.zeros((nf, nf))
    for i in range(nf):
        for j in range(nf):
            Vdt[i, j] = Smat[i, j] * (1 - np.exp(-(evals[i] + evals[j]) * dt)) / (evals[i] + evals[j])
            Vinf[i, j] = Smat[i, j] / (evals[i] + evals[j])

    Q = evecs @ Vdt @ evecs.T
    Q0 = evecs @ Vinf @ evecs.T

    # Initialization
    prevX = theta
    prevV = Q0
    Phi1 = expm(-kappa * dt)
    Phi0 = (np.eye(nf) - Phi1) @ theta
    log_likelihood_val = 0

    nobs = m_spot.shape[0]  # Number of observations

    for i in range(nobs):
        Xhat = Phi0 + Phi1 @ prevX
        Vhat = Phi1 @ prevV @ Phi1.T + Q

        y = m_spot[i, :]  # Observed yield
        y_implied = B @ Xhat + C  # Implied yield from the model
        er = y - y_implied  # Prediction error

        # Kalman gain
        ev = B @ Vhat @ B.T + H
        KG = Vhat @ B.T @ np.linalg.inv(ev)

        prevX = Xhat + KG @ er
        prevV = Vhat - KG @ B @ Vhat

        # Update log-likelihood
        log_likelihood_val += -0.5 * len(y) * np.log(2 * np.pi) - 0.5 * np.log(
            np.linalg.det(ev)) - 0.5 * er @ np.linalg.inv(ev) @ er

    return -log_likelihood_val

# Set the correct path to your Excel file
file_path = r'C:....'

# Read the spot rate data from Excel
df = pd.read_excel(file_path)

# Convert the DataFrame to a numpy matrix and remove the first column (dates)
m_spot = df.iloc[:, 1:].to_numpy() / 10000  # Assuming that the yields are in basis points

# Define the maturities based on the column names, convert to years
v_mat = np.array([float(mat.replace('Y', '')) if 'Y' in mat else float(mat.replace('M', ''))/12
                  for mat in df.columns[1:]])

nobs, nmat = m_spot.shape

# Define the number of observations and number of maturities
nobs, nmat = m_spot.shape

# Set the dt and nf based on the frequency of observations and number of factors
dt = 1/52  # weekly data
nf = 3  # Number of factors

# Initial guess for unconstrained parameters (para_un)
# Make sure that the number of initial parameters matches nmat + nf + 1 (for lambda)
init_para_un = np.array([
    1.226637, 0.840692, 0.603496,  # kappa
    0.006327, -0.005464, 0.003441,
    -0.000707, -0.003399, 0.010891, # sigma
    0.032577, -0.012536, -0.002748, # theta
    0.5,                             # lambda
] + [0.000524] * nobs)  # Measurement error

# You may need to adjust the bounds and constraints for your parameters
bounds = [(None, None) if i != 0 else (0, None) for i in range(len(init_para_un))]

# Define the optimization problem
opt_result = minimize(
    fun=log_likelihood,
    x0=init_para_un,
    args=(m_spot, v_mat, dt, nf, nobs),
    method='Powell',   # 'Powell', 'CG', 'TNC', 'Nelder-Mead'
    bounds=bounds,
    options={'maxiter': 5000, 'disp': True}
)

# Check the result
if opt_result.success:
    fitted_params = opt_result.x
    print("Optimization was successful.")
    print(f"Fitted parameters: {fitted_params}")
else:
    print("Optimization failed.")
    print(opt_result.message)
```

I use the following data for testing purpose:

## Answer by Achrbot (score 1, accepted)

https://quant.stackexchange.com/a/78936

I think your code is largely fine, but I have a couple of issues with it:

- You appear to be calibrating the independent-factor AFNS (since $K$ is diagonal). This means that $\Sigma$ must also be diagonal.

- I don't understand the purpose of your variable transformations, but removing it appears to help the calibration.

- Your dates are sorted as ascending, but are accessed in descending order. I'm not sure how much this matters for the in-sample fit though.

After removing the variable transformation, and making $\Sigma$ diagonal, I think the fit looks reasonable.



```
import numpy as np
from scipy.linalg import expm
from scipy.optimize import minimize
from scipy.special import expit  # Sigmoid function for the transformation
from scipy.stats import norm
import pandas as pd
import matplotlib.pyplot as plt

# DNS factor loading matrix
def NS_B(lambda_, tau):
    col1 = np.ones_like(tau)
    col2 = (1 - np.exp(-lambda_ * tau)) / (lambda_ * tau)
    col3 = col2 - np.exp(-lambda_ * tau)
    return np.column_stack((col1, col2, col3))

# yield adjustment term in AFNS
def AFNS_C(sigma, lambda_, tau):
    s = sigma
    A = s[0, 0] ** 2 + s[0, 1] ** 2 + s[0, 2] ** 2
    B = s[1, 0] ** 2 + s[1, 1] ** 2 + s[1, 2] ** 2
    C = s[2, 0] ** 2 + s[2, 1] ** 2 + s[2, 2] ** 2
    D = s[0, 0] * s[1, 0] + s[0, 1] * s[1, 1] + s[0, 2] * s[1, 2]
    E = s[0, 0] * s[2, 0] + s[0, 1] * s[2, 1] + s[0, 2] * s[2, 2]
    F = s[1, 0] * s[2, 0] + s[1, 1] * s[2, 1] + s[1, 2] * s[2, 2]

    t = tau
    la = lambda_
    r1 = -A * t ** 2 / 6
    r2 = -B * (1 / (2 * la ** 2) - (1 - np.exp(-la * t)) / (la ** 3 * t) + (1 - np.exp(-2 * la * t)) / (
                4 * la ** 3 * t))
    r3 = -C * (1 / (2 * la ** 2) + np.exp(-la * t) / (la ** 2) - t * np.exp(-2 * la * t) / (4 * la) -
               3 * np.exp(-2 * la * t) / (4 * la ** 2) - 2 * (1 - np.exp(-la * t)) / (la ** 3 * t) +
               5 * (1 - np.exp(-2 * la * t)) / (8 * la ** 3 * t))
    r4 = -D * (t / (2 * la) + np.exp(-la * t) / (la ** 2) - (1 - np.exp(-la * t)) / (la ** 3 * t))
    r5 = -E * (3 * np.exp(-la * t) / (la ** 2) + t / (2 * la) + t * np.exp(-la * t) / la -
               3 * (1 - np.exp(-la * t)) / (la ** 3 * t))
    r6 = -F * (1 / (la ** 2) + np.exp(-la * t) / (la ** 2) - np.exp(-2 * la * t) / (2 * la ** 2) -
               3 * (1 - np.exp(-la * t)) / (la ** 3 * t) + 3 * (1 - np.exp(-2 * la * t)) / (4 * la ** 3 * t))
    return r1 + r2 + r3 + r4 + r5 + r6

# Log-likelihood function
def log_likelihood(para, m_spot, v_mat, dt, nf, nobs):
    kappa = np.diag(para[:3])
    sigma = np.diag(para[3:6])
    theta = para[6:9]
    lambda_ = para[9]
    H = np.diag(para[10:10 + len(v_mat)]) 
    B = NS_B(lambda_, v_mat)
    C = AFNS_C(sigma, lambda_, v_mat)

    # Eigen decomposition for kappa
    evals, evecs = np.linalg.eig(kappa)
    # Avoiding numerical issues
    # if np.linalg.det(evecs) < 1e-10:
    #     evecs = evecs + np.eye(len(evecs)) * 1e-10

    Smat = np.linalg.inv(evecs) @ sigma @ sigma.T @ np.linalg.inv(evecs.T)

    Vdt = np.zeros((nf, nf))
    Vinf = np.zeros((nf, nf))
    for i in range(nf):
        for j in range(nf):
            Vdt[i, j] = Smat[i, j] * (1 - np.exp(-(evals[i] + evals[j]) * dt)) / (evals[i] + evals[j])
            Vinf[i, j] = Smat[i, j] / (evals[i] + evals[j])

    Q = evecs @ Vdt @ evecs.T
    Q0 = evecs @ Vinf @ evecs.T

    # Initialization
    prevX = theta
    prevV = Q0
    Phi1 = expm(-kappa * dt)
    Phi0 = (np.eye(nf) - Phi1) @ theta
    log_likelihood_val = 0
    
    X_estimate = np.zeros((nobs, nf))  # Store the estimated factors
    nobs = m_spot.shape[0]  # Number of observations

    for i in range(nobs):
        Xhat = Phi0 + Phi1 @ prevX
        Vhat = Phi1 @ prevV @ Phi1.T + Q

        y = m_spot[i, :]  # Observed yield
        y_implied = B @ Xhat + C  # Implied yield from the model
        er = y - y_implied  # Prediction error

        # Kalman gain
        ev = B @ Vhat @ B.T + H
        # # print(ev)
        # Avoiding numerical issues
        if np.abs(np.linalg.det(ev)) < 1e-10:
            ev = ev + np.eye(len(ev)) * 1e-10

        KG = Vhat @ B.T @ np.linalg.inv(ev)

        prevX = Xhat + KG @ er
        prevV = Vhat - KG @ B @ Vhat
        X_estimate[i, :] = prevX

        # Update log-likelihood
        log_likelihood_val += -0.5 * len(y) * np.log(2 * np.pi) - 0.5 * np.log(
            np.linalg.det(ev)) - 0.5 * er @ np.linalg.inv(ev) @ er

    return log_likelihood_val, X_estimate

# Create negative log-likelihood function for optimization
def neg_log_likelihood(para, m_spot, v_mat, dt, nf, nobs):
    return -log_likelihood(para, m_spot, v_mat, dt, nf, nobs)[0]

# Set the correct path to your Excel file
file_path = r'C:...'

# Read the spot rate data from Excel
df = pd.read_excel(file_path)
# Convert Date column to datetime
df['Date'] = pd.to_datetime(df['Date'])
# Sort the DataFrame by Date
df.sort_values(by='Date', inplace=True)
df = df.reset_index(drop=True)
# Convert yields to decimal assuming they are in Bps
df.iloc[:, 1:] = df.iloc[:, 1:] / 10**4
# Convert the DataFrame to a numpy matrix and remove the first column (dates)
m_spot = df.iloc[:, 1:].to_numpy()   

# Define the maturities based on the column names, convert to years
v_mat = np.array([float(mat.replace('Y', '')) if 'Y' in mat else float(mat.replace('M', ''))/12
                  for mat in df.columns[1:]])

# Define the number of observations and number of maturities
nobs, nmat = m_spot.shape

# Set the dt and nf based on the frequency of observations and number of factors
dt = 1/52  # weekly data
nf = 3  # Number of factors

# Initial guess for unconstrained parameters (para_un)
# Make sure that the number of initial parameters matches nmat + nf + 1 (for lambda)
init_para = np.array([
    1.226637, 0.840692, 0.603496,  # kappa
    0.006327, 0.003441, 0.010891, # sigma
    0.032577, -0.012536, -0.002748, # theta
    0.5,                             # lambda
] + [0.000524] * nmat)  # Measurement error

# You may need to adjust the bounds and constraints for your parameters
bounds = [(0, None) if i  >= 9 or i in [3,4,5] else (None,None) for i in range(len(init_para))]

# Test the log-likelihood function
log_likelihood(init_para, m_spot, v_mat, dt, nf, nobs)

# Create

# Define the optimization problem
opt_result = minimize(
    fun=neg_log_likelihood,
    x0=init_para,
    args=(m_spot, v_mat, dt, nf, nobs),
    method='L-BFGS-B',   # 'Powell', 'CG', 'TNC', 'Nelder-Mead'
    bounds=bounds,
    options={'maxiter': 10**4, 'disp': True}
)

# Check the result
if opt_result.success:
    fitted_params = opt_result.x
    print("Optimization was successful.")
    print(f"Fitted parameters: {fitted_params}")
else:
    print("Optimization failed.")
    print(opt_result.message)
    

# Compute factor loading matrix B and yield adjustment term C
kappa = np.diag(fitted_params[:3])
sigma = np.diag(fitted_params[3:6])#np.array([[para[3], 0, 0], [para[5], para[6], 0], [para[7], para[8], para[9]]])
theta = fitted_params[6:9]
lambda_ = fitted_params[9]
H = np.diag(fitted_params[10:10 + len(v_mat)])  # Änderung hier: Verwende die Länge von v_mat statt nobs.
B = NS_B(lambda_, v_mat)
C = AFNS_C(sigma, lambda_, v_mat)

X = log_likelihood(fitted_params, m_spot, v_mat, dt, nf, nobs)[1]

implied_yields = B @ X.T + C[:, np.newaxis]

# Plot n-latest yield curve along with the estimated one
plt.new()
n = 1
plt.plot(v_mat, m_spot[-n, :]*10**4, 'o-', label='Current Yield Curve')
plt.plot(v_mat, implied_yields[:,-n]*10**4, 'o-', label='Estimated Yield Curve')
plt.xlabel('Maturity (Years)')
plt.ylabel('Yield (Bps)')
plt.title('Current and Estimated Yield Curve (' + df.iloc[-n, 0].strftime('%Y-%m-%d') + ')')
plt.gca().spines['top'].set_visible(False)
plt.gca().spines['right'].set_visible(False)
plt.legend()
plt.grid(True)
plt.show()

# Estimate yields on a grid [0.5, 1, ... 30] years
grid_maturities = np.arange(0.5, 30.5, 0.5)
B_grid = NS_B(lambda_, grid_maturities)
C_grid = AFNS_C(sigma, lambda_, grid_maturities)

implied_yields_grid = B_grid @ X.T + C_grid[:, np.newaxis]
# Plot the estimated yield curve
plt.new()
plt.plot(grid_maturities, implied_yields_grid[:,-n]*10**4,'o-', label='Estimated Yield Curve on Grid')
plt.xlabel('Maturity (Years)')
plt.ylabel('Yield (Bps)')
plt.title('Estimated Yield Curve (' + df.iloc[-n, 0].strftime('%Y-%m-%d') + ')')
plt.gca().spines['top'].set_visible(False)
plt.gca().spines['right'].set_visible(False)
plt.grid(True)
plt.show()
```

Edit: For completeness, I have included the dataset as a table (with yields converted to Bps).

| Date | 6M | 1Y | 2Y | 3Y | 4Y | 5Y | 6Y | 7Y | 8Y | 9Y | 10Y | 12Y | 15Y | 20Y | 25Y | 30Y |
| 29/01/2024 | 386.6 | 349.3 | 290 | 268.5 | 259.9 | 256.7 | 256.2 | 257 | 258.4 | 260.1 | 262.2 | 266.4 | 269.5 | 264.6 | 254.3 | 245.5 |
| 22/01/2024 | 392.5 | 360.4 | 301.9 | 279.6 | 270.3 | 266.5 | 265 | 264.8 | 265.5 | 266.6 | 268 | 271.2 | 273.1 | 266.2 | 255 | 244.5 |
| 15/01/2024 | 387.4 | 351.6 | 292 | 269.3 | 260.5 | 257.5 | 256.8 | 257.4 | 258.7 | 260.4 | 262.4 | 266.4 | 269.3 | 264.1 | 254.1 | 245.3 |
| 08/01/2024 | 392.9 | 355.1 | 292.4 | 268.1 | 257.9 | 253.9 | 253 | 253.5 | 254.4 | 256.7 | 258.2 | 262 | 264.2 | 258.3 | 250 | 239.2 |
| 01/01/2024 | 386.1 | 346.2 | 280.8 | 256.7 | 246.2 | 243.6 | 243.4 | 244.2 | 245.9 | 247.7 | 249.4 | 253.2 | 255.8 | 251.5 | 241.2 | 232.6 |
| 25/12/2023 | 389.5 | 349.8 | 283.9 | 257.6 | 245.4 | 241.6 | 241.2 | 241.1 | 242.1 | 243.6 | 245 | 248.8 | 251.4 | 246.5 | 237.3 | 230.1 |
| 18/12/2023 | 390.4 | 357.7 | 296.3 | 270.7 | 258.9 | 253.7 | 251.9 | 251.4 | 252.1 | 253.1 | 254.7 | 257.9 | 260.1 | 254.3 | 243.8 | 235.4 |
| 11/12/2023 | 395.5 | 369.2 | 314.2 | 289.4 | 278.5 | 273.6 | 271.9 | 271.6 | 272.1 | 273.2 | 274.5 | 277.5 | 279.2 | 272.2 | 259.6 | 251.3 |
| 04/12/2023 | 394.5 | 368.8 | 314.9 | 292.1 | 282.4 | 278.9 | 277.9 | 278.4 | 279.6 | 281.2 | 283.2 | 287 | 289.6 | 283.9 | 273.1 | 264.9 |
| 27/11/2023 | 405 | 394.1 | 348 | 322.5 | 310.1 | 304.2 | 301.7 | 300.9 | 301.2 | 302.1 | 303.5 | 306.9 | 309.1 | 302.1 | 289.9 | 280.2 |
| 20/11/2023 | 405.8 | 396.2 | 350 | 325.3 | 313.5 | 308.1 | 305.6 | 304.9 | 305.2 | 306.3 | 307.8 | 311.1 | 313.2 | 305.7 | 292.1 | 282.5 |
| 13/11/2023 | 408.5 | 400.9 | 358.4 | 335.2 | 323.9 | 318.7 | 316.4 | 315.7 | 316.1 | 317 | 318.5 | 321.7 | 323.7 | 315.7 | 302.7 | 291.6 |
| 06/11/2023 | 406.6 | 397.8 | 353.1 | 329.9 | 319.8 | 316.1 | 315.4 | 316.2 | 317.9 | 320 | 322.5 | 327 | 330.2 | 323.4 | 312.4 | 300.6 |
| 30/10/2023 | 410.9 | 402.9 | 361.1 | 340.3 | 330.8 | 328 | 328 | 329.3 | 331.2 | 333.5 | 337.1 | 341.2 | 344.3 | 337.1 | 325.5 | 314.1 |
| 23/10/2023 | 409.6 | 410.4 | 376.2 | 355.9 | 345.9 | 341.5 | 339.8 | 339.5 | 340.3 | 341.6 | 343.4 | 345.4 | 347.8 | 337.3 | 323.6 | 312.6 |
| 16/10/2023 | 410.3 | 410.5 | 375.2 | 353.4 | 341.7 | 336 | 333.4 | 332.5 | 332.8 | 333.7 | 335.2 | 338.2 | 339.8 | 331.1 | 318.6 | 306.6 |
| 09/10/2023 | 412.9 | 408.3 | 367.6 | 345.8 | 335.4 | 331.1 | 329.7 | 330 | 331.3 | 332.9 | 334.9 | 338.8 | 341.1 | 332.6 | 318.5 | 307.9 |
| 02/10/2023 | 413.8 | 415.8 | 383.1 | 361.6 | 351.3 | 346.4 | 344.3 | 343.4 | 344.5 | 345.6 | 347.4 | 350.3 | 351.8 | 342.2 | 329.4 | 318.5 |
| 25/09/2023 | 412.2 | 415 | 380.9 | 357.7 | 344.9 | 338.2 | 334.6 | 333 | 332.6 | 333.1 | 334.1 | 336.8 | 338.2 | 330.5 | 319.6 | 307.9 |
| 18/09/2023 | 406.6 | 414.3 | 382.5 | 358.1 | 342.9 | 334.1 | 328.9 | 325.8 | 324.4 | 323.7 | 324 | 325.3 | 325.6 | 315.6 | 303.1 | 292.4 |