Nss

#Packages
import os
import pyodbc
import numpy as np
import pandas as pd
import itertools
from matplotlib import pyplot as plt
from numpy import linalg as LA

#Date of calculation
date = '10.01.2020'
actual_date = pd.Timestamp('{}'.format(date))

#Uploading data
tenors = pd.read_excel('Tenors.xlsx',header=None) #Tenors for Bank's yield curve
tenors_values = tenors[[0]]

russian_eurobonds_data = pd.read_excel('BI1168 10.01.2020.xlsx',header=None) #Russian Eurobonds data
tenors_real = russian_eurobonds_data[[0]] #Tenors for Russian Eurobonds
Y = russian_eurobonds_data[[1]] #Russian Eurobonds yields - like a real Y that is used for modelling

tenors_real.describe()

#Russian Eurobonds yields
plt.scatter(tenors_real,Y)

#%%time
# Nelson–Siegel model (NS model) for Russian eurobonds
variance = np.asscalar(Y.var(axis=0))

start_lam1 = 30 #Parameters for lambda1
max_to_lam1 = 100
step1 = 0.1

start_lam2 = 30 #Parameters for lambda2c
max_to_lam2 = 100
step2 = 0.1

lamb1 = [] #Generating list of lambdas1
for k in range(max_to_lam1):
    lamk1 = start_lam1 + k*step1
    lamb1.append(lamk1)

lamb2 = [] #Generating list of lambdas2
for p in range(max_to_lam2):
    lamk2 = start_lam2 + p*step2
    lamb2.append(lamk2)

ones1 = pd.DataFrame(np.ones(tenors_real.shape[0])) #Column of ones for X's massive
TNP = np.array(tenors_real) #Column of Russian Eurobonds' tenors
L = tenors_real.shape[0] #Getting an exact number of rows in Bank's tenors

X1 = []
X2 = []
X3 = []
R2 = []
rmse = []
b1 = []
b2 = []
b3 = []
b4 = []

num = len(lamb1)*len(lamb2)
index = []
for i in range(num):
    ind = i
    index.append(ind)
ind = index.copy()

z = pd.DataFrame(index = index, columns = ['lam1', 'lam2', 'R2', 'rmse', 'b1', 'b2', 'b3', 'b4'])

for lam1 in lamb1:
    for lam2 in lamb2:
        for j in range(tenors_real.shape[0]): #Generating independent variables for each lambda1 and lambda2
            x1j = (1-np.exp(-TNP[j]/lam1))/(TNP[j]/lam1)
            x2j = x1j-np.exp(-TNP[j]/lam1)
            x3j = ((1-np.exp(-TNP[j]/lam2))/(TNP[j]/lam2))-np.exp(-TNP[j]/lam2)
            X1.append(x1j)
            X2.append(x2j)
            X3.append(x3j)
        X1 = pd.DataFrame(X1) #Combining all rows of independent variables into one massive of X
        X2 = pd.DataFrame(X2)
        X3 = pd.DataFrame(X3)
        XX = pd.concat([ones1,X1,X2,X3], axis =1)
        X = XX.to_numpy()
        beta = LA.lstsq(X,Y)[0] #Getting coefficients from multiplying matrices

        residuals = np.array(Y)-np.matmul(X,beta)
        SSEi = np.asscalar(np.matmul(residuals.T,residuals))
        SSEi_aver = SSEi/L
        rmsei = (SSEi/L)**0.5
        R2i = 1 - SSEi_aver/variance
        b1i = np.asscalar(np.array(beta)[0])
        b2i = np.asscalar(np.array(beta)[1])
        b3i = np.asscalar(np.array(beta)[2])
        b4i = np.asscalar(np.array(beta)[3])
        zi = [lam1, lam2, R2i, rmsei, b1i, b2i, b3i, b4i]
        z.loc[какой-то индекс] = zi??????
z.index = range(z.shape[0])
rmse = list(z['rmse'])
index_optimal = rmse.index(min(rmse))
b1 = z['b1'][index_optimal]
b2 = z['b2'][index_optimal]
b3 = z['b3'][index_optimal]
b4 = z['b4'][index_optimal]
lbo1 = z['lm1'][index_optimal]
lbo2 = z['lm2'][index_optimal]
rmse = z['rmse'][index_optimal]

#Lambda function that will be used for bank's tenors
yc = lambda tau: b1 + b2*((1-np.exp(-tau/lbo1))/(tau/lbo1)) + b3*((1-np.exp(-tau/lbo1))/(tau/lbo1) - np.exp(-tau/lbo1)) + b4*((1-np.exp(-tau/lbo2))/(tau/lbo2) - np.exp(-tau/lbo2))

rus_eurobonds_nss = []
for t in range(tenors.shape[0]):
    yct = yc(tau = tenors[0][t])
    rus_eurobonds_nss.append(yct)
rus_eurobonds_nss = pd.DataFrame(rus_eurobonds_nss)

beta = [b1,b2,b3,b4]

X4 = []
X5 = []
X6 = []

for j in range(tenors_values.shape[0]):
    x4j = (1-np.exp(-TNP[j]/lbo1))/(TNP[j]/lbo1)
    x5j = x4j-np.exp(-TNP[j]/lbo1)
    x6j = (1-np.exp(-TNP[j]/lbo2))/(TNP[j]/lbo2)-np.exp(-TNP[j]/lbo2)
    X4.append(x4j)
    X5.append(x5j)
    X6.append(x6j)
X4 = pd.DataFrame(X4)
X5 = pd.DataFrame(X5)
X6 = pd.DataFrame(X6)
ones2 = pd.DataFrame(np.ones(tenors_values.shape[0]))
XX = pd.concat([ones2,X4,X5,X6], axis =1)
XX.columns = ['0','1','2','3']

y_hat = np.matmul(XX,beta)

fig= plt.figure(figsize=(20,8))
plt.plot(Y)
plt.plot(y_hat)