Calculating a Company’s Probability of Default with Python

For this analysis, we use several Python-based scientific computing technologies along with the AlphaWave Data Stock Analysis API. Jupyter Notebooks detailing this analysis are also available on Google Colab and Github.

import time
import requests
import selenium
import numpy as np
import pandas as pd
from sympy import *
from datetime import date
from datetime import timedelta
from selenium import webdriver
import pandas_datareader as dr
import matplotlib.pyplot as plt
from datetime import datetime as dt
from bs4 import BeautifulSoup as bs
from selenium.webdriver.common.by import By
from selenium.webdriver.support.ui import Select
from selenium.webdriver.chrome.options import Options
from selenium.webdriver.support.ui import WebDriverWait
from selenium.webdriver.support import expected_conditions as EC

We will use SM Energy (SM) in the subsequent code; however, this analysis works with the ticker of any large publicly traded company.

# Required
company_ticker = 'SM'# Optional
company_name = 'SM Energy'# Optional Input Choices:
# ALL, Annual, Anytime, Bi-Monthly, Monthly, N/A, None, Pays At Maturity, Quarterly, Semi-Annual, Variable
coupon_frequency = 'Semi-Annual'

The probability of default is the probability that a bond issuer will not meet its contractual obligations on schedule. Although the most common event of default is nonpayment leading to bankruptcy proceedings, the bond prospectus might identify other events of default, such as the failure to meet a different obligation or the violation of a financial covenant.

In the following example, we will determine the probability of default given corporate bond prices. The default probabilities that are reached in this exercise are called market-implied default probabilities. Historically, practitioners have focused on the one-year probability of default calculation. We know that over shorter horizons of one or two years firms are exposed to the credit cycle effect, while over longer horizons, the business cycle effect tends to have a lesser impact and the company’s capital structure becomes more important. This effect has made long-run risk levels less cyclical and more stable. Intuitively, default risk over a longer time period is less sensitive to the instantaneous default rates in the economy. For this reason, we will focus on corporate bonds with one or two years until maturity to calculate the market-implied default probabilities.

We will verify the accuracy of the market-implied default probabilities with the Standard & Poor’s Average One-Year Transition Rates For Global Corporates using historical data from 1981–2019. This transition matrix shows the observed historical probabilities of a particular rating transitioning to another rating, including default, over the course of one year.

In order to calculate the market-implied default probabilities, we must first acquire the company’s current bond prices. Using a short Selenium script that emulates a user’s keystrokes and clicks in a browser as a means of navigating to FINRA’s TRACE bond data, we can access the data needed to calculate the market-implied default probabilities.

The following is an example script. In case you do not have Selenium installed, you can visit their respective links and download them using pip in your terminal. We will also need a chromedriver (the simulated chrome browser Selenium controls) and to download it using Python you can use the webdriver-manager package also found in PyPi.

You will need to insert your own path to your chromedriver in the code block below.

