Jump Models with delay -- option pricing and logarithmic Euler-Maruyama scheme

In this paper, we obtain the existence, uniqueness and positivity of the solution to delayed stochastic differential equations with jumps. This equation is then applied to model the price movement of the risky asset in a financial market and the Black-Scholes formula for the price of European options is obtained together with the hedging portfolios. The option price is evaluated analytically at the last delayed period by using the Fourier transformation technique. But in general, there is no analytical expression for the option price. To evaluate the price numerically we then use the Monte-Carlo method. To this end, we need to simulate the delayed stochastic differential equations with jumps. We propose a logarithmic Euler-Maruyama scheme to approximate the equation and prove that all the approximations remain positive and the rate of convergence of the scheme is proved to be $0.5$.


Introduction
The risky asset in the classical Black-Scholes market is described by the geometric Brownian motion given by the stochastic differential equation driven by standard Brownian motion: Let N (dt, dz) be the Poisson random measure associated with a jump process which includes the hyper-exponential jump process as a special case and letÑ (dt, dz) denote its compensated Poisson random measure. Then the above equation with σ = 0 is a special case of the following equation zÑ (dz, ds) (1.2) and it has been argued in (eg. [4,10,8]) that the equation (1.2) is a better model for stock prices than (1.1).
In this paper, we propose a new model to describe the risky asset by combining the hyper-exponential process with delay. More precisely, we propose the following stochastic differential equation as a model for the price movement of the risky asset: where f and g are two given functions, and Z(t) is a Lévy process which include the hyper-exponential jump processes as a special case. The above model along with the Brownian motion component can be found in [14], where the coefficient of Brownian motion cannot be allowed to be zero. In this work, we let the coefficient of the Brownian motion to be zero and we use the Girsanov formula for the jump process to address the issue of completeness of the market and hedging portfolio missed in [14].
With the introduction of this new market model, the first question is that whether the equation has a unique solution or not and if the unique solution exists whether the solution is positive or not (since the price of an asset is always positive). We shall first answer these questions in Section 2, where we prove the existence, uniqueness and positivity of the solutions to a larger class of equations than (1.3). To guarantee that the solution is positive, we need to assume that the jump part g(t, S(t − b))dZ(t) of the equation is bounded from below by some constant (see the assumption (A3) in the next section for the precise meaning). The class of the equations our results can be applied is larger in the following two aspects: The first one is that Z(t) can be replaced by a more general Lévy process or more general Poisson random measure and the second one is that the equation can be multi-dimensional.
Following the Black-Scholes-Merton's principle we then obtain a formula for the fair price for the European option and the corresponding replica hedging portfolio is also given. To evaluate this formula during the last delay period, we propose a Fourier transformation method. This method appears more explicit than the partial differential equation method in the literature and is more closed to the original Black-Scholes formula in spirit. This is done in Section 4.
Due to the involvement of f (S(t − b)) and g(S(t − b)) the above analytical expression for the fair option price formula is only valid in the last delay period. Then how do we perform the evaluation by using this option price formula? We propose to use Monte-Carlo method to get the numerical value approximately. For this reason we need to simulate the equation (1.3) numerically. We observe that there have been a lot of works (eg. [20,11,25]) on Euler-Maruyama convergence scheme for SDDE models. There has already been study on the Euler-Maruyama scheme for SDDE models with jumps (e.g. [15]). However, in general the Euler-Maruyama scheme cannot preserve the positivity of the solution. Since the solution to the equation (1.3) is positive (when the initial condition is positive), we wish all of our approximations of the solution is also positive. To this end and motivated by the similar work in the Brownian motion case (see e.g. [13]) we introduce a logarithmic Euler-Maruyama scheme, a variant of the Euler-Maruyama scheme for (1.3). With this scheme all the approximate solutions are positive and the rate of the convergence of this scheme is also 0.5. This rate is optimal even in the Brownian motion case (e.g. [7]). Let us point out that the 0.5 rate of the usual Euler-Maruyama scheme for SDDE with jumps studied in [15] is only obtained in the L 2 sense. Not only our logarithmic Euler-Maruyama scheme preserves the positivity, its rate is 0.5 in L p for any p ≥ 2. This is done in Section 3.
Finally in Section 5 we present some numerical attempts and compared that with the classical Black-Scholes price formula against the market price for some famous call options in the real financial market.

