Maximum Entropy Evaluation of Asymptotic Hedging Error under a Generalised Jump-Diffusion Model

Abstract

In this paper we propose a maximum entropy estimator for the asymptotic distribution of the hedging error for options. Perfect replication of financial derivatives is not possible, due to market incompleteness and discrete-time hedging. We derive the asymptotic hedging error for options under a generalised jump-diffusion model with kernel bias, which nests a number of very important processes in finance. We then obtain an estimation for the distribution of hedging error by maximising Shannon’s entropy subject to a set of moment constraints, which in turn yields the value-at-risk and expected shortfall of the hedging error. The significance of this approach lies in the fact that the maximum entropy estimator allows us to obtain a consistent estimate of the asymptotic distribution of hedging error, despite the non-normality of the underlying distribution of returns.

1. Introduction

The theory of pricing and hedging options has been the centre of attention in modern mathematical finance since the seminal Black–Scholes model. It provides a theoretical value and hedging strategy for European options, under the key assumption that there exists a trading strategy that constructs a portfolio that perfectly replicates the pay-off of an option. Furthermore, Black and Scholes assume the underlying stock price follows a geometric Brownian motion, and trading may take place in continuous time. With these assumptions, they show that the initial value of the replicating portfolio provides the initial price of the option. Moreover, the Black–Scholes analysis demonstrates that an option can be created synthetically by dynamically trading in the underlying asset. Nevertheless, it is well accepted that the perfect replication of options by any self-financing strategy is impossible, due to market incompleteness as well as discrete-time hedging. These two sources of error are termed the jump error and the gamma error.

Many researchers, (see for example Hubalek et al. (2006), Cont et al. (2007), and Kennedy et al. (2009)) have studied the problem of hedging an option in an incomplete market, particularly where stock prices may jump. It is well understood that, except in very special cases, martingales with respect to the filtration of discontinuous processes cannot be represented in the form of a unique self-financing strategy, which leads to market incompleteness. At the jump time, both the model price of the option and value of the hedging portfolio jump. The former is a non-linear function of the stock price but the latter is a linear function of the stock price. Therefore, the jump induces a discrepancy between the value of the option and its replicating portfolio, and thus leads to a jump error. Furthermore, these researchers generally assume that the hedging portfolio can be continuously rebalanced, which is only possible in the absence of transcription costs. In practice, this level of liquidity is not possible and market practitioners rebalance their hedging portfolio using discrete-time observations, just a few times per trading day. The discrete hedging of derivatives securities leads to the gamma error. This error is not easy to measure because the stochastic analysis techniques are not available in discrete time. Wang et al. (2015) provides a recent analysis of literature regarding hedging in a discrete time incomplete market. We base our research on the seminal work of Bertsimas et al. (2000), Hayashi and Mykland (2005) developed a methodology to analyse the discrete hedging error in a continuous-time framework using an asymptotic approach. This methodology was further developed by Tankov and Voltchkova (2009) who investigate the gamma risk via establishing a limit theorem for the renormalised error when the discretisation step tends to zero. Additionally, Rosenbaum and Tankov (2014) discuss the optimality conditions of discretised hedging strategies in the presence of jump.

In this paper we contribute to the literature by approaching the problem from a different angle. We characterise the risk in a dynamic hedge of options through the asymptotic distribution of hedging error. In particular, we investigate the case of conventional delta hedging for a European call option, although other types of options may be treated in a similar manner. Furthermore, we obtain an estimation for the distribution of hedging error by maximising Shannon’s entropy (Shannon (1948)) subject to a set of moment constraints, which in turn yield the value-at-risk (VaR) and expected shortfall (ES) of the hedging error, two widely-used risk metrics in finance. In the literature there exist two dominant approaches for constructing the distribution of hedging error, namely, the parametric and non-parametric approaches. The new approach that we propose in this paper chooses the probability distribution with the most uncertainty, or maximum entropy, subject to what is known. This allows us to obtain a consistent estimate of the asymptotic distribution of hedging error, despite the non-normality of the underlying distribution of returns. As a result, we can drive a very generalised modelling framework, which can be applied in different areas of derivatives pricing.

We first extend the methodology introduced in Hayashi and Mykland (2005) to model the asymptotic hedging error for vanilla call options when the underlying asset is governed by a generalised jump-diffusion model with kernel bias. The class of kernel biased completely random measures is a wide class of jump-type processes that can nicely be represented by a generalised kernel-biased mixture of Poisson random measures. The importance of using a kernel biased completely random measure is to derive a variety of forms of distortion of jump sizes through the kernel link function. For example, Fard and Siu (2013) discuss the importance of this representation as it provides great flexibility in modelling different types of finite and infinite jump activities compared with existing studies.

Next, we estimate the probability density function of the hedging error using the maximum entropy (ME) methodology. It allows the agent to effectively combine the aforementioned risk factors dynamically over time to update its belief about the possible distribution governing the terminal hedging error. Maximum entropy is widely used in estimation and information theory, in which beliefs are updated so that the posterior coincides with the prior as closely as possible. Furthermore, ME methodology only updates those aspects of beliefs for which new evidence was gained (cf., Cover and Thomas (2012), Saghafian and Tomlin (2016) and references therein). Despite its widespread use in fields such as estimation theory, physics, statistical mechanics, and information theory among others, it is only recently that researchers have begun to appreciate its usefulness in econophysics. For instance, Xi et al. (2014) use the maximum entropy model to study business cycle synchronisation of the G7 economic system. Gzyl and Mayoral (2016) uses the maximum entropy principle to develop a non-parametric method of determining the prices of the zero coupon bonds, when the only information available consists of the prices of a few coupon bonds. Chan (2009) proposes a general modelling framework for the EM algorithm in approximating the distribution of financial returns in order to develop entropy-based risk metrics; Geman et al. (2014) and Xu et al. (2014) expand on Chan (2009) and develop a robust optimisation framework for portfolio construction; Mistrulli (2011) and Zhou et al. (2013) applies the entropy maximisation principle to measure financial contagion and systemic risk. Geman et al. (2014), who use ME density in the VaR context, make an interesting observation that the real world is mostly ignorant about the importance of true probability distributions. They further point out that historically, finance theory has had a preference for parametric, less robust, methods. An approach that is based on distributional and parametric certainties may be useful for research purposes but does not accommodate responsible risk taking. Their study shows the importance of the use of true probability distributions in VaR calculations.

The remainder of this paper is structured as follows. In Section 2 we present the calculations for pricing European options under a generalised jump-diffusion model with kernel bias. Furthermore, we generalise the Hayashi and Mykland (2005) framework to derive the asymptotic hedging error, stemmed from the market incompleteness and discrete hedging. In Section 3, we obtain an estimation for the distribution of hedging error by maximising Shannon’s entropy subject to a set of moment constraints, which in turn yield the value-at-risk and expected shortfall of the hedging error. Section 4 provides a numerical analysis to highlight the applicability of the method. Section 5 concludes the paper.

2. Modelling Framework

The inadequacy of the constant volatility of the Black–Scholes (henceforth, BS) option valuation model to replicate the characteristics of observed option prices is empirically evidenced in the finance literature. For example, the volatility smile for equity options exhibits a consistent pattern illustrating that implied volatilities are lower for options with higher strike prices. This is particularly the case particularly for short maturity options. Several methods are suggested to accommodate the volatility smile in option pricing, but, stochastic volatility models (see e.g., Heston 1993; Hull and White 1987) and jump diffusion models (see e.g., Bates 1991; Merton 1976; Naik and Lee 1990) have become popular alternatives to BS’s constant volatility.

Many studies provide evidence that stochastic volatility models are a significant improvement (cf. Dotsis et al. (2007)). Nevertheless, their applicability in empirical studies suffers from two major drawbacks, namely, implausible correlation structure between returns and volatility, and excessively high “volatility of volatility” (cf. Bates (2000); Diavatopoulos et al. (2012); Dotsis et al. (2007)). Amongst these, the seminal work of Bakshi et al. (1997) reviews a range of alternative option pricing models, and conclude that stochastic volatility is of “first-order importance in improving upon the [Black–Scholes] formula”, but adding jumps may provide further improvement based on an out-of-sample test.

As a result, a large body of the recent literature has focused on augmenting stochastic volatility models with jumps. For example, Bakshi et al. (2003), Eraker et al. (2003), and Bakshi et al. (2012) conduct a set of empirical investigations for the case of index options, evidencing the importance of adding the jump component in the returns process. Chang et al. (2013) finds that joint time series data of the underlying S&P 500 index and options on it strongly reject a stochastic volatility model without jumps. Shanahan et al. (2016) price long-maturity equity linked products, as a composition of three embedded options, using the Meixner process (to capture jumps) and a diffusion process to capture stochastic volatility.

