Actualmente estoy tratando de implementar el modelo Nelson Siegel sin Arbitraje (AFNS) en Python. Sin embargo, me encuentro con el problema de que mis resultados no coinciden en absoluto con la curva de rendimiento actual. ¿Hay alguien aquí que tenga experiencia con el AFNS y/o pueda ayudar con el filtro de Kalman para determinar las variables?
Sigo el trabajo de Christensen, Diebold y Rudebusch (https://www.frbsf.org/wp-content/uploads/wp07-20bk.pdf)
¡Muchas gracias de antemano! Aquí está mi código
import numpy as np
from scipy.linalg import expm
from scipy.optimize import minimize
from scipy.special import expit # Función sigmoide para la transformación
from scipy.stats import norm
import pandas as pd
# Matriz de carga de factores DNS
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))
# Término de ajuste del rendimiento en 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
# Restricciones de parámetros
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 # error de medición
return bb
# Función de log-verosimilitud
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)
# Descomposición en valores propios para 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
# Inicialización
prevX = theta
prevV = Q0
Phi1 = expm(-kappa * dt)
Phi0 = (np.eye(nf) - Phi1) @ theta
log_likelihood_val = 0
nobs = m_spot.shape[0] # Número de observaciones
for i in range(nobs):
Xhat = Phi0 + Phi1 @ prevX
Vhat = Phi1 @ prevV @ Phi1.T + Q
y = m_spot[i, :] # Rendimiento observado
y_implied = B @ Xhat + C # Rendimiento implícito del modelo
er = y - y_implied # error de predicción
# Ganancia de Kalman
ev = B @ Vhat @ B.T + H
KG = Vhat @ B.T @ np.linalg.inv(ev)
prevX = Xhat + KG @ er
prevV = Vhat - KG @ B @ Vhat
# Actualizar log-verosimilitud
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
# Establecer la ruta correcta a tu archivo de Excel
file_path = r'C:....'
# Leer los datos de tasas de interés del archivo de Excel
df = pd.read_excel(file_path)
# Convertir el DataFrame a una matriz numpy y eliminar la primera columna (fechas)
m_spot = df.iloc[:, 1:].to_numpy() / 10000 # Suponiendo que los rendimientos están en puntos básicos
# Definir los vencimientos basados en los nombres de las columnas, convertir a años
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
# Definir el número de observaciones y número de vencimientos
nobs, nmat = m_spot.shape
# Establecer el dt y nf basados en la frecuencia de observaciones y número de factores
dt = 1/52 # datos semanales
nf = 3 # Número de factores
# Suposición inicial para parámetros no restringidos (para_un)
# Asegúrate de que el número de parámetros iniciales coincida con nmat + nf + 1 (para 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) # Error de medición
# Es posible que necesites ajustar los límites y restricciones para tus parámetros
bounds = [(None, None) if i != 0 else (0, None) for i in range(len(init_para_un))]
# Definir el problema de optimización
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}
)
# Revisar el resultado
if opt_result.success:
fitted_params = opt_result.x
print("La optimización fue exitosa.")
print(f"Parámetros ajustados: {fitted_params}")
else:
print("La optimización falló.")
print(opt_result.message)Utilizo los siguientes datos con fines de prueba: 