# Selenium script
chrome_options = webdriver.ChromeOptions()
chrome_options.add_argument('--headless')
driver = webdriver.Chrome(options=chrome_options, executable_path=r'/PATH/TO/YOUR/chromedriver')# store starting time
begin = time.time()# FINRA's TRACE Bond Center
driver.get('http://finra-markets.morningstar.com/BondCenter/Results.jsp')# click agree
WebDriverWait(driver, 10).until(EC.element_to_be_clickable(
(By.CSS_SELECTOR, ".button_blue.agree"))).click()# click edit search
WebDriverWait(driver, 10).until(EC.element_to_be_clickable(
(By.CSS_SELECTOR, 'a.qs-ui-btn.blue'))).click()# input Issuer Name
WebDriverWait(driver, 10).until(EC.presence_of_element_located(
(By.CSS_SELECTOR, 'input[id=firscreener-issuer]')))
inputElement = driver.find_element_by_id("firscreener-issuer")
inputElement.send_keys(company_name)# input Symbol
WebDriverWait(driver, 10).until(EC.presence_of_element_located(
(By.CSS_SELECTOR, 'input[id=firscreener-cusip]')))
inputElement = driver.find_element_by_id("firscreener-cusip")
inputElement.send_keys(company_ticker)# click advanced search
WebDriverWait(driver, 10).until(EC.element_to_be_clickable(
(By.CSS_SELECTOR, 'a.ms-display-switcher.hide'))).click()# input Coupon Frequency
WebDriverWait(driver, 10).until(EC.presence_of_element_located(
(By.CSS_SELECTOR, 'select[name=interestFrequency]')))
Select((driver.find_elements_by_css_selector(
'select[name=interestFrequency]'))[0]).select_by_visible_text(coupon_frequency)# click show results
WebDriverWait(driver, 10).until(EC.element_to_be_clickable(
(By.CSS_SELECTOR, 'input.button_blue[type=submit]'))).click()# wait for results
WebDriverWait(driver, 10).until(EC.presence_of_element_located(
(By.CSS_SELECTOR, '.rtq-grid-row.rtq-grid-rzrow .rtq-grid-cell-ctn')))# create dataframe from scrape
frames = []
for page in range(1, 11):
bonds = []
WebDriverWait(driver, 10).until(EC.presence_of_element_located(
(By.CSS_SELECTOR, (f"a.qs-pageutil-btn[value='{str(page)}']"))))  # wait for page marker to be on expected page
time.sleep(2)headers = [title.text for title in driver.find_elements_by_css_selector(
'.rtq-grid-row.rtq-grid-rzrow .rtq-grid-cell-ctn')[1:]]tablerows = driver.find_elements_by_css_selector(
'div.rtq-grid-bd > div.rtq-grid-row')
for tablerow in tablerows:
tablerowdata = tablerow.find_elements_by_css_selector(
'div.rtq-grid-cell')
bond = [item.text for item in tablerowdata[1:]]
bonds.append(bond)# Convert to Dataframe
df = pd.DataFrame(bonds, columns=headers)frames.append(df)try:
driver.find_element_by_css_selector('a.qs-pageutil-next').click()
except:
breakbond_prices_df = pd.concat(frames)# store end time 
end = time.time()# total time taken 
print(f"Total runtime of the program is {end - begin} seconds")bond_prices_df

def bond_dataframe_filter(df):
# Drop bonds with missing yields and missing credit ratings
df['Yield'].replace('', np.nan, inplace=True)
df["Moody's®"].replace({'WR':np.nan,'':np.nan}, inplace=True)
df["S&P"].replace({'NR':np.nan,'':np.nan}, inplace=True)
df = df.dropna(subset=['Yield'])
df = df.dropna(subset=["Moody's®"])
df = df.dropna(subset=["S&P"])# Create Maturity Years column that aligns with Semi-Annual Payments from corporate bonds
df['Yield'] = df['Yield'].astype(float)
df['Coupon'] = df['Coupon'].astype(float)
df['Price'] = df['Price'].astype(float)
now = dt.strptime(date.today().strftime('%m/%d/%Y'), '%m/%d/%Y')
df['Maturity'] = pd.to_datetime(df['Maturity']).dt.strftime('%m/%d/%Y')
daystillmaturity = []
yearstillmaturity = []
for maturity in df['Maturity']:
daystillmaturity.append(
(dt.strptime(maturity, '%m/%d/%Y') - now).days)
yearstillmaturity.append(
(dt.strptime(maturity, '%m/%d/%Y') - now).days/360)
df = df.reset_index(drop=True)
df['Maturity'] = pd.Series(daystillmaturity)
#         df['Maturity Years'] = pd.Series(yearstillmaturity).round() # Better for Annual Payments
df['Maturity Years'] = round(pd.Series(yearstillmaturity)/0.5)*0.5 # Better for Semi-Annual Payments# Target bonds with short-term maturities 
df['Maturity'] = df['Maturity'].astype(float)
df = df.loc[df['Maturity'] >= 0]
years_mask = (df['Maturity Years'] > 0) & (df['Maturity Years'] <= 2)
df = df.loc[years_mask]
return dfbond_df_result = bond_dataframe_filter(bond_prices_df)
bond_df_result

