yield_curve_functions.py

#%% Package loading
import numpy as np
import math
import pandas as pd
import statsmodels.api as sm
import time
from datetime import timedelta, date
from numpy import linalg as LA
from statsmodels.tsa.api import VAR
from statsmodels.tsa.api import AR
import matplotlib.pyplot as plt
#from scipy.stats import t
from pykalman import KalmanFilter



#%% SECTION 0 : FUNCTION DEFINITIONS

#Computes Lambda parameter from the NS model
#INPUT: y = full data used for the estimation. tau = maturities vector
#OUTPUT: vector of lambdas computed for each cross-section. Recall the chosen lambda is the mean of this.
def lamb(y, tau):
    best_lamb=np.zeros(shape=(1,y.shape[0]))
    for i in range(0,y.shape[0]):
        current_best=1000;
        y_present = y[i]
        for j in range(1,2000):
            current_lamb = j/1000;
            fit = OLS_DNS_Error(y_present,current_lamb,tau)
            if (fit<current_best):
                current_best = fit;
                best_lamb[0,i] = current_lamb;
    return best_lamb

#Computes the error of the NS model on a cross section
#INPUT: data_i = vector of yields on one cross-section. lamb_i = lambda used in the NS model. tau_in = vector of maturities
#OUTPUT: square error of the NS model estimated by OLS with lambda = lamb_i, fitted to data_i.
def OLS_DNS_Error(data_i,lamb_i,tau_in):
    tau = tau_in.transpose()
    dummy= np.array(lamb_i*tau,dtype=float)
    col2 = (np.ones(shape=(tau.size,1))-np.exp(-1*dummy))/dummy
    col3 = ((np.ones(shape=(tau.size,1))-np.exp(-1*dummy))/dummy)-np.exp(-1*dummy)
    X = np.hstack((np.ones((tau.size,1)),col2,col3))
    est=sm.OLS(data_i, X)
    est = est.fit()
    f = est.params
    
    squ_error=0
    for i in range(tau.shape[0]):
        squ_error = squ_error + (DNS_formula(tau[i],f,lamb_i)-data_i[i])**2
    
    return squ_error
    
#Computes the NS value
#INPUT: x = vector of maturities. f = factors [l,s,c] of the NS model. lambb = lambda parameter of the NS model.
#OUTPUT: yield following the NS model from the input.
def DNS_formula(x,f,lambb):
    [l1,s1,c1]=f
    y=l1+s1*((1-np.exp(-lambb*x))/(lambb*x))+c1*((1-np.exp(-lambb*x))/(lambb*x)-np.exp(-lambb*x))
    return y

#Computes table of the OLS fitted time series [l,s,c]
#INPUT: data = yield data that needs to be fitted. tau_in = maturities vector. dat = dates of the time series. lamb_i = lambda parameter of the NS model
#OUTPUT: table with rows indexed by dat and 3 columns ([l,s,c] respectively).
def DNS_OLS(data,tau_in,lamb_i): #IMPORTANT : data is a numpy WITHOUT index on row or cols, tau_in is a horizontal vector, dat is a horizontal vector   
    tau = tau_in.transpose()
    dummy = np.array(lamb_i*tau, dtype=float)
    
    #matrix which will support the values of f=[l,s,c] at each time t
    f_concat = np.array(np.zeros(shape = (data.shape[0],3)));
    
    #Static Nelson Siegel for each time t
    for i in range(0,data.shape[0]):
        #selection of the current data
        y_i = np.array([data[i]]).transpose()

        #Computation of the least squares estimator for f    
        col2 = (np.ones(shape=(tau.size,1))-np.exp(-1*dummy))/dummy
        col3 = ((np.ones(shape=(tau.size,1))-np.exp(-1*dummy))/dummy)-np.exp(-1*dummy)
        
        X = np.hstack((np.ones((tau.size,1)),col2,col3))

        #OLS solution :      
        est=sm.OLS(y_i, X)
        est = est.fit()
        
        f = est.params
        f_concat[i] = f

    return f_concat

#Chosen loss function
#INPUT : n=0 : squared error / n=1 : absolute error
#OUTPUT : g(e_it)
def g(e_it,n=0):
    v=0
    if (n==0):
        v=e_it**2
    if (n==1):
        v=abs(e_it)
    return v

#Forecasting with last known value
#INPUT : da = sample of data used for the estimation. pred =  number of forecasts ahead (in days)
#OUTPUT : table with pred rows. Row i is the yield curve forecast i days ahead of the last date in da
def forecast_RW_fct(da,pred=1):
    forecast = np.zeros(shape=(pred,da.shape[1]))
    for i in range(0,da.shape[1]):
        forecast[:,i] = np.ones(shape=(pred))*da[da.shape[0]-1,i]
    return forecast