Delayed stochastic differential equations with jumps
Let (Ω, F , P) be a probability space with a filtration (F t ) {t≥0} satisfying the usual conditions. On (Ω, F , P) let Z(t) be a Lévy process adapted to the filtration F t . We shall consider the following delayed stochastic differential equation driven by the Lévy process Z(t): where (i) f, g : R → R are some given bounded measurable functions; (ii) b > 0 is a given number representing the delay of the equation; To study the above stochastic differential equation, it is common to introduce the Poisson random measure associated with this Lévy process Z(t) (see e.g. [2,8,9,23] and references therein). First, we write the jump of the process Z at time t by Denote R 0 := R\{0} and let B(R 0 ) be the Borel σ-algebra generated by the family of all Borel subsets U ⊂ R, such thatŪ ⊂ R 0 . For any t > 0 and for any U ∈ B(R 0 ) we define the Poisson random measure, N : [0, T ] × B(R 0 ) × Ω → R, associated with the Lévy process Z by where χ U is the indicator function of U . The associated Lévy measure ν of the Lévy process Z is given by and the compensated Poisson random measureÑ associated with the Lévy process Z(t) is defined byÑ For some technical reason, we shall assume that the process Z(t) has only bounded negative jumps to guarantee that the solution S(t) to (2.1) is positive. This means that there is an interval J = [−R, ∞) bounded from the left such that ∆Z(t) ∈ J for all t > 0.
With these notations, we can write ) are known given functions of t (and z). Thus, (2.5) is a linear equation driven by Poisson random measure. The standard theory (see e.g. [2,23]) can be used to show that the equation has a unique solution. Moreover, it is also well-known (see the above mentioned books or [1]) that by Itô's formula the solution to (2.5) can be written as From this formula we see that if φ(0) > 0, then the random variable X(t) > 0 almost surely for every t ∈ [0, b].
In similar way, we can consider the equation (2.5) on t ∈ [kb, (k + 1)b] recursively for k = 1, 2, 3, · · · , and obtain the same statements on this interval from previous results on the interval t ∈ [−b, kb].
Since (2.1) is a special case of (2.5), we can write down a corresponding result of the above theorem for (2.1).
Corollary 2.2. Let the Lévy process Z(t) have bounded negative jumps (e.g. ∆Z(t) ∈ J ⊆ [−R, ∞)). Suppose that f, g : R → R are bounded measurable functions such that there is a constant α 0 > 1 satisfying g(x) ≤ α0 R for all x ∈ R. Then, the stochastic differential delay equation (2.1) admits a unique pathwise solution with the property that if φ(0) > 0, then for all t > 0 the random variable X(t) > 0 almost surely.
Proof Equation (2.1) is a special case of (2.5) with g(z, x) = zg(x). The condition g(x) ≤ α0 R implies g(z, x) ≥ α 0 > −1 for all z ∈ J and for all x ∈ R. Thus, Theorem 2.1 can be applied. Example 2.3. One example of the Lévy process Z(t) we have in mind which is used in finance is the hyper-exponential jump process, which we explain below. Let Y i , i = 1, 2, · · · be independent and identically distributed random variables with the probability distribution given by where η i > 0, p i ≥ 0, θ j > 0, q j ≥ 0 , i = 1, · · · , m, j = 1, · · · , n with m i=1 p i + n j=1 q j = 1. Let N t be a Poisson process with intensity λ. Then is a Lévy process. If m = 1, n = 1 then Z(t) is called a double exponential process. The assumption on the boundedness of the negative jumps can be made possible by requiring that q j = 0 for all j = 1, · · · , n or by replacing the negative exponential distribution by truncated negative exponential distributions, namely, For this truncated hyper-exponential process, we can take J = [−R, ∞) with R = max{R 1 , · · · , R n }.
Although this paper will mainly concern with the one dimensional delayed stochastic differential equation (2.5) or (2.1) it is interesting to extend Theorem 2.1 to more than one dimension.
LetÑ j (ds, dz), j = 1, · · · , d be independent compensated Poisson random measures. Consider the following system of delayed stochastic differential equations driven by Poisson random measures: and φ i (0) ≥ 0 , i = 1, · · · , d, then, the stochastic differential delay equation (2.6) admits a unique pathwise solution with the property that for all i = 1, · · · , d and for all t > 0, the random variable S i (t) ≥ 0 almost surely.
Proof We can follow the argument as in the proof of Theorem 2.1 to show that the system of delayed stochastic differential equations (2.6) has a unique solution S(t) = (S 1 (t), · · · , S d (t)) T . We shall modify slightly the method of [12] to show the positivity of the solution. Denoteg ij (t, z) = g ij (z, S(t− b)). Let Y i (t) be the solution to the stochastic differential equation Jg where ν j is the associated Lévy measure forÑ j (ds, dz). Denotef ij (t) = f ij (S(t − b)) and let p i (t) be the solution to the following system of equations By the assumption on f we have that when i = j,f ij (t) ≥ 0 almost surely. By a theorem in [5, p.173] we see that p i (t) ≥ 0 for all t ≥ 0 almost surely. Now it is easy to check by the Itô formula thatS i (t) = p i (t)Y i (t) is the solution to (2.6) which satisfies that S i (t) ≥ 0 almost surely. By the uniqueness of the solution we see that S i (t) =S i (t) for i = 1, · · · , d. The theorem is then proved.

