import numpy as np import pandas as pd import matplotlib.pyplot as plt import seaborn as sns from scipy.stats import norm import random from scipy.optimize import minimize, rosen, rosen_der from scipy import optimize

v0 = 100 mu = .03 sigma_v = 0.2 T = 5 dt = 1/12 N = 85 random.seed(7) paths_p1 = np.zeros(shape = (7,60)) paths_p1[:,0] = v0 for path in range(0,7): for column in range(1, 60): paths_p1[path, column] = paths_p1[path, column-1] + paths_p1[path, column-1]*mu*dt + sigma_v*paths_p1[path, column-1]*np.random.standard_normal(1)*np.sqrt(dt)

for i in range(0,7): plt.plot(np.arange(0,60,1), paths_p1[i,:]) plt.legend(['path 1', 'path 2', 'path 3', 'path 4', 'path 5', 'path 6', 'path 7'], bbox_to_anchor = (1.01,1), borderaxespad = 0) plt.axhline(y = N, linestyle = '--') plt.title('Asset Value Paths Over a 5 year Horizon') plt.ylabel('Asset Value') plt.xlabel('Time (Months)')

V_0=100 miu=0.03 sigma_v=0.2 T=5 dt=1/12 N=85 def Yield(V,Mu,sigma,T,t,dt,N,r): d1=1/(sigma_v*(np.sqrt(T-t)))*(np.log(V/N)+(r+0.5*sigma**2)*(T-t)) d2=d1-sigma*np.sqrt(T-t) D=N*np.exp(-r*(T-t))-(N*np.exp(-r*(T-t))*norm.cdf(-d2)-V*norm.cdf(-d1)) Yield=1/(T-t)*np.log(N/D) Spread=Yield-r return Spread YieldSpread_1b=np.zeros(61) #YieldSpread_1b[0]=0 for i in range(0,60): month=i/12 YieldSpread_1b[i]=Yield(V=V_0,Mu=miu,sigma=sigma_v,T=5,t=month,dt=dt,N=N,r=0.03) plt.figure() ax= sns.lineplot(data=YieldSpread_1b) ax.set_xlabel('Time To Maturity') ax.set_ylabel('Spread') plt.title('Term Structure of Spreads')

YieldSpread_1c=np.zeros(61) YieldSpread_1c[0]=0 for i in range(0,60): month=i/12 YieldSpread_1c[i]=Yield(V=100,Mu=miu,sigma=0.07,T=5,t=month,dt=dt,N=50,r=0.03) plt.figure() ax= sns.lineplot(data=YieldSpread_1c) ax.set_xlabel('Time To Maturity') ax.set_ylabel('Spread') plt.title('Term Structure of Spreads W/ N=50 and Sigma=7%')

Return=np.array([-0.33333,-0.103448,-0.038462,-0.1112,0.278128,-.176056,-.059829,.090909,-.020833,-.256596,-.029765,-.262537]) Price=np.array([2.9,2.6,2.5,2.222,2.84,2.34,2.2,2.4,2.35,1.747,1.695,1.25]) NumOfShares=np.array([20778000,20778000,20808000,20808000,20808000,20808000,20808000,20808000,20856000,20856000,20856000,23056000]) EquityValue=Price*NumOfShares Volatility_2a=np.std(Return) print("the volatility is:",Volatility_2a)

rf=0.0025 STD=0 LTD=33456000 N=STD+0.5*LTD T = 1 def Merton_Model (V,Mu,Sigma,T,dt,N,r): d1=np.log(V/N)+(r+0.5*sigma**2)*(T)/(sigma*np.sqrt(T)) d2=d1-np.sqrt(T) Equity=V*norm.cdf(d1)-N*np.exp(-r*T)*norm.cdf(d2) return Equity def Black_scholes(V,N=N,Sigma=Volatility_2a,T=T,r=rf): d1=(np.log(V/N)+(r+0.5*Sigma**2)*(T))/(Sigma*np.sqrt(T)) d2=d1-np.sqrt(T) Equity=V*norm.cdf(d1)-N*np.exp(-r*T)*norm.cdf(d2) return Equity def estimates_func(V_est,E=np.mean(EquityValue)): return Black_scholes(V=V_est)-E #