def params_VAR(da):
    model = VAR(da)
    lag_order=1
    model_fitted = model.fit(lag_order)
    trans_matrix = (model_fitted.params).to_numpy()
    return trans_matrix

#Forecasting with two-step DNS, VAR(1).
#INPUT : ts = table of OLS fitted factors (l,s,c) per cross-section (they are generated by ts_DNS_creator_daily.py). pred =  number of forecasts ahead (in days)
#OUTPUT : table with pred rows. Row i is the factor forecast i days ahead of the last date in ts.
def forecast_DNS_VAR(ts, pred): #IMPORTANT : ts has undergone the DNS_OLS function previously. pred is the date pred months after the last entry of the time series ts

    model = VAR(ts)
    model_fitted = model.fit(1, method='mle') #See Diebold and Rudebusch. All use VAR(1)
    
    lag_order=model_fitted.k_ar    
    
    return model_fitted.forecast(ts.values[-lag_order:],pred)


#Forecasting with two-step DNS, VAR(1).
# the fitting is done by the method of moments method assuming stationarity
#INPUT : ts = table of OLS fitted factors (l,s,c) per cross-section (they are generated by ts_DNS_creator_daily.py). pred =  number of forecasts ahead (in days)
#OUTPUT : table with pred rows. Row i is the factor forecast i days ahead of the last date in ts.
def forecast_DNS_VAR_yw(ts, pred):
    
    T, n_ts = ts.shape
    mu = np.mean(ts, axis=0)
    
    cov0 = 1/T *  (ts - mu).transpose() @ (ts - mu)
    cov1 = 1/T *  (ts[1:,:] - mu).transpose() @ (ts[0:-1,:] - mu)
    A = cov1 @ np.linalg.inv(cov0)
    
    # forecasts
    step_two_VAR_new = np.zeros([pred, n_ts])
    step_two_VAR_new[0,:] = mu + A@(ts[-1,:]-mu)
    for t in range(1,pred):
        step_two_VAR_new[t,:] = mu + A @ (step_two_VAR_new[t-1,:]-mu)
        
    return(step_two_VAR_new)

#Forecasting with two-step DNS, VAR(p).
#INPUT : ts = table of OLS fitted factors (l,s,c) per cross-section (they are generated by ts_DNS_creator_daily.py). pred =  number of forecasts ahead (in days)
#OUTPUT : table with pred rows. Row i is the factor forecast i days ahead of the last date in ts.
def forecast_DNS_VARm(ts,pred):
    model = VAR(ts)
    x = model.select_order(maxlags=3)
    lag_order = x.selected_orders["bic"] #we select best model based on the BIC criterion
    if lag_order==0:    #constrains not turning into a random walk
        lag_order=1
    model_fitted = model.fit(lag_order)
        
    
    return model_fitted.forecast(ts.values[-lag_order:],pred)

#Forecasting with VAR(1) using the OLS fit.
#INPUT : da = sample of data used for the estimation. pred =  number of forecasts ahead (in days)
#OUTPUT : table with pred rows. Row i is the yield curve forecast i days ahead of the last date in da.
def forecast_VAR(da, pred):
    model = VAR(da)
    lag_order = 1
    model_fitted = model.fit(lag_order)
    return model_fitted.forecast(da.values[-lag_order:],pred)

#Forecasting with VAR(1) using the Yule Walker fitting method.
#INPUT : da = sample of data used for the estimation. pred =  number of forecasts ahead (in days)
#OUTPUT : table with pred rows. Row i is the yield curve forecast i days ahead of the last date in da.
def forecast_VAR_yw(da, pred):
    
    # extract size parameters
    T, n_matu = da.shape
    mu = np.mean(da, axis = 0)
    
    # yule-walker estimation of the transition matrix
    cov0 = 1/T *  (da - mu).transpose() @ (da - mu)
    cov1 = 1/T *  (da[1:,:] - mu).transpose() @ (da[0:-1,:] - mu)
    A = cov1 @ np.linalg.inv(cov0)
    
    # forecasting
    forecast = np.zeros([pred, n_matu])
    forecast[0,:] = mu + A@(da[-1,:]-mu)
    for t in range(1,pred):
        forecast[t,:] = mu + A @ (forecast[t-1,:]-mu)
        
    return(forecast)