Logarithmic Euler-Maruyama scheme
The equation (2.1) or (2.5) is used in Section 4 to model the price of a risky asset in a financial market and its the solution is proved to be positive as in Theorem 2.1. As it is well-known the usual Euler-Maruyama scheme cannot preserve the positivity of the solution (e.g. [13] and references therein). Motivated by the work [13], we propose in this section a variant of the Euler-Maruyama scheme (which we call logarithmic Euler-Maruyama scheme) to approximate the solution so that all approximations are always non-negative. For the convenience of the future simulation, we shall consider only the equation (2.1), which we rewrite here: Here N t is a Poisson process with intensity λ and Y 1 , Y 2 , · · · , are iid random variables.
The solution to the above equation can be written as We shall consider a finite time interval [0, T ] for some fixed T > 0. Let ∆ = T n > 0 be a time step size for some positive integer n ∈ N. For any nonnegative integer k ≥ 0, denote t k = k∆. We consider the partition π of the time interval [0, T ]: On the subinterval [t k , t k+1 ] the solution (3.2) can also be written as Motivated by the formula (3.3), we propose a logarithmic Euler-Maruyama scheme to approximate (2.1) as follows.
surely for all k = 0, 1, 2, ..., n. Then our approximations S π (t k ) are always positive. Notice that the approximations from usual Euler-Maruyama scheme is always not positive preserving (see e.g. [13] and references therein).
We shall prove the convergence and find the rate of convergence for the above scheme. For the convergence of the usual Euler-Maruyama scheme of jump equation with delay, we refer to [15]. To study the convergence of the above logarithmic Euler-Maruyama scheme, we make the following assumptions.
(A2) f is bounded. f and g are global Lipschitz. This means that there exists a constant ρ > 0 such that For notational simplicity we introduce two step processes Define the continuous interpolation of the logarithmic Euler-Maruyama approximate solution on the whole interval [−b, T ] (not only on t k , k = 0, · · · , n) as follows: (3.7) With this interpolation, we see that S π (t) > 0 almost surely for all t ≥ 0.
Lemma 3.1. Let the assumptions (A1)-(A4) be satisfied. Then for any q ≥ 1 there exists K q , independent of the partition π, such that where and throughout the remaining part of this paper, we denote T = [0, t] × J. Now we are going to handle the factor Let h = ((1 + zg(v 2 (u)) 2q − 1))/z. Then where we used boundedness of g and the assumption (A4). Now an application of the Cauchy-Schwartz inequality yields Inserting this estimate of I into (3.8) proves E sup 0≤t≤T |S π (t)| q ≤ K q < ∞. In the same way we can show E sup 0≤t≤T |S(t)| q ≤ K q < ∞. This completes the proof of the lemma.
An application of the Hölder inequality yields that for any p > 1, Now we want to bound .
(we use the same notation I to denote different quantities in different occasions and this will not cause ambiguity). We write the above sum as an integral: By the Burkholder-Davis-Gundy inequality, we have Thus, we have Inserting this bound into (3.9) yields the lemma. Our next objective is to obtain the rate of convergence of our logarithmic Euler-Maruyama approximation S π (t) to the true solution S(t). Theorem 3.3. Assume (A1)-(A4). Let S π (t) be the solution to (3.4) and let S(t) be the solution to (3.1). Then there is a constant K p,T , independent of π such that (3.10) Proof We write S(t) = φ(0) exp (X(t)) and S π (t) = φ(0) exp (p(t)). Then

