Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
104,386 characters · 20 sections · 68 citation commands
Closed-form approximations of moments and densities of continuous--time Markov models
Continuous-time jump-diffusion processes are used in economics and finance to model the dynamics of state variables (see, e.g., bjork2008). They lead to a simple and elegant analysis of problems such as the pricing of financial assets, portfolio management and other dynamic phenomena. This comes at a big computational cost though: Many relevant characteristics, such as moments and densities, of such processes cannot be expressed in closed-form except in a few special cases. This hampers their practical use and implementation. This has led researchers to develop numerical methods for the computation of these. Broadly speaking, these methods fall in three categories: Finite--difference methods Ames1992, simulation--based methods elerian2001, brandt2002, durham2002,beskos2009,kristensen2012,sermaidis2013 and series expansions aitsahalia2002, bakshi2006, yu2007, aitsahalia2008, Filipovic2013, li2013 . This paper focuses on the latter category.
Most existing expansions proposed in the literature are application specific: Depending on the particular features of the chosen moment and model of the underlying stochastic process, different methods have been developed. One exception is kristensen2011 who developed power series expansions that covered a general class of moment functions and the transition density of multivariate diffusion processes. Their focus was on applications to option pricing but the class of expansions applies more generally. The current paper makes four contributions:
First, we demonstrate that the class of series expansions of kristensen2011 are easily extended to cover fully general continuous--time Markov models, including any jump--diffusion process. Thus, the proposal of kristensen2011 can in principle be applied to any moment of any Markov process. As part of this extension, we present a novel derivation and representation of the series expansion of kristensen2011. This new representation highlights important features of the original expansion that was perhaps not obvious from the analysis of kristensen2011.
Second, we revisit the recent work of Yangetal2019 and Wan2021 and demonstrate that in fact their proposed expansions of transition densities and option prices are special cases of kristensen2011. Thus, at a theoretical level the expansions in Yangetal2019 and Wan2021 are not new. At the same time, it should be emphasized that Yangetal2019 and Wan2021 make important contributions in terms of the practical implementation of the proposal in kristensen2011. They develop numerical algorithms that allow for fast implementation of the general method of kristensen2011 when applied to transition densities and option prices of diffusion processes and a limited set of jump--diffusion processes. As such, the current paper should hopefully clarify the relationship between these three existing papers and their relative contributions to the literature.
Third, we propose a novel numerical implementation of our series expansions when applied to general jump--diffusion models. The algorithms of kristensen2011 and Yangetal2019 are restricted to pure diffusions while the extension found in Wan2021 requires the jump component to be fully independent of the diffusive component. That is, the jump intensity and the jump sizes are not allowed to be state--dependent. Our numerical implementation allows for both to be state--dependent. We demonstrate through a series of numerical studies that our numerical method works well in practice.
Fourth, we provide a novel theory for the validity of power series expansions of moment functions of continuous--time Markov processes used here and elsewhere in the literature, including all above references to papers employing series--based approximations. Most existing theoretical results for these expansions only show that a given moment expansion converges as the time interval over which the conditional moment is defined shrinks to zero. As such existing results provide no guarantees that the approximation error will get smaller as more terms are added to expansion; in fact, nothing rules out that the approximation error may actually explode as more terms are added. For the power series expansion to be reliable, it is desirable with conditions under which the expansions converge not only over shrinking time intervals but also over a fixed time interval. We here provide guarantees for the approximations to be numerically stable as the order of the approximation grows. Our theoretical results rely on semi--group theory as also used by, e.g., Scheinkman1995 to analyze the properties of continuous--time Markov processes.
Our theoretical results demonstrate that power series expansions of Markov moments may very well not converge: The chosen moment and model has to satisfy certain regularity conditions for this to hold. In particular, we demonstrate that the expansions of transition densities and option prices proposed by kristensen2011, Yangetal2019 and Wan2021 do not converge. That is, these methods are bound to fail as the number of series terms grows. As such, the expansions proposed in these papers and the extension to general jump--diffusions developed here should be used with care. In particular, researchers may not wish to add more than, say, 4--5 terms to the expansion in order to avoid the numerical error to blow up.
The remains of the paper are organized as follows. Section (ref) presents series expansions of a broad class of moments and densities of basically any continuous--time Markov process. In section (ref), we propose a numerical implementation of the general method when applied to general jump--diffusion models. Section (ref) analyzes the theoretical properties of the power series expansion over both shrinking and fixed time distances. Section (ref) examines the numerical performance of our numerical algorithm. Section (ref) concludes. Appendix (ref) gathers all proofs.
We first provide a motivating example of a jump--diffusion model and some of the moments researchers often are interested in computing. We then proceed to consider more general framework and develop a general moment expansion method in this setting.
Consider a $d$-dimensional process, $x_{t}\in \mathcal{X\subseteq }\mathbb{R} ^{d}$ that solves the following stochastic differential equation (SDE):
where $\mu \left( x\right) $ and $\sigma \left( x\right) $ are the so-called drift and diffusion functions, respectively, $W_{t}$ is a $d$-dimensional standard Brownian motion, $N_{t}$ is a Poisson process with jump intensity $ \lambda \left( x_{t}\right) $, and $J_{t}$ captures the jump-sizes and has conditional density $\nu \left( \cdot |x_{t}\right) $. The precise form of $ \mu \left( x\right) $, $\sigma \left( x\right) $, $\lambda \left( x\right) $ and $\nu \left( \cdot |x\right) $ are chosen by the researcher according to the dynamic problem that is being considered and so are known to us. To keep notation simple, we restrict ourselves to the time--homogenous case meaning that none of the functions entering the model depend on $t$; the extension to the time--inhomogenous case can be found in Appendix (ref).
We are interested in computing conditional moments on the form
where
is a conditional moment operator. This family of operators, indexed by the time variable $t\geq 0$, constitutes a so--called semi--group of linear operators; for an overview of the general theory of semi--groups with applications to Markov processes we refer to ethier1986; for applications of semi--group theory in econometrics and finance, see AitSahaliaetal2010.\footnote{ Note that we here opt for the so--called Musiela parameterization where $t$ measures the time distance between the current and some future calendar time point. One could alternatively have defined the function of interest as, for some given $T<\infty $,
where now $\tau \leq T$ is a calendar time point. In the current time--homogenous case, it is easily seen that $\tilde{u}_{\tau }\left( x\right) =u_{T-\tau }\left( x\right) $, where $u_{t}$ was defined in ((ref)).}
The functions $r\left( x\right) $ and $f\left( x\right) $ entering ((ref))--((ref)) are chosen by the researcher according to the problem of interest. For example, with $r\left( x\right) =0$ and $f\left( x\right) =\delta \left( y-x\right) $ for some fixed $y\in \mathcal{X}$, where $\delta \left( x\right) $ is Dirac's Delta function, $u_{t}\left( x\right) =p_{t}\left( y|x\right) $, where $p_{t}$ is the transition density of $x_{t}$,
If instead we choose $r_{t}\left( x\right) =r>0$ and $f\left( x\right) =\left( \exp \left( x_{1}\right) -K\right) ^{+}$ then $u_{t}\left( x\right) $ becomes the price of a European call option with time to maturity $t$ when the state variables $x_{t}$ satisfy ((ref)) under the risk--neutral measure with the first component, $x_{1,t}$, being the log-price of the underlying asset and the short-term interest rate equals the constant $r$. When $r\left( x\right) =x_{2}$, $u_{t}\left( x\right) $ is the price of the same option but now allowing for a stochastic short--term interest rate, which is the second component of $x_{t}$.
In most cases, an analytic expression of ((ref)) is not available and $u_{t}\left( x\right) $ has to be computed using numerical approximations. To motivate our proposed approximation of $u_{t}\left( x\right) $, observe that an equivalent representation of it is the solution to a partial integro-differential equation (PIDE). An important component of this PIDE is the so--called (infinitesimal) generator $A$ of $x_{t}$ which fully characterizes the dynamics. The generator is given by, for any sufficiently regular function $f\left( x\right) $,
where, with $\sigma ^{2}\left( x\right) :=\sigma \left( x\right) \sigma \left( x\right) ^{\top }\in \mathbb{R}^{d\times d}$,
and
are the generators of the diffusive and jump component of $x_{t}$, respectively. Here, $\partial _{x_{i}}f\left( x\right) =\partial f\left( x\right) /\left( \partial x_{i}\right) $, $\partial _{x_{i},x_{j}}^{2}f\left( x\right) =\partial ^{2}f\left( x\right) /\left( \partial x_{i}\partial x_{j}\right) $ and similar for other partial derivatives.
It can then be shown, c.f. Section (ref), that $u_{t}\left( x\right) $ solves the following PIDE:
with initial condition $u_{0}\left( x\right) =f\left( x\right) $ for all $ x\in \mathcal{X}$. In the case of pure diffusions ($A_{J}=0$), the reader may recognize ((ref)) as the celebrated Feynman--Kac representation of the solution to ((ref)) which also holds for the general case of jump--diffusions. The solution to this PIDE can be represented in the following abstract manner: $u_{t}\left( x\right) =e^{\left( A-r\right) t}f\left( x\right) $, where $e^{\left( A-r\right) t}$ is the exponential of the operator $A-r$ in the sense that
We are now interested in obtaining an approximation of $u_{t}\left( x\right) $ based on a series expansion w.r.t. time $t$. A simple version of this would be a Taylor series expansion around $t=0$ on the form
for some $M\geq 1$, where the second equality uses ((ref) ). This type of moment approximations have found widespread use in the literature; see, e.g., aitsahalia2002, bakshi2006, yu2007, aitsahalia2008, Filipovic2013, li2013 . However, this expansion is not valid (well-defined) when, for example, $ f\left( x\right) $ is a non--smooth function since the domain of the operator $A_{D}$ is restricted to smooth functions, c.f. ((ref) ). The transition density and option pricing examples provided above fall in this category. We will now present a generalized version of above expansion that circumvents this issue; this is done for fully general semi--groups/Markov processes.
We take as given some semi--group\footnote{ A family of linear operators $\left\{ E_{t}:t\geq 0\right\} $ is said to be a semi--group if it satisfies (i) $E_{0}f=f$ and and (ii) $ E_{s+t}f=E_{s}E_{t}f$ for all $s,t\geq 0$..} $\left( t,f\right) \mapsto E_{t}f\left( x\right) $ of interest. It could, for example, be on the form ( (ref)) for some continuous--time Markov process $x_{t}$, not necessarily a jump--diffusion process. But we do not restrict ourselves to this case.
Suppose that $E_{t}f\left( x\right) $ is not available on closed form for a given choice of $f$; we here show how this can be approximated through a Taylor series expansion of $E_{t}f\left( x\right) $ w.r.t. $t$ when either $ f $ is sufficiently regular, where the notion of "regular" will be made clear below, or it can be expressed as the limit of a regular function. The proposal is a generalisation of the one of kristensen2011, but the derivation will be carried out using semi--group theory which simplifies the derivations substantially compared to kristensen2011 and provides new insights into the expansion.
Let $\mathcal{D}\left( E\right) $ denote the domain of $E$ and let $B$ be the infinitesimal operator of the semi--group defined as
and let $\mathcal{D}\left( B\right) $ denote the domain of $B$; that is, the set of functions $f\in \mathcal{D}\left( E\right) $ for which the above limit exists. In the motivating jump--diffusion example, $B=A-r$. For any $ f\in \mathcal{D}\left( E\right) $, we write
where as before this should be interpreted as
If the chosen $f$ satisfies $f\in \mathcal{D}\left( B^{M}\right) $, ((ref)) also holds at $t=0$ and so the following Taylor series expansion of $E_{t}f\left( x\right) $ w.r.t. $t$ around $t=0$ is valid,
where under weak conditions $\tilde{E}_{t}f\left( x\right) =E_{t}f\left( x\right) +O\left( t^{M}\right) $. This is a generalised version of ((ref)).
We are now interested in generalising this series expansion to also work when $f\notin \mathcal{D}\left( B^{M}\right) $. An important ingredient of this is to first identify/construct a smoothed version of $f\left( x\right) $ , denoted $u_{0,s}\left( x\right) $, $s\geq 0$ and $x\in \mathcal{X}$, which we require to satisfy the following two conditions:
The function $u_{0,s}\left( x\right) $ is chosen by the researcher and needs to be available on closed form for the subsequent approximation to be operational. The choice of $u_{0,s}\left( x\right) $ is application specific in the sense that Assumption A.0 has to be satisfied: Part (i) requires $ u_{0,s}\left( x\right) $ to converge towards the irregular function of interest $f\left( x\right) \notin \mathcal{D}\left( B\right) $ as $ s\rightarrow 0^{+}$. Part (ii) says that, for some $s>0$, $u_{0,s}\left( x\right) $ is sufficiently regular in the sense that it is $M$ times continuously differentiable in $s$ and each of these derivatives belongs to $ \mathcal{D}\left( B^{M}\right) $.
Assumption A.0 allows for a broad range of smoothers. One choice of $ u_{0,s}\left( x\right) $ which under great generality will satisfy A.0 is $ u_{0,s}\left( x\right) =E_{0,s}f\left( x\right) $ where $E_{0}$ is another semi--group chosen such that $u_{0,s}\left( x\right) $ is available on closed form. This choice clearly satisfies part (i) and if $E_{0}$ has similar properties as the one of interest, $E$, so that their respective generators have shared domain, then part (ii) will also hold. A simple choice of $E_{0}$, as proposed by kristensen2011, is $E_{0,s}f\left( x\right) =E\left[ f\left( x_{0,s}\right) |x_{0}=x\right] ,$where $x_{0,s}$ is another stochastic process specified by the researcher. The process $ x_{0,t}$ could, for example, be chosen as a random walk type stochastic process with transition density $p_{0,s}\left( y|x\right) =K\left( \left( y-x\right) /\sqrt{s}\right) /\sqrt{s}$, for some kernel density $K:\mathbb{R} ^{d}\mapsto \mathbb{R}^{d}$, in which case
This choice satisfies (i) and if $K$ is $M_{1}$ times differentiable then $ u_{0,s}\left( x\right) $ has the same property. The final requirement, $ u_{0,s}\in \mathcal{D}\left( B^{M_{2}}\right) $, has to be checked on a case by case basis.
Under A.0, the following identity holds:
where the second equality simply states that $u_{0,0}\left( x\right) =u_{0,s}\left( x\right) +\int_{0}^{s}\left( -\partial _{\tau }\right) u_{0,\tau }\left( x\right) d\tau $. Substituting this into ((ref)) yields
where the last equality uses the following fundamental result: If two infinitesimal operators, say, $B_{1}$ and $B_{2}$, commute in the sense that $B_{1}B_{2}f=B_{2}B_{1}f$ then $ e^{B_{1}s}e^{B_{2}t}f=e^{B_{1}s+B_{2}t}f=e^{B_{2}t}e^{B_{1}s}f$. This applies to the case of $B$ and $\partial _{s}$, $\partial _{s}B=B\partial _{s}$, since $B$ acts on $x$ while $\partial _{s}$ acts on $s$.
Finally, carry out a Taylor series expansion w.r.t. $\left( s,t\right) $ to obtain
where the order of $\partial _{s}$ and $B$ can be exchanged since $\partial _{s}^{m_{1}}B^{m_{2}}u_{0,s}\left( x\right) =B^{m_{2}}\partial _{s}^{m_{1}}u_{0,s}\left( x\right) $. The resulting approximation error is of order $O\left( s^{M_{1}}\right) +O\left( t^{M_{2}}\right) $. In particular, the above expansion will generally be more precise as $s$ gets smaller. Thus, we ideally want to choose $s$ as small as possible to reduce the approximation error. However, for the chosen value of $s\geq 0$ A.0(ii) has to be satisfied. This rules out, for example, $s=0$ when $f$ is irregular since $u_{0,0}\left( x\right) =f\left( x\right) $.
However, if the approximation error is not a major concern (which is, for example, the case if the order of approximation can be chosen sufficiently large) then one can choose $s=t$ in which case $E_{t}f\left( x\right) =e^{\left( B-\partial _{t}\right) t}u_{0,t}\left( x\right) $ and the following special case of ((ref)) can be employed,
We show in Appendix (ref) that ((ref) ) is a generalized version of the proposal of kristensen2011 which in turn contains as special cases the expansions of Yangetal2019 and Wan2021.
This section provides details regarding the practical implementation of the proposed approximation in the jump--diffusion case. We here focus on the special case of $r\left( x\right) =0$ and $s=t$, in which case $u_{t}\left( x\right) =E_{t}f\left( x\right) =\mathbb{E}\left[ f\left( x_{t}\right) |x_{0}=x\right] $ and
This is done to avoid overly complicated notation. Most of the ideas and arguments extend to the general case.
Following kristensen2011, a simple choice of $u_{0,s}\left( x\right) $ that satisfies A.1 is $u_{0,s}\left( x\right) =E_{0,s}f\left( x\right) = \mathbb{E}\left[ f\left( x_{0,s}\right) |x_{0,0}=x\right] $ where $x_{0,s}$ is chosen as the solution to an auxiliary jump--diffusion model,
where $N_{0,t}$ is a Poisson process with jump intensity $\lambda _{0}\left( x\right) $ and $J_{0,t}$ has density $\nu _{0}\left( \cdot |x\right) $. The auxiliary model should be chosen so that $u_{0,t}\left( x\right) $ is available on closed form. One such model is the multivariate Brownian motion with drift model,
where $\mu _{0}\in \mathbb{R}^{d}$ and $\sigma _{0}\in \mathbb{R}^{d\times d} $ are constants, or the multivariate Vasicek (Ornstein--Uhlenbeck) model,
both of which have a Gaussian transition density on known form. In either case,
where $p_{0,s}\left( y|x\right) $ is the transition density of the auxiliary model. For example, in the case of ((ref)),
Note that with $\mu _{0}=0$ above specification corresponds to ((ref)) with $K$ chosen as the Gaussian kernel.
Recall the two motivating examples of transition density and option price approximation. In the case of $f\left( x\right) =\delta \left( y-x\right) $, we get $u_{0,s}\left( x\right) =p_{0,s}\left( y|x\right) $. If $f\left( x\right) =\left( \exp \left( x_{1}\right) -K\right) ^{+}$, and we set $\mu _{0,1}=r-\sigma _{0,11}^{2}/2$ to ensure risk--neutrality in the auxiliary model, then $u_{0,s}\left( x\right) $ takes the form of the well-known formula for the risk--neutral expected pay-off of a call option in the Black--Scholes model,
where $d_{\pm }\left( x,s\right) =\left( x-\log \left( K\right) +\left( r\pm \frac{1}{2}\sigma _{0,11}^{2}\right) s\right) /\left( \sigma _{0,11}\sqrt{s} \right) $ and $\Phi \left( \cdot \right) $ denotes the cdf of the $N\left( 0,1\right) $ distribution.
In the pure diffusion case, where no jump component is present so that $ A_{J}=0$, analytical expressions of $\left( A_{D}-\partial _{t}\right) ^{m}u_{0,t}\left( x\right) $ are in principal straightforward to obtain relying on symbolic software packages, such as Mathematica, since $A_{D}$ is a differential operator. We refer to kristensen2011, Yangetal2019 and Wan2021 for more details on this for the two leading examples of density and option price approximations and with $u_{0,t} $ chosen as the corresponding solution under ((ref)).
Next, consider jump--diffusion models where either the diffusive component or the jump component of $x_{t}$ are state--independent; the latter case corresponds to the class of jump--diffusions considered in Wan2021.
These two cases correspond to (i) $\mu \left( x\right) =\mu $ and $\sigma ^{2}\left( x\right) =\sigma ^{2}$ are constant or (ii) $\lambda \left( x\right) =\lambda $ and $\nu \left( \cdot |x\right) =\nu \left( \cdot \right) $ are independent of $x$, respectively. In either case, we can write $x_{t}=x_{D,t}+x_{J,t}$ where the diffusive component, $x_{D,t}$, and the jump component, $x_{J,t}$, are now mutually independent. As a consequence, the two generators $A_{D}$ and $A_{J}$ commute, $A_{D}A_{J}=A_{J}A_{D}$, in which case
where
Now, consider first the case where (ii) is satisfied. In this scenario, $ x_{J,t}|x_{J,0}=x$ has density
where $\nu _{k}\left( y\right) $ is the density of the sum of $k$ independent jumps, $\sum_{i=1}^{k}J_{i}$, $J_{i}\sim \nu \left( \cdot \right) $. Since $p_{J,t}\left( y|x\right) $ is a smooth function then $ x\mapsto B_{J,t}\left( \Delta x_{D,T}+x\right) $ is also a smooth function even if $f\left( x\right) $ is irregular. Thus, if $B_{J,t}\left( x\right) $ is available on closed form then the smoothing device is not needed and we can approximate $u_{t}$ by
Similar, if (i) is satisfied then $x_{D,t}|x_{D,0}$ is a Brownian motion with drift and has Gaussian density as given in ((ref)). Because of its simple dynamics, $B_{D,t}\left( x\right) $ is available on closed form in many cases and will again be a smooth function; if so, we propose to approximate $u_{t}$ by
If closed form expressions of neither $B_{D,t}$ nor $B_{J,t}$ are available, it is still possible to simplify the computation using, for example,
assuming that closed form expressions of $e^{A_{J}\left( T-t\right) }\left( A_{D}-\partial _{t}\right) ^{m}u_{0,t}\left( x\right) $ can be computed. This last version is the one proposed by Wan2021 for jump--diffusions with state--independent jumps.
Finally, consider the general case where $A_{J}\neq 0$ and both the diffusion and jump component are state--dependent. First observe that when the jumps are state--dependent, or have a complex distribution, $ A_{J}f\left( x\right) $ cannot be evaluated analytically for a given function $f$ in general. We propose to resolve this issue by approximating the integral part of $A_{J}f\left( x\right) $, $A_{J1}f\left( x\right) =\lambda \left( x\right) \int_{\mathbb{R}^{d}}f\left( x+c\right) \nu \left( c\right) dc$, by
where $\omega _{s}$ and $c_{s}$, $s=1,...,S$, are integration weights and nodes, respectively. For example, in the case of Monte Carlo integration with $S$ random draws from $\nu $, $\omega _{s}=1/S$ and $c_{s}$ is the $s$ th draw from $\nu \left( \cdot \right) $. The resulting approximate operator $\hat{A}_{J}f\left( x\right) =\hat{A}_{J1}f\left( x\right) -\lambda \left( x\right) f\left( x\right) $ is on closed form and so we can now continue as in the pure diffusion case. Also note that $\hat{A}_{J1}f\left( x\right) \rightarrow A_{J1}f\left( x\right) $ as $S\rightarrow \infty $ which ensures that the added numerical error can be controlled by choosing $S$ large enough.
In the case that $v\left( c\right) $ belongs to the exponential family, the generator of jump component, $A_{J1}$, is well--approximated using Gauss-Hermite or Gauss-Laguerre quadrature. For example, when $J_{t}$ is i.i.d. scalar with double exponential distribution with mean zero and standard deviation $\sigma _{J}$, it follows from a change of variables that
Then, given the nodes and weights, $c_{s}^{GL}$ and $\omega _{s}^{GL}$, for the Gauss-Laguerre quadrature, the approximation takes the following form:
We use this approximation method in our numerical studies when we cannot obtain an exact expression of the integral (as discussed with standard packages such as Mathematica, or it may not be evaluated through equally standard packages such as Matlab. We find that Gaussian quadrature is more accurate and easier to implement than Monte Carlo methods with low computational cost.
With $\hat{A}_{J1}$ replacing $A_{J1}$, we can now use a symbolic software package to obtain expressions of $\left( A_{D}+\hat{A}_{J}-\partial _{t}\right) ^{m}u_{0,t}$, $m=1,2,\ldots $. For example,
where the evaluation of $\left( A_{D}-\lambda -\partial _{t}\right) ^{2}f$ and $\left( A_{D}-\lambda -\partial _{t}\right) f$ can done using symbolic methods while (here in the univariate case for simplicity)
and
We first present a general theory of series expansions on the form ((ref)) when the function $f$ is regular in the sense that $f\in D\left( B\right) $. We provide two sets of results: First, we derive an error bound for any given value $M$ of the order of the expansion. Second, we provide conditions under which the error bound vanishes as $M\rightarrow \infty $ at a given value of the time horizon $t>0$. The conditions for the second set of results come in two forms: We first provide conditions under which the proposed power series expansion converges globally, i.e., over the whole domain of $x_{t}$. These conditions are somewhat restrictive though and rule out certain models and functions of interest. We therefore proceed to examine how the approximation behaves on a given compact subset of the full domain, and show that the power series expansion is consistent over compact subsets under weak regularity conditions that most known models satisfy. We then apply the theory to moments of jump--diffusions on the form ((ref)) and provide primitive conditions under which the expansion is valid. Some of the results presented here rely on the important insights found in the unpublished work of schaumburg2004 which we are indebted to.
Next, we then proceed to analyze the "smoothed" expansion ((ref)). As in the regular case, we are able to derive an error bound for a given choice of $M$. But at the same time, this expansion is generally not consistent in the sense that it will not converge as $M\rightarrow \infty $ for a fixed value of $t>0$. This is an important result since this shows that the approximation error will eventually blow up as we increase $M$. Thus, researchers should use the generalized version with caution.
We take as given a semi--group $E_{t}:\mathcal{F\mapsto F}$ where $\mathcal{F }$ is equipped with some function norm $\left\Vert \cdot \right\Vert _{ \mathcal{F}}$. In the leading case of $E_{t}f\left( x\right) =\mathbb{E} \left[ f\left( x_{t}\right) |x_{0}=x\right] $, two standard choices of $ \left( \mathcal{F},\left\Vert \cdot \right\Vert _{\mathcal{F}}\right) $ are the following: The first is the space of bounded functions equipped with a $ \sup $ norm, $\left\Vert f\right\Vert _{\mathcal{F}}=\sup_{x\in \mathcal{X} }\left\vert f\left( x\right) \right\vert $. The second is the space of functions with second moments equipped with the following $L_{2}$ norm,$\ \left\Vert f\right\Vert _{\mathcal{F}}^{2}=\int_{\mathcal{X}}^{\infty }f^{2}\left( x\right) \pi \left( x\right) dx$ for some weighting function $ \pi \left( x\right) $. In case of $x_{t}$ being stationary, a natural choice for $\pi $ is the stationary marginal distribution in which case $\left\Vert f\right\Vert _{\mathcal{F}}^{2}=\mathbb{E}\left[ f^{2}\left( x_{t}\right) \right] $; this norm was, for example, used by Scheinkman1995.
We now formally introduce the so--called generator associated with $E_{t}$. We will here work with the so--called extended generator which is defined as follows (see, e.g., Meyn1993):
For a given function $f\in \mathcal{F}$, we will in the following frequently use $u_{t}\left( x\right) $ to denote
to economize on notation. As a first step, we show that $u_{t}\left( x\right) $ solves ((ref)) if $f\in \mathcal{D}\left( B\right) $:
or, equivalently,
The continuity condition in the second part of the theorem is satisfied under great generality when $E_{t}$ is on the form ((ref)). A sufficient condition is that the mapping $\left\{ x_{t}:t\geq 0\right\} $ is Borel measurable w.r.t. the product sigma algebra, c.f. p. 771 in Scheinkman1995. The above result, and many subsequent ones, requires the function $f$ defining $u_{t}\left( x\right) $ to satisfy $f\in \mathcal{D} \left( B\right) $. Unfortunately, it rarely easy to give an explicit characterization of $\mathcal{D}\left( B\right) $. Instead, we will often work in a smaller subspace, say, $\mathcal{D}_{0}\left( B\right) \subseteq \mathcal{D}\left( B\right) $ which is known to us; see Section (ref) for an example. One says that $\mathcal{D}_{0}\left( B\right) $ is a core of $\mathcal{D}\left( B\right) $ if it is a dense subset of the latter.
We recognize ((ref)) as a generalized version of the celebrated Kolmogorov's backward equation for jump-diffusion models. In particular, it implies that $\lim_{t\rightarrow 0^{+}}\partial _{t}u_{t}\left( x\right) =Bf\left( x\right) $. More generally, under suitable regularity conditions, $ t\mapsto u_{t}\left( x\right) $ will be $M\geq 1$ times differentiable with
in which case the following Taylor series approximation is valid,
In order for $\hat{u}_{t}\rightarrow u_{t}$ as $M\rightarrow \infty $, we need $t\mapsto u_{t}$ to be analytic:
The definition of $B$ and the convergence result ((ref)) are stated w.r.t. the chosen function norm $\left\Vert \cdot \right\Vert _{\mathcal{F}}$ introduced earlier. As we shall see, different assumptions regarding the model and the chosen function $f$ defining $u$ motivate different spaces and norms. Ideally, we would like the convergence to take place uniformly over all values of $x\in \mathcal{X}$, but this will only hold for a small set of functions $f$ and models, and so in some applications it is necessary to work with the weaker $L_{2}$ norm.
In order for $u_{t}$ to be analytic, we need as a minimum that $u_{t}$ is infinitely differentiable so that ((ref)) holds for all $m\geq 1$ . This in turn requires $B^{m}f\left( x\right) $, $m\geq 1$, to be well--defined. That is, $f\in \mathcal{D}\left( B^{m}\right) $, $m\geq 1$, where the domains are defined recursively as
The following result shows that the Taylor series $\hat{u}_{t}\left( x\right) $ is a valid approximation for any $f\in \mathcal{D}\left( B^{M+1}\right) $ and also provide an error bound for it:
We recognize the error bound as a generalized version of the one that holds for a Taylor series approximation of a $M+1$ times differentiable function. The error bound can be used to show convergence of our expansion of the transition density with $M\geq 1$ fixed as the time distance between observations, corresponding to $t$, shrinks to zero. This is the standard result found in the existing literature on expansions of moments of continuous-time processes. But, based on this result alone, the corresponding approximate moment is then only guaranteed to converge towards the exact one when high-frequency data is available. That is, when $t$ shrinks to zero as the number of observations diverge. For a fixed $t$, there is no reason why the error bound provided in the theorem will not blow up as $M\rightarrow \infty $.
We will therefore now derive conditions that guarantee convergence for a given fixed $t>0$. From Theorem (ref) we see that convergence of $\hat{u}_{t}\left( x\right) $ requires the following two conditions to be satisfied: $f\in \mathcal{D}\left( B^{\infty }\right) $ and $\left\Vert \frac{t^{m}}{m!}B^{m}f\right\Vert _{\mathcal{F}}\rightarrow 0$ as $m\rightarrow \infty $. The convergence result will generally not hold for all $t>0$. Formally, the radius of convergence is given by
Often the exact value of $T_{0}$ cannot be derived, but it may still be possible to identify a lower bound for it. Similarly, it is in many applications difficult to provide a precise characterization of $\mathcal{D}\left( B^{\infty }\right) =\bigcap\nolimits_{m=1}^{\infty }\mathcal{D}\left( B^{m}\right) $. One partial characterization is that it constitutes a core of $\mathcal{D}\left( B\right) $, c.f. Theorem 7.4.1 of davies2007, so that most functions in $\mathcal{D}\left( B\right) $ also belongs to $\mathcal{D}\left( B^{\infty }\right) $. But this provides no guarantees for that a given function in $\mathcal{D}\left( B\right) $ belongs to $\mathcal{D}\left( B^{\infty }\right) $.
Instead one may seek to identify a subset $\mathcal{F}_{0}\subseteq \mathcal{ F}$ so that (i) $\mathcal{F}_{0}\subseteq \mathcal{D}\left( B\right) $ and (ii) the image $B\left( \mathcal{F}_{0}\right) =\left\{ Bf|f\in \mathcal{F} _{0}\right\} \subseteq \mathcal{F}_{0}$. For a given $f\in \mathcal{F}_{0}$, part (i) ensures that $Bf$ is well-defined while part (ii) implies that $ Bf\in \mathcal{F}_{0}$. In particular, (i)--(ii) guarantee that $\mathcal{F} _{0}\subseteq \mathcal{D}\left( B^{m}\right) $ for all $m\geq 1$. As a consequence, $\mathcal{F}_{0}\subseteq \mathcal{D}\left( B^{\infty }\right) $ thereby providing us with a partial characterization of $\mathcal{D}\left( B^{\infty }\right) $. In particular, for any given $f\in \mathcal{F}_{0}$, we have that $t\mapsto u_{t}$ is infinitely differentiable. The following theorem states the formal result of the above analysis:
The last part of the theorem provides one sufficient condition for $u_{t}$ to be analytic. There are two tensions when seeking such a suitable set $ \mathcal{F}_{0}$: First, we would like to choose $\mathcal{F}_{0}$ as large as possible in order to guarantee convergence of $\hat{u}_{t}$ over a large set of functions. But at the same time we need to restrict $\mathcal{F}_{0}$ so that it satisfies $B\left( \mathcal{F}_{0}\right) \subseteq \mathcal{F} _{0}$. Second, to ensure a strong convergence result, we would like to choose the norm $\left\Vert \cdot \right\Vert _{\mathcal{F}}$ as "strong" as possible, e.g., as the $\sup $ norm. But establishing $T_{0}>0$ then proves more difficult.
One way of designing the function class $\mathcal{F}_{0}$ is to build it from the so--called eigenfunctions of $B$. Eigenfunctions are defined in terms of the so--called spectrum of $B$,
In particular, for any given eigenvalue $\xi \in \sigma \left( B\right) $ there exists a corresponding eigenfunction $\phi \in \mathcal{D}\left( B\right) $ so that $\left( \xi I-B\right) \phi =0\Leftrightarrow B\phi =\xi \phi $. This in turn implies that $\phi \in \mathcal{D}\left( B^{\infty }\right) $ with $B^{m}\phi =\xi ^{m}\phi $. Thus,
which is clearly analytic and so our power series expansion will converge for any eigenfunction. The following corollary shows that in principle $ \mathcal{F}_{0}$ can be chosen as the span of any given countable set of eigenfunctions:
This particular choice of $\mathcal{F}_{0}$ is in some cases somewhat restrictive in the sense that it may be only a small subset of $\mathcal{D} \left( B^{\infty }\right) $. However, in the special case of a given semi--group's spectrum being countable, we generally have that $\mathcal{F} _{0}=\mathcal{D}\left( B^{\infty }\right) $. One example of this is so--called time reversible Markov processes whose spectra are countable with the corresponding eigenfunctions forming an orthnormal basis of $\mathcal{F}$ ; see, e.g., Hansen1998. But many Markov processes are irreversible and have an uncountable spectrum in which case $\mathcal{F}_{0}$ is a proper subset of $\mathcal{D}\left( B^{\infty }\right) $.
The corollary does not guarantee that for any $f\in \mathcal{F}_{0}$ the corresponding $u_{t}\left( x\right) $ is analytic -- only that it is infinitely differentiable. To see the complications of ensuring analyticity, observe that, for any given $f\in \mathcal{F}_{0}$ with $\mathcal{F}_{0}$ defined above, $B^{m}f=\sum_{i=1}^{\infty }\alpha _{i}\xi _{i}^{m}\phi _{i}$ , $m\geq 1$, so that
Thus,
and so we need at a minimum $\sup_{i\geq 1}\left\vert \sum_{m=0}^{M}\frac{ \left( \xi _{i}t\right) ^{m}}{m!}-e^{-\xi _{i}t}\right\vert \rightarrow 0$, $ M\rightarrow \infty $. But this convergence result will generally not hold; for example, if $\xi _{i}\in \mathbb{R}$ and $\xi _{i}\rightarrow \infty $ as $i\rightarrow \infty $ then convergence will fail.
In conclusion, to ensure convergence, we need to impose restrictions on the eigenvalues/the spectrum. We will now present such a set of conditions. These will involve the so--called resolvent of the generator defined as
The first part of the theorem states necessary and sufficient conditions for $E_{t}f\left( x\right) $ to be analytic at any given $t>0$ and for any $f\in \mathcal{F}$. The conditions ((ref))--((ref)) ensure that the spectrum of $B$ is such that the convergence problem discussed before the theorem does not occur. This is a strong result but at the same time ((ref))--((ref)) are rather strong conditions. Moreover, they tend to be difficult to verify in practice since this requires knowledge of the spectrum $\sigma \left( B\right) $. Primitive sufficient conditions for them to hold are provided in the next section. Both the conditions and the results are relative to the chosen function space and norm $\left( \mathcal{F},\left\Vert \cdot \right\Vert _{\mathcal{F}}\right) $. By choosing $\mathcal{F}$ suitably small, we expect that ((ref))--((ref)) will hold in great generality. We give an example of this in Section (ref).
The second part then shows that for the subclass of functions $f$ that satisfy $f\left( x\right) =E_{\tau _{0}}g\left( x\right) $, for some $\tau _{0}>0$ and $g$, analyticity of $E_{t}f\left( x\right) $ extends to $t=0$. This part follows as a direct consequence of the first part since this implies that $E_{t}f\left( x\right) =E_{t+\tau _{0}}g\left( x\right) $ is analytic at $t=0$. The lower bound of the radius of convergence $T_{0}$ depends on the degree of smoothness of $f$, as measured by $\tau _{0}$, and the properties of the model, specifically the bound $C_{A}$ on its resolvent.
The requirement $f\in E_{\tau _{0}}\left( \mathcal{F}\right) $ is difficult to verify in a given application. In the leading case of $E_{t}f\left( x\right) =\mathbb{E}\left[ f\left( x_{t}\right) |x_{0}=x\right] $, the condition amounts to showing that there exists a solution $g\left( x\right) $ to the following integral equation $f\left( x\right) =\int g\left( y\right) p_{\tau _{0}}\left( y|x\right) dy$ for some $\tau _{0}>0$, assuming that $ x_{t}$ has a transition density $p_{t}\left( y|x\right) $. This is a so--called Fredholm equation of the first kind; conditions for a solution to this to exist are available but not easily verified in a given application. However, it can be shown that, for any given $\tau _{0}$, $E_{\tau _{0}}\left( \mathcal{F}\right) $ is dense in $\mathcal{F}$, see, e.g., Theorem 7.4.4 in davies2007, and so the result will hold for "almost every" $f\in \mathcal{F}$.
We now apply the general theory to our jump--diffusion model. In the following, let $x_{t}$ be a weak solution to ((ref)) for a given specification of $\left( \mu ,\sigma ^{2},\lambda ,\nu \right) $ with generator $A$ given in ((ref))--((ref)) and $ E_{t}f\left( x\right) =\mathbb{E}\left[ f\left( x_{t}\right) |x_{0}=x\right] $.
We first need to get a handle on the generator of the process and its domain $\mathcal{D}\left( A\right) $. A complete characterization of $\mathcal{D} \left( A\right) $ is unfortunately not possible and we will instead only work with a subset of $\mathcal{D}\left( A\right) $ where the generator takes the form ((ref)). Let $\mathcal{C}^{m}=\mathcal{C}^{m}\left( \mathcal{X}\right) $ denote the space of functions $f\left( x\right) $ with domain $\mathcal{X}$ that are $m\geq 0$ times continuously differentiable w.r.t. $x$. If $f\in \mathcal{C}^{2}$ then Ito's Lemma for jump--diffusions (see, e.g., Cont2003, Proposition 8.14) yields
where $A_{D}$ is defined in ((ref)), $\sigma _{i}\left( x\right) =\left[ \sigma _{i1}\left( x\right) ,...,\sigma _{id}\left( x\right) \right] $ while $\tau _{i}$ and and $\Delta x_{i}$ denote the time and the size, respectively, of the $i$th jump. Assuming $E_{t}\left\vert f\right\vert \left( x\right) <\infty $ and $E_{t}(\frac{\partial f}{\partial x_{i}}\sigma _{i})^{2}\left( x\right) <\infty $, $i=1,...,d$, we can take conditional expectations w.r.t. the natural filtration on both sides of the above to obtain ((ref)) with $A$ given in ((ref)). Thus, the following is a subset of the domain of the generator,
In the following we will only consider functions situated in $\mathcal{D} _{0}\left( A\right) $ and so not distinguish between the general generator and the one restricted to $\mathcal{D}_{0}\left( A\right) $. Under the assumption that $\mu $, $\sigma ^{2}$ and $\lambda $ and $f$ all belong to $ \mathcal{C}^{2m}$, we can apply Ito's Lemma repeatedly and it follows straightforwardly that
Implicit in this definition is the requirement that $\int_{\mathbb{R} ^{d}}\left\vert A^{k}f\left( x+c\right) \right\vert \nu _{t}\left( c\right) dc<\infty $ for $k=0,...,m$. Thus, a given $f\in \mathcal{C}^{2m}$ belongs to $\mathcal{D}_{0}\left( A^{m}\right) $ if relevant moments w.r.t the jump measure $\nu $ and the probability measure of $x_{t}$ exist. For example, if $f$ and all its derivatives are bounded, $\mu $, $\sigma ^{2}$ and $\lambda $ and all their derivatives are bounded by some function $V\left( x\right) \geq 0$ with $\mathbb{E}\left[ V\left( x_{t}\right) \right] <\infty $, and $ \nu $ has bounded support then $f\in \mathcal{D}_{0}\left( A^{\infty }\right) $. Similarly, if $f$ is a polynomial of order $q$, $\mu $, $\sigma ^{2}$ and $\lambda $ are linear w.r.t. $x$, $\nu _{t}$ has all polynomial moments, and $\mathbb{E}\left[ \left\Vert x_{t}\right\Vert ^{q}\right] <\infty $, $0\leq t\leq T$, then $f\in \mathcal{D}\left( A^{\infty }\right) $ .
Since $f\in \mathcal{D}_{0}\left( A^{\infty }\right) $ is necessary for our expansion to work, we will maintain the following assumption on the model:
Part (i) ensures that, under suitable moment conditions as described above, if $f\in \mathcal{C}^{\infty }$ then $f\in \mathcal{D}_{0}\left( A^{\infty }\right) $. Part (ii) is imposed to simplify subsequent arguments since it entails the following result (see Pazy, 1983, Theorem 3.2.1):
Thus, for a given jump--diffusion model satisfying A.1(ii), or any other conditions ensuring $A_{J}$ is bounded, we only need to ensure that the diffusive component is analytic. In the following, we will implicitly assume that indeed $A_{J}$ is bounded and derive conditions under which $E_{t}$ for pure diffusion processes ($A_{J}=0$) is analytic.
Ideally we would now provide primitive conditions for general jump--diffusion processes to satisfy the high-level conditions found in the theorems and corollaries stated in the previous section. This is unfortunately not possible since the spectral properties of jump--diffusions are still not fully understood. We will therefore only state results for special cases for which results do exist. At the same time we would like to emphasise that we expect the results to hold more broadly.
We first develop conditions under which polynomial moment functions are analytic. We start out with a few definitions: For a given multi-index $ \alpha =\left( \alpha _{1},...\alpha _{d}\right) \in \mathbb{N}_{0}^{d}$ and $x=\left( x_{1},...,x_{d}\right) ^{\prime }\in \mathbb{R}^{d}$ let $ \left\vert \alpha \right\vert =\alpha _{1}+\cdots +\alpha _{d}$ and $\left( x\right) ^{\alpha }=x_{1}^{\alpha _{1}}\cdots x_{d}^{\alpha _{d}}$. We then let
denote the family of polynomials of order $k$ and $\mathcal{P}_{k|\mathcal{X} }$ be these polynomials restricted to the domain of $x_{t}$. Observe here that $\mathcal{P}_{k|\mathcal{X}}$ is a finite-dimensional function space. In particular, we can choose a set of basis functions $e=\left( e_{1},...,e_{N}\right) \in \mathcal{P}_{k|\mathcal{X}}$, where $N=\dim \mathcal{P}_{k|\mathcal{X}}$, so for any $p\in \mathcal{P}_{k|\mathcal{X}}$ there exists $c=\left( c_{1},...c_{N}\right) $ so that
If $\mathcal{P}_{k|\mathcal{X}}$ satisfies the two conditions of Theorem (ref) then analyticity follows automatically from the fact that when we restrict the domain of $A$ to $\mathcal{P}_{k|\mathcal{X}}$ then it becomes a finite--dimensional operator and therefore bounded:
The second result provides primitive conditions for the high--level assumptions ((ref))--((ref)) to hold in the context of jump diffusions, where $E_{t}^{\ast }$ and $A^{\ast }$ denotes the so--called adjoint operators of $E_{t}$ and $A$, respectively:
The time--reversibility condition implies that $A$'s spectrum is discrete and contained in the negative half--line which suffices for ((ref))--((ref)) to hold. The three examples referred to in the second part of the last theorem are time--homogenous scalar diffusions, multivariate factor diffusion models, and a restricted class of multivariate diffusions; see Scheinkman1995 for the precise details.
Note here that the corollary imposes no smoothness conditions on $\mu $, $ \sigma ^{2}$ and $f$. This is because that $Af$ may still be well--defined even without smoothness, c.f. above discussion of $\mathcal{D}\left( A\right) $. However, its particular form in these cases is generally unknown to us. Thus, in order to compute $Af$ in practice we restrict ourselves to smooth models, as in Assumption A.1(i), and smooth choices of $f$, as in $ \mathcal{D}_{0}\left( A\right) $.
Our third result again uses Theorem (ref) but focuses on a different class of "test functions" to obtain results for such models. We restrict the function set to
where $\partial _{x}^{\alpha }f=\partial ^{\alpha }f/\left( \partial x^{\alpha }\right) $, which we equip with the norm $\left\Vert f\right\Vert _{\mathcal{F}_{0}}=\sup_{\left\vert \alpha \right\vert \geq 0}\left\Vert \partial _{x}^{\alpha }f\right\Vert _{\mathcal{F}}.$Importantly, if $f\in \mathcal{F}_{0}$ then, for any $\alpha \in \mathbb{N}_{0}^{\infty }$, $ \partial _{x}^{\alpha }f\in \mathcal{F}_{0}$ with $\left\Vert \partial _{x}^{\alpha }f\right\Vert _{\mathcal{F}_{0}}\leq \left\Vert f\right\Vert _{ \mathcal{F}_{0}}.$This property of $\mathcal{F}_{0}$ ensures that if $\mu $ and $\sigma $ in $\mathcal{F}_{0}$ then $Af\in \mathcal{F}_{0}$ for all $ f\in \mathcal{F}_{0}$ and so $\mathcal{F}_{0}\subseteq \mathcal{D}\left( A^{\infty }\right) $. Moreover, the generator, when restricted to $\mathcal{F }_{0}$, is bounded and so the radius of convergence is infinite:
Note here that convergence holds for all $t>0$ and that the convergence rate is super-geometric. Moreover, the result allows for a broad class of non-linear multivariate diffusion models. On the other hand, it rules out unbounded drift and diffusion terms.
The above results are strong in the sense that they guarantee convergence w.r.t a function norm over the full state space $\mathcal{X}$. But at the same time they are restrictive in that they do not apply to general multivariate jump--diffusion models. One way of allowing for a broader class of models and functions is to restrict attention to solutions defined on a bounded subset of $\mathcal{X}$ leading to the following class of so--called localized Cauchy problems. We here focus on the case of pure diffusions since for this class of models results exist on analytic solutions on bounded sets.
Let $\mathcal{X}_{0}\subseteq \mathcal{X}$ be a bounded open set and let $ u_{t}^{\ast }\left( x\right) $ be a function chosen by the researcher which satisfies $u_{0}^{\ast }\left( x\right) =f\left( x\right) $. We then consider the following "trimmed" version of the Cauchy problem for diffusion models:
with initial condition $w_{0}\left( x\right) =f\left( x\right) $ for $x\in \mathcal{X}_{0}$. We now only require the solution $w_{t}\left( x\right) $ to solve the Cauchy problem on a bounded open subset $\mathcal{X}_{0}$ of the full domain $\mathcal{X}$ and then pin down its behaviour outside of $ \mathcal{X}_{0}$ through the pre-specified function $u^{\ast }$. The class of problems on the form ((ref))--((ref)) can be described by a semi--group $E_{t}$ so that $w_{t}=E_{t}f$. By choosing $ \mathcal{X}_{0}$ as a bounded set, the requirements for the semi--group to be analytic becomes a lot less restrictive and essentially requires $\mu $ and $\sigma ^{2}$ to be sufficiently smooth; see, e.g., Chapter 3 in Lunardi1995. The following theorem states the precise conditions:
This provides simple and relatively weak conditions under which a series expansion of $w_{t}$ will converge. But will such series expansion be a good approximation to $u_{t}$? By eq. ((ref)) together with the initial condition
Thus, under the conditions of the theorem, our proposed power series approximation shares derivatives with $u_{t}$ on $\mathcal{X}_{0}$. At the same time, the solution $w_{t}$ will generally differ from the global solution $u_{t}$. However, if we restrict $f\in \mathcal{D}^{\infty }\left( A_{D}\right) $ then $\left. \partial _{t}^{m}u_{t}\left( x\right) \right\vert _{t=0^{+}}=A_{D}^{m}f\left( x\right) $ and so $w_{t}\left( x\right) =u_{t}\left( x\right) $, $x\in \mathcal{X}_{0}$, and the power series will be consistent on $\mathcal{X}_{0}$. In particular, if we can show that $w_{t}\left( x\right) $ is analytic on $\mathcal{X}_{0}$ then the same will hold for $u_{t}\left( x\right) $ when considered as a function with domain $\mathcal{X}_{0}$. This result combined with Lemma (ref) shows that our power series expansions converges for a very broad class of diffusion models over bounded subsets of their domains.
Finally, we provide an analysis of smoothed expansions on the form ((ref)). First, by following the same arguments as in Theorem (ref), it is easily shown using Taylor's Theorem that if $ u_{0,s}$ satisfies A.0 then $\hat{u}_{t}\left( x\right) =\hat{E}_{t}f\left( x\right) $ given in ((ref)) with $M_{1}+M_{2}=M$ satisfies
One could now hope for that as long as $u_{0,s}\left( x\right) $ is sufficiently regular then the expansion would converge under conditions similar to the ones in the "regular" case analyzed in the previous section. This is unfortunately not the case. To see this, observe that in order for the expansion to be asymptotically valid $s\mapsto u_{0,s}\left( x\right) $ has to be analytic so that
where $r\left( x\right) $ is the remainder term from a $M_{1}$th order Taylor expansion of $s\mapsto u_{0,s}$ around $s=0$. If $f\notin \mathcal{D} \left( B\right) $ and, for some $M_{1}\geq 1$, $\hat{f}\left( x\right) \in \mathcal{D}\left( B\right) $ then obviously $r\left( x\right) \notin \mathcal{D}\left( B\right) $. Thus, as $M$ grows large enough, we must have $ \hat{f}\notin \mathcal{D}\left( B\right) $ in which case $\hat{E}_{t}f\left( x\right) =\sum_{m_{2}=0}^{M_{2}}\frac{t^{m_{2}}}{m_{2}}B^{m_{2}}\hat{f} \left( x\right) $ is not well-defined. In practice, we expect $\hat{E} _{t}f\left( x\right) $ in ((ref)) to become numerically unstable as $M_{2}\rightarrow \infty $. That is, the numerical error will start blowing up.
This demonstrates that the proposed series expansions of irregular functions such as densities and option prices should be used with caution: As more terms are added to the expansions, they will most eventually become numerically unstable and produce unreliable estimates. However, as we shall see in the next section, the expansions still work well when a reasonably small number of terms are used.
We assess the performance of our approximations when applied to the problem of option pricing when the underlying asset's dynamics are described by a stochastic volatility model with jumps under the risk--neutral measure. We consider the following class of asset pricing models where the log--price $ s_{t}$ of a given asset exhibits both stochastic volatility and jumps,
where the volatility process $v_{t}$ is solution to either
or
Here, $\mu =r-\delta $ where $r$ and $\delta $ are the risk-free rate and the constant dividend, respectively. To ensure that the model has a well-defined solution, $\kappa _{V}$, $\alpha _{V}$, $\sigma _{V}$ are restricted to be positive and $1/2\leq \beta \leq 1$.
The jump component consists of a Cox process $N_{t}$ with a jump intensity function given by $\lambda \left( v\right) =\lambda _{0}+\lambda _{1}v$, and a random variable $J_{t}$ with support $[-1,\infty )$, and expectation $\bar{ J}$. We include $-\lambda \left( v_{t}\right) \bar{J}$ in the drift as a compensator such that the jump part is a martingale. For example, if $J+1$ is chosen to be log-normally distributed with parameters $\mu _{J}$ and $ \sigma _{J}$, then $\bar{J}=\exp \left( \mu _{J}+\sigma _{J}^{2}/2\right) -1$ . Special cases of this model include Merton1976, where both volatility and jump intensity are constant, $v_{t}=\sigma _{0}$ and $\lambda \left( v\right) =\lambda _{0}$. Eq. ((ref)) together with either ( (ref)) or ((ref)) is a special case of ((ref)) with $x_{t}=\left( s_{t},v_{t}\right) $.
This class of models subsumes models in andersen2002 and Wan2021 as well as a number of other special cases. Compared to andersen2002, our specification allows the variance process to be the non-affine continuous-time GARCH model ($\beta =1$) and the CEV model ($ 1/2<\beta <1$). Also, compared to Wan2021, we allow for state-dependent jump intensity ($\lambda _{1}\neq 0$) which they rule out.
We consider a European call option with payoff $f\left( s_{T}\right) \equiv \max \left\{ \exp \left( s_{T}\right) -K,0\right\} $ at maturity time $T>0$, where $K=100$ is the strike price. With the above model formulated under the so--called risk--neutral measure, let $u_{\Delta }\left( s,v\right) =E\left[ f\left( s_{T}\right) |s_{T-\Delta }=s,v_{T-\Delta }=v\right] $ be the expected risk--neutral pay-off the option expires in $\Delta $ time units and the current log stock price and volatility is $s$ and $v$, respectively. Within the above class of models for $s_{t}$, no closed-form formula for the option price is available. We here implement our proposed series expansion of the unknown price, $\hat{u}_{t}\left( s,v\right) $, as given in ((ref)), where we choose $u_{0,\Delta }$ as the pay-off under the Black- Scholes model as given in ((ref)).
In the case of state--dependent jumps, we need to compute the integration part of $A_{J}$ using numerical methods. Since $\log \left( J_{t}+1\right) $ is i.i.d. with normal distribution with mean $m_{J}$ and standard deviation $ \sigma _{J}$ for all models in this section, we use the Gauss-Hermite quadrature with different numbers of nodes and weights, whose values are fixed after choosing the number of nodes and weights, c.f. Section (ref).
To assess the numerical performance of our expansion, we will use as benchmark the option price obtained via Monte Carlo methods, where the total number of simulation trials $10,000,000$ and the time-step is $10,000$ per year, see Chapter 3 in Giesecke2018 for details. We measure the accuracy of the approximations by the maximum absolute error and the absolute percentage error defined as follows: $\max_{S\in \left[ 90,110 \right] }\lvert \hat{u}_{\Delta }\left( s,v\right) -u_{\Delta }^{MC}\left( s,v\right) \rvert $ and $\max_{S\in \left[ 90,110\right] }\lvert \hat{u} _{\Delta }\left( s,v\right) -u_{\Delta }^{MC}\left( s,v\right) \rvert /u_{\Delta }^{MC}\left( s,v\right) $, respectively, where $\hat{u}_{\Delta }\left( s,v\right) $ and $u_{\Delta }^{MC}\left( s,v\right) $ are the series expansion and the Monte Carlo version of the option price, respectively.
We consider increasingly challenging experiments, aiming to assess the resilience of our method to the approximation of option prices under increasingly complex models.
In this subsection, we explore the performance our method when jumps are state-independent ($\lambda _{1}=0$).
In Figure (ref), we depict the approximation errors resulting from our method for ((ref))--((ref)) with $\beta =0.5$ across different levels of the current asset price. As in Wan2021, the parameter values used in this experiment are chosen as the estimates reported in Eraker2004, which are displayed in the figure legend. From left to right, the time to maturity ranges from $\Delta :=T-t=1/52$, 1/12, and 1/4, respectively. In the top three panels, the maximum absolute error has been plotted for 1st, 2nd, 3rd, and 4th order approximation, respectively; whereas the horizontal axis denotes the number of nodes and weights for the Gauss-Hermite quadrature. For the bottom three panels, the vertical axis denotes the absolute percentage error; whereas the horizontal axis denotes the stock price.
We make the following observations: First, for all maturities, as $M$ increases, the approximation error decreases. Second, for a given order of approximation, our method is more accurate as the time to maturity decreases. Third, one can achieve accurate approximations with small number of nodes and weights used in the quadrature approximation of the jump component. For small time to maturity ($t=\Delta =1/52$), it is sufficient to use the Gauss-Hermite quadrature with 10 nodes and weights, but, for larger time to maturity ($t=\Delta =1/12$, 1/4), only 4 or 5 nodes and weights. It indicates that the error in computing the integration part of the jump component is smaller than the error of our Taylor series approximation as the number of nodes and weights for the Gauss-Hermite quadrature increases.
Figure (ref) and (ref) investigate the numerical performance of our approximation for the call option under the stochastic volatility model ((ref))--((ref)) with the CEV ($\beta =0.8$) and GARCH ($\beta =1$) specifications of variance, respectively. The parameters for the CEV and GARCH specification are from aitsahalia2007 and Yang2017, respectively, but we added or changed the parameters for the jump part, which is the same as in Wan2021.
For $\Delta =1/52$ and $\Delta =1/12$, the performances of the approximation for both two models share three patterns arose in the outcome in Figure (ref). However, for longer time-to-maturity, $\Delta =1/4$, the higher order of approximation does not guarantee smaller approximation error. In general, the performance of the approximation error is good with shorter maturities and/or $\beta $ takes on a relatively small value.
In this subsection, we provide results for the case where the option price is computed under models with state-dependent jump intensities ($\lambda _{1}\neq 0$).
In Figures (ref)--(ref), we depict the relative error of the approximation for the same three models considered in the previous subsection, except that now $\lambda _{1}=1,10,30$, when time-to-maturity equal to one month, $\Delta =1/12$. In each figure, from left to right, the state dependency of jump intensity ranges $\lambda _{1}=1,10,30$. Overall, the approximation errors for each of the three models are comparable to that of the same model with state-independent jumps ($\lambda _{1}=0$). Furthermore, the absolute percentage error is smaller for all orders of approximation for larger $\lambda _{1}$. It indicates that the magnitude of $\lambda _{1}$ affects the level of option prices but does not affect the approximation errors. That is, the performance of our approximation is not very sensitive to the degree of state dependence of the jumps as measured by the value of $\lambda _{1}$.
Next, we consider the performance when $v_{t}$ solves the log--volatility model ((ref)) with parameters chosen as $\left( r,\delta ,\kappa _{V},\alpha _{V},\sigma _{V},\rho \right) =\left( 0.0304,0,0.0145,-0.8276,0.1153,-0.6125\right) $ and $\left( \lambda _{0},\mu _{J},\sigma _{J}\right) =\left( 0.0137,-0.000125,0.015\right) $; these are the estimates reported in andersen2002. Figures (ref) and (ref) display the relative error of the approximation for the call option under this model for different values of $\Delta $ and $\lambda _{1}$ with $v_{0}=\alpha _{V}=-0.8276$. We see that even for the 2nd order approximation, the approximation error is quite small for all choices of time--to--maturity and $\lambda _{1}$. The plotted errors are now more ragged which we conjecture is due to bigger numerical errors in the Monte Carlo benchmark that we use for comparison.
We here wish to examine the robustness of our method when applied to more complex models that go beyond one-factor volatility. We consider the stochastic volatility model with two factors for the volatility used in filipovic2016. In their specification, the dynamics of $s_{t}$ under the risk--neutral measure are given by
where $W_{1}$, $W_{2}$, and $W_{3}$ are mutually independent standard Brownian motions. Compared to the models of the previous subsection, there is a second variance factor $m_{t}$, which represents a stochastic level around which $v_{t}$ reverts. The jump component consists of: (i) $N_{t}$, a Cox process with a bounded intensity function given by $\lambda \left( v,m\right) =\lambda _{0}+\lambda _{1}v+\lambda _{2}m$, and the variance jump size $J^{V}$ is exponentially distributed with parameter $\mu _{J}^{V}= \mathbb{E}\left[ J_{t}^{V}\right] $.
Figure (ref) reports the performance of our approximation for different times to maturity, and with different numbers of nodes for Gauss-Hermite quadrature. The parameter values we used are estimates in aitsahalia2020. The performance of the approximation shares the same patterns that we found in Figure (ref). It indicates that the performance of the approximation is still very good when we add more factors to the volatility specification.
This paper provides a general framework for developing and analyzing series expansions of moments of continuous-time Markov processses, including jump-diffusions. The expansions come in two versions depending on the features of the moment. For "regular" moments, we provide conditions under which the corresponding expansion will converge towards the actual moments as more terms are added. For the "smoothed" expansion, no such theoretical guarantees exist: The expansion will eventually become imprecise as the number of terms grows. A numerical study shows that the smoothed expansions still work well in practice when a relatively small number of terms are used in its implementation.