def bonds_probability_of_default(coupon, maturity_years, bond_price, principal_payment, risk_adjusted_discount_rate):
price = bond_price
prob_default_exp = 0
#     times = np.arange(1, maturity_years+1) # For Annual Cashflows
#     annual_coupon = coupon # For Annual Cashflows
times = np.arange(0.5, (maturity_years-0.5)+1,0.5) # For Semi-Annual Cashflows
semi_annual_coupon = coupon/2 # For Semi-Annual Cashflows
# Calculation of Expected Cash Flow
cashflows = np.array([])
for i in times[:-1]:
#         cashflows = np.append(cashflows, annual_coupon) # For Annual Cashflows
#     cashflows = np.append(cashflows, annual_coupon+principal_payment)#  For Annual Cashflows
cashflows = np.append(cashflows, semi_annual_coupon) # For Semi-Annual Cashflows
cashflows = np.append(cashflows, semi_annual_coupon+principal_payment) # For Semi-Annual Cashflows
for i in range(len(times)):
#         This code block is used if there is only one payment remaining
if len(times) == 1:
prob_default_exp += ((cashflows[i]*(1-P) + cashflows[i]*recovery_rate*P) / 
np.power((1 + risk_adjusted_discount_rate), times[i]))
#         This code block is used if there are multiple payments remaining
else:
#             For Annual Cashflows
#             if times[i] == 1:
#                 prob_default_exp += ((cashflows[i]*(1-P) + principal_payment*recovery_rate*P) / 
#                                     np.power((1 + risk_adjusted_discount_rate), times[i]))
#             For Semi-Annual Cashflows
if times[i] == 0.5:
prob_default_exp += ((cashflows[i]*(1-P) + principal_payment*recovery_rate*P) / 
np.power((1 + risk_adjusted_discount_rate), times[i]))
#             Used for either Annual or Semi-Annual Cashflows
else:
prob_default_exp += (np.power((1-P), times[i-1])*(cashflows[i]*(1-P) + principal_payment*recovery_rate*P)) / 
np.power((1 + risk_adjusted_discount_rate), times[i])prob_default_exp = prob_default_exp - price    
implied_prob_default = solve(prob_default_exp, P)
implied_prob_default = round(float(implied_prob_default[0]) * 100, 2)
if implied_prob_default < 0:
return 0.0
else:
return implied_prob_default

# Ten-Year Risk-free Rate
timespan = 100
current_date = date.today()
past_date = current_date-timedelta(days=timespan)
ten_year_risk_free_rate_df = dr.DataReader('^TNX', 'yahoo', past_date, current_date) 
ten_year_risk_free_rate = (ten_year_risk_free_rate_df.iloc[len(ten_year_risk_free_rate_df)-1,5])/100
ten_year_risk_free_rate

# Market Risk Premium
market_risk_premium = 0.0472

# Market Equity Beta
stock_market_beta = 1

# Market Rate of Return
market_rate_of_return = ten_year_risk_free_rate + (stock_market_beta*market_risk_premium)
market_rate_of_return

# One-Year Risk-free Rate
one_year_risk_free_rate = (1 + ten_year_risk_free_rate)**(1/10) - 1
one_year_risk_free_rate

# Vanguard Short-Term Corporate Bond Index Fund ETF Shares
bond_fund_ticker = 'VCSH'url = "https://stock-analysis.p.rapidapi.com/api/v1/resources/key-stats"querystring = {"ticker":bond_fund_ticker}headers = {
'x-rapidapi-host': "YOUR_X-RAPIDAPI-HOST_WILL_COPY_DIRECTLY_FROM_RAPIDAPI_PYTHON_CODE_SNIPPETS",
'x-rapidapi-key': "YOUR_X-RAPIDAPI-KEY_WILL_COPY_DIRECTLY_FROM_RAPIDAPI_PYTHON_CODE_SNIPPETS"
}bond_fund_key_stats_response = requests.request("GET", url, headers=headers, params=querystring)# Create Key Statistics DataFrame
bond_fund_key_stats_df = pd.DataFrame.from_dict(bond_fund_key_stats_response.json())
bond_fund_key_stats_df = bond_fund_key_stats_df.transpose()
bond_fund_key_stats_df