#Forecasting with AR(1) using the OLS fit.
#INPUT : da = sample of data used for the estimation. pred =  number of forecasts ahead (in days)
#OUTPUT : table with pred rows. Row i is the yield curve forecast i days ahead of the last date in da.
def forecast_AR(da, pred):
    forecast = np.zeros(shape=(pred,da.shape[1]))
    #forecasting per maturity (suppose no dependance between maturities)
    for i in range(0,da.shape[1]):
        #fit to time series y_{tau} for fixed tau
        model = AR(da[:,i])
        model_fitted = model.fit(1)
        
        #assume y_{tau,t+1} = b0 + b1*y_{tau,t} + \epsilon
        b0=model_fitted.params[0]
        b1=model_fitted.params[1]
        
        #print(model_fitted.params)
        
        # compute predictions
        mu = b0 / (1-b1) # compute the mean in the reparametrisation: x_t - mu = b_1 ( x - mu_{t-1}) + e_t
        
        for j in range(0,pred):
            forecast[j,i] = mu + b1**(j+1) * (da[-1,i]-mu)
        
#        #compute predictions by rolling out the recursion : \hat{y}_{tau,t+j} = b0*(1+...+b1^{j-1}) + b1^j * \hat{y}_{tau,t} 
#        for j in range(0,pred):
#            for k in range(j,pred-1):
#                forecast[k,i]=forecast[k,i]+b0*b1**k
#            forecast[j,i]=forecast[j,i] + b1**j*da[da.shape[0]-1,i] # heere's was bug
    return forecast

#%% Forecasting with AR(1) with Yule Walker equations.
#INPUT : da = sample of data used for the estimation. pred =  number of forecasts ahead (in days)
#OUTPUT : table with pred rows. Row i is the yield curve forecast i days ahead of the last date in da.
def forecast_AR_yw(da, pred):
    
    # extract size parameters
    T, n_matu = da.shape
    mu = np.mean(da, axis = 0)
    
    # yule-walker estimation of the transition matrix
    cov0 = 1/T *  (da - mu).transpose() @ (da - mu)
    cov1 = 1/T *  (da[1:,:] - mu).transpose() @ (da[0:-1,:] - mu)
    
    # force the matrices to be diagonal, therefore it is going to be AR at each maturity
    cov0 = np.diag(np.diag(cov0))
    cov1 = np.diag(np.diag(cov1))
    
    # transition matrix
    A = cov1 @ np.linalg.inv(cov0)
    
    # forecasting
    forecast = np.zeros([pred, n_matu])
    forecast[0,:] = mu + A@(da[-1,:]-mu)
    for t in range(1,pred):
        forecast[t,:] = mu + A @ (forecast[t-1,:]-mu)
        
    return(forecast)

