import pandas as pd
def cubicSpineInterpolation(ZeroRates,ZeroRatesTenors,RequestTenor):
    import numpy as np
    RequestZeroRates = RequestTenor.astype(np.float)
    h_original = ZeroRatesTenors[1:]-ZeroRatesTenors[0:-1]
    h = np.delete(h_original,-1)
    A = np.zeros(shape=(len(ZeroRates)-2,len(ZeroRates)-2))
    A[0,0] = 4
    A[0,1] = 1
    A[-1,-1] = 4
    A[-1,-2] = 1
    for i in range(1,(A.shape[0]-1)):
        A[i,i-1]=1
        A[i,i]=4
        A[i,i+1]=1

    Right = ZeroRates[2:]-2*ZeroRates[1:-1]+ZeroRates[0:-2]
    Right = np.transpose(np.multiply(Right,pow(1/h,2))*6)
    #print(np.dot(np.linalg.inv(A),b))
    sigma = np.dot(np.linalg.inv(A),Right)
    sigma= np.insert(sigma,0,0)
    sigma = np.append(sigma,0)

    #print(sigma.shape)
    #print(h.shape)
    a = (sigma[1:]-sigma[0:-1])/(6*h_oroginal)
    b = sigma[0:-1]/2
    c = (ZeroRates[1:]-ZeroRates[0:-1])/h_original - h_original*(2*sigma[0:-1]+sigma[1:])/6
    d = ZeroRates[0:-1]
    #print(c)
    #print(b)
    for i in range(0,len(RequestTenor)):
        j = 0
        while RequestTenor[i]>ZeroRatesTenors[j+1]:
            j=j+1
        #print(a[j]*pow(RequestTenor[i]-ZeroRatesTenors[j],3)+b[j]*pow(RequestTenor[i]-ZeroRatesTenors[j],2)+c[j]*(RequestTenor[i]-ZeroRatesTenors[j])+d[j])
        RequestZeroRates[i] = a[j]*pow(RequestTenor[i]-ZeroRatesTenors[j],3)+b[j]*pow(RequestTenor[i]-ZeroRatesTenors[j],2)+c[j]*(RequestTenor[i]-ZeroRatesTenors[j])+d[j]
    return RequestZeroRates

RequestTenor = np.arange(2,361,1)
ZeroRatesTenors = np.array([1,3,6,12,24,36,60,84,120,240,360])
ZeroRates = np.array([1.68,1.85,2.05,2.26,2.50,2.66,2.82,2.93,2.97,3.03,3.31])

test = cubicSpineInterpolation(ZeroRates,ZeroRatesTenors,RequestTenor)
import matplotlib.pyplot as plt
plt.plot(RequestTenor,test)
plt.show()