Let’s reset the simulation and run 100 simulations.

pathos.helpers.shutdown()
configs.clear()

sim_config_dict = {
    'T': range(360),
    'N': 100,
"M": {'a': [.5], 'b':[.04], 'sigma':[.02], 't':[1/12]}
}
c = config_sim(sim_config_dict)

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

#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)

With seaborn’s distplot we get a nice estimate of the shape of the distribution from our 100 simulations.

#Plot terminal cash distribution
import seaborn as sns

sns.distplot(result[result['timestep'] == 360]['cash'], kde=True)
plt.show()

If we assume that the bonds cost a total of $30,000 we find the yield the normal way and then can plot the distribution of yields from this portfolio. This is useful when we want to be able to understand bounds around what yield we can get from an investment.

#Plot yields distribution

yields = result[(result['timestep'] == 360) & (result['substep'] == 5)]['cash']
yields = (yields / 30000) ** (1/30) - 1
yields.name = 'Yield'
sns.distplot(yields, kde=True)
plt.show()

We could also find a 95% confidence interval around the cash position.

#Find the bounds for two standard deviations around cash

cash = result[result['substep'] == 5].pivot("timestep", "run", "cash")
mu = cash.mean(axis=1)
std = cash.std(axis=1)
lower_bound = mu - 2 * std
upper_bound = mu + 2 * std

mu.plot(kind='line', color='black')
lower_bound.plot(kind='line', color='red', linestyle='--')
upper_bound.plot(kind='line', color='red', linestyle='--')
plt.show()

For our final analysis, let’s compare two investments which cost $30,000. One with higher quality, lower coupon bonds and one with lower quality higher coupon bounds. Below are the definitions of the portfolios involved. The question we want to answer is what are the differences between investing in these two portfolios from a yield point of view not only in terms of the expected yield but also the possible variations of yield.

import pandas as pd
import numpy as np

#Create the portfolio
portfolio1 = pd.DataFrame(["AAA"]*10+["AA"]*10+["A"]*10, columns = ['Rating'])
portfolio2 = pd.DataFrame(["BBB"]*10+["BB"]*10+["B"]*10, columns = ['Rating'])


portfolio1['Coupon Rate'] = [.05] * 10 + [.06] * 10 + [.07] * 10
portfolio1['Payment Frequency'] = [6,12] * 15
portfolio1['Face Value'] = 1000

portfolio2['Coupon Rate'] = [.08] * 10 + [.09] * 10 + [.1] * 10
portfolio2['Payment Frequency'] = [6,12] * 15
portfolio2['Face Value'] = 1000

#Set initial portfolios
initial_state1 = {
    'portfolio': portfolio1,
    'default': 0,
    'r': .06,
    'cash': 0,
    'payment': 0
}

initial_state2 = {
    'portfolio': portfolio2,
    'default': 0,
    'r': .06,
    'cash': 0,
    'payment': 0
}

Now we have two defined initial states. We are going to create two experiments, one for each initial setting. Then we will run the simulation with these two initial states.

pathos.helpers.shutdown()
configs.clear()

sim_config_dict1 = {
    'T': range(360),
    'N': 100,
"M": {'a': [.5], 'b':[.04], 'sigma':[.02], 't':[1/12]}
}
c1 = config_sim(sim_config_dict1)

sim_config_dict2 = {
    'T': range(360),
    'N': 100,
"M": {'a': [.5], 'b':[.04], 'sigma':[.02], 't':[1/12]}
}
c2 = config_sim(sim_config_dict2)

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

exp = Experiment()

exp.append_configs(
    initial_state = initial_state2,
    partial_state_update_blocks = PSUBs,
    sim_configs = c2
)

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)

By first indexing only where the timestep is 360 (the final state) and substate is 5, we get 200 rows, 100 for each simulation’s ending cash. Next we convert it to a yield, and finally, we plot the two distribution plots. It should be obvious that the spread of simulation 2’s yield is much larger because the low quality bonds have much more randomness/chance to default.

yields = result[(result['timestep'] == 360) & (result['substep'] == 5)].copy()
yields['Yield'] = (yields['cash'] / 30000) ** (1/30) - 1


sns.distplot(yields[yields['simulation'] == 0]['Yield'], kde=True, hist=False, label='Simulation 1')
sns.distplot(yields[yields['simulation'] == 1]['Yield'], kde=True, hist=False, label='Simulation 2')
plt.legend()

plt.show()