On the other hand, option pricing with both jumps and stochastic volatility is rather cumbersome, as analytical solutions are rarely achievable. Even in the relatively simple context of vanilla American options, the studies of Bakshi et al. (2003) and Bakshi et al. (2012) uses European option calculations on a restricted set of out-of-the-money options.

Recently, some theoretical improvement have been offered in the literature, such as the polynomial approximation of stochastic volatility in Shanahan et al. (2016). However, due to the absence of supportive empirical evidence, we chose a more practical alternative, proposed by Andersen and Andreasen (2000). By so doing, we aim to extend the analysis of Dupire (1994) to the case of jumps, whereby a model combining jumps with a deterministic volatility is developed. This type of model captures the observed behaviour of implied volatilities (see e.g., Andersen and Andreasen 2000). In addition, it is much easier to use for contracts where no analytic solution is readily available. Moreover, note that adopting a deterministic local volatility function approach has an advantage in calibrating European option prices. Indeed, instead of solving a set of backward equations, one for each option of a different strike and maturity, one only needs to solve a single one-dimensional forward equation for all strikes and maturities.

2.1. Financial Markets

We fix a complete probability space $(\Gamma ,\mathcal{F},\mathbb{P})$, where $\mathbb{P}$ is the real-world probability measure. Let $\mathcal{T}$ denote the time index set $[0,T]$ of the economy. Let ${\left\{r\left(t\right)\right\}}_{t\in \mathcal{T}}$ be the instantaneous market interest rate of a money market account. Then, the dynamics of the value of the risk-free asset, ${\left\{B_t\right\}}_{t\in \mathcal{T}}$ would be:

$$\frac{d{B}_t}{B_t}=r\left(t\right)dt,B_0=1.$$

To model different types of finite and infinite jump activities, we adopt the kernel biased representation of completely random measures of James (2002, 2005).

To proceed, consider the measurable space $(\mathcal{T},\mathcal{B}(\mathcal{T}\left)\right)$, where $\mathcal{B}\left(\mathcal{T}\right)$ is the Borel $\sigma$-field (generated by the open subsets of $\mathcal{T}$). We denote by ${\mathcal{B}}_0$ the family of Borel sets $U\in {\mathcal{R}}^{+}$, with closure $\overline{U}$ not containing 0. Define $\mathcal{X}=\mathcal{T}\times {\mathcal{R}}^{+}$. Under these definitions, the measurable space $(\mathcal{X},\mathcal{B}(\mathcal{X}\left)\right)$ is explicitly given by $(\mathcal{T}\times {\mathcal{R}}^{+},\mathcal{B}\left(\mathcal{T}\right)\otimes {\mathcal{B}}_0)$.

Let $N(.,U)$ define a Poisson random measure for all $U\in {\mathcal{B}}_0$. We denote by $N(dt,dz)$ the differential form of measure $N(t,U)$. In addition, define $\rho \left(dz\right|t)$ as a Lévy measure that depends on t. Let $\eta$ be a $\sigma$-finite (nonatomic) measure on $\mathcal{T}$. Following James (2005), we assume that there exists an arbitrary positive function $h\left(z\right)$ on ${\mathcal{R}}^{+},$ which along with $\rho$ and $\eta$ are chosen such that:

$${\displaystyle \sum _{i=1}^N}{\int }_{\mathcal{B}}{\int }_{{\mathcal{R}}^{+}}min(h\left(z\right),1)\rho \left(dz\right|t)\eta \left(dt\right)<\infty .$$

and $h^2\left(z\right)\le {ϝ}_t\rho \left(z\right)$, where $ϝ$ is a cádlág ${\mathcal{F}}_t$-adapted process. Define the intensity measure:

$$\nu (dt,dz):=\rho \left(dz\right|t\left)\eta \right(dt),$$

as well as the kernel biased completely random measure:

$$\mu \left(dt\right):={\int }_{{\mathcal{R}}^{+}}h\left(z\right)N(dt,dz).$$

The latter is a kernel-biased Poisson random measure $N(dt,dz)$ over the state space of the jump size ${\mathcal{R}}^{+}$ with the mixing kernel function $h\left(z\right)$. We can replace the Poisson random measure with any random measure and choose some quite exotic functions for $h\left(z\right)$ to generate different types of finite and infinite jump activities. Let ${\left\{W_t\right\}}_{t\in \mathcal{T}}$ denote a standard Brownian motion on $(\Omega ,\mathcal{F},\mathbb{P})$ with respect to the $\mathbb{P}$-augmentation of its natural filtration ${\mathcal{F}}^W:={\left\{{\mathcal{F}}_t^W\right\}}_{t\in \mathcal{T}}$. Let ${\tilde{N}}_{X_t}(dt,dz)$ denote the compensated Poisson random measure defined by:

$$\tilde{N}(dt,dz)=N(dt,dz)-\rho \left(dz\right|t)\eta \left(dt\right).$$

Let ${\mu }_t$ and ${\sigma }_t$ denote the drift and volatility of the market value of the underlying asset, respectively. Consider a random jump process $A:=\left\{A\right(t\left)\right|t\in \mathcal{T}\}$, such that:

$$d{A}_t=A_{t^{-}}\left[{\mu }_{t}dt+{\sigma }_{t}d{W}_t+{\int }_{{\mathcal{R}}^{+}}h\left(z\right)\tilde{N}(dt,dz)\right],$$

(1)

where $A_0=0$. We assume under $\mathbb{P}$ the price process ${\left\{S_t\right\}}_{t\in \mathcal{T}}$ is defined as $S_t:=exp\left(A_t\right)$ so that:

$$d{S}_t=\left({\mu }_t+\frac{1}{2}{\sigma }_t^2\right)dt+{\sigma }_{t}d{W}_t-{\int }_{{\mathcal{R}}^{+}}\left\{h\left(z\right)-e^{h\left(z\right)}+1\right\}\rho d\left(z\right|t)\eta \left(dt\right)+{\int }_{{\mathcal{R}}^{+}}\left(e^{h\left(z\right)}-1\right)\tilde{N}(dt,dz),$$

(2)

with $S_0=1$.

2.2. Esscher Transform

It is well known that no arbitrage opportunities are necessary for the determination of a unique equivalent risk neutral martingale measure, thus ensuring the fair valuation of the option (Pilska 1997). In this paper, we emphasise on the possibility of having incomplete markets, and as such there may be more than one equivalent martingale measure, i.e., more than one no-arbitrage price. There are different methods to price and hedge derivative securities in incomplete financial markets. For example, one can choose an equivalent martingale measure by minimising the quadratic utility of the terminal hedging errors (see e.g., Follmer and Sondermann 1986; Follmer and Schweizer 1991; Schweizer 1995). One can also adopt an economic approach based on the marginal rate of substitution to select a pricing measure via a utility maximisation problem (see e.g., Davis 1997). Finally, one may employ the minimum entropy martingale measure method to determine the equivalent martingale measure (see e.g., Avellaneda 1998; Fard and Siu 2012; Frittelli 2000).

In this study we employ the Esscher transform to determine an equivalent martingale measure for the valuation of the option (see Gerber and Shiu 1994). The method provides market practitioners with a convenient and flexible way to value options. For example, it has been shown in Elliott et al. (2005) that for exponential Lévy models, the Esscher martingale transform for the linear process is also the minimal entropy martingale measure, i.e., the equivalent martingale measure which minimises the relative entropy, and this measure also has the property of preserving the Lévy structure of the model. In the framework of exponential Lévy models, the study of equivalent martingale measures the relationships and optimality properties, which has been developed in several directions, see Esche and Schweizer (2005), Hubalek and Sgarra (2006), and Tankov (2003) and the references therein.

Let ${\mathcal{F}}^A:={\left\{{\mathcal{F}}_t^A\right\}}_{t\in \mathcal{T}}$ and ${\mathcal{F}}^S:={\left\{{\mathcal{F}}_t^S\right\}}_{t\in \mathcal{T}}$ denote the $\mathbb{P}$-augmentation of the natural filtration generated by A and S, respectively. Since, ${\mathcal{F}}^A$ and ${\mathcal{F}}^S$ are equivalent, we can use either one as an observed information structure. Write $B\left(\mathcal{T}\right)$ for the Borel $\sigma$-field of $\mathcal{T}$ and let $BM\left(\mathcal{T}\right)$ denote the collection of $B\left(\mathcal{T}\right)$-measurable and nonnegative functions with compact support on $\mathcal{T}$. For each process $\theta \in BM\left(\mathcal{T}\right)$, write:

$${(\theta .A)}_t:={\int }_0^t\theta \left(u\right)dA\left(u\right),t\in \mathcal{T},$$

such that $\theta$ is integrable with respect to the return process.

Let ${\left\{{\Lambda }_t\right\}}_{t\in \mathcal{T}}$ denote a $\mathcal{G}$-adapted stochastic process:

$${\Lambda }_t:=\frac{e^{{(\theta .A)}_t}}{\mathcal{M}{\left(\theta \right)}_t},t\in \mathcal{T},$$

where $\mathcal{M}{\left(\theta \right)}_t:=E\left[e^{{(\theta .A)}_t}\right|{\mathcal{F}}_t^A]$ is a Laplace cumulant process and takes the following form:

$$\mathcal{M}{\left(\theta \right)}_t=exp\left[{\int }_0^t{\theta }_s\left({\mu }_s-\frac{1}{2}{\sigma }_s^2\right)ds+\frac{1}{2}{\int }_0^t{\theta }_s^2{\sigma }_s^{2}ds+{\int }_0^t{\int }_{{\mathcal{R}}^{+}}\left(e^{{\theta }_{s}h\left(z\right)}-1-\theta h\left(z\right)\right)\rho \left(dz\right|s)\eta \left(ds\right)\right].$$

Therefore,

$$\begin{array}{cc}\hfill {\Lambda }_t& =exp[{\int }_0^t{\theta }_s{\sigma }_{s}d{W}_s-\frac{1}{2}{\int }_0^t{\theta }_s^2{\sigma }_s^{2}ds+{\int }_0^t{\int }_{{\mathcal{R}}^{+}}{\theta }_{s}h\left(z\right)\tilde{N}(dz,ds)\hfill \\ \hfill & -{\int }_0^t{\int }_{{\mathcal{R}}^{+}}\left(e^{{\theta }_{s}h\left(z\right)}-1+{\theta }_{s}h\left(z\right)\right)\rho \left(dz\right|s)\eta \left(ds\right)].\hfill \end{array}$$

(3)

The goal is to use ${\Lambda }_t$ in (3) as the Radon–Nikodym derivative to change the historical probability measure to the risk-neutral measure. Therefore, (3) is an essential part of our pricing formulation. An important characteristic of risk-neutral measure is that every discounted price process is a martingale under this measure. To establish this key property, it is paramount to show that (3) is ${\mathcal{G}}_t$-martingale.

Lemma 1.

${\Lambda }_t$ is $\mathbb{P}$ martingale w.r.t ${\mathcal{G}}_t$.

See Appendix A for proof.

For each $\theta \in L\left(A\right)$ define a new probability measure ${\mathbb{P}}^{\theta }\sim \mathbb{P}$ on $\mathcal{G}\left(T\right)$ by the Radon–Nikodym derivative:

$$\frac{d{\mathbb{P}}^{\theta }}{d\mathbb{P}}{|}_{\mathcal{G}\left(T\right)}:={\Lambda }_T.$$

(4)

This new measure $d{\mathbb{P}}^{\theta }$ is defined by the Esscher transform ${\Lambda }_T$ associated with $\theta \in L\left(A\right)$.

The local-martingale condition, i.e., there exists an equivalent martingale measure under which discounted asset prices are local-martingales in the absence of arbitrage, is the foundation of asset pricing theory. Below, we state a necessary and sufficient condition for the local martingale condition in our framework.

Proposition 1.

For each $t\in \mathcal{T}$, let the discounted price of the risky asset at time t be:

$$\tilde{S}\left(t\right):=e^{-rt}S\left(t\right).$$

Then the discounted price process $\tilde{S}:=\left\{\tilde{S}\left(t\right)\right|t\in \mathcal{T}\}$ is an ${\mathbb{P}}^{\theta }$-local-martingale if and only if ${\theta }_t:=\langle \theta ,X_t\rangle ,t\in \mathcal{T},$ is such that $\theta :=({\theta }_1,{\theta }_2,\dots ,{\theta }_N)\in {\mathcal{R}}^N$ satisfies the following equation:

$${\theta }_t{\sigma }_t^2+{\int }_{{\mathcal{R}}^{+}}\left\{e^{\theta h\left(z\right)}(e^{h\left(z\right)}-1)-h\left(z\right)\right\}\rho d\left(z\right|t){\eta }^{\prime }\left(t\right)=r_t-{\mu }_t.$$

(5)

See Appendix A for proof.

The results from the Lemma 1, Equation (4), and Proposition 1, allow us to use (3) to drive the risk-neutral dynamics of the return process.

Proposition 2.

Suppose ${\tilde{W}}_t=W_t-{\int }_0^t{\sigma }_s\theta ds$ is a ${\mathbb{P}}^{\theta }$-Brownian motion, ${\rho }^{\theta }\left(dz\right|t):=e^{\theta h\left(z\right)}\rho \left(dz\right|t)$ is the ${\mathbb{P}}^{\theta }$ compensator of $N^{\theta }(dz,dt)$ then:

$$\begin{array}{c}\hfill d{A}_t=\left({\mu }_t+2\theta {\sigma }_s^2-\frac{1}{2}{\sigma }_s^2\right)dt+{\sigma }_{t}d{\tilde{W}}_t+{\int }_{{\mathcal{R}}^{+}}h\left(z\right)\left(1-e^{-\theta h\left(z\right)}\right){\rho }^{\theta }\left(dz\right|t)\eta \left(dt\right)+{\int }_{{\mathcal{R}}^{+}}h\left(z\right){\tilde{N}}^{\theta }(dz,dt).\end{array}$$

(6)

See Appendix B for proof.

Similarly, we can derive the risk-neutral price process of the reference portfolio.

Proposition 3.

The price process of the reference portfolio S under ${\mathbb{P}}^{\theta }$ is:

$$\begin{array}{cc}\hfill d{S}_t& =(r_t-\frac{1}{2}{\sigma }_t^2)dt+{\sigma }_{t}d{\tilde{W}}_t+{\int }_{{\mathcal{R}}^{+}}\left(e^{\theta h\left(z\right)}(e^{h\left(z\right)}-1)-h\left(z\right)\right){\tilde{N}}^{\theta }(dt,dz).\hfill \end{array}$$

(7)

Proof. 

Proof of Proposition 3. Recall $S_t:=exp\left(A_t\right)$. Then the proof can easily follow by applying Ito’s Lemma and the martingale condition (5) to (6). □

We study the hedging of a European option with pay-off function G using the popular delta hedging strategy. The option price is given by:

$$\begin{array}{c}\hfill C(t,S)=E^{\theta }\left[G\left(S_T\right)\right|S_t=S],\end{array}$$

(8)

where we assume $C\in C^{\infty }([0,T)\times \mathcal{R})$. Furthermore, the delta hedging strategy is $H_t:=\frac{\partial C(t,S)}{\partial S}$, which is the most widely-used hedging strategy with a mathematically tractable structure. Detailed discussion about the hedging strategy is provided in the next subsection.

2.3. Continuous Hedging Strategy

We assume the existence of a continuous-time trading strategy H. If continuous-time hedging was possible, agents in the market would like to follow this strategy. This strategy may be chosen in several ways and may not lead to a perfect replication when markets are incomplete. In this study, we do not address the relative advantages of different choices of H. Rather, we suppose its existence from another generalised jump-diffusion process that satisfies the same set of assumptions as S. As such, we assume that:

$$\begin{array}{cc}\hfill d{H}_t& =a_{t}dt+b_{t}d{W}_t+{\int }_{{\mathcal{R}}^{+}}{\gamma }_t\tilde{N}(dt,dz),\hfill \end{array}$$

(9)

under $\mathbb{P}$, where $a_t$, $b_t$, and ${\gamma }_t$ are the parameters of the process that will be determined below.

By applying the Itô Lemma to the definition of $H_t$, we can show the following decomposition:

$$\begin{array}{cc}\hfill d{H}_t:=d\frac{\partial C(t,S)}{\partial S}=& \left\{\frac{{\partial }^{2}C}{\partial t\partial S}(t,S)+\left({\mu }_t+\frac{1}{2}{\sigma }_t^2+{\int }_{{\mathcal{R}}^{+}}h\left(z\right)\rho \left(dz\right|t){\eta }^{\prime }\right)\frac{{\partial }^{2}C}{\partial S^2}(t,S)+\frac{1}{2}{\sigma }_t^2\frac{{\partial }^{3}C}{\partial S^3}(t,S)\right\}dt\hfill \\ & +{\sigma }_t\frac{{\partial }^{2}C}{\partial S^2}(t,S)d{W}_t+{\int }_{{\mathcal{R}}^{+}}\left\{\frac{\partial C}{\partial S}(t,S+e^{h\left(z\right)}-1)-\frac{\partial C}{\partial S}(t,S)\right\}\rho \left(dz\right|t){\eta }^{\prime }dt\hfill \\ & +{\int }_{{\mathcal{R}}^{+}}\left\{\frac{\partial C}{\partial S}(t,S+e^{h\left(z\right)}-1)-\frac{\partial C}{\partial S}(t,S)\right\}\tilde{N}(dt,dz).\hfill \end{array}$$

Then,

