CIR Model - Square-root Diffusion

Reference: Python for Finance by Yves Hilpisch, O'Reilly Media, Chapter 10, Stochastics

Multidimensional Stochastic Processes

import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt
import pandas as pd

Do also:

• Vasicek model https://en.wikipedia.org/wiki/Vasicek_model
• Hull–White model

Used to model short interest rates, bond prices, or volatility. It fixes the possibility of negative interest rates allowed by the Vasicek model.

Properties:

• Mean reversion; $\theta$ is the long-term mean
• Level dependent volatility
• $x_t$ are chi-squared distributed. Asymptotic distribution is gamma.
• For given positive $r_0$ the process will never touch zero, if $2k\theta \geq \sigma^2$

Square-root Diffusion SDE:

$dx_t = k(\theta - x_t)dt + \sigma \sqrt{x_t}\sqrt{dt} z_t $ where all $x_t$ must be positive values, or zero if else.

Positivity cannot be guaranteed in the discretization, so $x_t$ is replaced by $x^+_t = max(x_t,0)$.

#Simulate I CIR paths: 
T = 50
I = 100
dt = 1/50

#Real data: ?
sigma = 0.1
x0 = 0.05
kappa = 3.0
theta= 0.02

# z=npr.standard_normal((T, I))

def discretize():
    x = np.zeros((T+1, I))
    x[0, :] = x0

    for t in range(1, T+1):
        clamped = np.maximum(x[t-1,:], 0)
        z=npr.standard_normal(I)
        x[t, :]=x[t-1,:] +kappa*(theta-clamped)*dt+sigma*np.sqrt(clamped*dt)*z

    return x                                                                                                           
                                                                                                               
x = discretize()
xplus = np.maximum(x,0)
  

df = pd.DataFrame(xplus)
df.describe()
# Google total re

	0	1	2	3	4	5	6	7	8	9	...	90	91	92	93	94	95	96	97	98	99
count	51.000000	51.000000	51.000000	51.000000	51.000000	51.000000	51.000000	51.000000	51.000000	51.000000	...	51.000000	51.000000	51.000000	51.000000	51.000000	51.000000	51.000000	51.000000	51.000000	51.000000
mean	0.031736	0.022249	0.029834	0.022945	0.028403	0.033074	0.026718	0.027869	0.029797	0.028638	...	0.024937	0.033960	0.027133	0.026626	0.022302	0.036980	0.036411	0.032499	0.029943	0.023069
std	0.011440	0.008807	0.010369	0.006771	0.009802	0.008398	0.010294	0.010324	0.008604	0.010872	...	0.012102	0.006467	0.008781	0.010018	0.013579	0.008649	0.008207	0.006215	0.008651	0.009157
min	0.017584	0.013395	0.016429	0.015194	0.011142	0.020937	0.014882	0.014746	0.013700	0.014933	...	0.009291	0.023790	0.016950	0.012480	0.010877	0.024166	0.019854	0.023093	0.018221	0.009199
25%	0.021142	0.016715	0.023227	0.018691	0.019324	0.026558	0.018805	0.019103	0.022258	0.021062	...	0.014750	0.029718	0.021050	0.016018	0.013246	0.027674	0.031756	0.027814	0.024403	0.016866
50%	0.027772	0.019524	0.026500	0.020264	0.030133	0.030689	0.023885	0.023491	0.030308	0.025061	...	0.023301	0.032174	0.023196	0.028553	0.015341	0.038455	0.035399	0.031464	0.028018	0.021315
75%	0.043208	0.023234	0.031572	0.026314	0.036255	0.040322	0.034751	0.033694	0.034665	0.032991	...	0.032547	0.037501	0.030657	0.032317	0.029070	0.044235	0.042234	0.036161	0.031172	0.025624
max	0.053148	0.050000	0.055284	0.050000	0.050000	0.050489	0.050129	0.050000	0.050000	0.054442	...	0.051944	0.050000	0.051101	0.050000	0.050738	0.050256	0.052650	0.050000	0.051338	0.050000

8 rows × 100 columns

ST = xplus[-1]

fig, ax1 = plt.subplots()#
ax1.hist(ST, bins=50); ax1.set_title('Terminal Values - Square-root Diffusion')

Text(0.5, 1.0, 'Terminal Values - Square-root Diffusion')

df = pd.DataFrame(ST)
df.describe()
# Google total return over this period is 569.07%.  

	0
count	100.000000
mean	0.021164
std	0.005762
min	0.010452
25%	0.016419
50%	0.021217
75%	0.025562
max	0.036022

plt.plot(xplus[:, :20], lw = 1.5)
plt.xlabel('# periods'); plt.ylabel('Index Level'); plt.title("Square-Root Diffusion Process. The process converges to the long-term mean Theta")
plt.grid(True)

np.cov(x[:,0], x[:,1])

array([[1.30880184e-04, 7.97155787e-05],
       [7.97155787e-05, 7.75650550e-05]])

#xx = pd.DataFrame(xplus)
## save to xlsx file
#filepath = './xplus.xlsx'
#xx.to_excel(filepath, index=False)

Application to Bond Values: TO DO

import pandas_datareader as pr
from pandas_datareader import data

#start_date='6/30/2010'; end_date='07/31/2020'
#google = data.DataReader('bond_name', 'yahoo', start_date, end_date)