Iterative weak approximation and hard bounds for switching diffusion

Owing to the increasing demands on regime-switching diffusion from emerging applications in financial engineering and wireless communications, much attention has been given to switching diffusion processes. A salient feature of such systems is that they include both continuous dynamic and discrete events.

States of the environment process was referred to in [20] as economic circumstances or political regime-switching. It is therefore theoretically appealing to modify the classical diffusion process so that both expected return and the variance reflect the changes in the external environmental factors. The modelling framework that is advocated in this paper achieves this. The rationale is that in the different modes or regimes, the volatility and return rates can be very different. Another example is from a wireless communication networks. Consider the performance analysis of an adaptive, linear, multi-user detector in a cellular, direct-sequence, code-division, multiple-access wireless network with changing user activity due to an admission or access controller at the base station. Under certain conditions, an optimization problem for the aforementioned application leads to certain issues, which in turn lead to switching diffusion limit [26].

One early attempt at such hybrid models for financial applications can be traced back to the 1960s in which both the appreciation rate and the volatility of a stock depend on a continuous-time Markov chain. The introduction of such models makes it possible to describe stochastic volatility in a relatively simple manner. To illustrate, in the simplest case, a stock market may be considered to have two “modes” or “regimes”, up and down, resulting from the state of the underlying economy, the general mood of investors in the market, and so on. The primary motivation for the generalization is to enhance the flexibility of the model parameter settings for the classical diffusion process. The examples usually given are weather conditions and epidemic outbreaks, even though seasonality would play a role and can probably not be modelled by a Markov regime-switching model.

A variety of numerical methods extend those that already exist in the single-regime setting, but the extensions are rarely straightforward due to the inter-dependence nature of regime switching systems. Lattice methods in regime-switching settings have been attempted; first in [5] in a two-regime setting and later in the form of the adaptive mesh model [11], the $2k$-branch lattice for a $k$-regime model [1], and the trinomial method [28]. Mathematical programming is employed for obtaining hard bounds on regime-switching diffusion in [4, 19]. Approaches based on penalty methods for pricing American options have also been extended in a non-trivial way to regime-switching frameworks. For instance, in [14, 17, 18], exponential time difference schemes are applied in combination with a penalty method to approximate a solution to the complex system of partial differential equations with coupled free boundary conditions that arise in the presence of a regime-switching diffusion. The Fourier space time-stepping method [15] overcomes some of the problems encountered by lattice methods and finite-difference schemes when it comes to option pricing in models in higher dimensions or models containing jumps. This approach is applicable to a wide class of path-dependent options and options on multiple assets in regime-switching diffusion containing jumps.

A Monte Carlo method is presented in [13] for pricing barrier option in regime-switching diffusion, where the author demonstrates the non-trivial modifications that must be made in order to extend standard simulation-based techniques from the single-regime setting to multiple regimes. Esscher transform methods are employed in [9] and [10] for pricing European and American options in regime-switching diffusion and jump-diffusion settings. The Canadization method was extended to regime-switching settings to price American options [7] and double barrier options [6]. We mention the radial basis collocation method [3], a mesh-free method for pricing American options under regime-switching jump diffusion models. A MATLAB toolbox specially designed for the estimation, simulation and forecasting of a general Markov regime-switching diffusion is developed in [23]. In all the aforementioned studies, a great deal of special attention must be paid to the switching nature.

In this paper, we establish a novel convergent iteration framework on a generalized Feynman-Kac formula under regime-switching, which are the solutions to a set of weakly coupled partial differential equations with initial and boundary value conditions. We construct our recursive approximation mechanism by restricting the maximum number of switching occurring before the terminal time in the probabilistic representation in a similar manner to [24, 25]. The significance of the present approach is that it succeeds to untangle and compensate a challenging set of weakly coupled partial differential equations to a set of independent partial differential equations. Such untangling of the inter-dependency is quite remarkable, as undoubtedly a variety of accurate and efficient numerical methods are available for the latter untangled partial differential equations.

