1
$\begingroup$

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: Test_data

$\endgroup$

1 Answer 1

1
$\begingroup$

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

  1. You appear to be calibrating the independent-factor AFNS (since $K$ is diagonal). This means that $\Sigma$ must also be diagonal.
  2. I don't understand the purpose of your variable transformations, but removing it appears to help the calibration.
  3. 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.

  • Note: I had to change the solver to Nelder-Mead, as well as a few other tweaks to get it to work. The altered code is included below.
    • After playing around with it a bit, I get the best results with the 'L-BFGS-B' solver.

Yield Curve

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
$\endgroup$
9
  • $\begingroup$ Thanks a lot for your answer and the help!! One point, I added the following: ev_reg = ev + np.eye(ev.shape[0]) * 1e-5 KG = Vhat @ B.T @ np.linalg.inv(ev_reg) as I got the error that I try to calculate the inverse of a singular matrix ("numpy.linalg.LinAlgError: Singular matrix"). When I do this I end with a curve which looks like your curve but multiplied with -1 till 20Y and fits afterward the current curve (and your estimate) again. Did you change something at the uploaded code compared to the one you used to generate the graph? $\endgroup$
    – Marc157
    Apr 10 at 21:17
  • $\begingroup$ The pasted code is identical to the one used the generate the figure (except for the path ofc). I didn’t get any trouble with a singular ev. Perhaps double check that your data is loaded correctly? I initially had trouble with the dates, and yields not being in bps. I have added a table with the data I used, to my answer. $\endgroup$
    – Achrbot
    Apr 11 at 8:11
  • $\begingroup$ Thanks a lot! Maybo two last things: I need the model to obtain interest rates for every maturity in np.arange(0.5, 30, 0.5). When I try I always run into errors...Could you help a last time on this? $\endgroup$
    – Marc157
    Apr 11 at 18:57
  • $\begingroup$ Second: When I add more data (more observations)in negative interest environment, I end up with a pretty bad fit. Is there a way to improve this again? I use a shift of +3% to all values, as I want to price caps with Black76 based on this model estimates. At some point I even end with "numpy.linalg.LinAlgError: Singular matrix" $\endgroup$
    – Marc157
    Apr 11 at 19:31
  • $\begingroup$ 1): To obtain yields at some maturity grid, you can construct appropriate $B$ and $C$ matrices, and compute yields, using the appropriate factors. 2): Running if np.linalg.det(ev) < 1e-10: ev = ev + np.eye(len(ev)) * 1e-10 alleviates the singular matrix error, and switching solver to 'L-BFGS-B' helps the calibration. I haven't tried it on additional data, but I have no problem calibrating to your yield curve, shifted down 300 bps. $\endgroup$
    – Achrbot
    Apr 12 at 7:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.