Now that we have the results, we can first find the mean rate at each time step for each simulation.

#Find the mean rates over time
mu = result.groupby(['timestep', 'simulation'])['r'].mean().unstack()
print(mu)

Let’s also convert the parameters to a dataframe.

#Turn the parameters into a dataframe
params = pd.DataFrame(params, columns = ['sigma', 'b', 't', 'a'])

The following code will graph all the results in a way that we can see two things. One, we see how the a parameter causes the convergence faster and second, we see that the simulation is not quite as stable with higher levels of sigma.

Another interesting metric is the standard deviation of the cross-section (meaning at each time step) for the rates. We will end up seeing how it starts at 0 (since all rates are the same to begin with), rises for a time and then levels off to a constant level at equilibrium.

#Find the standard deviation of rates between each simulation
std_cross_section = result.groupby(['timestep', 'simulation'])['r'].std().unstack()
print(std_cross_section)

The graphs created by the code below will allow us to compare the effect of a on the cross-sectional volatility. What you will see is how even with the same volatility term, the cross-sectional volatility is negatively related to the a term because it shrinks rates towards the mean more quickly.

fig, ax = plt.subplots(nrows=len(params['sigma'].unique()),
                      ncols=len(params['b'].unique()),
                      figsize=(15,15))

for i1, sigma in enumerate(params['sigma'].unique()):
    for i2, b in enumerate(params['b'].unique()):
        ax[i1,i2].set_title("Standard Deviation Path-b={} sigma={}".format(b,sigma))
        ts = std_cross_section[params[(params['sigma'] == sigma) & (params['b'] == b)].index]
        ts.columns = params[(params['sigma'] == sigma) & (params['b'] == b)]['a']
        ts.columns = "a="+ts.columns.astype(str)
        ts.plot(ax=ax[i1,i2])
plt.show()

We could also of course see the effect of volatilities in the graphs created by the code below as well.

fig, ax = plt.subplots(nrows=len(params['a'].unique()),
                      ncols=len(params['b'].unique()),
                      figsize=(15,15))

for i1, a in enumerate(params['a'].unique()):
    for i2, b in enumerate(params['b'].unique()):
        ax[i1,i2].set_title("Standard Deviation Path-b={} a={}".format(b,a))
        ts = std_cross_section[params[(params['a'] == a) & (params['b'] == b)].index]
        ts.columns = params[(params['a'] == a) & (params['b'] == b)]['sigma']
        ts.columns = "sigma="+ts.columns.astype(str)
        ts.plot(ax=ax[i1,i2])
plt.show()