<jupyter_start><jupyter_text>Pricing a Basket CDS Author: <NAME>  IntroductionBasket CDS is an exotic structured credit product, which is linked to a variety of underlying reference assets. The protection seller provides protection against a specified reference entity suffering a credit event, usually the 'k-th to default', if the k-th default occurs, the protection buyer receives a contingent payment of only 1 name in the basket, and the contract terminates. For a protection buyer that wants to buy protection to hedge default risk on credit exposure of several securities, buying CDS for each reference entity in the portfolio will cost a large amount of spreads, a basket CDS is then a reasonable choice to consider. In this project, 5 names are equally weighted, each of them has 2 million notional with maturity  years, let  be the total notional, we have:$$LGD = (1-R)\times\frac{1}{5}\times N$$Let $PL$ be the premium leg and $s$ be the spread basis points, the protection seller receives $PL = s \times 10$ million per year. If we consider the basket CDS as a credit derivative that has a structure that is similar to synthetic CDO, the 1st to default name is like the equity tranche, and the 2nd to default name is like the Mezzanine tranche, and so on. Thus, like the methodology of pricing synthetic CDOs, we need to cope with the joint default distribution using copula method.  Import DataDue to limited resources, I have only got CDS spread quotations of 5 reference names with tenors of 1Y, 2Y, 3Y, 4Y and 5Y on 15 December 2020. As for estimating correlation matrix, I choose standarized stock returns using stock price data form [Yahoo Finance](https://finance.yahoo.com). The risk free rate is set to be 10 year US Treasury rate on 15 December 2020, which is 0.92.<jupyter_code>import numpy as np
import pandas as pd
from numpy.linalg import inv, det, cholesky
from scipy.special import gamma
import math
import matplotlib.pyplot as plt
from scipy.stats import norm, t, chi2
import openpyxl
%matplotlib inline
import cufflinks as cf
cf.set_config_file(offline=True)<jupyter_output><empty_output><jupyter_text>Snapshot of Credit spreadsThe choosen entities are: Google(GOOG), Amazon(AMZN), Microsoft(MSFT), Apple(AAPL), Netflix(NFLX), which are top firms worldwide.The table below shows spread term structures of 5 reference names, we can see that spread increases according to the increase of maturity.<jupyter_code>spread = pd.read_csv('Data\Reference Credit Spreads\spread.csv', header=None)
spread.columns=['1Y','2Y','3Y','4Y','5Y']
name = ['GOOG', 'AMZN', 'MSFT', 'AAPL', 'NFLX']
spread.index = name
print(spread)
# spread2 = pd.DataFrame(spread.values.T, index=spread.columns, columns=name)
# print(spread2)
spread2 = pd.DataFrame(spread.values.T, index=spread.columns, columns=spread.index)
spread2.iplot(
    xTitle = 'Maturity',
    yTitle = 'Spread',
    title = 'CDS Term Structure'
)<jupyter_output><empty_output><jupyter_text>Historical Stock PriceThe data contains the historical stock prices of the 5 companies form Jan 5 2016 to Dec 14 2020, they are converted to standardized returns.<jupyter_code>def import_stock_returns(names):
    stocks=pd.DataFrame()
    for name in names:
        df=pd.read_csv('Data\Asset Prices for Correlation\{}.S.csv'.format(name),parse_dates=[0],usecols=['Date','Adj Close'])
        df.columns=['Date',name]
        if name == names[0]:
            stocks=df.copy()
        else:
            stocks=stocks.merge(df,on='Date')
   
    data=stocks.drop('Date',axis=1)
    # Calculate stock returns
    returns=(data-data.shift(1))/data.shift(1)

    # Convert returns to standardised returns
    returns=(returns-returns.mean())/returns.std()
    returns.dropna(inplace=True)
    # returns=returns.values
    
    return returns

standardized_returns = import_stock_returns(name) 
standardized_returns.to_csv('Data\Asset Prices for Correlation\st_returns.csv') # This file is for pseudo sample calculation on MATLAB 
s_returns = standardized_returns.values
# print(s_returns)<jupyter_output><empty_output><jupyter_text>Discount FactorGiven risk-free rate $r$, discount factor between time $i$ and $j$ ($j>i$) is:$$D(i,j) = e^{-r(j-i)}$$This is very useful when calculating payments of k-th to default CDS.<jupyter_code>risk_free_rate = 0.0092

def DF(i,j,rate):
    return np.exp(-rate*(j-i))

discount_factor = []
for i in range(1, 6):
    discount_factor.append(DF(0,i,risk_free_rate))
print(discount_factor)<jupyter_output>[0.9908421905166154, 0.9817682465077646, 0.9727773999494099, 0.9638688898509309, 0.9550419621907147]<jupyter_text>Bootstrap Hazard Rate MethodologyGiven the number of default events occured by time $t$ follows an Inhomogenous Poisson Process, the survival probability can be derived as:$$P(T_n) = e^{-\int_{0}^{T_n}\lambda_{s}ds}$$where $\lambda_{s}$ is the hazard rate of this time interval. Given a set of CDS with different tenors and spreads, we can calibrate hazard rate term structure for every one of reference names.If we descretize the time the maturity into number of $\Delta t$ (1 year), the corresponding end of period maturities are then expressed as $T_n = n\Delta t$, probability of survival given discretized time and (piecewise constant) hazard rates will then be represented as$$P(T_n)=e^{-\sum_{k=1}^{n}\lambda_k\Delta t}$$The fair quots of the spread $S_N$ is such that the expected present value of the payments made by buyer and seller are equal:$$PL_N = DL_N$$Thus, we have:$$S_N=\frac{(1-R)\sum_{n=1}^{N}D(0,T_n)(P(T_n-1)-P(T_n))}{\sum_{n=1}^{N}P(T_n)(\Delta t_n)}$$where $R$ is recovery rate. Given the obligor is solvent at time 0, we have:$$P(T_0)=P(0)=1$$assume $S_1$ is the spread of year 1 ($T_1$), form the above equation we have:$$S_1 = \frac{(1-R)(1-P(T_1))}{P(T_1)\Delta t_1}$$ So the first step of bootstraping is:$$P(T_1)=\frac{1-R}{(1-R)+\Delta t_1 S_1}$$Let $L = 1 - R$, we can derive the second step:$$P(T_2)=\frac{D(0,T_1)[L·1-(L+\Delta t_1 S_2)P(T_1)]}{D(0,T_2)(L+\Delta t_2 S_2)}+\frac{P(T_1)L}{L+\Delta t_2 S_2}$$In general, the expression of $P(T_n)$ is:$$P(T_n)=\frac{\sum_{n=1}^{N-1}D(0,T_n)[L·P(T_{n-1})-(L+\Delta t_n S_N)P(T_n)]}{D(0,T_N)(L+\Delta t_n S_N)}+\frac{P(T_{N-1})L}{L+\Delta t_n S_N}$$Given survival probability, the piecewise constant hazard rate is given by:$$\lambda_m = -\frac{1}{\Delta t}log\frac{P(T_m)}{P(T_{m-1})}$$<jupyter_code># Bootstrap default probability and derive hazard rate term structure
# Note: spreads need to be converted !!! 
def bootstrap(spread, dt, lgd, discount_rate):
    surv_prob=np.zeros(len(spread)+1,dtype=np.float64) # include T = 0
    lambdas=np.zeros(len(surv_prob),dtype=np.float64)
    
    #Computing the survival probability of the first year by the formula: RR/(RR+spread*delta_t) 
    surv_prob[0]=1
    surv_prob[1]=lgd/(lgd+dt*spread[0])
    lambdas[0]=0
    
    for i in range(2,len(surv_prob)):
        t1=0
        for j in range(1,i):
            t1+=DF(0,j,discount_rate)*(lgd*surv_prob[j-1]-(lgd+dt*spread[i-1])*surv_prob[j])
    
        surv_prob[i]=(t1/(DF(0,i,discount_rate)*(lgd+dt*spread[i-1])))+(surv_prob[i-1]*lgd)/(lgd+dt*spread[i-1])
        
    for i in range(1,len(lambdas)): # Calibrate hazard rates
        lambdas[i]=(-1/dt)*np.log(surv_prob[i]/surv_prob[i-1])
    
    return lambdas[1:len(lambdas)],surv_prob[1:len(surv_prob)] # exclude T = 0
recovery = 0.4 # recovery rate
lgd = 1 - recovery
dt = 1

hazard_rate_term_structure = np.zeros((5,5), dtype = np.float64)
survival_probabilities = np.zeros((5,5), dtype = np.float64)

# Get bootstraped hazard rates and survival probabilities for all reference names
for i in range(0, len(name)):
    l, s = bootstrap(spread.loc[i,:]/10000, dt, lgd, risk_free_rate)
    hazard_rate_term_structure[i,:] = l
    survival_probabilities[i,:] = s

print(hazard_rate_term_structure)
print(survival_probabilities)<jupyter_output>[[0.00169523 0.00314736 0.00612442 0.00690061 0.00863557]
 [0.0022957  0.00375074 0.00582665 0.00791905 0.0101803 ]
 [0.00103113 0.00193483 0.00335872 0.00611197 0.00802148]
 [0.00133744 0.00231128 0.00357359 0.00559771 0.00913865]
 [0.00685313 0.01256349 0.01679794 0.02094858 0.03881454]]
[[0.99830621 0.99516912 0.9890929  0.98229106 0.97384493]
 [0.99770694 0.99397181 0.98819712 0.98040244 0.97047229]
 [0.9989694  0.99703843 0.99369527 0.98764035 0.97974971]
 [0.99866346 0.99635793 0.99280371 0.98726181 0.97828067]
 [0.9931703  0.98077067 0.96443334 0.94443998 0.90848429]]<jupyter_text>Term StructuresThe tables and pictures below shows term structures of hazard rate and survival probabilities. As maturity increases, hazard rate goes up while survival probability goes down.<jupyter_code>term_structure = pd.DataFrame(hazard_rate_term_structure)
term_structure.columns = ['1Y','2Y','3Y','4Y','5Y']
term_structure.index = name
term_structure2 = pd.DataFrame(term_structure.values.T, index=term_structure.columns, columns=name)
print(term_structure2)
term_structure2.iplot(title='Hazard Rate Term Structure', 
                                          xTitle='CDS Maturity', 
                                          yTitle='Hazard Rate' 
                                          )
term_structure3 = pd.DataFrame(survival_probabilities)
term_structure3.columns = ['1Y','2Y','3Y','4Y','5Y']
term_structure3.index = name
term_structure4 = pd.DataFrame(term_structure3.values.T, index=term_structure3.columns, columns=name)
print(term_structure4)
term_structure4.iplot(title='Survival Probability Term Structure', 
                                          xTitle='CDS Maturity', 
                                          yTitle='Survival Probability' 
                                          )<jupyter_output>GOOG      AMZN      MSFT      AAPL      NFLX
1Y  0.998306  0.997707  0.998969  0.998663  0.993170
2Y  0.995169  0.993972  0.997038  0.996358  0.980771
3Y  0.989093  0.988197  0.993695  0.992804  0.964433
4Y  0.982291  0.980402  0.987640  0.987262  0.944440
5Y  0.973845  0.970472  0.979750  0.978281  0.908484<jupyter_text>Calculation of Exact Default TimeAfter sampling from Gaussian and Student-t copulae, we get a set of realizations of correlated uniform random variables $U = [u_1, u_2, u_3, u_4, u_5]$, given a set of bootstraped hazard rates, the exact future default time for each reference name can be derived by converting $u_i$ to $\tau_i$ ($i=1,2,...,5$ for the index of reference names), where $\tau$ stands for actual default time. From the properties of Poisson Process, $\tau_i$ follows an exponential distribution whoose intensity has a term structure:$$\tau_{i} \ \sim \ Exp(\lambda_{1Y},\lambda_{2Y},...,\lambda_{5Y})$$ From the intuitive of Probability Integral Transformation we have:$$u_i = 1-e^{-\lambda_{\tau_i}\tau_i}$$$$log(1-u_i)=-\lambda_{\tau_i}\tau_i$$Hazard rates are piecewise constant, size of each piece is equally 1 year, so at first the year that the default happens is needed, next, a year fraction is needed to present the actual default time in that year.Suppose the default occurs between year $t_m$ and $t_{m-1}$, we have:$$\tau = t_{m-1}+\delta t$$where $\delta t$ stands for year fraction.As the hazard rates are piecewise constant, we can get $t_{m-1}$ through an iterative process:$$\tau = inf\{t>0 \ : \ log(1-u) \geq -\sum^{t}\lambda_m \}$$where default occurs if inequality holds and $t{m-1}\leq \tau \leq t_m$, if the inequality holds after adding $\lambda_m$ then the default occurs.As for year fraction, we have:$$1-u = e^{-\int_{0}^{t_{m-1}+\delta t}\lambda_s ds}=P(0,t_{m-1})e^{-\int_{t_{m-1}}^{t_{m-1}+\delta t}\lambda_s ds}$$$$log(\frac{1-u}{P(0,t_{m-1})})=-\delta t \lambda_m$$$$\delta t = -\frac{1}{\lambda_m}log(\frac{1-u}{P(0,t_{m-1})})$$<jupyter_code># Compute actual default time
# usim: simulated vector of uniform variables
def default_times(lambdas,survival_probabilities,usim):

    default_year=[]
    delta_t=[]
    # Find the year that default occurs
    for i in range(0,len(lambdas)):
        flag=0
        for j in range(0,len(lambdas.values[i])):
            if np.sum(lambdas.values[i][0:j]) >= np.abs(np.log(1-usim[i])): # for number i reference name
                default_year.append(j-1)
                flag=1
                break
        if flag==0: # no default occurs
            default_year.append(999)

    for i in range(0,len(default_year)):
        if default_year[i] == 999: # at that year no default occurs
            delta_t.append(999)
        else:
            delta_t.append((-1/lambdas.values[i][default_year[i]])*np.log((1-usim[i])/(survival_probabilities.values[i][default_year[i]]))) # year fraction        
        
    default_time=np.array(default_year)+np.array(delta_t)
    
    return default_time,default_year 
# A test case
sim = [0.002, 0.008, 0.008, 0.4397, 0.1235]
d1, d2 = default_times(term_structure, term_structure3, sim)
print(d1)
print(d2)<jupyter_output>[9.74701147e-02 1.34080229e+00 2.27936703e+00 1.99800000e+03
 1.99800000e+03]
[1, 2, 3, 999, 999]<jupyter_text>Calculation of CorrelationsCorrelation matrix is essential for copula fitting. Gaussian requires linear correlation under the assumption of normality, while Student-t needs to calculate rank correlation at first and then convert it to linear form.  Linear Correlation for Gaussian CopulaFirstly, the normality of the data must be checked though it has been standardized, because with wrong correlation we may not sample an uniform vector from Gaussiam Copula. Kolmogorov-Smirnov Test results shows the standardized returns of all 5 companies are significantly different from normal distribution, so we need extra measures to calculate linear correlation.First, we convert the data to uniform distribution $U$ with empirical CDFs of 5 standardized returns, then we apply inverse normal CDF on $U$ to obtain a normal sample, now we can calculate Pearson correlation matrix of this normal sample, this procedure is done using MATLAB and result is shown below. The pseudo sample $U$ is also used in estimating degree of freedom of t-Copula.<jupyter_code>from scipy.stats import kstest
kstest(standardized_returns.GOOG, cdf='norm')
kstest(standardized_returns.AMZN, cdf='norm')
kstest(standardized_returns.AAPL, cdf='norm')
kstest(standardized_returns.MSFT, cdf='norm')
kstest(standardized_returns.NFLX, cdf='norm')
c_p = pd.read_csv('Pearson_corr.csv', header=None)
Corr_Pearson = c_p.values
Corr_Pearson # Pearson Correlation<jupyter_output><empty_output><jupyter_text>Kendall's Tau for t-CopulaThe Kendall's Tau correlation is a kind of rank correlation, formally defined as the difference between probability of concordance and dis-concordance. It measures the degree to which the large values of one variable associate with the large values of another variable. Kendall's tau is calculated on standardized returns rather than pseudo samples already used, and it needs to be convert to linear correlation by $\rho = sin(\frac{\pi}{2}\rho_{\tau})$.The non-converted Kendall's tau is calculated by MATLAB, the converted matrix is shown below.<jupyter_code>c_k = pd.read_csv('Corr_Kendall.csv', header=None) # the original Kendall's tau
Corr_Kendall = np.sin(np.pi/2*c_k.values) # convert to linear
Corr_Kendall<jupyter_output><empty_output><jupyter_text>Maximum Likelihood for t-Copula's Degree of FreedomIn the flexible Canonical Maximum Likelihood approach, we first obtain the default correlation matrix and then use the same underlying data (pseudo-samples $U^{Hist}$, which is $U$ mentioned above), in order to calculate the log-likelihood function (of $\nu$) as sum of densities.$$argmax_{\nu} \ \sum_{t=1}^{T}log c_t(U_{t}^{Hist}; \nu, \hat{\Sigma}) $$The above representation means given a set of candidate $\nu$, for each day of our dataset, we aggregate the logarithm of t-Copula density using candidate $\nu$ and the estimated correlation matrix, then the $\nu$ that maximazes the above likelihood function can be derived, this procedure is completed using MATLAB. The optimal $\nu$ is 4.7304.     Payment CalculationWhen simulating the premium legs and default legs, in each iteration, we have 5 simulated default times, if there are no defaults, or the shortest default time is greater than or equal to 5 years, or the default time of this sample is less than k, all 5 premiums should be payed, and the default leg equals 0.For simplicity, the reference entities that defaulted before the k-th default are not removed, thus the spread of the k-th to default basket CDS is:$$s=\frac{(1-R)D(0,\tau_k)\times\frac{1}{5}}{D(0,\tau_k)\tau_k}$$<jupyter_code># Fair PL and DL for k-th to default CDS
# k: k-th to default
def payments(default_time,default_year,k,recovery,risk_free_rate):
    
# Disregarding defaults which happen just after the contract is signed. To get appropriate results
    for i in range(len(default_time)):
        if default_time[i] <= 0.25:
            default_time[i] = 999+999

    order_of_default = np.argsort(default_time)

    defaults = len(default_year)-default_year.count(999) # compute actual defaults in this iteration
    
    pl = 0
    dl = 0

    if defaults == 0 or np.min(default_time) >= 5:
        pl = np.sum(discount_factor) # All 5 names are equally weighted
        dl = 0
    else:
        if defaults >= k: # Compute PL, DL of k-th to default
            tau = default_time[order_of_default[k-1]]
            pl = tau*DF(0, tau, risk_free_rate)
            dl = (1-recovery)/5*DF(0, tau, risk_free_rate)
        else:
            pl = np.sum(discount_factor)    
            dl = 0

    return pl*10000,dl*10000 # convert to bps<jupyter_output><empty_output><jupyter_text>Sampling from Copulae Sampling default time from Gaussian Copula- Use Cholesky decomposition to get $\Sigma_{Pearson}=AA^T$- For each simulation, draw a 5-dimentional vector of independent standard normal variables $Z = (z_1, z_2,...z_5)^T$- Compute a vector of correlated variables by $X = AZ$- Use normal CDF to map to a uniform vector $U=\Phi(X)$, where $U = (u_1, u_2,...u_5)^T$- Convert $u_i \ (i=1,2,...,5)$ to $\tau_i$ to get default time vector Sampling default time from t-Copula- Use Cholesky decomposition to get $\Sigma_{Kendall}=AA^T$- For each simulation, draw a 5-dimentional vector of independent standard normal variables $Z = (z_1, z_2,...z_5)^T$- Draw an independent chi-square random variables $s \sim \chi^2_{\nu}$ by $s = z_1^2+z_2^2+\cdots+z_{\nu}^2$, where $\nu$ is the estimated degree of freedom- Compute 5-dimentional t vector $Y=Z/\sqrt{\frac{s}{\nu}}$- Impose correlation by $X=AY$- Map to a correlated uniform vector using t CDF by $U=T_{\nu}(X)$- Convert $u_i \ (i=1,2,...,5)$ to $\tau_i$ to get default time vector<jupyter_code># Simulate a uniform vector from Gaussian Copula

def sampling_from_Gaussian(corr):
    zsim = np.random.standard_normal(5) # draw z vector
    A = cholesky(corr)
    usim = norm.cdf(np.dot(A,zsim)) 
    return usim
# Simulate a uniform vector from t Copula
# df: degree of freedom
def sampling_from_t(corr, df):          
    A = cholesky(corr)
    z = np.random.standard_normal(5)
    z_chi = np.random.standard_normal(df) # draw another set of standard normal variables to aggregate chi-square
    chisquare = np.sum(z_chi**2)
    y = z/((chisquare/df)**0.5)
    usim = t.cdf(np.dot(A,y),df)
    return usim
u = sampling_from_Gaussian(Corr_Pearson)
u
v = sampling_from_t(Corr_Kendall, 4)
v<jupyter_output><empty_output><jupyter_text>Monte Carlo MC procedureIn each iteration:- Sample a uniform vector $u$ from Gaussian or t coopula- Convert $u$ to actual default time- Compute PL and DLThen we have a set of PL and a set of DL, average them and compute fair spread<jupyter_code># Calculate Fair Spreads of k-th to default basket CDS by sampling from Gaussian-Copula
# k: k-th to default
# corr: correlation
# lambdas: hazard rate term structure
def fair_spread_MC_Gaussian(k, corr, recovery, risk_free_rate, lambdas, survival_probabilities, iteration):
    running_pl = 0
    running_dl = 0
    
    i = 0

    for i in range(0, iteration+1):
        u = sampling_from_Gaussian(corr) # sampling from Gaussian copula
        default_time, default_year = default_times(lambdas,survival_probabilities,u) # calibrate actual default time
        pl, dl = payments(default_time, default_year, k, recovery, risk_free_rate) # compute PL, DL
        running_pl += pl
        running_dl += dl
        i += 1

    avg_pl = running_pl/iteration
    avg_dl = running_dl/iteration

    s = avg_dl/avg_pl
    return s*10000 # Convert to bps
# Calculate Fair Spreads of k-th to default basket CDS by sampling from t-Copula
# df: degree of freedom
def fair_spread_MC_t(k, df, corr, recovery, risk_free_rate, lambdas, survival_probabilities, iteration):
    running_pl = 0
    running_dl = 0
    
    i = 0

    for i in range(0, iteration+1):
        u = sampling_from_t(corr, df)
        default_time, default_year = default_times(lambdas,survival_probabilities,u)
        pl, dl = payments(default_time, default_year, k, recovery, risk_free_rate)
        running_pl += pl
        running_dl += dl
        i += 1

    avg_pl = running_pl/iteration
    avg_dl = running_dl/iteration

    s = avg_dl/avg_pl
    return s*10000 # Convert to bps
spread_t = []
for k in range(1, 6): # 1st, 2nd, 3rd, 4th, 5th to default
    spread_t.append(fair_spread_MC_t(k, 4, Corr_Kendall, recovery, risk_free_rate, term_structure, term_structure3, 10000))

df_t = {
    "1st": [spread_t[0]],
    "2nd": [spread_t[1]],
    "3rd": [spread_t[2]],
    "4th": [spread_t[3]],
    "5th": [spread_t[4]]
}     

fair_spread_t0 = pd.DataFrame(df_t)
fair_spread_t0.index = ['fair spread']
fair_spread_t = pd.DataFrame(fair_spread_t0.values.T, columns=fair_spread_t0.index, index=fair_spread_t0.columns)
print(fair_spread_t)
fair_spread_t.iplot(
    title = 'Fair Spread of Basket CDS by Sampling form t Copula',
    xTitle = 'N-th to default',
    yTitle = 'Spread'
)
spread_g = []
for k in range(1, 6):
    spread_g.append(fair_spread_MC_Gaussian(k, Corr_Pearson, recovery, risk_free_rate, term_structure, term_structure3, 10000))

df_t = {
    "1st": [spread_g[0]],
    "2nd": [spread_g[1]],
    "3rd": [spread_g[2]],
    "4th": [spread_g[3]],
    "5th": [spread_g[4]]
}     

fair_spread_g0 = pd.DataFrame(df_t)
fair_spread_g0.index = ['fair spread']
fair_spread_g = pd.DataFrame(fair_spread_g0.values.T, columns=fair_spread_g0.index, index=fair_spread_g0.columns)
print(fair_spread_g)
fair_spread_g.iplot(
    title = 'Fair Spread of Basket CDS by Sampling form Gaussian Copula',
    xTitle = 'N-th to default',
    yTitle = 'Spread'
)<jupyter_output><empty_output><jupyter_text>Result interpretationGiven 4 degree of freedom, recovery rate equals 40% and iterate 10000 times, the above tables and plots can be delivered. We can see that though the 2 structures of fair spreads are similar, when $k > 1$, the result derived by t-Copula is slightly higher than the result of Gaussian, that is because t-Copula implies stronger co-movement, so multiple defaults together raise the spread of k-th to default compared to Gaussian results.  Sensitivity Analysis<jupyter_code># Vary the Recovery Rate
# 10%
spread_g10 = []
for k in range(1, 6):
    spread_g10.append(fair_spread_MC_Gaussian(k, Corr_Pearson, 0.1, risk_free_rate, term_structure, term_structure3, 10000))

df_t10 = {
    "1st": [spread_g10[0]],
    "2nd": [spread_g10[1]],
    "3rd": [spread_g10[2]],
    "4th": [spread_g10[3]],
    "5th": [spread_g10[4]]
}     

fair_spread_g10 = pd.DataFrame(df_t10)
fair_spread_g10.index = ['fair spread']
fair_spread_g_10t = pd.DataFrame(fair_spread_g10.values.T, columns=fair_spread_g10.index, index=fair_spread_g10.columns)
print(fair_spread_g_10t)
# 20%
spread_g20 = []
for k in range(1, 6):
    spread_g20.append(fair_spread_MC_Gaussian(k, Corr_Pearson, 0.2, risk_free_rate, term_structure, term_structure3, 10000))

df_t20 = {
    "1st": [spread_g20[0]],
    "2nd": [spread_g20[1]],
    "3rd": [spread_g20[2]],
    "4th": [spread_g20[3]],
    "5th": [spread_g20[4]]
}     

fair_spread_g20 = pd.DataFrame(df_t20)
fair_spread_g20.index = ['fair spread']
fair_spread_g_20t = pd.DataFrame(fair_spread_g20.values.T, columns=fair_spread_g20.index, index=fair_spread_g20.columns)
print(fair_spread_g_20t)
# 30%
spread_g30 = []
for k in range(1, 6):
    spread_g30.append(fair_spread_MC_Gaussian(k, Corr_Pearson, 0.3, risk_free_rate, term_structure, term_structure3, 10000))

df_t30 = {
    "1st": [spread_g30[0]],
    "2nd": [spread_g30[1]],
    "3rd": [spread_g30[2]],
    "4th": [spread_g30[3]],
    "5th": [spread_g30[4]]
}     

fair_spread_g30 = pd.DataFrame(df_t30)
fair_spread_g30.index = ['fair spread']
fair_spread_g_30t = pd.DataFrame(fair_spread_g30.values.T, columns=fair_spread_g30.index, index=fair_spread_g30.columns)
print(fair_spread_g_30t)
# 50%
spread_g50 = []
for k in range(1, 6):
    spread_g50.append(fair_spread_MC_Gaussian(k, Corr_Pearson, 0.5, risk_free_rate, term_structure, term_structure3, 10000))

df_t50 = {
    "1st": [spread_g50[0]],
    "2nd": [spread_g50[1]],
    "3rd": [spread_g50[2]],
    "4th": [spread_g50[3]],
    "5th": [spread_g50[4]]
}     

fair_spread_g50 = pd.DataFrame(df_t50)
fair_spread_g50.index = ['fair spread']
fair_spread_g_50t = pd.DataFrame(fair_spread_g50.values.T, columns=fair_spread_g50.index, index=fair_spread_g50.columns)
print(fair_spread_g_50t)
# 60%
spread_g60 = []
for k in range(1, 6):
    spread_g60.append(fair_spread_MC_Gaussian(k, Corr_Pearson, 0.6, risk_free_rate, term_structure, term_structure3, 10000))

df_t60 = {
    "1st": [spread_g60[0]],
    "2nd": [spread_g60[1]],
    "3rd": [spread_g60[2]],
    "4th": [spread_g60[3]],
    "5th": [spread_g60[4]]
}     

fair_spread_g60 = pd.DataFrame(df_t60)
fair_spread_g60.index = ['fair spread']
fair_spread_g_60t = pd.DataFrame(fair_spread_g60.values.T, columns=fair_spread_g60.index, index=fair_spread_g60.columns)
print(fair_spread_g_60t)<jupyter_output>fair spread
1st    13.159240
2nd     2.484906
3rd     0.729050
4th     0.080666
5th     0.016034<jupyter_text>Sensitivity Against Recovery RateCredit quality projects an important influence on credit spreads. We use the Gaussian Copula, keep every other variables constant, vary the recovery rate from 10% to 60% to have the following table, we can see that according to the increase of recovery rate, fair spreads decreases significantly, and vice versa.Moreover, the impact of credit quality can be seen more significantly on k-th to default Basket CDS with higher k.<jupyter_code>vary_recovery = pd.concat([fair_spread_g_10t, fair_spread_g_20t, fair_spread_g_30t, fair_spread_g, fair_spread_g_50t, fair_spread_g_60t], axis=1, ignore_index=True)
vary_recovery.columns = ['10%', '20%', '30%', '40%', '50%', '60%']
print(vary_recovery)<jupyter_output>10%        20%        30%        40%        50%        60%
1st  30.701916  24.709645  22.162621  19.742750  16.166147  13.159240
2nd   6.154454   5.503114   4.434660   4.096124   2.987338   2.484906
3rd   1.458160   1.392620   1.275903   0.948035   0.870211   0.729050
4th   0.217806   0.225347   0.282211   0.168970   0.080446   0.080666
5th   0.036051   0.032048   0.028177   0.024342   0.020068   0.016034<jupyter_text>Sensitivity Against CorrelationWe can shift the correlation between any 2 names by a percentage, or by a level, say 0.05, then:$$New \ Cov(i,j) = Cov(i,j) \times (1+level), \ i \neq j$$while the elements on the diagonal remains 1. Set levels to be $\pm 0.05$, $\pm 0.1$, set all other variables constant, use shifted Kendall's tau correlation matrix to do sensitivity analysis.<jupyter_code># Shift Correlation Matrix
def correlation_shift(Corr, level):
    for i in range(len(Corr)):
        for j in range(len(Corr)):
            if i != j:
                Corr[i,j] = Corr[i,j]*(1+level)
    return Corr            
# Shifted Correlation Matrix
cor_k_5p = correlation_shift(Corr_Kendall, 0.05)
cor_k_5m = correlation_shift(Corr_Kendall, -0.05)
cor_k_10p = correlation_shift(Corr_Kendall, 0.1)
cor_k_10m = correlation_shift(Corr_Kendall, -0.1)
# +5%
spread_t5p = []
for k in range(1, 6):
    spread_t5p.append(fair_spread_MC_t(k, 4, cor_k_5p, 0.4, risk_free_rate, term_structure, term_structure3, 10000))

df_t5p = {
    "1st": [spread_t5p[0]],
    "2nd": [spread_t5p[1]],
    "3rd": [spread_t5p[2]],
    "4th": [spread_t5p[3]],
    "5th": [spread_t5p[4]]
}     

fair_spread_t5p = pd.DataFrame(df_t5p)
fair_spread_t5p.index = ['fair spread']
fair_spread_t5p_t = pd.DataFrame(fair_spread_t5p.values.T, columns=fair_spread_t5p.index, index=fair_spread_t5p.columns)
print(fair_spread_t5p_t)
# level = -5%
spread_t5m = []
for k in range(1, 6):
    spread_t5m.append(fair_spread_MC_t(k, 4, cor_k_5m, 0.4, risk_free_rate, term_structure, term_structure3, 10000))

df_t5m = {
    "1st": [spread_t5m[0]],
    "2nd": [spread_t5m[1]],
    "3rd": [spread_t5m[2]],
    "4th": [spread_t5m[3]],
    "5th": [spread_t5m[4]]
}     

fair_spread_t5m = pd.DataFrame(df_t5m)
fair_spread_t5m.index = ['fair spread']
fair_spread_t5m_t = pd.DataFrame(fair_spread_t5m.values.T, columns=fair_spread_t5m.index, index=fair_spread_t5m.columns)
print(fair_spread_t5m_t)
# level = +10%
spread_t10p = []
for k in range(1, 6):
    spread_t10p.append(fair_spread_MC_t(k, 4, cor_k_10p, 0.4, risk_free_rate, term_structure, term_structure3, 10000))

df_t10p = {
    "1st": [spread_t10p[0]],
    "2nd": [spread_t10p[1]],
    "3rd": [spread_t10p[2]],
    "4th": [spread_t10p[3]],
    "5th": [spread_t10p[4]]
}     

fair_spread_t10p = pd.DataFrame(df_t10p)
fair_spread_t10p.index = ['fair spread']
fair_spread_t10p_t = pd.DataFrame(fair_spread_t10p.values.T, columns=fair_spread_t10p.index, index=fair_spread_t10p.columns)
print(fair_spread_t10p_t)
# level = -10%
spread_t10m = []
for k in range(1, 6):
    spread_t10m.append(fair_spread_MC_t(k, 4, cor_k_10m, 0.4, risk_free_rate, term_structure, term_structure3, 10000))

df_t10m = {
    "1st": [spread_t10m[0]],
    "2nd": [spread_t10m[1]],
    "3rd": [spread_t10m[2]],
    "4th": [spread_t10m[3]],
    "5th": [spread_t10m[4]]
}     

fair_spread_t10m = pd.DataFrame(df_t10m)
fair_spread_t10m.index = ['fair spread']
fair_spread_t10m_t = pd.DataFrame(fair_spread_t10m.values.T, columns=fair_spread_t10m.index, index=fair_spread_t10m.columns)
print(fair_spread_t10m_t)<jupyter_output>fair spread
1st    17.669866
2nd     4.525264
3rd     1.581462
4th     0.364116
5th     0.024241<jupyter_text>Result InterpretationFrom the table below, we can observe that the spread of 1st default decreases when correlation is shifted by -5% and -10%, it slightly goes up when level equals +5%, then goes down when level reaches 10%. Other spreads moves similarly corresponding to the change of shift level except 4-th to default basket CDS, its spread goes up as level increases. We can see that the variation of spreads is not large, if we consider the change of correlation as the change of original data, then this result corresponds to the property of rank correlation that rank correlation matrix is unlikely to vary as we transform data, they are robust estimates of dependence structure.<jupyter_code>vary_corr = pd.concat([fair_spread_t10m_t, fair_spread_t5m_t, fair_spread_t, fair_spread_t5p_t, fair_spread_t10p_t], axis=1, ignore_index=True)
vary_corr.columns = ['-10%', '-5%', 'original', '5%', '10%']
print(vary_corr)<jupyter_output>-10%        -5%   original         5%        10%
1st  17.669866  17.710303  17.988717  18.310303  18.082763
2nd   4.525264   3.998723   4.694135   5.160640   4.053011
3rd   1.581462   1.338160   1.366226   1.926123   1.390401
4th   0.364116   0.533812   0.363950   0.388945   0.582249
5th   0.024241   0.048383   0.072506   0.048460   0.072430
{'fraction_non_alphanumeric': 0.11562135033616364, 'fraction_numerical': 0.07114674410529971, 'mean_word_length': 4.322916666666667, 'pattern_counts': {'":': 55, '<': 59, '<?xml version=': 0, '>': 64, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '270902', 'license': 'MIT', 'n_tokens_mistral': 12705, 'n_tokens_neox': 11006, 'n_words': 3143, 'path': 'Basket CDS.ipynb', 'repo_name': 'Newtype07/Basket-CDS-Pricing', 'size': '31684'}