European Option Pricing of electricity under exponential functional of L\'evy processes with Price-Cap principle

TL;DR

This paper introduces a novel electricity pricing model using exponential functionals of jump Lévy processes, capturing market features like mean reversion and jumps, and provides a numerical scheme for option valuation.

Contribution

The paper develops a new Lévy process-based model for electricity prices incorporating the price cap principle and offers a finite difference method for option valuation.

Findings

The option price is characterized as a viscosity solution of a PIDE.

A finite difference scheme is proposed with proven stability and convergence.

Numerical simulations demonstrate the model's applicability.

Abstract

We propose a new model for electricity pricing based on the price cap principle. The particularity of the model is that the asset price is an exponential functional of a jump L\'evy process. This model can capture both mean reversion and jumps which are observed in electricity market. It is shown that the value of an European option of this asset is the unique viscosity solution of a partial integro-differential equation (PIDE). A numerical approximation of this solution by the finite differences method is provided. The consistency, stability and convergence results of the scheme are given. Numerical simulations are performed under a smooth initial condition.

Tables

Table 1. Table 1: scheme parameters

$T$	$M$	$N$	$A_l$	$A_r$
$1$	$100$	$175$	$-0.096$	$0.079$

Table 2. Table 2: Model parameters

figure	model	$r$	Strike	Product
1, 2, 3, 4	$\alpha =0.015$ $\beta =0.4$ ${\sigma }_J=0.5$ $\ell =1.5$ $\sigma =0.5$ $S_0=50$	$0.04$	$K=45$	Call

Equations

$$\text{mbox}d{S}_t=\left(\alpha (t)S_t-\beta (t)\right)\text{mbox}dt+{\sigma }_t{S}_t\text{mbox}d{W}_t+(J-1)S_t\text{mbox}d{q}_t,$$

$$S_t=S_0{e}^{X_t}-{\int }_0^t\beta (s)e^{X_t-X_s}ds,$$

$$X_t={\int }_0^t\alpha (s)-\frac{1}{2}\sigma (s{)}^{2}ds+{\int }_0^t\sigma (s)d{W}_s+{\int }_0^t\ln Jd{q}_s.$$

$$C_t=\mathbb{E}[e^{-r(T-t)}{H}_T|{\mathcal{F}}_t],$$

$$C(t,S)=\mathbb{E}[e^{-r(T-t)}{H}_T|S_t=S].$$

$$C:(0,T)\times (0,\infty )$$

$$(t,S)$$

$$\frac{\partial C}{\partial t}(t,S)+(\alpha (t)S-\beta (t))\frac{\partial C}{\partial S}(t,S)+\frac{\sigma (t{)}^2{S}^2}{2}\frac{{\partial }^{2}C}{\partial S^2}(t,S)-rC(t,S)$$

$$+\ell {\int }_{\mathbb{R}}\nu (dx)[C(t,xS)-C(t,S)]=0$$

$$d\widetilde{C_t}$$

$$\ell {\int }_{\mathbb{R}}\nu (dx)(C(t,S_{t}x)-C(t,S_t))dt,$$

$$d\widetilde{C_t}=a(t)dt+d{M}_t,$$

$$a(t)$$

$$d{M}_t$$

$$|C(t,x)-C(t,y)|$$

$$|C(t,x)-C(t,y)|\le c|x-y|e^{{\int }_t^T\alpha (s)ds}\le c_1|x-y|,$$

$$\displaystyle\mathbb{E}\big{[}\int_{0}^{T}\int_{\mathbb{R}}\nu(dx)|\varphi(t,x)|^{2}dt\big{]}$$

$${\int }_0^T{\int }_{\mathbb{R}}{e}^{-rt}[C(t,x{S}_{t^{-}})-C(t,S_{t^{-}})]d{\tilde{q}}_t$$

$$\mathbb{E}\big{[}\int_{0}^{T}S_{t}^{2}|\dfrac{\partial C}{\partial S}(t,S_{t})|^{2}dt\big{]}\leq c_{2}^{2}\mathbb{E}[\int_{0}^{T}S_{t}^{2}dt]<\infty.$$

$$u(\tau ,x)$$