The rest of the paper is organized as follows. We present the necessary background materials and problem formulation regarding regime-switching diffusion in Section 2. We prove in Sections 3.1 and 3.2 that such an approximation can be obtained as the solution to a set of independent partial differential equations with additional heat source (the non-homogeneous terms) and additional discounting rates or, equivalently, the approximation admits a probabilistic representation under a probability measure without regime-switching. Further, in Section 3.3, we derive hard bounding functions for the true solution under regime-switching based on the approximation at each iteration as well as their convergence to the unknown solution. After briefly discussing the proposed approach in the context the initial-boundary value problem in Section 4, we examine our theoretical findings through numerical results in Section 5. Finally, Section 6 concludes the present study and highlights future research directions.

We begin with the notation that will be used throughout the paper. We denote by $\mathbb{R}$ the set of real numbers, and by $\mathbb{N}:=\{1,2,\cdots\}$ the set of natural numbers excluding $0$, with $\mathbb{N}_{0}:=\mathbb{N}\cup\{0\}$. We reserve $p$, $n$ and $d$ for fixed positive integers. We denote by $\{R_{t}:t\geq 0\}$ a one-dimensional Markov process with right-continuous sample paths taking values in the finite-state space $\mathcal{M}:=\{1,\cdots,p\}$. We define an $n$-dimensional stochastic process $\{X_{t}:t\geq 0\}$ taking values in a suitable domain $D(\subseteq\mathbb{R}^{n})$ by the stochastic differential equation:

where $\{W_{t}:t\geq 0\}$ is a $d$-dimensional standard Brownian motion and the drift coefficient $b:[0,\infty)\times D\times\mathcal{M}\to\mathbb{R}^{n}$ and the diffusion coefficient $\sigma:[0,\infty)\times D\times\mathcal{M}\to\mathbb{R}^{n\times d}$ are deterministic functions satisfying the usual conditions [12, 21]: for each $i\in\mathcal{M}$, $\sigma(\cdot,\cdot,i)$ and $b(\cdot,\cdot,i)$ are at most of linear growth and Lipschitz continuous on every compact subset of $[0,\infty)\times D$. Here, the domain $D(\subseteq\mathbb{R}^{n})$ depends on the problem setting under consideration. For example, we may have $D=(0,\infty)$ when the underlying process $\{X_{t}:t\geq 0\}$ is a regime-switching geometric Brownian motion starting from a positive level, without an attainable boundary. For problems with attainable boundaries (Section 4), the domain $D$ is a bounded. We define the first hitting time of the underlying process on the boundary by

with $\eta_{t}=+\infty$, if the boundary is never attained. As usual, we denote by $\partial D$ the boundary of the domain $D$, and by $\overline{D}=D\cup\partial D$ the closure.

In this paper, the transition rate of the regime process $\{R_{t}:t\geq 0\}$ is Markov and can be state dependent, governed by the so-called generator matrix $Q:D\to\mathbb{R}^{p\times p}$, where we denote its $(i,j)$-entry by $q_{ij}({\bf x}):=(Q({\bf x}))_{i,j}$. To ensure that the stochastic differential equation (2.1) admits a unique solution, we impose the Assumption 2.1, the so-called $q$-property [2, 27], on top of the aforementioned usual conditions on the coefficients $b$ and $\sigma$.

We write $\mathbb{P}_{Q}$ for the probability measure under which the transition rate of the regime process $\{R_{t}:t\geq 0\}$ is given by the generator matrix $Q$. Throughout, we assume that the paired stochastic process $\{(X_{t},R_{t}):t\geq 0\}$ is defined on the probability space $(\Omega,\mathcal{F},\mathbb{P}_{Q})$ endowed with the right-continuous filtration $\mathcal{F}=(\mathcal{F}_{t})_{t\geq 0}$. For $i,j\in\mathcal{M}$, $(t,{\bf x})\in[0,\infty)\times D$ and an infinitesimal $\Delta$, we formally define

Hereafter, we denote by $\mathbb{E}_{Q}^{t,{\bf x},i}$ the expectation taken under the probability measure $\mathbb{P}_{Q}^{t,{\bf x},i}$. Next, for $t\geq 0$, we define the time of the first switching occur after time $t$ as:

Since all entries of the generator matrix $Q$ are bounded, the probability that the Markov process $\{R_{s}:s\geq 0\}$ switches more than once within an infinitesimal-sized interval $(t,t+\Delta]$ is as negligibly small as $o(\Delta)$, that is,

due to (2.3). For $m\in\mathbb{N}$, the time of $(m+1)$-st switching occur after time $t$ can be defined recursively as the first switching after the $m$-th switching after time $t$, that is,