Table_2=pd.DataFrame({'return':Return,"price":Price, "NumOfShares":NumOfShares,"EquityValue":EquityValue}) Table_2 Volatility_old=Volatility_2a Volatility_new=0 for j in range(0,100): #Asset Value #print("test") Asset_value_list=np.zeros(12) for i in range(0,12): def estimates_func(V_est,E=Table_2.EquityValue[i]): return Black_scholes(V=V_est,Sigma=Volatility_old)-E Asset_value_list[i] = optimize.newton(estimates_func, 8000000,maxiter=1000000) #Asset_Value_Return=np.zeros(12) #Volatility_old=Volatility_new Volatility_new=np.std((np.log(Asset_value_list[1:]/Asset_value_list[:-1]))) if(Volatility_old==Volatility_new): # print("Volatility_new",Volatility_new) print("The Volatility is",Volatility_old) Table_2["AssetValue"]=Asset_value_list break Volatility_old=Volatility_new #Asset_value_list #Volatility_old=Volatility_new Table_2

#Average Return Avg_Return_2_3=np.mean(np.log(Table_2.AssetValue.shift(1)/Table_2.AssetValue)) Avg_Return_2_3 #DD Distance_Default=(np.log(1+Avg_Return_2_3)+(0-0.5*Volatility_new**2)*(T))/(Volatility_new*np.sqrt(T)) PD=1-norm.cdf(Distance_Default) print("The distance to Default is:",Distance_Default) print("Probability of default is:",round(PD*100,2),"%")

Boeing = pd.read_csv("Boeing.csv")

Boeing

def A2(stock_vol, stock_price, gdr, debt_share, lam, t): return np.sqrt((stock_vol*(stock_price/(stock_price+gdr*debt_share))**2)*t + lam**2) def d(init_price, gdr, debt_share, L, lam): return ((init_price + gdr*debt_share)/(L*debt_share))*np.exp(lam**2) def survival_prob(stock_vol, stock_price, init_price, gdr, debt_share, lam, L, t): At = A2(stock_vol, stock_price, gdr, debt_share, lam, t) d1 = d(init_price, gdr, debt_share, L, lam) return norm.cdf(-(At/2)+(np.log(d1)/At))-d1*norm.cdf(-(At/2)-(np.log(d1)/At)) def credit_spread(r, R, stock_vol, stock_price, init_price, gdr, debt_share, lam, L, t): z = np.sqrt(0.25 + 2*(r/stock_vol**2)) d1 = d(init_price, gdr, debt_share, L, lam) s = (lam**2)*(stock_vol)**2 PS0 = survival_prob(stock_vol, init_price, init_price, gdr, debt_share, lam, L, t) PST = survival_prob(stock_vol, stock_price, init_price, gdr, debt_share, lam, L, t) def G(u, z, d1, stock_vol): term1 = (np.power(d1,z+0.5))*norm.cdf(-np.log(d1)/(stock_vol*np.sqrt(u))-z*stock_vol*np.sqrt(u)) term2 = (np.power(d1,-z+0.5))*norm.cdf(-np.log(d1)/(stock_vol*np.sqrt(u))-z*stock_vol*np.sqrt(u)) return term1 + term2 numerator = 1-PS0+(np.exp(-r*s)*(G(t+s, z, d1, stock_vol)-G(s, z, d1, stock_vol))) denominator = PS0 - PST*np.exp(-r*t) - (np.exp(r*s)*(G(t+s, z, d1, stock_vol)-G(s, z, d1, stock_vol))) return r*(1-R)*(numerator/denominator) credit_spread(Boeing['r'], Boeing['R'], Boeing['sigma_stock'], Boeing['S'], 62, 0.75, Boeing['D'], 0.3, Boeing['L'],5)

survival_prob(Boeing['sigma_stock'], Boeing['S'], 62, 0.75, Boeing['D'], 0.3, Boeing['L'],5)