Hence by Lemma 3.1 we have for any
Thus we need only to bound the above expectation I, which is given by the following. .
By the Lipschitz conditions we have (3.13) By Lemma 3.2 and by the assumption (A1) about the Hölder continuity of the initial data φ we have (3.14) We write the above sum I 3 with jumps as a stochastic integral: [ln(1 + zg(S(u − b))) − ln(1 + zg(v 2 (u)))]Ñ (du, dz) Using the Lipschitz condition on g and (A3), we have Using the Burkholder-Davis-Gundy inequality we have Similar to the bound for I 32 , we have Combining the estimates for I 31 and 32 , we see It is easy to verify Inserting the bounds obtained in (3.14)-(3.17) into (3.13), we see that Combining this estimate with (3.11), we see for any p ≥ 2 and for any r ∈ [0, T ]. Now we shall use (3.18) to prove the theorem on the interval [0, kb] recursively for k = 1, 2, · · · , for any p ≥ 2. Now taking r = 2b in (3.18), we have Continuing this way we obtain for any positive integer k ∈ N, Now since T is finite, we can choose a k such that (k − 1)b < T ≤ kb. This completes the proof of the theorem.

Option Pricing in Delayed Black-Scholes market with jumps
In this section we consider the problem of option pricing in a delayed Black-Scholes market which consists of two assets. One is risk free, whose price is described by Another asset is a risky one, whose price is described by the delayed equation (2.1) or (3.1), namely, where Z(t) = Nt i=1 Y i is a Lévy process, N t is a Poisson process with intensity λ, and Y 1 , Y 2 , · · · , are iid random variables. As in Section 2, we introduce the Poisson random measure N (dt, dz) and its compensatorÑ (dt, dz). The above delayed equation can be written as where f Y is the probability density of Y i (whose support is J). Then J zν(dz) = λL .

SetS (t) = S(t) B(t) .
Then by Itô's formula we have . We shall keep the assumptions (A1)-(A4) made in previous section and we need to make an additional assumption: [0, ∞) To find the risk neutral probability measure we apply Girsanov theorem for Lévy process (see [9,Theorem 12.21]). The θ(t) is predictable for t ∈ [0, T ]. From the assumptions above we also have that 0 < θ(s) ≤ 1 α1 . Thus, In order for us to obtain an equivalent martingale measure we need to verify the following Novikov condition: This is a consequence of our assumption (A5). In fact, we have first Hence we have Thus, we have (4.5). Now since we have verified the Novikov condition (4.5) we have then E[S θ (T )] = 1. Define an equivalent probability measure Q on F T by dQ := S θ (T )dP . (4.6) On the new probability space (Ω, F T , Q) (new probability Q) the random measurẽ is a compensated Poisson random measure. The corresponding Lévy measure is denoted by ν Q . With this new Poisson random measure we can write (4.4) as The following result gives the fair price formula for the European call option as well as the corresponding hedging portfolio.
where Q is the martingale measure on (Ω, F T ) given by (4.6).
there is an adapted and and the hedging strategy is given by where U (t) = E Q (e −rT (S T − K) + |F t ).
Proof Applying the Itô formula to (4.8) we get Denote X = (S T − K) + and consider In order to apply martingale representation theorem for Lévy process (see e.g. [2, Theorem 5.3.5]) we shall first show that U t ∈ L 2 , which is implied by t − b)). Then we can writẽ Applying the Hölder inequality we have From the definition ofh, we have zh = (1 + zh) 4 − 1. Then Thus, which is finite by the assumptions of the theorem. From the martingale representation theorem (see e.g. [2, theorem 5.3.5]) there exists a square integrable predictable mapping ψ : T × Ω → R such that Consider the strategy {(π B (t), π S (t)) : t ∈ [0, T ]} to invest π B (t) units in the riskyless asset B(t) and π S (t) units in the risky asset S(t) at time t. Then the value of the portfolio at time t is given by By the definition of the strategy we see that Hence the strategy is self-financing. Moreover, we have Hence the claim (referring to the European call option) is attainable stand therefore the market {S(t), B(t) : t ∈ [0, T ]} is complete. The pricing formula (4.9) is hard to evaluate analytically and we shall use a general Monte-Carlo method to find the approximate values. But when the time fall in the last delay period, namely, when t ∈ [T − b, T ] we have the following analytic expression for the price.
where w = ln(K/A) − rT and (4.14) Proof By (4.9) for any time t ∈ [0, T ] we have First, let us compute V 1 (t) and V 2 (t) can be computed similarly. The solutionS(t) is given by (4.11), which we rewrite here: Hence while computing the conditional expectation of h(S(T )) with respect to F t , we can consider the integrands ln(1 + zg(S(u − b))) and ln(1 + zg(S(u − b))) − zg(S(u − b)) as "deterministic" functions. Thus, the analytic expression for the conditional expectation is possible. But it is still complicated. To find the exact expression and to simplify the presentation, let us use the notation (4.14) and introduce Y = T t J ln (1 + zg(S(u − b)))Ñ Q (dz, du) .
With these notation we haveS To calculate E Q e Y I {v≥Y ≥w} we first express I [w,v] as the (inverse) Fourier transform of exponential function because E(e iξY ) is computable. Since the Fourier transform of we can write Therefore we have Taking w = ln(K/A) − rT , v → ∞ in the above formula we can evaluate (4.15) as follows.
Exactly in the same way (and now without the factor e Y ), we have This gives (4.13).