$$\{\begin{array}{l}\frac{\partial u}{\partial \tau }=\mathcal{L}u,on(0,T]\times O\\ u(0,x)=H(S_0{e}^x),x\in O,u(\tau ,x)=0,x\in O^c,\end{array}$$

$$\mathcal{L}u(\tau ,x)$$

$$\{\begin{array}{l}\frac{\partial u}{\partial \tau }=({\mathcal{L}}_D+{\mathcal{L}}_J)u,on(0,T]\times O\\ u(0,x)=H(S_0{e}^x),x\in O,u(\tau ,x)=0,x\in O^c,\end{array}$$

$${\mathcal{L}}_{D}u(\tau ,x)$$

$${\mathcal{L}}_{J}u(\tau ,x)$$

$$\frac{u^{n+1}-u^n}{\Delta t}=L_D{u}^{n+1}+L_J{u}^n.$$

$$\{\begin{array}{l}\frac{\partial u_{l,r}}{\partial \tau }=({\mathcal{L}}_D+{\mathcal{L}}_J)u_{l,r},on(0,T]\times (-A_l,A_r)\\ u_{l,r}(0,x)=H(S_0{e}^x),x\in (-Al,Ar)u_{l,r}(\tau ,x)=0,x\notin (-A_l,A_r).\end{array}$$

$$|u(\tau ,x)-u_{l,r}(\tau ,x)|\le C_{\tau }\Vert H{\Vert }_{\infty }{e}^{-\max (A_l,A_r)+|x|},\forall x\in (-A_l,A_r)$$

$$u_{l,r}(\tau ,x)=\mathbb{E}\left[\text{mathds}{1}_{\{M_{\tau }^x<\max (Al,Ar)\}}H(Y_{\tau }^x)\right]$$

$$u(\tau ,x)=\mathbb{E}\left[H(Y_{\tau }^x)\right].$$

$$|u-u_{l,r}|$$

Full text

**Keywords: ** Mean reverting jump-diffusion, option pricing, price-cap, integro-differential equation, Viscosity solution

**subclass MSC: ** 60H10, 60H35, 65M06.

1 Introduction

In finance, options are tools that help to guard against risks. But, it is difficult to know the value of an option before the maturity date since, for this end, one has to estimate the value of the underlying in future. In the early of 1970’s, Black and Scholes [4] brought a major contribution in the evaluation of options. In the case where the underlying is a share that does not pay dividends, they construct a risk-neutral portofolio that replicates the winning profile of an option, which allows to perform the theoretical value of a European option under a closed formula. In the case of the exponential Lévy Black and Scholes model, this formula is derived from some classical results in discounting, statistics, stochastic and differential calculus. On the contrary, the valuation of options remains an open topic in the case of jump-diffusion processes due to the additional jump term that complicates calculation of option prices. This has been investigated by several authors as [5], [24] and [30]. One may distinguish two main approaches used by these authors. The first one, which relies on Fourier transform-based methods, has been introduced by Carr and Madan [9] to price and analyse European option prices. Many other authors follow the same idea to evaluate option prices. Among them one find [10], [14], [22] and [31]. The advantage of this approach relies on its high computational efficiency when the characteristic function is available. The idea here is to apply the direct discounted expectation method to evaluate the integral of the discounted payoff and risk neutral density function of the underlying process. Continuous and dicrete Fourier transforms were sucessfully applied to pricing options of three exponential Lévy models: the classical Black-Scoles model (which is a continuous exponential Lévy model), the Merton jump diffusion model (which is an exponential Lévy model with finite arrival rate of jumps), and the variance Gamma model (which is an exponential Lévy model with infinite arrival rate of jumps). However, for many other jump diffusion processes, it is not possible even by applying the Lévy-Khinchin representation to determine an analytical expression of the characteristic function. This may be due to the complexity of diffusion processes or to the fact that two or more processes are combined in the same model. Therefore one expects to evaluate option prices with rather an approximation of the characteristic function. The debate is then on the choice of a good approximation. Indeed, when the characteristic function can not be expressed explicitly, the Fourier transform method requires two levels of approximation which increases the error level in option valuation. In order to reduce approximation errors, a second approach, namely the Feynman-Kac formula-based approach, has been introduced by some authors such as [2], [6] and [11]. This approach relates the risk-neutral valuation formula to either the solution of a partial differential equation (PDE), when the model is a continuous exponential Lévy one, or to the solution of a partial integro-differential equation (PIDE) in the case of continuous Lévy process. The PDE, respectively PIDE obtained may be complex and its theoretical analysis requires new mathematical tools. However, after a numerical approximation step, the solution leads directly to the value of the option. So, we have one level of approximation. Nevertheless, theoretical analysis of the approximation (consistency, stability and convergence of the scheme) remains a challenge. Several finite difference schemes are used in the literature (see e.g. [2], [12] and [28]). Most of these works focus on pricing option when the underlying is driven by either a Lévy processes or an exponential Lévy process. Some of the main difficulties are related to the local integral term due to the fact that, on the one hand, the approximation of the risk-neutral density can often involve an infinite summation, and on the other hand, a local integral term requires a specific treatment at both theoretical and numerical levels. We may have other difficulties as the smoothness of option prices and even the degeneracy of the diffusion coefficient. To overcome such difficulties, the notion of viscosity solution was introduced by Grandall and Lions [13] for PDEs and, more generally, for PIDEs (see e.g. [2] and [6]). Precisely, one can split the integro-differential operator into a local and a local part, and then treat the local term using an implicit step, and the local term using an explicit step. This idea was applied by Cont and Voltchkova [12] to obtain a better approximation of option prices than the previous ones.

In this work, we apply the Feynman-Kac formula-based approach to evaluate option prices in a jump-diffusion model which represents the electricity spot price in a regulated market under the Price-Cap principle. Indeed, in this case, the characteristic function of the jump diffusion process is unavailable. In the next section, we present our jump diffusion model and state the main assumptions. We then taking advantage of the Markov property of the electricity spot price to prove in Section 3 that the European option price solves a linear PIDE. In Section 4, we apply an explicit-implicit scheme to compute numerically the viscosity solution of the PIDE provided in the previous section. The consistency, the stability and the convergence of this scheme are studied. The approximation method is similar to [12] except that contrary to their work (i) the market price is the sum of the exponential Lévy process and Lévy process, (ii) the presence of this additional term (see equation 2 ) has created an additional difficulty at all levels of the analysis since it still depends on the process, (iii) the parameters used in numerical simulations are derived from the Cameroonian context and (iv) it should also be noted that we worked in the case of time-dependent parameters. Finally, some numerical simulations are performed under a smooth initial condition.

2 Electricity Spot Price Model

We suppose that the underlying asset risk-neutral dynamics is in the form

[TABLE]

with $\alpha(t)=I(t)-EF(t)$ and $\beta(t):=CP(t)+SQ(t)-FU(t)$, where:

• $S_{t}$ represents the electricity spot prices process;
• $I(t)$ represents the inflation rate;
• $EF(t)$ represents the efficiency factor;
• $\sigma(t)$ represents the volatility;
• $W$ represents the standard Brownian motion;
• $CP(t)$ represents the subvention;
• $SQ(t)$ represents the service quality penalties, if any;
• $FU(t)$ represents the uncontrollable cost;
• $J$ represents the proportional random jump size;
• $\mbox{d}q$ a Poisson process such that; \mbox{d}q_{t}=\left\{\begin{array}[]{rl}0&\text{with probability}\,\,1-\ell\mbox{d}t\\ 1&\text{with probability}\,\,\ell\mbox{d}t\end{array}\right..

with the following assumptions:

• Assumption 1: The proportional random jump size $J$ is log-Normally distributed, with $\mathbb{E}\left[J\right]=1$. Hence, $\ln J\sim\mathcal{N}\left(-\frac{\sigma_{J}^{2}}{2},\sigma_{J}^{2}\right)$.
• Assumption 2: The random jump size $J$, the Poisson process $\mbox{d}q_{t}$ and the diffusion $\mbox{d}W_{t}$ are independent.

This model has been partly inspired from the economic principle price cap proposed by Littlechild in [23] for regulated telecommunication market in UK. This principle has been recently applied to electricity market [19].

We will further assume that the electricity parameters $\alpha(.)$, $\beta(.)$ and $\sigma(.)$ are bounded. It follows from the Itô formula for jump diffusion process that the exact solution of (1) is given by

[TABLE]

where

[TABLE]

3 Partial Integro-Differential Equation for Option Prices

This part aims to evaluate the price of an option (Put or Call) in the regulated electricity market under risk-neutral probability, $\mathbb{Q}$, with the terminal payoff, $H_{T}$, which is given by:

[TABLE]

where $r$ represents a free risk discounting rate, $T$ denotes the maturity, and $K$ represents the strike price. Let $S_{T}$ be the solution of (1) at T, which is the equation of the underlying. $H_{T}=H(S_{T})$, with $H(S)=(S-K)^{+}$ for European Call or $H(S)=(K-S)^{+}$ for European Put. From the Markov property, $C_{t}$ becomes

[TABLE]

Proposition 3.1**.**

Assume that the European option C given by

[TABLE]

is $C^{1,2}$, with $\partial C/\partial S$ and $\partial^{2}C/\partial S^{2}$ bounded, then $C$ satisfies the partial integro-differential equation:

[TABLE]

on $(0,T)\times(0,\infty)$ with the terminal condition $C(T,S)=H(S)$, $\;\forall S>0$, where $\ell$ represents the intensity of the Poisson process under risk-neutral measure, and the measure $\nu(dx)$ is the jump size distribution.

Proof.

The proof consists of applying Itô formula with jump as in [15] to the martingale $\tilde{C}(t,S_{t})=e^{-rt}C(t,S_{t})$. Then the result follows from the fact that the drift term is equal to zero.

By construction $\tilde{C}$ is a martingale. Applying the Itô formula to $\tilde{C}$ we obtain:

[TABLE]

From (1) this equation is equivalent to

[TABLE]

Adding and subtracting

[TABLE]

one has

[TABLE]

where

[TABLE]

and

[TABLE]

with $\tilde{q}_{t}=q_{t}-\ell t$. We now show that $M_{t}$ is a martingale. Since the payoff function $H$ is Lipschitz. Then, $C$ is also Lipschitz with respect to the second variable $S$. Indeed:

[TABLE]

Since $e^{-\int_{t}^{T}(1/2)\sigma(s)^{2}ds+\int_{t}^{T}\sigma(s)dW_{s}}$ is a martingale and we also have from assumption 1 that $\mathbb{E}[e^{\int_{t}^{T}\ln Jdq_{s}}]=1$, then we get

[TABLE]

with $c_{1}=ce^{\int_{t}^{T}\alpha(s)ds}$.

Therefore the predictable random function $\varphi(t,x)=C(t,xS_{t^{-}})-C(t,S_{t^{-}})$ satisfies:

[TABLE]

where the last inequality holds because the distribution $\nu(dx)$ of the jump sizes is assumed $\log$-normal. Indeed, we have $\int_{\mathbb{R}}x^{2}\nu(dx)<\infty$, hence $\mathbb{E}[\int_{0}^{T}S_{t}^{2}dt]<\infty$. Therefore, the compensated Poisson integral

[TABLE]

is a square integrable martingale. Since $C$ is Lipschitz, $\dfrac{\partial C}{\partial S}(t,.)\in L^{\infty}$ and \big{\|}\dfrac{\partial C}{\partial S}(t,.)\big{\|}_{L^{\infty}}\leq c_{2}. Thus,

[TABLE]

Furthermore, using the isometry formula and the preceding result, it follows that $\int_{0}^{T}\dfrac{\partial C}{\partial S}(t,S_{t})S_{t}\sigma(t)dW_{t}$ is a square integrable martingale. Therefore, $M_{t}$ is also a square integrable martingale, implying $\tilde{C}_{t}-M_{t}$ is a square integrable martingale. But $\tilde{C}_{t}-M_{t}=\int_{0}^{t}a(t)dt$ is also a continuous process with finite variation, so, from Theorem 4.13-450 in [18], one must have $a(t)=0$ $\mathbb{Q}$-almost surely, leading to the PIDE (3.1). ∎

Note that the smoothness (particularly the uniform boundedness of derivatives) assumption made on the European call option is not generally verified as discussed in [11]. In this case, option prices should be considered as a viscosity solution of the PIDE obtained in Proposition 3.1. The following proposition gives the link between option prices and the viscosity solution of the PIDE.

Proposition 3.2**.**

*(Option prices as viscosity solutions)

The forward value of the European option defined by (4) is the (unique) viscosity solution of the Cauchy problem (3.1).*

Proof.

Existence and uniqueness of viscosity solutions for such parabolic integro-differential equations are discussed in [2] in the case (the one considered here) where $\nu$ is the finite measure. In what follows, we propose a numerical solution to the PIDE which converges to the viscosity solution as proven in [12]. ∎

4 An Explicit-Implicit Difference scheme

In this section we present a numerical procedure for solving the PIDE (3.1) obtained in Proposition 3.1. Introducing the change of variable $x=\ln\dfrac{S}{S_{0}}\quad\text{and}\quad\tau=T-t$ and defining: $u(\tau,x)=e^{r\tau}C(T-\tau,S_{0}e^{x})$, we obtain:

[TABLE]

where $Y^{x}_{\tau}=S_{0}e^{x+X_{T}-X_{T-\tau}}-\int_{T-\tau}^{T}\beta(s)e^{X_{T}-X_{s}}ds$. We then obtain a PIDE in terms of $u$, given by:

[TABLE]

where $O\subset\mathbb{R}$ is an open interval which is not necessarily bounded,

[TABLE]

with $g_{\ln J}$ denoting the density function of $\ln J$.

The main idea in this method is to split the operator $\mathcal{L}$ into two parts as in [11]. We replace the differential part with a finite difference approximation, and the integral part with a trapezoidal quadrature approximation. We treat the integral part with an explicit time stepping in order to avoid the inversion problem of the dense matrix $L_{J}$ associated to the discretization of the integral term. We then rewrite the PIDE (8) as follow:

[TABLE]

where

[TABLE]

Hence, we obtain the approximate problem using the following explicit-implicit time stepping scheme:

[TABLE]

Before showing the stability of this scheme and applying discretization, the equation must be localized to a bounded domain.

4.1 Localisation to a bounded domain

To numerically solve the Cauchy problem (9), we first truncate the space domain to a bounded interval $(-A_{l},A_{r})$. Usually this leads to defining some boundary conditions at $x=-A_{l}$ and $x=A_{r}$. But here, we are in an elliptic local PIDE due to the presence of an integral term. Thus, we need to extend the function $u(\tau,.)$ to a subset $\left\{x+y:x\in(-A_{l},A_{r}),y\in\text{supp}\,g_{\ln J}\right\}$, where $\text{supp}\,g_{\ln J}=\mathbb{R}_{+}$, is the support of $g_{\ln J}$. Let $u_{A}(\tau,x)$ be the solution of the following localization problem:

[TABLE]

We will show in the following proposition that the localization error decays exponentially with the domain size $A$.

Proposition 4.1**.**

Assume $\|H\|_{\infty}<\infty$ and $C_{\tau}=\mathbb{E}\left[e^{\sup\limits_{\eta\in[0,\tau]}|X_{T}-X_{T-\eta}|}\right]<\infty$. Let $u_{l,r}(\tau,x)$ and $u(\tau,x)$ be respectively the solutions of the Cauchy problems (9) and (13). Then

[TABLE]

where the constant $C_{\tau}$ does not depend on $A_{r}$ and $A_{l}$.

Proof.

Let $M^{x}_{\tau}=\sup\limits_{\eta\in[0,\tau]}|x+X_{T}-X_{T-\eta}|$. Then

[TABLE]

Hence

[TABLE]

Theorem 25.18 in [29] and the fact that $\int_{\mathbb{R}}e^{|x|}\nu(dx)<\infty$ imply

[TABLE]

But

[TABLE]

and, since $\sup\limits_{\eta\in[0,\tau]}|x+X_{T}-X_{T-\eta}|\leq\sup\limits_{\eta\in[0,\tau]}|X_{T}-X_{T-\eta}|+|x|$, then

$\left\{M^{x}_{\tau}\geqslant\max(Al,Ar)\right\}\subset\left\{\sup\limits_{\eta\in[0,\tau]}|X_{T}-X_{T-\eta}|\geq\max(Al,Ar)-|x|\right\}$. Hence

[TABLE]

Now, using Markov’s inequality we obtain

[TABLE]

Comparing these last inequalities with ((16)) gives the desired result. ∎

4.2 Truncation of the integral

To compute numerically the integral term of the PIDE (9), we need to reduce the region of integration to a bounded interval which leads to the truncation of large jumps. We then estimate the error resulting from this approximation. Precisely, suppose a new process, $\tilde{S}_{t}$, is characterized by the fact that logarithm of the jump size, $\ln\tilde{J}$, is bounded in $[B_{l},B_{r}]$, with the associated measure $\mathds{1}_{\{y\in[Bl,Br]\}}\nu$, where $B_{l}$ and $B_{r}$ are real. We further suppose, without loss generality, that $B_{l}<0$ and $B_{r}>0$. In this case the corresponding solution to the associated PIDE is denoted by $\tilde{u}(\tau,x)$. We analyse, in the following proposition, the difference $|u-\tilde{u}|$.

Proposition 4.2**.**

One has:

[TABLE]

where $C^{\tau}=C\left[\tau\ell S_{0}+\ell\int_{T-\tau}^{T}\beta(s)(T-s)ds\right]$.

Proof.

Firstly, let $\tilde{X_{t}}$ be a new Lévy process defined by:

[TABLE]

and let:

[TABLE]

where $\tilde{Y}^{x}_{\tau}=S_{0}e^{x+\tilde{X}_{T}-\tilde{X}_{T-\tau}}-\int_{T-\tau}^{T}\beta(s)e^{\tilde{X}_{T}-\tilde{X}_{s}}ds$. Setting $R_{\tau}=X_{T}-\tilde{X}_{T}-(X_{\tau}-\tilde{X}_{\tau})$ and using the Lipschitz property on $H$, we obtain:

[TABLE]

Since $R_{\tau}$ and $\tilde{X}_{T}-\tilde{X}_{\tau}$ are independent, we have

[TABLE]

Moreover, $\left(e^{\tilde{X}_{T}-\tilde{X}_{T-u}},u\in[0,T]\right)$ being a martingale, $\mathbb{E}[e^{\tilde{X}_{T}-\tilde{X}_{T-\tau}}]=1$ and $\mathbb{E}[e^{\tilde{X}_{T}-\tilde{X}_{s}}]=1$. As a consequence,

[TABLE]

Since, for every $a\in\mathbb{R}$, $|e^{a}-1|=(e^{a}-1)+2(e^{a}-1)^{+}$ and $(e^{a}-1)^{+}\leq|a|$, then

[TABLE]

But:

[TABLE]

where $C_{1}=-\int_{-\infty}^{Bl}ye^{|y|}g_{\ln J}(y)dy$ and $C_{2}=\int_{Br}^{+\infty}ye^{|y|}g_{\ln J}(y)dy$. Replacing (29) into (25), we get:

[TABLE]

where $C=c_{1}e^{x}$ ∎

From Proposition4.1 and 4.2, $\tilde{u}$ converges to $u$ when $|B_{l}|$ and $|B_{r}|$ grow to infinity.

4.3 Explicit-implicit scheme

Define a uniform grid on $(0,T]\times(-A_{l},A_{r})$ by $\tau_{n}=n\Delta t$, $n=0$,…..,$M$, $x_{i}=i\Delta x-A_{l}$, $i=0,..,N$, with $\Delta t=T/M$ and $\Delta x=\dfrac{A_{r}+A_{l}}{N}$. Let $(u^{n}_{i})$ be the solution of the numerical scheme which we define below: Firstly, to approximate the integral terms, we use the trapezoidal quadrature rule with the same resolution $\Delta x$. Let $K_{l}$ and $K_{r}$ be such that $[B_{l},B_{r}]\subset[(K_{l}-1/2)\Delta x,(K_{r}+1/2)\Delta x]$, then:

[TABLE]

where $\nu_{j}=\int_{(j-1/2)\Delta x}^{(j+1/2)\Delta x}g_{\ln J}(y)dy$. Notice that to compute the integral term, we need to extend the solution to $[-A_{l}+B_{l},Ar+B_{r}]$. Hence, we assume that this solution is zero except in $[-A_{l},A_{r}]$. The derivatives are discretized using the finite difference method thus:

[TABLE]

where $f(\tau,x)=\alpha(T-\tau)-\dfrac{1}{2}\sigma^{2}(T-\tau)-\dfrac{\beta(T-\tau)}{S_{0}}e^{x}$.

Using (31) and (32), and supposing $f(\tau,x)<0$, we obtain the following relation:

[TABLE]

where

[TABLE]

Finally, we replace the problem (9) with the following time-stepping numerical scheme:

[TABLE]

After defining the numerical scheme, we study some of its important properties, particularly, consistency, monotonicity, stability and convergence.

4.4 Consistency

The follow proposition shows that (35) is consistent with (9).

Proposition 4.3**.**

*(Consistency)

The finite difference scheme (35) is locally consistent with equation (9): That is, $\forall\,v\in C^{\infty}_{0}([0,T]\times(A_{l},A_{r}))$, and $\forall\,(\tau_{n},x_{i})\in[0,T]\times\mathbb{R}$, one has:*

[TABLE]

as $(\Delta t,\Delta x)\rightarrow(0,0)$. In other words, $\exists c>0$ such that: $\;\;|r^{n}_{i}(\Delta t,\Delta x)|\leq c(\Delta t+\Delta x).$

Proof.

Let

[TABLE]

Using the second order Taylor expansion with respect to $\tau$, we obtain

[TABLE]

Plugging this relation in the first equation in ((37)) we get:

[TABLE]

where $M=\left\|\dfrac{\partial^{2}v}{\partial\tau^{2}}\right\|_{\infty}$. We now show that $|a_{2}|$ is bounded. By the mean-value theorem, there exists $\theta\in]\tau_{n},\tau_{n+1}[$  such that:

[TABLE]

If one replaces this relation in the second equation in ((37)), then:

[TABLE]

Next, taking Taylor expansion of $v$ of order $4$ gives:

[TABLE]

[TABLE]

hence

[TABLE]

Putting this last result in ((39)) gives:

[TABLE]

since $\alpha$, $\sigma$, $\beta$ and $f$ are bounded functions. Also, since the derivatives \partial^{m+n}v\big{/}\partial\tau^{n}\partial x^{m} are bounded, it implies:

[TABLE]

One also has:

[TABLE]

Since $[B_{l},B_{r}]\subset[(K_{l}-1/2)\Delta x,(K_{r}+1/2)\Delta x]$, then we have:

[TABLE]

and using Taylor’s expansion of order one, we get:

[TABLE]

From the scheme, we have $\Delta x(j-1/2)\leq y\leq\Delta x(j+1/2)$, which leads to $-\dfrac{\Delta x}{2}\leq y_{j}-y\leq\dfrac{\Delta x}{2}$. Hence,

[TABLE]

where $M_{3}=\left\|\dfrac{\partial v}{\partial x}\right\|_{\infty}\left|\sum\limits_{j=K_{l}}^{K_{r}}\int_{(j-1/2)\Delta x}^{(j+1/2)\Delta x}g_{\ln J}(y)dy\right|$. Finally, (38), (41) and (42) imply

[TABLE]

as $(\Delta t,\Delta x)\rightarrow(0,0)$. ∎

4.5 Stability and monotonicity

Two properties are important to show convergence to viscosity solutions: stability and monotonicity of scheme.

Definition 4.4**.**

**Stability

The scheme (32) is stable if, and only if, for some bounded initial conditions, its solution exists and is bounded independently of $\Delta t$ and $\Delta x$, and uniformly bounded on $[0,T]\times\mathbb{R}$. That is to say,**

[TABLE]

We will say that a given vector $v$ is positive if all its elements are positive. We write $u\geq v$ if $u-v\geq 0$. In this part we show the stability property of the scheme, which in turn implies the discrete comparison principle, a property which has an important interpretation in finance. This property makes possible the fact that the options values computed using our numerical scheme will check arbitrage inequalities: Inequality between payoffs leading to inequality between options values.

Proposition 4.5**.**

*(Stability and the discrete comparison principle)

If \Delta t\leq 1\big{/}\sum_{j=K_{l}}^{K_{r}}\nu_{j}, the scheme (32) is stable, and hence verifies the discrete comparison principle:*

[TABLE]

Proof.

Firstly, consider (32) in the form:

[TABLE]

where

[TABLE]

The positivity of $a$ and $b$ arises from $g$ being positive. If $g$ is not positive, we change the approximation of the first-order derivatives in the scheme used. In either case, one has:

[TABLE]

It follows that $|1+a\Delta t|>|-c\Delta t|+|-b\Delta t|$, implying the matrix of linear system on $(u^{n+1}_{0},...,u^{n+1}_{N})$ has a strict dominant diagonal, hence invertible. Therefore, the solution of the linear system exists and is unique. We now show by mathematical induction that this solution is bounded. That is, if $\|H\|_{\infty}\leq\infty$ is the bounded initial condition, then, $\forall n\in\mathbb{N}$,

[TABLE]

By definition of $u^{0}$, we have $\|u^{0}\|_{\infty}\leq\|H\|_{\infty}$. Assume (45) holds for $n$. To show that it holds for $n+1$, we suppose on the contrary that $\|u^{n+1}\|_{\infty}>\|H\|_{\infty}$. By the definition of $\|.\|_{\infty}$, $\exists i_{0}\in\{0,...,n\}$ such that $|u^{n+1}_{i_{0}}|=\|u^{n+1}\|_{\infty}$, and $\forall i\in\mathbb{Z},$ $\;|u^{n+1}_{i}|<|u^{n+1}_{i_{0}}|$.

Since $a=b+c$, we can write,

[TABLE]

Moreover, as $|u^{n+1}_{i_{0}-1}|<|u^{n+1}_{i_{0}}|$ and $|u^{n+1}_{i_{0}+1}|<|u^{n+1}_{i_{0}}|$ we have

[TABLE]

Using (43) and (47), and the fact that \Delta t\leq 1\big{/}\sum_{j=K_{l}}^{K_{r}}\nu_{j}, we obtain:

[TABLE]

which contradicts our assumption. Hence $\|u^{n}\|_{\infty}\leq\|H\|_{\infty}$. ∎

Proposition 4.6**.**

*(Monotonicity)

Let $u^{n}$ and $v^{n}$ be two solutions to (32) corresponding to some initial conditions $f$ and $h$ respectively, satisfying $f(x)\geq h(x)$ $\;\forall x\in\mathbb{R}$. If \Delta t\leq 1\big{/}\sum_{j=K_{l}}^{K_{r}}\nu_{j}, then $u^{n}\geq v^{n}$, $\forall n\in\mathbb{N}$.*

Proof.

Define $w^{n}=u^{n}-v^{n}$. We show that $w^{n}\geq 0$ $\;\forall n\in\mathbb{N}$. As in Proposition 4.5, we proceed by induction. By construction, we have $w^{0}_{i}=f(x_{i})-b(x_{i})\geq 0$, $\forall i\in\mathbb{Z}$. Let $w^{n}\geq 0$, and suppose that: $\inf\limits_{i\in\mathbb{Z}}w^{n+1}_{i}<0$. Since $\forall i\in\mathbb{Z}\backslash\{0,...,N\}$, $w^{n+1}_{i}=0$, this implies that $\exists i_{0}\in\{0,...,N\}$ such that $w^{n+1}_{i_{0}}=\inf\limits_{i\in\mathbb{Z}}w^{n+1}_{i}$. Using (43) and $\Delta t\leq\dfrac{1}{\sum_{j=K_{l}}^{K_{r}}\nu_{j}}$, we have that

[TABLE]

which is a contradiction. Therefore, $\inf\limits_{i\in\mathbb{Z}}w^{n+1}_{i}\geq 0$, and hence $w^{n+1}\geq 0$. ∎

4.6 Convergence

As proved above, our scheme (35) is locally consistent, stable, monotone and verifies the discrete comparison principle. In the usual approach to the convergence of finite difference schemes for PDE’s, consistency and stability ensure convergence under regularity assumptions on the solution. These conditions are not sufficient here because the solution may not be smooth, and higher order derivatives may not exist. This is where the notion of viscosity solutions are introduced. In the second order parabolic/elliptic PDEs Barles and Souganidis [7] showed that for elliptic (or parabolic) PDEs, any locally consistent, stable and monotone finite difference scheme converge uniformly, on each compact subset $[0,T]\times\mathbb{R}$, to the unique continuous viscosity solution. Rama Cont and Ekaterina Voltchkova in [11] showed that the solution of a numerical scheme converges uniformly on each compact subset of $[0,T]\times\mathbb{R}$ to the unique viscosity solution even when the subsolution and the supersolution constructed using a numerical scheme may not have uniform continuity properties. The PIDE studied in this paper relies on the same assumptions as in [11], except that here, we are in the case of a finite activity measure since the sizes of the jumps is log-normal. Then we used the same technics to showed the convergence of the explicit-implicit scheme (35) to the viscosity solution of problem (9).

Proposition 4.7**.**

*(Convergence of the explicit-implicit scheme)

Let $H$ be a bounded piecewise continuous initial condition, then solution $u^{(\Delta t,\Delta x)}$ of the numerical scheme converges uniformly on each compact subset of $[0,T]\times\mathbb{R}$ to the solution $u$ of continuous problem (9).*

Proof.

Define

[TABLE]

The aim of this proof consists to show the following equalities $\underline{u}(\tau,x)=\overline{u}(\tau,x)=u(\tau,x)$. Before showing this equalities some preparatory results are needed.

We start by giving an equivalent expression for (35).

[TABLE]

One can define super and subsolution of 4.6 by the following definition

Definition 4.8**.**

A function w is a supersolution of 4.6 if

[TABLE]

A function z is a subsolution of 4.6 if

[TABLE]

To Avoid the problem of uniform continuity and smoothness which may not hold for $\underline{u}$ and $\overline{u}$ define in (48), it is convenient to consider smooth super and subsolutions of 9 and super and subsolutions of 4.6, and to derive the link with $\underline{u}$ and $\overline{u}$. The following results extends the comparison principle to the super and subsolutions.

Lemma 4.9**.**

For any supersolution $w$ and subsolution $z$ of 4.6 we have $z\leq u\leq w$.

Proof.

For ($i\notin{0,...,N}$) or ($n=0$ and $i\in{0,...,N}$) the above inequalities are satisfied by definition. For $n=1,...,M$,   $i\in{0,...,N}$ from monotonicity of the scheme we have

[TABLE]

∎

Lemma 4.10**.**

*Let $w$ and $z$ be a smooth supersolution and subsolution of 9 respectively. Then for all $\epsilon$, there exists $\Delta>0$ such that

$\forall\Delta t,\Delta x,\leq\Delta$,   $\forall n\geq 0$,   $\forall i\in\mathbb{Z},$  $z(\tau_{n},x_{i})-\epsilon<u(\tau_{n},x_{i})<w(\tau_{n},x_{i})+\epsilon$*

Proof.

Choose q such that $0<q(T+1)<\epsilon$ and let $\tilde{w}(\tau,x)=w(\tau,x)+q(1+\tau)$, notice that a constant function is always a solution. In fact one can see from the definition that the scheme is linear.

If $i\notin{0,...,N}$, we have

[TABLE]

If $n=0$ and $i\in{0,...,N}$,

[TABLE]

If $n\geq 1$,   $i\in{0,...,N}$ from the consistency of the scheme we obtain

[TABLE]

as $\Delta t,\,\Delta x\longrightarrow(0,0)$,   $(\tau_{n},x_{i})\longrightarrow(\tau,x)$, uniformly on $(0,T[\times O$. Therefore for any sufficiently small $\Delta>0$, for all $\Delta t,\,\Delta x\leq\Delta$, we have

[TABLE]

Combining 52, 4.6 and 55, show that function $\tilde{w}$ is supersolution of 4.6. Indeed, Lemma 4.9 implies that

[TABLE]

which is the desired property. The lower bound $z(\tau_{n},x_{i})-\epsilon$ can be proved in the same manner and then completes the proof. ∎

Following Lemma 52 and Lemma 4.6, we have the following main Lemma

Lemma 4.11**.**

Let $\underline{u}$ and $\overline{u}$ be the function define by 48. For any smooth supersolution $w(\tau,x)$ and any subsolution $z(\tau,x)$ of the problem 9, we have for $(\tau,x)\in[0,T]\times O$,

[TABLE]

Proof.

By the definition of upper and lower limits, Lemma 4.10 implies desired property. ∎

After giving some properties needed we can start the proof of convergence (i.e. showed that $\underline{u}=\overline{u}=u$). If $\overline{H},\,\underline{H}$ are smoothness functions on $\mathbb{R}$ such that $\underline{H}\leq H\leq\overline{H}$, then $w(\tau,x)=\mathbb{E}[\overline{H}(Y^{x}_{\tau})]$ and $z(\tau,x)=\mathbb{E}[\underline{H}(Y^{x}_{\tau})]$ are respectively a supersolution and a subsolution of the Cauchy problem 9. From Lemma 4.11 we obtain 56. Notice that If $w(\tau,x)-u(\tau,x),\;u(\tau,x)-z(\tau,x)$ could be made small this would imply that $\lim\limits_{(\Delta t,\Delta x)\rightarrow(0,0)\!(\tau_{n},x_{i})\rightarrow(\tau,x)}u^{(\Delta t,\Delta x)}(\tau_{n},x_{i})=u(\tau,x)$. Indeed, it remains to construct appropriate smooth approximations $\overline{H}$ and $\underline{H}$.

Let $\zeta_{1},....,\zeta_{I}$ be the discontinuity points of $H$. We suppose that the jumps of $H$ are bounded by $\delta$. Given $\epsilon>0$ and $\overline{H},\,\underline{H}$ smooth functions that satisfied the following relations

[TABLE]

We have

[TABLE]

Noting that $\bigcap\limits_{\epsilon>0}\{Y_{\tau}^{x}\in\bigcup\limits_{j=1}^{I}(\zeta_{j}-\epsilon,\zeta_{j}+\epsilon)\}=\{Y_{\tau}^{x}\in\{\zeta_{1},....,\zeta_{I}\}\}$. Since $Y_{\tau}^{x}$ has an absolutely continuous distribution, so we have $\mathbb{Q}(\{Y_{\tau}^{x}\in\{\zeta_{1},....,\zeta_{I}\}\})=0$. Consequently $\mathbb{Q}(Y_{\tau}^{x}\in\bigcup\limits_{j=1}^{I}(\zeta_{j}-\epsilon,\zeta_{j}+\epsilon))\longrightarrow 0$ as $\epsilon\longrightarrow 0$.

Therefore $w(\tau,x)-z(\tau,x)\longrightarrow 0$ as $\epsilon\longrightarrow 0$ and the inequalities $z(\tau,x)\leq u(\tau,x)\leq w(\tau,x)$ together with Lemma 4.11 implies desired result which completes the proof. ∎

Remark 4.12**.**

For $\tau=0$ the scheme does not converge to the initial condition at the discontinuous points of $H$. This is due to the fact that $\mathbb{Q}(Y_{\tau}^{x}\in\bigcup\limits_{j=1}^{I}(\zeta_{j}-\epsilon,\zeta_{j}+\epsilon))\longrightarrow\mathbb{Q}(S_{0}e^{x}\in\{\zeta_{1},....,\zeta_{I}\})=\mathds{1}_{\{x\in\{\ln(\frac{\zeta_{1}}{S_{0}}),...,\ln(\frac{\zeta_{I}}{S_{0}})\}\}}$. However, this has no practical interest since it is not important to compute the solution numerically at $\tau=0$.

5 Numerical Results

In this section we discuss to the details of the implementation of our schemas and present numerical results and some interpretation.

Before started simulation parameters scheme are take as follow

the parameters used to implement the following results are taking independent to time as follow.

Figure 1 plots call option for different spot price of the underlying at maturity versus remaining time. Note that, as remaining time increase the mean-reverting and price-cap effect in all four case cause the call prices converge to zero under the set parameters gived above. Figure 2 give comparison between illustrative curve of call for four different electricity spot price at maturity. This two figures illustrate that for different initials conditions call option prices converge verse to the same value. From this we can thus, say that different properties shown theoretically are effective under the established conditions.

Figure 3 and 4 illustrate a reality enough close to those in the classical financials markets. In the sens that in financial market the values of call option before the maturity evolved in the form of a curve which towards to the payoff line progressively and as we approach maturity.

From figures 5, 6 and 7 one can observed that the values of call decrease when strike price increase. This behaviour of the call values is from a risk management point of view what it is wished.

Analysis plots of figures 8 and 9 which illustrate call option price as a function of remaining time and spot price one can said that for a large remaining time jumps effects are not perceptible and then can effect call. Whereas a small remaining time to maturity call prices increase suddenly which express the effect of jump. We must therefore say that the jump term which allows to take into account certain reality of the electricity market is not inconsiderable since that it impact on the call prices are quite noticeable.

6 Discussion and Conclusion

The characterization of european options prices in terms of the classical solution, or, in general, in the terms of the viscosity solution of a PIDE allows the use of numerical methods to obtain efficient approximations of option values. This has been a centre of research in recent times in the case of exponential Lévy models with finite arrival or infinite arrival rate of jumps. Some authors use the finite difference method to approximate the PIDE solution (see [2]), while others like [12] approximate viscosity solutions in the case of nonsmoothness of option prices. In both cases, success (in terms of efficient approximation) has been obtained. In this paper we used their approach to evaluate European call option when the underlying is electricity. The motivation behind our approach arose from the fact that the electricity prices model presented here, by hypothesis, possesses most of the properties (as in their case) of an exponential Lévy model, and the Markov process property. We focus on the pricing of call option because put option can be deduced from the put-call using parity formula. In the mathematical point of view, numerical results have confirmed the established theoretical results. In the finance point of view, numerical results present an interpretation which was coherent with some realities in the electricity market, when it is regulated under price cap.

A limitation of this work is that it can be adapted only to option pricing with short maturity in a regulated market. In future work, we plan to study option pricing on future option.

Acknowledgements

We thank the African Center of Excellence in Technologies, Information and Communication ( CETIC ) which placed at our disposal its infrastructures. This helped us to improve conditions of work.