kou's jump model

pfhedge/instruments/__init__.py

@@ -15,6 +15,7 @@
 from .primary.brownian import BrownianStock  # NOQA
 from .primary.cir import CIRRate  # NOQA
 from .primary.heston import HestonStock  # NOQA
+from .primary.kou_jump import KouJumpStock  # noqa: F401
 from .primary.local_volatility import LocalVolatilityStock  # NOQA
 from .primary.merton_jump import MertonJumpStock  # NOQA
 from .primary.rough_bergomi import RoughBergomiStock  # NOQA

pfhedge/instruments/primary/kou_jump.py

@@ -0,0 +1,191 @@
+from math import ceil
+from typing import Optional, Tuple, cast
+
+import torch
+from torch import Tensor
+
+from pfhedge._utils.doc import _set_attr_and_docstring, _set_docstring
+from pfhedge._utils.str import _format_float
+from pfhedge._utils.typing import TensorOrScalar
+from pfhedge.stochastic import generate_kou_jump
+
+from .base import BasePrimary
+
+
+class KouJumpStock(BasePrimary):
+    r"""A stock of which spot prices follow the Kou's jump diffusion.
+
+    .. seealso::
+        - :func:`pfhedge.stochastic.generate_kou_jump`:
+          The stochastic process.
+
+    Args:
+        sigma (float, default=0.2): The parameter :math:`\sigma`,
+            which stands for the volatility of the spot price.
+        mu (float, default=0.0): The parameter :math:`\mu`,
+            which stands for the drift of the spot price.
+        jump_per_year (float, optional): Jump poisson process annual
+            lambda: Average number of annual jumps. Defaults to 1.0.
+        jump_eta_up (float, optional): 1/Mu for the up jumps:
+            Instaneous value. Defaults to 1/0.02.
+            This has to be larger than 1.
+        jump_eta_down (float, optional): 1/Mu for the down jumps:
+            Instaneous value. Defaults to 1/0.05.
+            This has to be larger than 0.
+        jump_up_prob (float, optional): Given a jump occurs,
+            this is conditional prob for up jump.
+            Down jump occurs with prob 1-jump_up_prob.
+            Has to be in [0,1].
+        cost (float, default=0.0): The transaction cost rate.
+        dt (float, default=1/250): The intervals of the time steps.
+        dtype (torch.device, optional): Desired device of returned tensor.
+            Default: If None, uses a global default
+            (see :func:`torch.set_default_tensor_type()`).
+        device (torch.device, optional): Desired device of returned tensor.
+            Default: if None, uses the current device for the default tensor type
+            (see :func:`torch.set_default_tensor_type()`).
+            ``device`` will be the CPU for CPU tensor types and
+            the current CUDA device for CUDA tensor types.
+
+    Buffers:
+        - spot (:class:`torch.Tensor`): The spot prices of the instrument.
+          This attribute is set by a method :meth:`simulate()`.
+          The shape is :math:`(N, T)` where
+          :math:`N` is the number of simulated paths and
+          :math:`T` is the number of time steps.
+
+    Examples:
+        >>> from pfhedge.instruments import KouJumpStock
+        >>>
+        >>> _ = torch.manual_seed(42)
+        >>> stock = KouJumpStock()
+        >>> stock.simulate(n_paths=2, time_horizon=5 / 250)
+        >>> stock.spot
+        tensor([[1.0000, 1.0101, 1.0137, 1.0144, 1.0211, 1.0180],
+                [1.0000, 1.0067, 1.0065, 1.0158, 1.0175, 1.0287]])
+
+        Using custom ``dtype`` and ``device``.
+
+        >>> stock = KouJumpStock()
+        >>> stock.to(dtype=torch.float64, device="cuda:0")
+        KouJumpStock(sigma=0.2000, dt=0.0040, jump_per_year=1., j_mu=0.,
+        j_sigma=0.2000, dtype=torch.float64, device='cuda:0')
+    """
+
+    def __init__(
+        self,
+        sigma: float = 0.2,
+        mu: float = 0.0,
+        jump_per_year: float = 68.0,
+        jump_eta_up: float = 1 / 0.02,
+        jump_eta_down: float = 1 / 0.05,
+        jump_up_prob: float = 0.5,
+        cost: float = 0.0,
+        dt: float = 1 / 250,
+        dtype: Optional[torch.dtype] = None,
+        device: Optional[torch.device] = None,
+    ) -> None:
+        super().__init__()
+
+        self.sigma = sigma
+        self.mu = mu
+        self.jump_per_year = jump_per_year
+        self.jump_eta_up = jump_eta_up
+        self.jump_eta_down = jump_eta_down
+        self.jump_up_prob = jump_up_prob
+        self.cost = cost
+        self.dt = dt
+
+        self.to(dtype=dtype, device=device)
+
+    @property
+    def default_init_state(self) -> Tuple[float, ...]:
+        return (1.0,)
+
+    @property
+    def volatility(self) -> Tensor:
+        """Returns the volatility of self.
+
+        It is a tensor filled with ``self.sigma``.
+        """
+        return torch.full_like(self.get_buffer("spot"), self.sigma)
+
+    @property
+    def variance(self) -> Tensor:
+        """Returns the volatility of self.
+
+        It is a tensor filled with the square of ``self.sigma``.
+        """
+        return torch.full_like(self.get_buffer("spot"), self.sigma ** 2)
+
+    def simulate(
+        self,
+        n_paths: int = 1,
+        time_horizon: float = 20 / 250,
+        init_state: Optional[Tuple[TensorOrScalar]] = None,
+    ) -> None:
+        """Simulate the spot price and add it as a buffer named ``spot``.
+
+        The shape of the spot is :math:`(N, T)`, where :math:`N` is the number of
+        simulated paths and :math:`T` is the number of time steps.
+        The number of time steps is determinded from ``dt`` and ``time_horizon``.
+
+        Args:
+            n_paths (int, default=1): The number of paths to simulate.
+            time_horizon (float, default=20/250): The period of time to simulate
+                the price.
+            init_state (tuple[torch.Tensor | float], optional): The initial state of
+                the instrument.
+                This is specified by a tuple :math:`(S(0),)` where
+                :math:`S(0)` is the initial value of the spot price.
+                If ``None`` (default), it uses the default value
+                (See :attr:`default_init_state`).
+                It also accepts a :class:`float` or a :class:`torch.Tensor`.
+
+        Examples:
+            >>> from pfhedge.instruments import KouJumpStock
+            >>>
+            >>> _ = torch.manual_seed(42)
+            >>> stock = KouJumpStock()
+            >>> stock.simulate(n_paths=2, time_horizon=5 / 250, init_state=(2.0,))
+            >>> stock.spot
+            tensor([[2.0000, 2.0032, 2.0091, 2.0149, 1.9865, 1.9817],
+                    [2.0000, 1.9839, 1.9954, 2.0022, 2.0157, 2.0364]])
+        """
+        if init_state is None:
+            init_state = cast(Tuple[float], self.default_init_state)
+
+        spot = generate_kou_jump(
+            n_paths=n_paths,
+            n_steps=ceil(time_horizon / self.dt + 1),
+            init_state=init_state,
+            sigma=self.sigma,
+            mu=self.mu,
+            jump_per_year=self.jump_per_year,
+            jump_eta_up=self.jump_eta_up,
+            jump_eta_down=self.jump_eta_down,
+            jump_up_prob=self.jump_up_prob,
+            dt=self.dt,
+            dtype=self.dtype,
+            device=self.device,
+        )
+
+        self.register_buffer("spot", spot)
+
+    def extra_repr(self) -> str:
+        params = ["sigma=" + _format_float(self.sigma)]
+        if self.mu != 0.0:
+            params.append("mu=" + _format_float(self.mu))
+        if self.cost != 0.0:
+            params.append("cost=" + _format_float(self.cost))
+        params.append("dt=" + _format_float(self.dt))
+        params.append("jump_per_year=" + _format_float(self.jump_per_year))
+        params.append("jump_eta_up=" + _format_float(self.jump_eta_up))
+        params.append("jump_eta_down=" + _format_float(self.jump_eta_down))
+        params.append("jump_up_prob=" + _format_float(self.jump_up_prob))
+        return ", ".join(params)
+
+
+# Assign docstrings so they appear in Sphinx documentation
+_set_docstring(KouJumpStock, "default_init_state", BasePrimary.default_init_state)
+_set_attr_and_docstring(KouJumpStock, "to", BasePrimary.to)