# Bond's Beta
bond_beta = float(bond_fund_key_stats_df.loc[r'Beta (5Y monthly)'])
bond_beta

# Expected Risk Premium
expected_risk_premium = (market_rate_of_return - one_year_risk_free_rate) * bond_beta
expected_risk_premium

# One-Year Risk-free Rate
one_year_risk_free_rate = (1 + ten_year_risk_free_rate)**(1/10) - 1
one_year_risk_free_rate

# Risk-adjusted Discount Rate
risk_adjusted_discount_rate = one_year_risk_free_rate + expected_risk_premium
risk_adjusted_discount_rate

# Variables defined for bonds_probability_of_default function
principal_payment = 100
recovery_rate = .40
P = symbols('P')

# This calculation may take some time if there are many coupon payments
bond_df_result['Probability of Default %'] = bond_df_result.apply(lambda row : bonds_probability_of_default(row['Coupon'], row['Maturity Years'], 
row['Price'], principal_payment, risk_adjusted_discount_rate), axis = 1)bond_df_result

def prob_default_term_structure(df):
fig, (ax1, ax2) = plt.subplots(1, 2, clear=True)
fig.subplots_adjust(wspace=.5)
Mgroups = df.groupby("Moody's®")
ax1.clear()
ax1.margins(0.5)
ax1.set_xlabel('Days Until Maturity')
ax1.set_ylabel('Probability of Default %')
ax1.set_title("Moody's® Ratings")
for name, group in Mgroups:
ax1.plot(group['Maturity'], group['Probability of Default %'],
marker='o', linestyle='', ms=12, label=name)
ax1.legend(numpoints=1, loc='upper left')SPgroups = df.groupby("S&P")
ax2.clear()
ax2.margins(0.5)
ax2.set_xlabel('Days Until Maturity')
ax2.set_ylabel('Probability of Default %')
ax2.set_title("S&P Ratings")for name, group in SPgroups:
ax2.plot(group['Maturity'], group['Probability of Default %'],
marker='o', linestyle='', ms=12, label=name)
ax2.legend(numpoints=1, loc='upper left')
pltprob_default_term_structure(bond_df_result)

tgt_website = r'https://www.spglobal.com/ratings/en/research/articles/200429-default-transition-and-recovery-2019-annual-global-corporate-default-and-rating-transition-study-11444862'def get_transition_matrix(tgt_website):
df_list = pd.read_html(tgt_website)
matrix_result_df = df_list[22]
return matrix_result_df
transition_matrix_df = get_transition_matrix(tgt_website)sp_clean_result_df = pd.DataFrame(transition_matrix_df.iloc[:34,:19].dropna(axis=0))
sp_clean_result_df

# Will scrape the default probability for each ratingsp_rating_list = ['AAA','AA+','AA','AA-','A+','A','A-',
'BBB+','BBB','BBB-','BB+','BB','BB-',
'B+','B','B-']ccc_list = ['CCC+','CCC','CCC-','CC+','CC','CC-','C+','C','C-']sp_rating = Nonefor i in sp_rating_list:
if bond_df_result['S&P'].iloc[0] == i:
sp_rating = bond_df_result['S&P'].iloc[0]if sp_rating == None:
for i in ccc_list:
if bond_df_result['S&P'].iloc[0] == i:
sp_rating = 'CCC/C'sp_transition_dp = 0for i in range(33):
if transition_matrix_df.loc[i][0] == sp_rating:
sp_transition_dp += float(sp_clean_result_df.loc[i][18])sp_transition_dp

# Compare the nearest maturity Market-implied probability of default with
# the historical probability of default in the Standard & Poor’s 2019 transition matrix
print("Market-implied probability of default = %s"%(bond_df_result['Probability of Default %'].iloc[0])+"%")
print("Standard & Poor’s probability of default = %s"%(sp_transition_dp)+"%")