import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy.stats import norm from datetime import datetime def CreditGrades_CDS(df_t): """ Calculates the credit default spread using the CreditGrades model from JPMorgan :param df_t: a SINGLE dataframe row w/ the args needed to find survival proba under the CG model. :param N_sims: number of simulations to perform for each element. :return: CDS of the firm at a given time 't'. """ # These terms are fixed at 't': S_0, Lbar, D, Lambda, r, R = df_t[["S", "L", "D", "lambda", "r", "R"]] LD = Lbar * D """ We would like to obtain 'd' (which is constant for a given expectation of a stochastic default barrier). """ d = (S_0 + Lbar * D) / LD * np.exp(Lambda**2) """ Closed form solution for survival probability utilizing CDF of a normal distribution """ T, Sigma_S = df_t[["T", "sigma_stock"]] A_0 = np.sqrt(Sigma_S * (S_0 / (S_0 + Lbar * D))**2 * 0 + Lambda**2) A_T = np.sqrt(Sigma_S * (S_0 / (S_0 + Lbar * D))**2 * T + Lambda**2) P_0 = norm.cdf(-A_0/2 + np.log(d)/A_0) - d * norm.cdf(-A_0/2 - np.log(d)/A_0) P_T = norm.cdf(-A_T/2 + np.log(d)/A_T) - d * norm.cdf(-A_T/2 - np.log(d)/A_T) """ Abstracting away technical details on page 46 of technical document """ Sigma_V = Sigma_S * S_0 / (S_0 + LD) # Asset volatility relation documented on pg. 11 of technical paper z = np.sqrt(1/4 + 2*r / Sigma_V**2) Epsilon = Lambda**2 / Sigma_V**2 G_TEpsilon = \ d**(z+0.5) * \ norm.cdf(-np.log(d) / (Sigma_V * np.sqrt(T+Epsilon)) - \ z * Sigma_V * np.sqrt(T+Epsilon)) + \ d**(-z+0.5) * \ norm.cdf(-np.log(d) / (Sigma_V * np.sqrt(T+Epsilon)) + \ z * Sigma_V * np.sqrt(T+Epsilon)) G_Epsilon = \ d**(z+0.5) * \ norm.cdf(-np.log(d) / (Sigma_V * np.sqrt(Epsilon)) - \ z * Sigma_V * np.sqrt(Epsilon)) + \ d**(-z+0.5) * \ norm.cdf(-np.log(d) / (Sigma_V * np.sqrt(Epsilon)) + \ z * Sigma_V * np.sqrt(Epsilon)) """ Finally combine terms to find the credit default spread implied at time 't' """ CDS_t = \ r * (1-R) * \ (1 - P_0 + np.exp(r*Epsilon) * (G_TEpsilon - G_Epsilon)) /\ (P_0 - P_T * np.exp(-r*T) - np.exp(r*Epsilon) * (G_TEpsilon - G_Epsilon)) return CDS_t Boeing = pd.read_csv("Boeing.csv") Boeing["date"] = pd.to_datetime(Boeing["date"]) # Reformat dates to datetime objects Boeing["cds_model"] = Boeing.apply(CreditGrades_CDS, axis=1) Boeing["cds_model"] = Boeing["cds_model"] / 0.001 # Convert CDS model values to % """ Plotting results of CDS Model implied values against 5y Spreads from data Note: Our model would indicate we would want to buy CDS when the area is green and Sell CDS when area is red. """ plt.plot(Boeing["date"], Boeing["cds_model"], label="5y (Model)") plt.plot(Boeing["date"], Boeing["spread5y"], label="5y (Real)") plt.fill_between(Boeing["date"], Boeing["cds_model"], Boeing["spread5y"], where=Boeing["cds_model"] > Boeing["spread5y"], interpolate=True, color="green", alpha=0.35) plt.fill_between(Boeing["date"], Boeing["cds_model"], Boeing["spread5y"], where=Boeing["cds_model"] <= Boeing["spread5y"], interpolate=True, color="red", alpha=0.35) plt.legend() plt.title("Boeing 5y CDS Spread (CG Model vs. Real)") plt.xlabel("Date") plt.ylabel("Spread (%)") plt.show()