pfhedge/stochastic/__init__.py

@@ -2,6 +2,7 @@
 from .brownian import generate_geometric_brownian  # NOQA
 from .cir import generate_cir  # NOQA
 from .heston import generate_heston  # NOQA
+from .kou_jump import generate_kou_jump  # noqa: F401
 from .local_volatility import generate_local_volatility_process  # NOQA
 from .merton_jump import generate_merton_jump  # NOQA
 from .random import randn_antithetic  # NOQA

pfhedge/stochastic/kou_jump.py

@@ -0,0 +1,175 @@
+import math
+from typing import Callable, Optional, Tuple, Union
+
+import torch
+from torch import Tensor
+
+from pfhedge._utils.typing import TensorOrScalar
+
+from ._utils import cast_state
+
+
+def generate_kou_jump(
+    n_paths: int,
+    n_steps: int,
+    init_state: Union[Tuple[TensorOrScalar, ...], TensorOrScalar] = (1.0,),
+    sigma: float = 0.2,
+    mu: float = 0.0,
+    jump_per_year: float = 68.0,
+    jump_eta_up: float = 1 / 0.02,
+    jump_eta_down: float = 1 / 0.05,
+    jump_up_prob: float = 0.5,
+    dt: float = 1 / 250,
+    dtype: Optional[torch.dtype] = None,
+    device: Optional[torch.device] = None,
+    engine: Callable[..., Tensor] = torch.randn,
+) -> Tensor:
+    r"""Kou's Jump Diffusion Model for stock prices.
+       Assumes number of jumps to be poisson distribution
+       with ASYMMETRIC jump for up and down movement; the
+       lof of these jump follows exponential distribution
+       with mean 1/jump_eta_up and 1/jump_eta_down resp.
+
+       See Glasserman, Paul. Monte Carlo Methods in Financial
+       Engineering. New York: Springer-Verlag, 2004.for details.
+       Copy available at "https://www.bauer.uh.edu/spirrong/
+       Monte_Carlo_Methods_In_Financial_Enginee.pdf"
+
+       Combined with the original paper by Kou:
+       A Jump-Diffusion Model for Option Pricing
+       https://www.columbia.edu/~sk75/MagSci02.pdf
+
+    Args:
+        n_paths (int): The number of simulated paths.
+        n_steps (int): The number of time steps.
+        init_state (tuple[torch.Tensor | float], default=(0.0,)): The initial state of
+            the time series.
+            This is specified by a tuple :math:`(S(0),)`.
+            It also accepts a :class:`torch.Tensor` or a :class:`float`.
+            The shape of torch.Tensor must be (1,) or (n_paths,).
+        sigma (float, default=0.2): The parameter :math:`\sigma`,
+            which stands for the volatility of the time series.
+        mu (float, default=0.0): The parameter :math:`\mu`,
+            which stands for the drift of the time series.
+        jump_per_year (float, optional): Jump poisson process annual
+            lambda: Average number of annual jumps. Defaults to 1.0.
+        jump_eta_up (float, optional): 1/Mu for the up jumps:
+            Instaneous value. Defaults to 1/0.02.
+            This has to be larger than 1.
+        jump_eta_down (float, optional): 1/Mu for the down jumps:
+            Instaneous value. Defaults to 1/0.05.
+            This has to be larger than 0.
+        jump_up_prob (float, optional): Given a jump occurs,
+            this is conditional prob for up jump.
+            Down jump occurs with prob 1-jump_up_prob.
+            Has to be in [0,1].
+        dt (float, default=1/250): The intervals of the time steps.
+        dtype (torch.dtype, optional): The desired data type of returned tensor.
+            Default: If ``None``, uses a global default
+            (see :func:`torch.set_default_tensor_type()`).
+        device (torch.device, optional): The desired device of returned tensor.
+            Default: If ``None``, uses the current device for the default tensor type
+            (see :func:`torch.set_default_tensor_type()`).
+            ``device`` will be the CPU for CPU tensor types and the current CUDA device
+            for CUDA tensor types.
+        engine (callable, default=torch.randn): The desired generator of random numbers
+            from a standard normal distribution.
+            A function call ``engine(size, dtype=None, device=None)``
+            should return a tensor filled with random numbers
+            from a standard normal distribution.
+            Only to be used for the normal component,
+            jupms uses poisson distribution
+
+    Shape:
+        - Output: :math:`(N, T)` where
+          :math:`N` is the number of paths and
+          :math:`T` is the number of time steps.
+
+    Returns:
+        torch.Tensor
+
+    Examples:
+        >>> from pfhedge.stochastic import generate_kou_jump
+        >>>
+        >>> _ = torch.manual_seed(42)
+        >>> generate_kou_jump(2, 5)
+        tensor([[1.0000, 1.0053, 1.0119, 0.9272, 0.9174],
+                [1.0000, 1.0321, 1.0275, 1.0373, 1.0446]])
+    """
+    assert jump_eta_up > 1.0, "jump_eta_up must be larger than 1.0"
+    assert jump_eta_down > 0.0, "jump_eta_down must be larger than 0.0"
+    assert 0 <= jump_up_prob <= 1.0, "jump prob must be in [0,1]"
+
+    init_state = cast_state(init_state, dtype=dtype, device=device)
+
+    init_value = init_state[0]
+
+    t = dt * torch.arange(n_steps, device=device, dtype=dtype)[None, :]
+    returns = (
+        engine(*(n_paths, n_steps), dtype=dtype, device=device) * math.sqrt(dt) * sigma
+    )
+
+    returns[:, 0] = 0.0
+
+    # Generate jump components
+    n_jumps = torch.poisson(torch.full((n_paths, n_steps - 1), jump_per_year * dt)).to(
+        returns
+    )
+    # if n_steps is greater than 1
+    if (n_steps-1 ) > 0:
+        # max jumps used to aggregte jump in between dt time
+        max_jumps = int(n_jumps.max())
+        size_paths = (n_paths, n_steps - 1, max_jumps)
+
+        # up exp generator
+        up_exp_dist = torch.distributions.Exponential(jump_eta_up)
+
+        # down exp generator
+        down_exp_dist = torch.distributions.Exponential(jump_eta_down)
+
+        log_jump = torch.where(
+            torch.rand(size_paths) < jump_up_prob,
+            up_exp_dist.sample(size_paths),
+            -down_exp_dist.sample(size_paths),
+        ).to(returns)
+
+        # for no jump condition
+        log_jump = torch.cat(
+            (torch.zeros(n_paths, n_steps - 1, 1).to(log_jump), log_jump), dim=-1
+        )
+
+        exp_jump_ind = torch.exp(log_jump)
+
+        # filter out jump movements that did not occur in dt time
+        indices_expanded = n_jumps[..., None]
+        k_range = torch.arange(max_jumps + 1).to(returns)
+        mask = k_range > indices_expanded
+        # exp(0) as to no jump after n_jump
+        exp_jump_ind[mask] = 1.0
+
+        # aggregate jumps in time dt--> multiplication of exponent
+        exp_jump = torch.prod(exp_jump_ind, dim=-1)
+
+        # no jump at time 0--> exp(0.0)=1.0
+        exp_jump = torch.cat((torch.ones(n_paths, 1).to(exp_jump), exp_jump), dim=1)
+
+    else:
+        exp_jump = torch.ones(n_paths, 1).to(returns)
+
+    # aggregate jumps upto time t
+    exp_jump_agg = torch.cumprod(exp_jump, dim=-1)
+
+    # jump correction for drift: see the paper
+    m = (
+        (1 - jump_up_prob) * (jump_eta_down / (jump_eta_down + 1))
+        + (jump_up_prob) * (jump_eta_up / (jump_eta_up - 1))
+        - 1
+    )
+
+    prices = (
+        torch.exp((mu - jump_per_year * m) * t + returns.cumsum(1) - (sigma ** 2) * t / 2)
+        * init_value.view(-1, 1)
+        * exp_jump_agg
+    )
+
+    return prices