In this lesson, we focus on the idea of how credit ratings on bonds can transition. First, let’s look at a portfolio where there are 10 bonds rated AAA, 10 rated AA, and 10 rated A.

import pandas as pd
import numpy as np

portfolio = pd.Series(["AAA"]*10+["AA"]*10+["A"]*10)
print(portfolio)

Credit migration can be modeled as a markov chain which means that we can think of the only important factor being the current state for predicting the state in the future. Below is code to create a markov chain that has empirical estimations of transition probabilities. The way to think of the table is that each row shows the current state, and each column shows the probability of transitioning to a state. The row AAA, and column AAA has a probability of 90.81% meaning the odds of a bond going from state AAA in time t to state AAA in time t+1 is 90.81%. On the other hand, the row AAA and column BB only has .12% meaning the chance of transitioning is extremely low for one time period. It does not however, mean that after 10 time steps the probability will be .12%. To get the probabilities for state AAA, we can run the following:

#Probabilities for state AAA
print(transition_matrix.loc['AAA'])

The numpy function choice in the random module allows for picking a random choice based on a given probability distribution. We give the function a list of labels followed by setting p equal to the probabilities for each variable. We can see 10 random choices for the AAA state followed by 10 random choices for the A state.

#Get the random move
for _ in range(10):
    print(np.random.RandomState().choice(transition_matrix.index, p=transition_matrix.loc["AAA"].values))
print()
print()
for _ in range(10):
    print(np.random.RandomState().choice(transition_matrix.index, p=transition_matrix.loc["A"].values))

Let’s set the initial conditions. We are going to use the default variable to track the number of defaults.

#Set initial portfolio of bond's with credit ratings
initial_state = {
    'portfolio': portfolio,
    'default': 0
}

Let’s start by building out the functions for credit migration. We can use the map function on the series to map the same random choice function we used above. We are setting the portfolio variable each time.

#Function to transition the matrix randomly
def credit_migration(_params, substep, sH, s):
    new_ratings = s['portfolio'].map(lambda x: np.random.RandomState().choice(
  transition_matrix.index, 
  p=transition_matrix.loc[x].values))
    return {"ratings": new_ratings}

#Update the ratings    
def update_portfolio(_params, substep, sH, s, _input):
    return ('portfolio', _input['ratings'])

For this simulation we need a second set which will be used to find the number of defaults after the migrations. The functions should be easy to understand, it is just finding the proportion of bonds in the portfolio that have transitioned to the default state.

#Function to value the portfolio, for now just finding default %
def value_portfolio(_params, substep, sH, s):
    return {"default": (s['portfolio'] == "Default").sum() / len(s['portfolio'])}

#Update the value
def update_value(_params, substep, sH, s, _input):
    return ('default', _input['default'])

With these four functions defined, we are building out the update blocks by splitting into two pieces. The first piece will be the credit migration step, and the second one will be the valuation. They will run sequentially.

#Create partial updates
PSUBs = [
    {
        "policies": {
            "migration_policy": credit_migration,
        },
        "variables": {
            "portfolio": update_portfolio,
        }
    },
        {
        "policies": {
            "valuation_policy": value_portfolio,
        },
        "variables": {
            "default": update_value,
        }
    }
]

Create and run the experiment.

#Set parameters for 50 runs over 30 years
from cadCAD.configuration.utils import config_sim

sim_config_dict = {
    'T': range(30),
    'N': 50}
c = config_sim(sim_config_dict)

#Create experiment
from cadCAD.configuration import Experiment

exp = Experiment()
exp.append_configs(
    initial_state = initial_state,
    partial_state_update_blocks = PSUBs,
    sim_configs = c
)

from cadCAD.engine import ExecutionMode, ExecutionContext, Executor
from cadCAD import configs

#Execute!
exec_mode = ExecutionMode()
local_mode_ctx = ExecutionContext(context=exec_mode.single_mode)
simulation = Executor(exec_context=local_mode_ctx, configs=configs)
raw_result, field, sessions = simulation.execute()
result = pd.DataFrame(raw_result)

Since we now have two substeps, the substeps column has 1 and 2 depending on which step in the process is recorded. All steps get recorded when we do this. Let’s grab only records that have substep equal to 2 (the last of the substeps) and find the evolution of the default variable for the first 5 simulations.

import matplotlib.pyplot as plt

default_rates = result[result['substep'] == 2].pivot('timestep', 'run', 'default')
default_rates.iloc[:,:10].plot(kind='line')
plt.xlabel("Year")
plt.ylabel("Perfect Defaulted")
plt.title("Credit Migration")
plt.show()

Get the average default rate across the simulations.

#Find the average default rate across simulations
result[result['substep'] == 2].groupby("timestep")['default'].mean().plot(kind='line')
plt.xlabel("Year")
plt.ylabel("Percent Defaulted")
plt.title("Credit Migration (Average)")
plt.show()

One nice property of markov chains is that by starting with a vector and doing the dot product with a markov chain we can find the next period predicted states vector. Don’t worry too much about the code below, but it will find the expected states. From there it is easy to find the predicted default rate.

#Find the predicted default rate based on markov chain
transition_matrices = [pd.DataFrame(np.diag([1]*8),index=transition_matrix.index, columns=transition_matrix.columns)]
for _ in range(30):
    transition_matrices.append(transition_matrices[-1].dot(transition_matrix))
    
predicted = pd.Series([1/3, 1/3, 1/3, 0, 0, 0, 0, 0], index=transition_matrix.index)
predicted = pd.Series([predicted.dot(x)['Default'] for x in transition_matrices],index=list(range(31)))
print(predicted)

This will be good for comparing with our simulated values to make sure they are feasible.

result[result['substep'] == 2].groupby("timestep")['default'].mean().plot(kind='line')
predicted.plot(kind='line')
plt.xlabel("Year")
plt.ylabel("Percent Defaulted")
plt.title("Credit Migration (Average)")
plt.legend(['Simulated', 'Predicted'])
plt.show()