Source code for economic.BasicHullWhite.HullWhite

"""The main Space in the :mod:`~economic.BasicHullWhite` model.

Mathematical notations are defined consistent with those in Brigo and Mercurio (2001, 2nd Ed. 2006)

.. seealso::

   * Damiano Brigo, Fabio Mercurio (2001, 2nd Ed. 2006). Interest Rate Models — Theory and Practice with Smile, Inflation and Credit


Attributes:

    scen_size: Number of scenarios. 1000 by default.

    np: numpy module.

    step_size: Number of time steps. 360 by default.

    time_len: Simulation length in years. 30 by default.

    a: Parameter :math:`a` in the Hull-White stochastic differential equation. 0.1 by default.

    sigma: Parameter :math:`\sigma` in the Hull-White stochastic differential equation. 0.1 by default.

    seed1: Seed number to generate random numbers. 1234 by default. See :meth:`std_norm_rand`.

    seed2: Seed for the second random numbers used for :meth:`accum_short_rate2`.
"""

from modelx.serialize.jsonvalues import *

_formula = None

_bases = []

_allow_none = None

_spaces = []

# ---------------------------------------------------------------------------
# Cells

[docs]def A_t_T(i, j): r""":math:`A(t_i, t_j)` in :math:`P(t_i, t_j)` See :meth:`P_t_T`. """ t = t_(i) P_t = mkt_zcb(i) P_T = mkt_zcb(j) f_t = mkt_fwd(i) B = B_t_T(i, j) return P_T / P_t * np.exp( B * f_t - sigma**2 / (4*a) * (1 - np.exp(-2*a*t)) * B**2 )
[docs]def B_t_T(i, j): r""":math:`B(t_i, t_j)` in :math:`P(t_i, t_j)` See :meth:`P_t_T`. """ t, T = t_(i), t_(j) return (1 / a) * (1 - np.exp(-a * (T-t)))
[docs]def E_rt(): """The expected values of :math:`r(t_i)` at time 0 for all :math:`t_i`. Returns, in a numpy array, the expected values of :math:`r(t_i)` for all :math:`t_i`. Calculated as :math:`E\{r(t_i) | \mathcal{F}_{0}\}`. .. seealso:: * :meth:`E_rt_s` """ return np.array([E_rt_s(0, i)[0] for i in range(step_size + 1)])
[docs]def E_rt_s(i, j): r"""Conditional expected values of :math:`r(t_j)` Returns, in a numpy array, :math:`E\{r(t_j) | \mathcal{F}_{i}\}`, the expected values of :math:`r(t_j)` conditional on :math:`\mathcal{F}_{i}` for all scenarios. :math:`E\{r(t) | \mathcal{F}_{s}\}` is calculated as: .. math:: r(s)e^{-a(t-s)} + \alpha(t) - \alpha(s)e^{-a(t-s)} where :math:`\alpha(t)` is calculated by :meth:`alpha`. .. seealso: * :meth:`short_rate` * :meth:`alpha` """ s, t = t_(i), t_(j) r_s = short_rate(i) return r_s * np.exp(-a * (t-s)) + alpha(j) - alpha(i) * np.exp(-a * (t-s))
[docs]def P_t_T(i, j): r"""The price at :math:`t_i` of a zero-coupon bond paying off 1 at :math:`t_j` This formula corresponds to :math:`P(t, T)` in Brigo and Mercurio, which is defined as: .. math:: P(t, T)=A(t, T)e^{-B(t, T)r(t)} where :math:`t_i` and :math:`t_j` substitute for :math:`t` and :math:`T`. .. seealso:: * :meth:`A_t_T` for :math:`A(t, T)` * :meth:`B_t_T` for :math:`B(t, T)` * :meth:`short_rate` for :math:`r(t)` * Brigo and Mercurio (2001, 2nd Ed. 2006). Interest Rate Models — Theory and Practice with Smile, Inflation and Credit """ return A_t_T(i, j) * np.exp(-B_t_T(i, j) * short_rate(i))
[docs]def V_t_T(i, j): r"""The variance of :math:`\int_{t_j}^{t_i} r(u)du|\mathcal{F}_{t_i}` This formula corresponds to :math:`V(t, T)` in Brigo and Mercurio, which is defined as: .. math:: V(t, T)=\frac{\sigma^2}{a^2}\left[T-t+\frac{2}{a}e^{-a(T-t)}-\frac{1}{2a}e^{-2a(T-t)}-\frac{3}{2a}\right] where :math:`t_i` and :math:`t_j` substitute for :math:`t` and :math:`T`. .. seealso:: * :attr:`sigma` for :math:`\sigma` * :attr:`a` for :math:`a` * Brigo and Mercurio (2001, 2nd Ed. 2006). Interest Rate Models — Theory and Practice with Smile, Inflation and Credit """ dt = t_(j) - t_(i) return sigma**2 / a**2 * (dt + (2/a)*np.exp(-a*dt) - (1/(2*a))*np.exp(-2*a*dt) - (3/(2*a)))
[docs]def Var_rt(): r"""The variance of :math:`r(t_i)` at time 0 for all :math:`t_i`. Returns, in a numpy array, the variance of :math:`r(t_i)` for all :math:`t_i`. Calculated as :math:`Var\{r(t_i) | \mathcal{F}_{0}\}`. .. seealso:: * :meth:`Var_rt_s` """ return np.array([Var_rt_s(0, i) for i in range(step_size + 1)])
[docs]def Var_rt_s(i, j): r"""The variance of :math:`r(t_j)` conditional on :math:`\mathcal{F}_{t_i}` :math:`Var\{r(t_{j}) | \mathcal{F}_{t_i}\}`, the variance of :math:`r(t_j)` conditional on :math:`\mathcal{F}_{t_i}`, calculated as: .. math:: Var\{ r(t) | \mathcal{F}_s \} = \frac{\sigma^2}{2a} (1 - e^{-2a(t-s)}) .. seealso:: * :attr:`a` * :attr:`sigma` """ s, t = t_(i), t_(j) return sigma**2 / (2*a) * (1 - np.exp(-2 * a * (t-s)))
[docs]def accum_short_rate(i): r"""Accumulated short rates. a descrete approximation to the integral :math:`\int_0^{t_i}r(t)dt`, calculated as :math:`\sum_{j=1}^{i}r(t_{j-1})(t_j-t_{j-1})` .. seealso:: * :meth:`disc_factor` """ if i == 0: return np.full(scen_size, 0.0) else: dt = t_(i) - t_(i-1) return accum_short_rate(i-1) + short_rate(i-1) * dt
[docs]def accum_short_rate2(i): r"""Alternative implementation of accumulated short rates. An alternative approach to simulate :math:`Y(t_i)=\int_0^{t_i}r(t)dt` by using the fact that :math:`Y(t_i)` follows a normal distribution, and by simulating the joint distribution of :math:`(r(t_i), Y(t_i))`, as suggested in Glasserman (2003). .. seealso:: * :meth:`accum_short_rate` * :attr:`seed2` * Paul Glasserman (2003). Monte Carlo Methods in Financial Engineering """ if i == 0: return np.full(scen_size, 0.0) else: t, T = t_(i-1), t_(i) dt = T - t cov = sigma**2/(2*a**2)*(1 + np.exp(-2*a*dt) -2 * np.exp(-a*dt)) z1 = std_norm_rand(seed1)[:, i-1] z2 = std_norm_rand(seed2)[:, i-1] rho = cov / (Var_rt_s(i-1, i)**0.5 * V_t_T(i-1, i)**0.5) mean = B_t_T(i-1, i) * (short_rate(i-1) - alpha(i-1)) + np.log(mkt_zcb(i-1)/mkt_zcb(i)) + 0.5*(V_t_T(0, i)-V_t_T(0, i-1)) return accum_short_rate2(i-1) + mean + V_t_T(i-1, i)**0.5 * (rho*z1 + (1-rho**2)**0.5*z2)
[docs]def alpha(i): r""":math:`\alpha(t_i)` Returns, in a numpy array, :math:`\alpha(t_i)` for all scenarios. :math:`\alpha` appears in the expression of :math:`E\{r(t) | \mathcal{F}_{s}\}` and is defined as: .. math:: \alpha(t) = f^M(0, t) + \frac{\sigma^2} {2a^2}(1-e^{-at})^2 .. seealso:: * :meth:`E_rt_s` """ t = t_(i) return mkt_fwd(i) + 0.5 * sigma**2 / a**2 * (1 - np.exp(-a*t))**2
[docs]def disc_factor(i): """Discount factors Returns, in a numpy array, the discount factors for cashflows at :math:`t_i` for all scenarios. Defined as:: np.exp(-accum_short_rate(i)) .. seealso:: * accum_short_rate """ return np.exp(-accum_short_rate(i))
[docs]def disc_factor_paths(): """Discount factor scenarios. Returns, as a 2D numpy array, the simulated discount factors for all scenarios. .. seealso:: * :meth:`disc_factor` """ return np.array([disc_factor(i) for i in range(step_size + 1)]).transpose()
[docs]def mean_disc_factor(): """Discount factor means Returns, as a numpy array, the mean of discount factors of all scenarios for each :math:`t_i`. .. seealso:: * :meth:`disc_factor` """ return np.array([np.mean(disc_factor(i)) for i in range(step_size + 1)])
[docs]def mean_short_rate(): """The means of generated short rates Returns, as a numpy array, the means of short rates of all scenarios for all :math:`t_i`. This should converge to the theoretical variances calculated by :meth:`E_rt`. .. seealso:: * :meth:`short_rate` * :meth:`E_rt` """ return np.array([np.mean(short_rate(i)) for i in range(step_size + 1)])
[docs]def mkt_fwd(i): """The initial instantaneous forward rate for :attr:`t_(i)<t_>`. By default, returns 0.05 for all ``i``. """ return 0.05
[docs]def mkt_zcb(i): """The initial price of zero coupon bond The initial price of the unit zero coupon bond maturing at :attr:`t_(i)<t_>`. If ``i=0`` returns 1. Otherwise, defined as:: mkt_zcb(i-1) * np.exp(-mkt_fwd(i-1)*dt) where ``dt = t_(i) - t_(i-1)``. .. seealso:: * :attr:`t_` * :attr:`mkt_fwd` """ if i == 0: return 1.0 else: dt = t_(i) - t_(i-1) return mkt_zcb(i-1) * np.exp(-mkt_fwd(i-1)*dt)
[docs]def short_rate(i): r"""Stochastic short rates at :attr:`t_(i)<t_>` Returns, in a numpy array, simulated stochastic short rates at :attr:`t_(i)<t_>` for all scenarios. For ``i=0``, defined as :meth:`mkt_fwd(0)<mkt_fwd>`. For ``i>0``, defined as :math:`r(t_i) = E\{r(t_i) | \mathcal{F}_{i-1}\} + \sqrt{Var\{ r(t_i) | \mathcal{F}_{i-1} \}} * Z`, where :math:`E\{r(t_i) | \mathcal{F}_{i-1}\}`, the expected value of :math:`r(t_i)` conditional on :math:`\mathcal{F}_{i-1}` is calculated by :meth:`E_rt_s(i-1, i)<E_rt_s>`, :math:`Var\{ r(t_i) | \mathcal{F}_{i-1} \}` the variance of :math:`r(t_i)` conditional on :math:`\mathcal{F}_{i-1}` is calculated by :meth:`Var_rt_s(i-1, i)<Var_rt_s>`, and :math:`Z`, a random number drawn from :math:`\mathcal{N}(0, 1)` a standard normal distribution calculated by :meth:`std_norm_rand`. .. seealso:: * :attr:`scen_size` * :meth:`mkt_fwd` * :meth:`E_rt_s` * :meth:`Var_rt_s` * :meth:`std_norm_rand` * :attr:`seed1` """ if i == 0: return np.full(scen_size, mkt_fwd(0)) else: return E_rt_s(i-1, i) + Var_rt_s(i-1, i)**0.5 * std_norm_rand(seed1)[:, i-1]
[docs]def short_rate_paths(): """Short rate paths. Returns, as a 2D numpy array, the simulated short rate paths for all scenarios. .. seealso:: * :meth:`short_rate` """ return np.array([short_rate(i) for i in range(step_size + 1)]).transpose()
[docs]def std_norm_rand(seed=1234): """Random numbers from the standard normal distribution. Returns a numpy array shaped :attr:`scen_size` x :attr:`step_size`. The elements are random numbers drawn from the standard normal distribution. """ size = (scen_size, step_size) if hasattr(np.random, 'default_rng'): gen = np.random.default_rng(seed) return gen.standard_normal(size) else: np.random.seed(seed) return np.random.standard_normal(size)
t_ = lambda i: i * time_len / step_size """time index :math:`t_i`"""
[docs]def var_short_rate(): """Variance of generated short rates Returns, as a vector in a numpy array, the variances of the generated short rates for all :math:`t_i`. This should converge to the theoretical variances calculated by :meth:`Var_rt`. .. seealso:: * :meth:`short_rate` * :meth:`Var_rt` """ return np.array([np.var(short_rate(i)) for i in range(step_size + 1)])
# --------------------------------------------------------------------------- # References np = ("Module", "numpy") step_size = 360 time_len = 30 a = 0.1 sigma = 0.1 seed1 = 1234 seed2 = 5678 scen_size = 1000