import statsmodels.api as sm def CreditGrades_RPV01(df_t): """ :param df_t: a SINGLE row from the DataFrame at time 't' :return: the value of a risky annuity given args """ # These terms are fixed at 't': S_0, Lbar, D, Lambda, r, R = df_t[["S", "L", "D", "lambda", "r", "R"]] LD = Lbar * D """ We would like to obtain 'd' (which is constant for a given expectation of a stochastic default barrier). """ d = (S_0 + Lbar * D) / LD * np.exp(Lambda**2) """ Closed form solution for survival probability utilizing CDF of a normal distribution """ T, Sigma_S = df_t[["T", "sigma_stock"]] A_0 = np.sqrt(Sigma_S * (S_0 / (S_0 + Lbar * D))**2 * 0 + Lambda**2) A_T = np.sqrt(Sigma_S * (S_0 / (S_0 + Lbar * D))**2 * T + Lambda**2) P_0 = norm.cdf(-A_0/2 + np.log(d)/A_0) - d * norm.cdf(-A_0/2 - np.log(d)/A_0) P_T = norm.cdf(-A_T/2 + np.log(d)/A_T) - d * norm.cdf(-A_T/2 - np.log(d)/A_T) Sigma_V = Sigma_S * S_0 / (S_0 + LD) # Asset volatility relation documented on pg. 11 of technical paper z = np.sqrt(1/4 + 2*r / Sigma_V**2) Epsilon = Lambda**2 / Sigma_V**2 G_TEpsilon = \ d**(z+0.5) * \ norm.cdf(-np.log(d) / (Sigma_V * np.sqrt(T+Epsilon)) - \ z * Sigma_V * np.sqrt(T+Epsilon)) + \ d**(-z+0.5) * \ norm.cdf(-np.log(d) / (Sigma_V * np.sqrt(T+Epsilon)) + \ z * Sigma_V * np.sqrt(T+Epsilon)) G_Epsilon = \ d**(z+0.5) * \ norm.cdf(-np.log(d) / (Sigma_V * np.sqrt(Epsilon)) - \ z * Sigma_V * np.sqrt(Epsilon)) + \ d**(-z+0.5) * \ norm.cdf(-np.log(d) / (Sigma_V * np.sqrt(Epsilon)) + \ z * Sigma_V * np.sqrt(Epsilon)) rpv01_t = \ (P_0 - P_T * np.exp(-r*T) - np.exp(r*Epsilon) * (G_TEpsilon - G_Epsilon)) /\ (r) return rpv01_t def CreditGrades_sensitivity(df): """ Finds the partial deriv of CDS w.r.t. partial deriv. of stock price """ X = df["cds_model"] y = df["S"] model = sm.OLS(y, X).fit() sensitivity = model.params[0] # P. Deriv(CDS Model) / P. Deriv(Stock Price) return sensitivity """ Compute Boeing CDS hedge ratio """ Boeing_sensitivity = CreditGrades_sensitivity(Boeing) # <- This is the hedge ratio Boeing["$hedge"] = Boeing.apply(CreditGrades_RPV01, axis=1) * Boeing_sensitivity Boeing["$hedge"]

""" Testing a trading signal """ alpha = 3.0 # Trade if divergence >= 300% of 1y rolling avg spread Boeing["1y_AvgDiff"] = (Boeing["cds_model"]-Boeing["spread5y"]).rolling(250).mean() Boeing["signal"] = 0 Boeing["signal"] = np.where(Boeing["cds_model"]-Boeing["spread5y"] \ >= Boeing["1y_AvgDiff"]*alpha, 1, Boeing["signal"]) # Buy signal Boeing["signal"] = np.where(Boeing["cds_model"]-Boeing["spread5y"] \ >= Boeing["1y_AvgDiff"]/alpha, -1, Boeing["signal"]) # Sell signal # Backtesting out strategy - find our P&L (in %) if we held each signal position for 1 month. c = 0.0 Boeing["CDS Price"] = (Boeing["spread5y"] - c) / (Boeing["r"] + (Boeing["spread5y"] / (1 - Boeing["R"]))) * \ (1 - np.exp(-(Boeing["r"] + (Boeing["S"] / (1 - Boeing["R"]))) * Boeing["T"])) Position_Length = 30 # We hold each position for 1-month; Approx. 30 days Boeing["P&L"] = 0.0 Boeing["P&L"] = np.where(Boeing["signal"] == 1, (Boeing["CDS Price"].shift(-Position_Length)+Boeing["CDS Price"])/Boeing["CDS Price"] +\ (Boeing["$hedge"].shift(-Position_Length)+Boeing["$hedge"])/Boeing["$hedge"], Boeing["P&L"]) Boeing["P&L"] = np.where(Boeing["signal"] == -1, (-Boeing["CDS Price"].shift(-Position_Length)-Boeing["CDS Price"])/-Boeing["CDS Price"] +\ (-Boeing["$hedge"].shift(-Position_Length)-Boeing["$hedge"])/-Boeing["$hedge"], Boeing["P&L"]) Boeing_plotted = Boeing.dropna() plt.scatter(Boeing_plotted["date"], Boeing_plotted["P&L"], s=0.8, alpha=0.75, color="red") plt.title("P&L in % of Monthly CSA strategy (Alpha=3.0)") plt.xlabel("Date") plt.ylabel("P&L (%)") plt.show() # Used for analytics below avgMonthlyRet = np.mean(Boeing_plotted["P&L"]) avgAnnualRet = avgMonthlyRet*12

print("Average Monthly-Lagged Return on Traded Days: ", avgMonthlyRet, "%")