Numerical attempt
In this section we make an attempt to carry out some numerical computations of our formula (3.39) against the American call options Microsoft stock traded in Questrade platform. To apply our model in the financial market, we need to estimate all the parameters including the delay factor b from the real data. To the best of our knowledge the theory on the parameter estimation is still unavailable even in the case of the classical model of [3]. Motivated by the work of [19], we try our best guess of the parameters in the model (3.31)-(3.32).
The real market option prices we consider is for the American call option on Microsoft stock. The data we use is from Questrade trading/investment platform on October 5, 2020 at 12:25 PM (EDT). We take T to be one, three and six months active trading period respectively. The real prices of the options of different strike prices are listed in the last column of the three tables below.
The readers may wonder that since the option pricing formulas for both our model and the classical Black-Scholes model are for the European call option, why we use the market price for the American option. The reason is that we can only find the market price for the American option. On the other hand, as stated in [21, p.251] "There is no advantage to exercise an American call prematurely when the asset received upon early exercise does not pay dividends. The early exercise right is rendered worthless when the underlying asset does not pay dividends, so in this case the American call has the same value as that of its European counterpart". See also [16, p.61, Theorem 6.1]. This justifies our use of the market price for the American option.
Using Monte-Carlo simulation we calculate the prices of European option given by (4.9) and the analogous Black-Scholes formula obtained from the model: dS(t) = S(t)[αdt + σdW (t)]. We simulate 2000 paths of the solutions to both equations using the logarithmic Euler-Maruyama scheme [for Black-Scholes model the logarithmic Euler-Maruyama scheme is the same by replacing the jump process by Brownian motion]. In the simulations we take the time step ∆ to be the trading unit minute. So when T = 1 month, there are n = trading hours × 60 × trading days = 6.5 × 60 × 22 = 8580 minutes. So ∆ = 1 8580 . We do the same for T = 3 and T = 6. In our calculation for the delayed jump model we use the double exponential jump process as our Y i 's with parameters p = .60, q = 1 − p = .40, η = 12.8, θ = 8.40 with the intensity λ = .03. The interest rate r = .01 is the risk free rate. The delay factor was taken to be one day which is b = 6.5×60 8580 because there are trading 6.5 hours in a trading day. The function f (x) was taken to be a fixed constant f (x) = .1, g(x) = .15 * sin(x/209.11) and φ(x) = exp(αx/n) with α = .11. We choose α = .11 since the initial price we have taken is 209.11 and the predicted average price target of Microsoft stock for next one year (around 12 months from October 5, 2020) is 230 which is 11%.
For the simulation of the Black-Scholes model, based on stock prices for the year 2019 we take volatility of the Microsoft stock as σ = 15% to calculate Black-Scholes price. We have taken r = 1% since in the last one year the range of 10 year treasury rate has been between .52% to 1.92%.
The computations are summarized in the following tables. Notice an interesting phenomenon that the price we obtain by using our formula is comparable to the Black-Scholes price for shorter maturities and is more closer to the real market price for longer maturity. This may be because of our choice of the parameters by guessing.

Conclusion
In this paper we introduce and study a stochastic delay equation with jump and derive a formula for the fair price of the European call option. We assume that the jump is dictated by a compensated Lévy process, which includes the process like asymmetric double exponential, hyper-exponential jump process. In the numerical execution we consider the asymmetric double exponential process. Furthermore, we propose a logarithmic Euler-Maruyama scheme (a variant of Euler-Maruyama scheme) which preserve the positivity of the approximate solutions and show that the convergence rate of this scheme is 0.5 in any L p norm, the optimal rate for the classical Euler-Maruyama scheme for the stochastic differential equations driven by standard Brownian motion (see e.g. [7]). From the above tables we see that the parameters guessed here may not be the best possible values but our formula still gives a good fit to the real market prices compared to the Black-Scholes formula. We note further that potential research problem of parameter estimation is still open before we can come up with the best possible simulated results.
This research was funded by an NSERC discovery fund and a startup fund of University of Alberta.