#Forecasting with one-step DNS
#INPUT : da = sample of data used for the estimation. tau = vector of maturities. lamb = lambda factor estimated. state_init = initialisation of the states means for KF .pred =  number of forecasts ahead (in days)
#OUTPUT : table with pred rows. Row i is the yield curve forecast i days ahead of the last date in ts.
def forecast_DNS_KF(da,tau,lamb,state_init,offset_init,trans_init,pred):
    t = np.array(tau)
    t = t.transpose()
    
    dat = da.to_numpy()
    
    #STATE EQUATION
    #-- state initialisation
    initial_state_mean = state_init
    m = np.random.normal(0,1,[3,3])
    initial_state_covariance = m@m.transpose() # to achieve positive definite
    
    #-- transitions initialisation
    initial_transition_matrix = trans_init
    initial_transition_offset = offset_init
    m = np.random.normal(0,1,[3,3])
    initial_transition_covariance = m@m.transpose()
    

    #OBSERVATION EQUATION
    #Transition Matrix X is determined by Nelson Siegel dynamics
    dummy = np.array(lamb*t, dtype=float)
    
    col2 = (np.ones(shape=(tau.size,1))-np.exp(-1*dummy))/dummy
    col3 = ((np.ones(shape=(tau.size,1))-np.exp(-1*dummy))/dummy)-np.exp(-1*dummy)

    
    X = np.hstack((np.ones((tau.size,1)),col2,col3))
    
    true_observation_matrices = X
    true_observation_offsets = np.zeros(shape=(dat.shape[1]))
    
    m = np.random.normal(0,1,[dat.shape[1],dat.shape[1]])
    initial_observation_covariance = m@m.transpose()
    
    #-- Uncomment for a non-random initialisation of covariance matrices
    #random_init_observation_covariance = 0.1*np.eye(dat.shape[1])
    #initial_state_covariance = 0.1*np.eye(3)
    #initial_transition_covariance = 0.1*np.eye(3)
    
    
    #Kalman Filter estimation
    kf = KalmanFilter(transition_matrices=initial_transition_matrix,\
                      observation_matrices=true_observation_matrices,\
                      transition_covariance=initial_transition_covariance,\
                      observation_covariance=initial_observation_covariance,\
                      transition_offsets=initial_transition_offset,\
                      observation_offsets=true_observation_offsets,\
                      initial_state_mean = initial_state_mean,\
                      initial_state_covariance=initial_state_covariance,\
                      em_vars=['transition_matrices', 'transition_offsets','transition_covariance', 'observation_covariance',\
                               'initial_state_mean','initial_state_covariance'],\
                      n_dim_state=3,\
                      n_dim_obs=dat.shape[1])
    kf = kf.em(dat, n_iter=20) #EM estimation of quantities specified in em_vars
    smoother_x, smoother_var = kf.smooth(dat) #KF smoother estimation of the states
    
    #Uncomment for plots of [l,s,c] smoother estimations
    """
    fig = plt.figure()
    plt.plot( smoother_x[:,0] )
    plt.legend(['smoother x, l state'])
    plt.show()
    
    fig = plt.figure()
    plt.plot( smoother_x[:,1] )
    plt.legend(['smoother x, s state'])
    plt.show()
    
    fig = plt.figure()
    plt.plot( smoother_x[:,2] )
    plt.legend(['smoother x, c state'])
    plt.show()
    """
    # recall we had f_{t+1} = transition_matrices*f_{t} + transition_offsets + epsilon where epsilon ~ N (0, transition_covariance)
    # therefore forecasting h days ahead is done as is f*_{t+1} = (transition_matrices)*f_{t} + transition_offsets
    forecast=np.zeros(shape=(pred,3))
    last_value = np.array([smoother_x[-1]])
    offs = kf.transition_offsets
    for i in range(pred):
        forecast[i,:] = np.array(last_value @ kf.transition_matrices.transpose() + offs)
        last_value = forecast[i,:]
    return forecast

# Forecasting with one-step DNS with explosivity correction
# if the estimated transition matrix corresponds to an explosive VAR process, return the RW forecast instead
#INPUT : da = sample of data used for the estimation. tau = vector of maturities. lamb = lambda factor estimated. state_init = initialisation of the states means for KF .pred =  number of forecasts ahead (in days)
#OUTPUT : table with pred rows. Row i is the yield curve forecast i days ahead of the last date in ts.
def forecast_DNS_KF_explosivcor(da,tau,lamb,state_init,offset_init,trans_init,pred):
    t = np.array(tau)
    t = t.transpose()
    
    dat = da.to_numpy()
    
    #STATE EQUATION
    #-- state initialisation
    initial_state_mean = state_init
    m = np.random.normal(0,1,[3,3])
    initial_state_covariance = m@m.transpose() # to achieve positive definite
    
    #-- transitions initialisation
    initial_transition_matrix = trans_init
    initial_transition_offset = offset_init
    m = np.random.normal(0,1,[3,3])
    initial_transition_covariance = m@m.transpose()
    

    #OBSERVATION EQUATION
    #Transition Matrix X is determined by Nelson Siegel dynamics
    dummy = np.array(lamb*t, dtype=float)
    
    col2 = (np.ones(shape=(tau.size,1))-np.exp(-1*dummy))/dummy
    col3 = ((np.ones(shape=(tau.size,1))-np.exp(-1*dummy))/dummy)-np.exp(-1*dummy)

    
    X = np.hstack((np.ones((tau.size,1)),col2,col3))
    
    true_observation_matrices = X
    true_observation_offsets = np.zeros(shape=(dat.shape[1]))
    
    m = np.random.normal(0,1,[dat.shape[1],dat.shape[1]])
    initial_observation_covariance = m@m.transpose()
    
    #-- Uncomment for a non-random initialisation of covariance matrices
    #random_init_observation_covariance = 0.1*np.eye(dat.shape[1])
    #initial_state_covariance = 0.1*np.eye(3)
    #initial_transition_covariance = 0.1*np.eye(3)
    
    
    #Kalman Filter estimation
    kf = KalmanFilter(transition_matrices=initial_transition_matrix,\
                      observation_matrices=true_observation_matrices,\
                      transition_covariance=initial_transition_covariance,\
                      observation_covariance=initial_observation_covariance,\
                      transition_offsets=initial_transition_offset,\
                      observation_offsets=true_observation_offsets,\
                      initial_state_mean = initial_state_mean,\
                      initial_state_covariance=initial_state_covariance,\
                      em_vars=['transition_matrices', 'transition_offsets','transition_covariance', 'observation_covariance',\
                               'initial_state_mean','initial_state_covariance'],\
                      n_dim_state=3,\
                      n_dim_obs=dat.shape[1])
    kf = kf.em(dat, n_iter=20) #EM estimation of quantities specified in em_vars
    smoother_x, smoother_var = kf.smooth(dat) #KF smoother estimation of the states
    
    A = kf.transition_matrices
    offs = kf.transition_offsets
    if (np.linalg.norm( np.linalg.matrix_power(A,100000), ord=2)) > 1:
        print("explosivity corrected")
        A = np.eye(3)
        offs = np.zeros_like(offs)
    
    
    #Uncomment for plots of [l,s,c] smoother estimations
    """
    fig = plt.figure()
    plt.plot( smoother_x[:,0] )
    plt.legend(['smoother x, l state'])
    plt.show()
    
    fig = plt.figure()
    plt.plot( smoother_x[:,1] )
    plt.legend(['smoother x, s state'])
    plt.show()
    
    fig = plt.figure()
    plt.plot( smoother_x[:,2] )
    plt.legend(['smoother x, c state'])
    plt.show()
    """
    # recall we had f_{t+1} = transition_matrices*f_{t} + transition_offsets + epsilon where epsilon ~ N (0, transition_covariance)
    # therefore forecasting h days ahead is done as is f*_{t+1} = (transition_matrices)*f_{t} + transition_offsets
    forecast=np.zeros(shape=(pred,3))
    last_value = np.array([smoother_x[-1]])    
    for i in range(pred):
        forecast[i,:] = np.array(last_value @ A.transpose() + offs)
        last_value = forecast[i,:]
    return forecast