$$\begin{array}{cc}\hfill a_t=& \frac{{\partial }^{2}C}{\partial t\partial S}(t,S)+({\mu }_t+\frac{1}{2}{\sigma }_t^2+\frac{1}{2}{\sigma }_t^2\frac{{\partial }^{3}C}{\partial S^3}(t,S)\hfill \\ & +{\int }_{{\mathcal{R}}^{+}}\left\{\frac{\partial C}{\partial S}(t,S+e^{h\left(z\right)}-1)-\frac{\partial C}{\partial S}(t,S)+h\left(z\right)\frac{{\partial }^{2}C}{\partial S^2}(t,S)\right\}\rho \left(dz\right|t){\eta }^{\prime }dt\hfill \\ \hfill b_t=& {\sigma }_t\frac{{\partial }^{2}C}{\partial S^2}(t,S)\hfill \\ \hfill {\gamma }_t=& \frac{\partial C}{\partial S}(t,S+e^{h\left(z\right)}-1)-\frac{\partial C}{\partial S}(t,S).\hfill \end{array}$$

2.4. Asymptotic Hedging Error

In what follows we drive the asymptotic distribution of the hedging error, generalising the methodology proposed in Hayashi and Mykland (2005) and Tankov and Voltchkova (2009).

Let,

$$\begin{array}{ccccccc}\hfill {\mu }_t^n:=& (-n)\vee {\mu }_t\wedge n;\hfill & \hfill {\sigma }_t^n:=& (-n)\vee {\sigma }_t\wedge n;\hfill & \hfill h^n\left(z\right):=& -\sqrt{n\rho \left(n\right)}\vee h\left(z\right)\wedge \sqrt{n\rho \left(n\right)};\hfill & \\ \hfill a_t^n:=& (-n)\vee a_t\wedge n;\hfill & \hfill b_t^n:=& (-n)\vee b_t\wedge n;\hfill & \hfill {\gamma }_t^n\left(z\right):=& -\sqrt{n\rho \left(n\right)}\vee {\gamma }_t\left(z\right)\wedge \sqrt{n\rho \left(n\right)}.\hfill \end{array}$$

Then the processes $S^n$ and $H^n$ are Lévy–Itô processes with bounded coefficients and bounded jumps that coincide with S and H on the set:

$$\begin{array}{cc}\hfill {\Omega }_n:=& \{{\displaystyle \underset{0\le t\le T}{sup}}max\left(\right|{\mu }_t|,|{\sigma }_t|,|{ϝ}_t\left|\right)\le n;{\displaystyle \underset{0\le t\le T}{sup}}max\left(\right|a_t|,|b_t|,|{ϝ}_t\left|\right)\le n;N([0,T],((-\infty ,-n)\cup (n,\infty ))=0)\}.\hfill \end{array}$$

Since all processes are supported cádlág, $P\left[{\Omega }_n\right]\to 1$.

The continuous re-balancing of a portfolio is practically unfeasible. Typically, holders of a position in an option $\Delta$-hedge in discrete time intervals of $t_i=iT/n$. Therefore, the trading strategy is piecewise constant and given by $F_{{\varphi }_n\left(t\right)}$, where ${\varphi }_n\left(t\right)=sup\{t_i,t_i<t\}$. The value of the hedging portfolio at time t is $V_0+{\int }_0^t{H}_{s^{-}}d{S}_s$ with continuous hedging and $V_0+{\int }_0^t{H}_{{\varphi }_n\left(t\right)}d{S}_s$ with discrete hedging. Then the asymptotic distribution of the difference between discrete and continuous hedging is:

$$\begin{array}{c}\hfill U_t^n={\int }_0^t(H_s^{-}-H_{{\varphi }_n\left(t\right)})d{S}_s={\int }_0^t{H}_{s^{-}}^{n}d{S}_s\end{array}$$

(10)

where $n\to \infty$. For any process A we set $A_t^n:=A_t-A_{{\varphi }_n\left(t\right)}$. Under the above conditions, Hayashi and Mykland (2005) provides a thorough discussion on the stable convergence of the bounded processes to their respective original processes.

Furthermore, define the renormalised hedging error process by:

$$\begin{array}{c}\hfill Z_t^n=\sqrt{n}{U}_t^n=\sqrt{n}{\int }_0^t{H}_{s^{-}}^{n}d{S}_s.\end{array}$$

(11)

Let ${\left\{{\widehat{W}}_t\right\}}_{t\in \mathcal{T}}$ be a standard Brownian motion independent of W and N, and let ${\left({\xi }_k\right)}_{k\ge 1}$ and ${\left({\xi }_k^{\prime }\right)}_{k\ge 1}$ be two sequences of a standard normal random variable and ${\left({\zeta }_k\right)}_{k\ge 1}$ sequence of independent uniform random variables on [0,1], such that the three sequences are independent from each other and other random elements. Let ${\left(T_i\right)}_{i\ge 1}$ be an enumeration of the jump times of N, and define:

$$\Delta F_{T_i}=\frac{\partial C}{\partial S}(T_i,S_{T_i})-\frac{\partial C}{\partial S}(T_i,S_{T_i^{-}}),{\tilde{\sigma }}_{T_i^{-}}={\sigma }_{T_i}\frac{{\partial }^{2}C}{\partial S^{T_i}}(T_i,S_{T_i}),{\tilde{\sigma }}_s={\sigma }_s\frac{{\partial }^{2}C}{\partial S^2}(s,S_s).$$

Then, applying Theorem 1 in Tankov and Voltchkova (2009) to the renormalised discrete delta hedging error $Z_t^n$ gives the following asymptotical convergence result in finite-dimensional laws:

$$\begin{array}{cc}\hfill Z_t=& \sqrt{\frac{T}{2}}{\int }_0^t{\sigma }_s{\tilde{\sigma }}_{s}d{\widehat{W}}_s+\sqrt{T}{\displaystyle \sum _{i:T_i\le t:\Delta S_t\ne 0}}\Delta F_{T_i}\sqrt{{\zeta }_i}{\xi }_i{\sigma }_{T_i}+\sqrt{T}{\displaystyle \sum _{i:T_i\le t:\Delta S_t\ne 0}}\Delta S_{T_i}\sqrt{1-{\zeta }_i}{\xi }_i^{\prime }{\tilde{\sigma }}_{T_i^{-}}\hfill \\ \hfill \equiv & \sqrt{\frac{T}{2}}{\int }_0^t{\sigma }_s^2\frac{{\partial }^{2}C}{\partial S^2}d{\widehat{W}}_s+\sqrt{T}{\displaystyle \sum _{i:T_i\le t:\Delta S_t\ne 0}}\left(\frac{\partial C}{\partial S}(T_i,S_{T_i})-\frac{\partial C}{\partial S}(T_i,S_{T_i^{-}})\right)\sqrt{{\zeta }_i}{\xi }_i{\sigma }_{T_i}+\hfill \\ \hfill & \sqrt{T}{\displaystyle \sum _{i:T_i\le t:\Delta S_t\ne 0}}\Delta S_{T_i}\sqrt{1-{\zeta }_i}{\xi }_i^{\prime }{\sigma }_{T_i}\frac{{\partial }^{2}C}{\partial S^{T_i}}(T_i,S_{T_i}).\hfill \end{array}$$

(12)

3. Estimation of the Density of the Hedging Error

Let $Z_T$ be a realisation of the hedging error in (12) at date T (say the maturity date of the option). Suppose that we have a random sample $\left\{Z\left(j\right):j=1,\dots ,n\right\}$ of n $i.i.d.$ observations from $Z_T$, each with pdf $f_{{}_Z}$ and cdf $F_{{}_Z}$ on a support $\Omega \subseteq \mathbb{R}.$ Note from (12) that $Z_T$ always has a continuous density almost everywhere (a.e). In this section, the VaR and the ES associated to $Z_T$ are of inferential interest. The VaR associated to $Z_T$ at level $\alpha$ is defined as:

$$\begin{array}{ccc}\hfill {\mathrm{VaR}}_{\alpha }\left(Z_T\right)& =& inf\left\{z\in \mathbb{R}:F_{{}_Z}\left(z\right)\ge \alpha \right\}=F_{{}_Z}^{-1}\left(\alpha \right),\hfill \end{array}$$

(13)

while the expected shortfall is given as:

$$\begin{array}{ccc}\hfill {\mathrm{ES}}_{\alpha }\left(Z_T\right)& =& \frac{1}{\alpha }{\int }_0^{\alpha }{\mathrm{VaR}}_{\gamma }\left(Z_T\right)d\gamma \hfill \\ & =& -\frac{1}{\alpha }\left(E\left[Z_T{1}_{\left\{Z_T\le z_{\alpha }\right\}}\right]+z_{\alpha }[\alpha -F_{{}_Z}\left(z_{\alpha }\right)]\right),\hfill \end{array}$$

(14)

where $z_{\alpha }=inf\left\{z\in \mathbb{R}:F_{{}_Z}\left(z\right)\ge \alpha \right\}$ is the lowest $\alpha$-quantile and $1_A\left(z\right)=1$ if $x\in A$ and $1_A\left(z\right)=0$ else. The dual representation of (14) is given by:

$$\begin{array}{ccc}\hfill {\mathrm{ES}}_{\alpha }\left(Z_T\right)& =& {\displaystyle \underset{Q\in {\mathcal{Q}}_{\alpha }}{inf}}{E}^Q\left[Z_T\right],\hfill \end{array}$$

(15)

where Q is absolutely continuous such that $\frac{dQ}{dP}\le {\alpha }^{-1}$. If the density of $Z_T$ is continuous, then the expected shortfall is equivalent to the tail conditional expectation defined by ${\mathrm{TCE}}_{\alpha }\left(Z\right)=E[-Z_T|Z_T\le -{\mathrm{VaR}}_{\alpha }\left(Z_T\right)]$ . If $F_{{}_Z}$ (or $f_{{}_Z}$) were given, then the computation of ${\mathrm{VaR}}_{\alpha }\left(Z_T\right)$ and ${\mathrm{ES}}_{\alpha }\left(Z_T\right)$ is straightforward from (13)–(15). This is unfortunately not the case, and one has to resort to estimation techniques to approximate them. Suppose that we have a consistent estimate of $f_{{}_Z}$, say ${\widehat{f}}_{{}_Z}$. So, we can also estimate its cdf and compute the estimate of the value-at-risk and the expected shortfall as:

$$\begin{array}{ccc}\hfill {\widehat{\mathrm{VaR}}}_{\alpha }\left(Z_T\right)& =& {\widehat{F}}_{{}_Z}^{-1}\left(\alpha \right),{\widehat{\mathrm{ES}}}_{\alpha }\left(Z_T\right)=\frac{1}{\alpha }{\int }_0^{\alpha }{\widehat{\mathrm{VaR}}}_{\gamma }\left(Z_T\right)d\gamma .\hfill \end{array}$$

(16)

Our main goal is to find an estimator ${\widehat{f}}_{{}_Z}$ of $f_{{}_Z}$ that captures the maximum uncertainty in $Z_T$. To achieve this goal, we will use the information entropy approach.

3.1. Information Entropy and Density Estimation

The information entropy associated with $Z_T$ is defined as:

$$\begin{array}{ccc}\hfill I_E\left(Z_T\right)& =& -{\int }_{\Omega }{f}_{{}_Z}\left(z\right)ln{f}_{{}_Z}\left(z\right)dz,\hfill \end{array}$$

(17)

where, by convention, we assume that $0\times ln\left(0\right)=0$. $I_E\left(Z_T\right)$ is a measure of the information carried by $Z_T$. As data are communicated more, they are corrupted with more noise so that the entropy increases, therefore they carry less information.

Let $g_{{}_Z}^{\left(k\right)}:\Omega \to \mathbb{R}$ ($k=0,1,2,\dots$) be a moment function of $Z_T.$ Then, the moment of $Z_T$ with respect to $g_{{}_Z}^{\left(k\right)}$ is defined as:

$$\begin{array}{ccc}\hfill {\mu }_k& =& \mathbb{E}\left[g_{{}_Z}^{\left(k\right)}\left(z\right)\right]={\int }_{\Omega }{g}_{{}_Z}^{\left(k\right)}\left(z\right)f_{{}_Z}\left(z\right)dx,k=0,1,2,\dots .\hfill \end{array}$$

(18)

In practice, polynomial functions are often used for the moment function (for example, see Zellner and Highfield (1988)). $g_{{}_Z}^{\left(k\right)}$ ($k=0,1,2,\dots$), i.e., $g_{{}_Z}^{\left(k\right)}\left(z\right)=z^k,k=0,1,2,\dots ,$ with the normalisation $g_{{}_Z}^{\left(0\right)}\left(z\right)=1.$ In this case, we have:

$$\begin{array}{ccc}\hfill {\mu }_0& =& {\int }_{\Omega }{f}_{{}_Z}\left(z\right)dz=1,{\mu }_k={\int }_{\Omega }{z}^k{f}_{{}_Z}\left(z\right)dz,k=1,2,\dots .\hfill \end{array}$$

(19)

For the remainder of the paper, we assume the sequence ${\left({\mu }_k\right)}_{k\ge 1}$ satisfies the Carleman’s condition [see Akhiezer (1965)], i.e.,

$$\begin{array}{c}\hfill {\displaystyle \sum _{k=1}^{\infty }}{\mu }_{2k}^{-\frac{1}{2k}}=+\infty .\end{array}$$

(20)

Note that when $\Omega =\mathbb{R}$ and $g_{{}_Z}^{\left(k\right)}\left(z\right)=z^k$ for all $k=0,1,2,\dots ,$ condition (20) is sufficient for the determinacy of the Hamburger moment problem.

Due to issues associated with using many moments, we only use $k=0,\dots ,m$ moments in (19), and our goal is to find the density function $f_{{}_Z}\left(z\right)$ such that:

$$\begin{array}{ccc}& & {\displaystyle \underset{f_{{}_Z}\in {\mathbf{L}}^1(\Omega ):f_{{}_Z}\ge 0}{max}}\left\{I_E\left(Z_T\right)=-{\int }_{\Omega }{f}_{{}_Z}\left(z\right)ln{f}_{{}_Z}\left(z\right)dz\right\}\hfill \\ & & \hfill \\ & & \mathrm{subject}\mathrm{to}\hfill \end{array}$$

(21)

$$\begin{array}{ccc}& & {\overline{\mu }}_k={\int }_{\Omega }{g}_{{}_Z}^{\left(k\right)}\left(z\right)f_{{}_Z}\left(z\right)dz,k=0,1,2,\dots ,m\hfill \end{array}$$

(22)

where ${\overline{\mu }}_k=n^{-1}{\sum }_{j=1}^n{g}_{{}_Z}^{\left(k\right)}\left(Z\left(j\right)\right)$ is the empirical moment of $Z_T$ with respect to $g_{{}_Z}^{\left(k\right)}.$ For now, we assume that the number m of selected moments of the pdf in (19) that have been chosen to match the empirical moments is fixed. Later on, we discuss how m can be selected in a data-dependent manner based on a Bayesian- or Schwartz-type information criterion. Since $g_{{}_Z}^{\left(0\right)}=1,$ the zeroth moment implies that ${\mu }_0={\int }_{\Omega }{f}_{{}_Z}\left(z\right)dz={\overline{\mu }}_0=1$ in (22). For example if $g_{{}_Z}^{\left(k\right)}\left(z\right)=z^k,k=0,1,2,3,4$ with $g_{{}_Z}^{\left(0\right)}\left(z\right)=1,$ we have $m=4$ and problem (22) becomes:

$$\begin{array}{ccc}& & {\displaystyle \underset{f_{{}_Z}\in {\mathbf{L}}^1(\Omega ):f_{{}_Z}\ge 0}{max}}{I}_E\left(Z_T\right)\mathrm{subject}\mathrm{to}\hfill \\ & & \hfill \\ & & {\overline{\mu }}_0={\int }_{\Omega }{f}_{{}_Z}\left(z\right)dz=1,{\overline{\mu }}_1={\int }_{\Omega }z{f}_{{}_Z}\left(z\right)dz,{\overline{\mu }}_2={\int }_{\Omega }{z}^2{f}_{{}_Z}\left(z\right)dz\hfill \\ & & {\overline{\mu }}_3={\int }_{\Omega }{z}^3{f}_{{}_Z}\left(z\right)dz,{\overline{\mu }}_4={\int }_{\Omega }{z}^4{f}_{{}_Z}\left(z\right)dz.\hfill \end{array}$$

(23)

For the remainder of the paper, we call the solution ${\widehat{f}}_{{}_Z}\left(z\right)$ of problem (22) the maximum-entropy estimator of $f_{{}_Z}\left(z\right).$ We will now establish that such a solution exists and is unique. Let,

$$\begin{array}{ccc}\hfill L(f_{{}_Z},{\lambda }_0,\dots ,{\lambda }_m)& =& -I_E\left(Z\right)+{\displaystyle \sum _{k=0}^m}{\lambda }_k\left({\int }_{\Omega }{g}_{{}_Z}^{\left(k\right)}\left(z\right)f_{{}_Z}\left(z\right)dz-{\overline{\mu }}_k\right)\hfill \end{array}$$

(24)

denote the Lagrange function associated with (22), where ${\lambda }_0,\dots ,{\lambda }_m$ are the Lagrange multipliers. By noting that solving (22) is equivalent to static problem:

$$\begin{array}{c}\hfill {\displaystyle \underset{f_{{}_Z}\in {\mathbf{L}}^1(\Omega ):f_{{}_Z}\ge 0,{\lambda }_0,\dots ,{\lambda }_m}{min}}L(f_{{}_Z},{\lambda }_0,\dots ,{\lambda }_m),\end{array}$$

(25)

it is sufficient to show that (25) has a unique solution with respect to $f_{{}_Z}\in {\mathbf{L}}^1(\Omega ):f_{{}_Z}\ge 0.$ More often, the minimisation over $f_{{}_Z}$ is carried out via the first variation of $L(f_{{}_Z},{\lambda }_0,\dots ,{\lambda }_m)$ with respect to $f_{{}_Z}.$ This type of calculation is misleading because $L(f_{{}_Z},{\lambda }_0,\dots ,{\lambda }_m)$ has support on the set $\left\{f_{{}_Z}\in {\mathbf{L}}^1(\Omega ):f_{{}_Z}\ge 0\mathrm{a}.\mathrm{e}.\right\},$ and the complement of this set is dense in ${\mathbf{L}}^1$, meaning that not only is $L(f_{{}_Z},{\lambda }_0,\dots ,{\lambda }_m)$ nowhere differentiable in ${\mathbf{L}}^1$, it is also nowhere continuous. Due to these reasons, we follow the convex duality approach in Borwein et al. (2003). From (Borwein et al. 2003, Theorems 1–2), there is a unique solution $f_{{}_Z}$ of (25) satisfying:

$$\begin{array}{ccc}\hfill f_{{}_Z}(z;\lambda )& =& exp\left(-1-{\lambda }_0-{\displaystyle \sum _{k=1}^m}{\lambda }_k{g}_{{}_Z}^{\left(k\right)}\left(z\right)\right)\hfill \\ & =& \frac{1}{\theta }exp\left(-{\displaystyle \sum _{k=1}^m}{\lambda }_k{g}_{{}_Z}^{\left(k\right)}\left(z\right)\right)\hfill \end{array}$$

(26)

where $\lambda ={({\lambda }_0,{\lambda }_1,\dots ,{\lambda }_m)}^{\prime },$$\theta =exp(1+{\lambda }_0),$ and the zeroth-moment equality implies that:

$$\theta ={\int }_{\Omega }exp\left(-{\displaystyle \sum _{k=1}^m}{\lambda }_k{g}_{{}_Z}^{\left(k\right)}\left(z\right)\right)dz.$$

To complete the closed form of $f_{{}_Z}(z;\lambda )$, we must substitute $\lambda$ in (26) by an optimal value $\widehat{\lambda }$. Zellner and Highfield (1988) consider the case in which $g_Z^{\left(k\right)}\left(z\right)=z^k$ where $k=1,\dots ,m=4,$ and use an algorithm based on a Newton method to compute ${\lambda }_1,\dots ,{\lambda }_4$ from the restrictions in (22). This numerical approximation is cumbersome even for a moderate choice of $m=4.$ In this paper, we propose the maximum likelihood (ML) method to estimate ${\lambda }_1,\dots ,{\lambda }_m,$ and $\theta .$

Since $\left\{Z\left(j\right):j=1,\dots ,n\right\}$ are $i.i.d.$ with common pdf $\frac{1}{\theta }exp\left(-{\sum }_{k=1}^m{\lambda }_k{g}_{{}_Z}^{\left(k\right)}\left(z\left(j\right)\right)\right),$ the likelihood function of the sample can be written as:

$$\begin{array}{ccc}\hfill \mathcal{L}\left(Z\right(1),\dots ,Z(n);\lambda )& =& {\displaystyle \prod _{j=1}^n}\frac{1}{\theta }exp\left(-{\displaystyle \sum _{k=1}^m}{\lambda }_k{g}_{{}_Z}^{\left(k\right)}\left(z_j\right)\right)\hfill \\ & =& \frac{1}{{\theta }^n}exp\left[-{\displaystyle \sum _{j=1}^n}{\displaystyle \sum _{k=1}^m}{\lambda }_k{g}_{{}_Z}^{\left(k\right)}\left(z_j\right)\right].\hfill \end{array}$$

(27)

The log-likelihood function from (27) is then given by $\ell (Z\left(1\right),\dots ,Z\left(n\right);\lambda )=-nln\left(\theta \right)-{\sum }_{j=1}^n{\sum }_{k=1}^m{\lambda }_k{g}_{{}_Z}^{\left(k\right)}\left(z_j\right).$ The ML estimator of $\lambda ={({\lambda }_0,{\lambda }_1,\dots ,{\lambda }_m)}^{\prime }$ satisfies:

$$\begin{array}{c}\hfill {\displaystyle \underset{\lambda }{max}}\ell (Z\left(1\right),\dots ,Z\left(n\right);\lambda )\end{array}$$

(28)

and we can prove the following on the existence of unique solutions for both problems (28) and (22).

Proposition 4.

Suppose that ${\overline{\mu }}_k=n^{-1}{\sum }_{j=1}^n{g}_{{}_Z}^{\left(k\right)}\left(Z\left(j\right)\right)<\infty \forall k=1,\dots ,m.$ Then:

(a)

Problem (28) has a unique solution with respect to λ;

(b)

Problem (22) has a unique solution with respect to f_Z (·).

Proof. 

First, observe first that $\ell \left(Z\right(1),\dots ,Z(n);\lambda )$ is a concave function in $\lambda ,$ and twice differentiable with respect to $\lambda .$ The first order condition of (28) is given by:

$$\begin{array}{ccc}\hfill \frac{\partial \ell \left(Z\right(1),\dots ,Z(n);\lambda )}{\partial {\lambda }_k}& =& -n\frac{\partial ln\left(\theta \right)}{\partial {\lambda }_k}-{\displaystyle \sum _{j=1}^n}{g}_{{}_Z}^{\left(k\right)}\left(z_j\right)=0\hfill \\ & \iff & \frac{\partial ln\left(\theta \right)}{\partial {\lambda }_k}=-n^{-1}{\displaystyle \sum _{j=1}^n}{g}_{{}_Z}^{\left(k\right)}\left(z_j\right)=-{\overline{\mu }}_k,k=1,\dots ,m.\hfill \end{array}$$

(29)

Moreover, $\theta =exp(1+{\lambda }_0)={\int }_{\Omega }exp\left(-{\sum }_{k=1}^m{\lambda }_k{g}_{{}_Z}^{\left(k\right)}\left(z_j\right)\right)d{z}_j$ from the zeroth-moment equality. Thus,

$$\begin{array}{ccc}\hfill \frac{\partial ln\left(\theta \right)}{\partial {\lambda }_k}& =& \frac{1}{\theta }\frac{\partial }{\partial {\lambda }_k}{\int }_{\Omega }exp\left(-{\displaystyle \sum _{k=1}^m}{\lambda }_k{g}_{{}_Z}^{\left(k\right)}\left(z_j\right)\right)d{z}_j=-{\int }_{\Omega }{g}_{{}_Z}^{\left(k\right)}\left(z_j\right)\frac{1}{\theta }exp\left(-{\displaystyle \sum _{k=1}^m}{\lambda }_k{g}_{{}_Z}^{\left(k\right)}\left(z_j\right)\right)d{z}_j\hfill \\ & =& -{\int }_{\Omega }{g}_{{}_Z}^{\left(k\right)}\left(z_j\right)f_{{}_Z}\left(z_j\right)d{z}_j=-{\mu }_k\equiv {\mu }_k\left(\lambda \right),\hfill \end{array}$$

(30)

where ${\mu }_k\left(\lambda \right)$ is a continuous function of $\lambda .$ Therefore, (29) and (30) entail that:

$$\begin{array}{ccc}\hfill {\mu }_k\left(\lambda \right)& =& {\overline{\mu }}_k,k=1,\dots ,m.\hfill \end{array}$$

(31)

So, as long as ${\overline{\mu }}_k<\infty ,$ (31) has a solution with respect to $\lambda ,$ i.e., the likelihood function $\ell \left(Z\right(1),\dots ,Z(n);\lambda )$ has a critical point $\widehat{\lambda }$.

We will now show that this critical point, $\widehat{\lambda },$ is the unique point where $\ell \left(Z\right(1),\dots ,Z(n);\lambda )$ is maximised. Let $H\left(Z\right(1),\dots ,Z(n);\lambda )$ denote the hessian matrix of the log-likelihood function $\ell \left(Z\right(1),\dots ,Z(n);\lambda )$. The $(k,p)$ entree of $H\left(Z\right(1),\dots ,Z(n);\lambda )$ is:

$$\begin{array}{ccc}& & \frac{\partial }{\partial {\lambda }_p}\left(\frac{\partial \ell \left(Z\right(1),\dots ,Z(n);\lambda )}{\partial {\lambda }_k}\right)=\frac{\partial }{\partial {\lambda }_p}\left(-n\frac{\partial ln\left(\theta \right)}{\partial {\lambda }_k}-{\displaystyle \sum _{j=1}^n}{g}_{{}_Z}^{\left(k\right)}\left(z_j\right)\right)=-n\frac{\partial }{\partial {\lambda }_p}\frac{\partial ln\left(\theta \right)}{\partial {\lambda }_k}\hfill \\ & =& n\frac{\partial }{\partial {\lambda }_p}\left(\frac{1}{\theta }{\int }_{\Omega }{g}_{{}_Z}^{\left(k\right)}\left(z_j\right)exp\left(-{\displaystyle \sum _{k=1}^m}{\lambda }_k{g}_{{}_Z}^{\left(k\right)}\left(z_j\right)\right)d{z}_j\right)\hfill \\ & =& -n\frac{1}{\theta }{\int }_{\Omega }{g}_{{}_Z}^{\left(p\right)}\left(z_j\right)g_{{}_Z}^{\left(k\right)}\left(z_j\right)exp\left(-{\displaystyle \sum _{k=1}^m}{\lambda }_k{g}_{{}_Z}^{\left(k\right)}\left(z_j\right)\right)d{z}_j+\hfill \\ & & n\left(\frac{1}{{\theta }^2}{\int }_{\Omega }{g}_{{}_Z}^{\left(k\right)}\left(z_j\right)exp\left(-{\displaystyle \sum _{k=1}^m}{\lambda }_k{g}_{{}_Z}^{\left(k\right)}\left(z_j\right)\right)d{z}_j\right)\times {\int }_{\Omega }{g}_{{}_Z}^{\left(p\right)}\left(z_j\right)exp\left(-{\displaystyle \sum _{k=1}^m}{\lambda }_k{g}_{{}_Z}^{\left(k\right)}\left(z_j\right)\right)d{z}_j\hfill \\ & =& -n\left[{\int }_{\Omega }{g}_{{}_Z}^{\left(p\right)}\left(z_j\right)g_{{}_Z}^{\left(k\right)}\left(z_j\right)f_Z\left(z_j\right)d{z}_j-{\int }_{\Omega }{g}_{{}_Z}^{\left(p\right)}\left(z_j\right)f_Z\left(z_j\right)d{z}_j{\int }_{\Omega }{g}_{{}_Z}^{\left(k\right)}\left(z_j\right)f_Z\left(z_j\right)d{z}_j\right]\hfill \\ & =& -n\left[\mathbb{E}\left(g_{{}_Z}^{\left(p\right)}\left(z_j\right)g_{{}_Z}^{\left(k\right)}\left(z_j\right)\right)-\mathbb{E}\left(g_{{}_Z}^{\left(p\right)}\left(z_j\right)\right)\mathbb{E}\left(g_{{}_Z}^{\left(k\right)}\left(z_j\right)\right)\right]\hfill \\ & =& -nCoV\left(g_{{}_Z}^{\left(p\right)}\left(z_j\right),g_{{}_Z}^{\left(k\right)}\left(z_j\right)\right),\forall p,k=1\dots ,m,\hfill \end{array}$$

(32)

where $CoV\left(a,b\right)$ is the covariance of the two random variables a and b. From (32), it is clear that $H\left(Z\right(1),\dots ,Z(n);\lambda )$ is a symmetric and strictly negative definite for all values of the vector of Lagrange multipliers. Thus, the critical point $\widehat{\lambda }$ is a unique maximum, which completes the proof of Proposition 4(a). It follows immediately from (26) that ${\widehat{f}}_{{}_Z}(z;\widehat{\lambda })$ is also a unique maximum, thus establishing Proposition 4(b). □

3.2. Choice of m and the Moment Functions in the Density Estimation

To avoid the problems associated with using many moments, we suggest using a finite number m of moments in the estimation of the density function $f_{{}_Z}\left(z\right).$ In practice, the choice of $m,$ as well as that of the moment functions $g_{{}_Z}^{\left(k\right)}(·),k=1,\dots ,m,$ may not be obvious. In this section, we briefly discuss how both m and $g_{{}_Z}^{\left(k\right)}(·),k=1,\dots ,m$ can be approximated in a data dependent manner. First, note from (26) that the density function $f_{{}_Z}(z;\lambda )$ satisfies:

$$\begin{array}{ccc}\hfill f_{{}_Z}(z;\lambda )& =& \frac{1}{\theta }exp\left(-{\displaystyle \sum _{k=1}^m}{\lambda }_k{g}_{{}_Z}^{\left(k\right)}\left(z\right)\right).\hfill \end{array}$$

(33)

For the choice of moment functions, we replace each $g_{{}_Z}^{\left(k\right)}\left(z\right)$ in (33) with its truncated Taylor series expansion of around the expectation of z. From this, the problem of looking for the proper $g_{{}_Z}^{\left(k\right)}\left(z\right),k=0,1,\dots ,m,$ moment functions is the same as finding an optimal order of truncation for each of the expansion:

$$\begin{array}{ccc}\hfill \widehat{\lambda }& =& arg{\displaystyle \underset{\lambda }{max}}{\displaystyle \prod _{j=1}^n}\frac{1}{\theta }exp\left(-{\displaystyle \sum _{k=1}^m}{\lambda }_k[a_{0k}+a_{1k}{z}_j+a_{2k}{z}_j^2+a_{3k}{z}_j^3+\dots ]\right)\hfill \\ & =& arg{\displaystyle \underset{\lambda }{max}}{\displaystyle \prod _{j=1}^n}\frac{1}{\theta }exp\left({\beta }_0+{\beta }_1{z}_j+{\beta }_2{z}_j^2+{\beta }_3{z}_j^3+\dots \right),\hfill \end{array}$$

(34)

where ${\beta }_0,{\beta }_1,{\beta }_1,{\beta }_3,\dots$ are unknown coefficients to be estimated. We are interested in finding m and the truncation order $l_0$ of the power series in (34). So, the search for m and the moment functions is converted to a search for an optimal truncation orders m and $l_0$, that yields the best fit of ${\widehat{f}}_{{}_Z}(z;\widehat{\lambda })$ to the data. We suggest using Bayesian information criterion (BIC) or Schwartz information criterion (SIC) based on (34) to select the optimal m and $l_0.$ Let ${\widehat{\lambda }}_{MLE}$ denote the ML estimator of $\lambda$ using the optimal order of truncation $l_0.$ Then, the estimated pdf of $f_{{}_Z}\left(z\right)$ is ${\widehat{f}}_{{}_Z}(z;{\widehat{\lambda }}_{MLE}),$ which can then be used to compute ${\widehat{\mathrm{VaR}}}_{\alpha }\left(Z\right)$ and ${\widehat{\mathrm{ES}}}_{\alpha }\left(Z\right)$ in (16).

4. Numerical Analysis

In this section we conduct a numerical experiment to analyse the sensitivity of the hedging error with respect to model parameters. In the previous sections we have defined a general jump-diffusion process with the jump component specified by a kernel biased completely random measure. This generalised framework nests a number of very important models in mathematical finance, including, but not limited to, the jump diffusion model of Merton (1976), the generalised gamma process discussed in Lo and Weng (1989), the variance gamma process by Madan et al. (1998), and the CGMY model of Carr et al. (2002). Here, for simplicity, we only use the generalised gamma (GG) process. The analysis can be easily extended to other classes of models or even their Markovian regime switching versions discussed in Fard and Siu (2013).

The GG process generalises several famous models in finance. For example, the inverse gamma (IG) and weighted gamma (WG) processes are special forms of the GG process. Let $\alpha \le 1$ denote a constant shape parameter and $\delta \left(t\right)$ be the time-dependent scale parameter of the GG process. We can express the intensity process of the GG process as:

$$\rho \left(dz\right|t)\eta \left(dt\right)=\frac{1}{\Gamma \left(\alpha \right)\alpha e^{\delta \left(t\right)z}{z}^{\alpha }}dz\eta \left(dt\right),$$

where $\Gamma (·)$ is the gamma function. We note that this process simplifies to a WG process when $\alpha =1$ and to a IG process when $\alpha =\frac{1}{2}$.

Several observations are of order. First, to obtain the GG process, we must set the kernel function $h\left(z\right)$ to $c{z}^q$, for some constants c and q, and choose a particular parametric form of the compensator measure. Second, different versions of GG processes are obtained for different sets of the kernel function parameters values. For example, if $c>0$ and $q=1$ (linear kernel function), then we obtained the scale-distorted version of the GG process. If $c>1$, then the jump sizes are overstated, while they are understated if $0<c<1$. When $q>0$ and $c=1,$ we obtain the power-distorted version of the GG process. When $q>1$, small jump sizes (i.e., $0<z<1$) are understated and large jump sizes (i.e., $z>1$) are overstated. Finally, when $0<q<1$, small jump sizes are overstated and large jump sizes are understated.

To simulate the GG process, we use the Poisson weighted algorithm by Lee and Kim (2004). Throughout the simulations, we set at $\alpha =1/2.$ It is worth mentioning that the Poisson weighted algorithm applies to a wide class of completely random measures, which are very difficult to simulate directly, see Lee and Kim (2004).

To implement the algorithm, divide the time interval to maturity $[0,T]$ into $nT$ equally spaced subintervals. Then for each $j=0,1,...,n-1$, let $[t_j,t_{j+1}]$ be the $(j+1)$st subinterval. Let M denote the number of jumps of the completely random measure over the term to maturity, such that, M controls the accuracy of the approximation of the algorithm. To implement the Poisson weighted algorithm, we take the following steps:

Step 1.

Let ℵ denote a normalised density function defined by ${\aleph }^{-1}:=\frac{1}{{\eta }^{\prime }\left(t\right)}{\int }_0^T{\eta }^{\prime }\left(s\right)ds$;

Step 2.

Generate a set of $i.i.d$ random variables $T_i\ge 0$ for each $i=1,2...,MT$, from ℵ;

Step 3.

Generate the jump size ${\Omega }_i$ from the conditional density function $g_{T_i}$ (i.e., gamma);

Step 4.

Evaluate whether $T_i\in [t_j,t_{j+1})$. If yes, calculate:

$$\begin{array}{c}\hfill {\lambda }_i=\frac{{\rho }^{\theta }\left({\Omega }_i\right|T_i)}{MT{g}_{t_i}\left({\Omega }_i\right)}{\int }_0^T{\eta }^{\prime }\left(s\right)ds.\end{array}$$

${\lambda }_i$ is the intensity of the Poisson distribution, used to generate the Poisson weights $\mathcal{W}$.

We provide a comparative analysis of the density function of the hedging error obtained by our ME estimation of the hedging error with that obtained by Monte Carlo simulation. In Monte Carlo simulation, we calculate the hedging error at $t_n$, $\Delta {\Pi }_n$, $n=1,\dots ,N$, with hedging portfolio ${\Pi }_n:=\Pi (t_n,S)$ defined by:

$$\Pi (t_n,s)=S\left(t_n\right)H(t_n,S)-C(t_n,S)$$

where the option price and delta are computed by evaluating the risk neutral expectation in Equation (8).

We consider hedging a one-month at the money call option with $S\left(0\right)=K=1$. We analyse different hedging frequencies, namely, $N=4$ (weekly hedging), $N=20$ (daily hedging), $N=80$ (hedging every hour and a half), and $N=320$ (hedging every 20 min). For brevity, we suppose that the WG process with $c=q=1$, equal model parameters under $\mathbb{P}$ and ${\mathbb{P}}^{\theta }$, zero interest rate, and zero drift for the underlying asset. Furthermore, we assume constant volatility, $\sigma =18\%$. Then, the set of rebalancing times is $\left\{t_n\right\},t_n=n\delta t,n=1\dots ,N$ with $\delta =0.1.$ All computations are fully vectorised in ${\mathrm{JULIA}}^{\mathrm{TM}}$.

In this setting, the kernel-biased completely random measure on $\mathcal{T}$ is $\kappa \left(t\right)=$${\sum }_{i=1}^{MT}h\left({\Omega }_i{\mathcal{W}}_i\right)1_{(T_i\le t)}$. It can be shown that the process generated by the Poisson weighting algorithm converges in distribution to $\kappa \left(t\right)$ defined on $[0,T]$ with the Skorohod topology as $M\to \infty$ (see Lee and Kim 2004). Note that the true values of the truncation lag (m) and truncation order ($l_0$) in the entropy density are set at $m=l_0=4.$ The optimal choice of these parameters with the BIC criterion varies within the pseudo-samples, but the minimum and maximum choice of each parameter are 2 and 5, respectively.

Figure 1 compares the distribution of hedging error obtained from the Monte Carlo Simulation and estimated from the ME method. In

Table 1. Summary statistics: Moments (Mj for j = 1, 2, 3, 4) of the distribution of hedging error under the GG (generalised gamma) process. VaR: Value-at-risk.

	N	M1 (%)	M2 (%)	M3 (%)	M4 (%)	99% VaR	95% VaR
Max Entropy	4	0.173	43.21	−1.97	7.01	99.87	74.38
Max Entropy	20	0.154	28.44	−2.41	16.01	88.31	72.11
Max Entropy	80	−0.058	26.04	−3.36	19.74	77.59	51.6
Max Entropy	320	−0.041	22.14	−3.99	14.32	61.33	46.31
Monte Carlo	4	0.131	42.52	−1.03	6.51	101.73	75.12
Monte Carlo	20	0.143	29.37	−2.33	15.36	76.12	62.49
Monte Carlo	80	−0.072	25.41	−3.21	19.85	74.06	49.33
Monte Carlo	320	−0.074	23.33	−4.01	18.94	64.83	48.99

we report the summary statistics for the distribution of hedging error, as well as the 99% and 95% VaR, obtained by the above two methods. From the numerical analysis, it follows that the distribution of hedging error for a small number of trades has a high volatility and is negatively skewed, which indicates that it is more likely to have large losses than large wins. It is noteworthy that the ME estimation of the distribution fails to sufficiently capture the right tail of the empirical distribution, nevertheless it adequately describes the left tail.

As we increase the frequency of the delta hedging trades, the distribution becomes more symmetric, and our ME estimation performs better. For $N=20$, our distribution is strongly leptokurtic, indicating the option writers’ large exposure to jumps. When trading frequency increases to 80 and 320, the distribution peaks around zero and the volatility significantly decreases, however, the left tail is still noticeable. From our analysis, it appears that the volatility is negatively related to the frequency of delta hedging trades. However, further research is required to conclusively establish the relationship between the number of trades and moments of the distributions.

5. Conclusions

Perfect replication of financial derivatives is not possible, given market incompleteness and discrete-time hedging. We characterised the risk in dynamic hedge of European options through the asymptotic distribution of hedging error. Furthermore, we obtained an estimation for the distribution of hedging error by maximising Shannon’s entropy (Shannon (1948)) subject to a set of moment constraints, which in turn yielded the value-at-risk (VaR) and expected shortfall (ES) of the hedging error, two widely-used risk metrics in finance. This new approach chooses the probability distribution with the most uncertainty subject to what is known. Thus we obtained a consistent estimate of the asymptotic distribution of hedging error, despite the non-normality of the underlying distribution of returns. As a result, we derived a very generalised modelling framework, which can be applied in different areas of derivatives pricing. Some parametric specifications of this framework include, but are not limited to, the jump diffusion model of Merton (1976), the generalised gamma process discussed in Lo and Weng (1989), the variance gamma process by Madan et al. (1998), and the CGMY model of Carr et al. (2002).

Finally, we conducted a robust numerical simulation of the result, to highlight the practical applications of our model.