with $\tau^{(0)}_{t}:=t$. Let us remind that the stopping times $\tau_{t}^{(m)}$ here are defined for the regime process $\{R_{t}:\,t\geq 0\}$ alone, whereas the stopping time $\eta_{t}$ defined in (2.2) is defined for the entire stochastic process $\{X_{t}:t\geq 0\}$.

Hereafter, we reserve $T\in(0,+\infty)$ for the terminal time. For $i\in\mathcal{M}$, we define the differential operator $\mathcal{L}_{i}$ by:

provided that the right-hand side is finite valued, where the gradient and Hessian are taken with respect to the state variable. Moreover, we reserve $r$ for a bounded continuous real-valued function on $[0,T]\times D\times\mathcal{M}$, and define

for $0\leq t_{1}\leq t_{2}\leq T$ and $i\in\mathcal{M}$. Note that the notation $\Theta_{t_{1},t_{2}}^{i}$ is, strictly speaking, not necessary in some instances (for instance, $\Theta_{t_{1},t_{2}}=\Theta_{t_{1},t_{2}}^{i}$ for $t\leq t_{1}\leq t_{2}\leq T$ under $\mathbb{P}_{0}^{t,{\bf x},i}$), whereas we employ it for clearer presentation as well as some ease of coding. Then, for $(i,\gamma)\in\mathcal{M}\times[0,1]$, $f(\cdot,\cdot,i)\in\mathcal{C}^{1,2}_{0}([0,T)\times D;\mathbb{R})$ and a generator matrix $Q$, the infinitesimal generator of the paired process $\{(X_{s},R_{s}):s\geq 0\}$ is given by [2, Theorem 3]:

where we denote by $\mathcal{C}_{0}^{1,2}$ the class of functions, whose first derivative in the time variable and second derivatives in the state variables are continuous with compact support. Since $Q$ is a generator matrix satisfying the $q$-property (Assumption 2.1), so is the matrix $\gamma Q(\cdot)$ for any $\gamma\in[0,1]$. If $\gamma>0$, the infinitesimal generator (2.4) is said to be weakly coupled, in the sense that the value at $(t,{\bf x},i)$ depends on $f(t,{\bf x},j)$ for $j\neq i$. If $\gamma=0$, however, such inter-dependency gets untangled, that is, the infinitesimal generator (2.4) reduces to the usual one $((\partial/\partial t)+\mathcal{L}_{i}-r(t,{\bf x},i))f(t,{\bf x},i)$ in the absence of involvement of $f(\cdot,\cdot,j)$ for $j\neq i$. If $\gamma=0$, then the corresponding generator matrix is the zero matrix in $\mathbb{R}^{p\times p}$. Hence, the probability of a switching occur within an infinitesimal-sized interval $(t,t+\Delta]$ is as negligibly small as $o(\Delta)$ due to (2.3).

We are concerned with weak approximation of the following function written in terms of mathematical expectation:

with $g:D\times\mathcal{M}\to\mathbb{R}$ and $\phi:[0,T]\times D\times\mathcal{M}\to\mathbb{R}$. For instance, Monte Carlo approximation of the probabilistic representation (3.1) is not trivial in the presence of switching ($Q\neq 0$), clearly because one needs to keep track of the additional regime process $\{R_{s}:\,s\in[t,T]\}$ as well as modulate the diffusion process $\{X_{s}:\,s\in[t,T]\}$ according to the switching. It is known [2, 22, 27] that the probabilistic representation (3.1) is a unique solution to the following initial value problem:

under the following regularity conditions on the input data.

The translation of the probabilistic representation (3.1) into the initial value problem (3.2) does not allow us to avoid the additional complexity caused by the switching, because here $p(=|\mathcal{M}|)$ partial differential equations in (3.2) are not mutually independent but weakly coupled due to the last term $\sum_{j\in\mathcal{M}\setminus\{i\}}q_{ij}({\bf x})v(t,{\bf x},j)$. The primary contribution of this work is to construct an iterative weak approximation framework (Sections 3.1 and 3.2) and associated hard upper and lower bounding functions (Section 3.3) with suitable convergence towards the target solution (3.1), without dealing with such regime-switching or inter-dependency of multiple partial differential equations.