#Sample autocovariance at lag k
#INPUT:  k = lag, d = time series
#OUTPUT: sample autocovarance of time series d at lag k
def gamh(k,d):
    gam = 0
    db = np.mean(d)
    for tt in range((abs(k)+1),(d.shape[0])):
        gam = gam + (d[tt]-db)*(d[tt-abs(k)]-db)    
    
    gam = gam/(d.shape[0])
    return gam

#HLN test statistic
#INPUT: y1,y2 = vectors of forecasts of the two methods to compare (here y1 -> RW). y = vector of corresponding true value. h = number of steps ahead of the forecast
#OUTPUT: HLN test statistic
def DM(y1,y2,y,h):
    dm = 0
    if (y1.shape==y2.shape) and (y1.shape==y.shape):
        e_1 = y1 - y
        e_2 = y2 - y
        T = y.shape[0]
        d = g(e_1) - g(e_2)
        dbar = np.mean(d)
        fh0 = gamh(0,d)
        M = int(math.floor(math.pow(T,1/3))+1)
        for k in range(-M,M):
            fh0 = fh0 + 2*gamh(k,d)
        fh0 = fh0*(1/(2*math.pi))
        dm = dbar/(math.pow((2*math.pi*fh0)/T,1/2))
        hln=math.pow((T+1-2*h+h*(h-1))/T,1/2)*dm
        return hln
    else: 
        return -10000

#Student Test
#INPUT: tval = table with sample values of the test statistic. n = degrees of freedom.
#OUTPUT: two tables.
#Table 1: p value table ("-" in front of the p_value indicates that p-value is for RW significantly better, "+" sign for method better)
#Table 2: indicator table. entry = sign(test_stat). +1: method is (maybe) better, -1: RW is (maybe) better -> see Table 1 for p-values
def student_test(tval,n):
    
    ptab=np.zeros(shape=(tval.shape[0],tval.shape[1]))
    qtab=np.zeros(shape=(tval.shape[0],tval.shape[1]))
    
    for i in range(ptab.shape[0]):
        for j in range(ptab.shape[1]):

            if (tval[i][j]>0):  #a positive value means the method performs better than RW
                ptab[i][j] = t.cdf(-tval[i][j],n,0,1)
                if(t.cdf(-tval[i][j],n,0,1)<=0.05): qtab[i][j] = 1
                
            if (tval[i][j]<0): #a negative value means the method performs worse than the RW
                ptab[i][j] = -t.cdf(tval[i][j],n,0,1)
                if(t.cdf(tval[i][j],n,0,1)<=0.05): qtab[i,j] = -1
    return np.hstack((ptab,qtab))