# Pricing of American Carbon Emission Derivatives and Numerical Method under the Mixed Fractional Brownian Motion

Authors: Yuling Wang; Jing Wang
Journal: Discrete Dynamics in Nature and Society
 (2021)

---

## Abstract

This paper studies the pricing of American carbon emission derivatives and its numerical method under the mixed fractional Brownian motion. To capture the long memory properties such as self-similarity and long-range dependence in the price process, we proposed a model based on a fractional Black–Scholes, which is more in line with the actual characteristics of the option market. We have outlined a power penalty approach using parabolic variation inequality and linear complementarity (LCP) which arises from mixed fractional Brownian motion. In addition, we introduced a nonuniform grid-based modification of the fitted finite volume method for numerical solution. Numerically, we show the impact of Hurst exponent on the pricing and analyze the convergence rates of the proposed penalty method. In conclusion, since mfBm is a well-developed mathematical model of strongly correlated stochastic processes, this model will be an efficient model for pricing carbon financial derivative.


---

## Body

## 1. Introduction

The global harmony of action against carbon emission is becoming a vital challenge for human beings. It has gradually attracted too much focus in various researcher fields. Evidently, energy saving and emission reduction is the most effective measure to deal with carbon emissions. [1] Carbon trading not only could control the total amount of carbon emissions but also could minimize the cost of emission reduction through market-oriented approaches. Because the carbon emission derivatives can be applied as a tool for arbitrage and hedging against adverse movements in the energy saving and emission reduction, carbon derivative trading has timely appeared in the international financial market. [2] Since the EU officially launched the EU carbon emissions trading system in 2005 (EU ETS), it has been running for 11 years. In recent years, a large number of academic research topics has been emerged. Among of them, the rational pricing of carbon emission derivatives is essentially important.The pricing is critically dependent on assumptions about the distribution of underlying asset. A large amount of empirical evidence indicates that the logarithmic returns of underlying in general indicate self-similarity, heavy tails, and long-range dependence (see [3–9]). In the literature, Guy Jumarie [10] derived two fractional Black–Scholes equations and endeavored to obtain the solutions of these pricing equations [11–14]. Cheridito [15] proved that the model with mixed fractional Brownian motion may be arbitrage free under certain condition. Elliott [16] presented a closed solution of perpetual American options with fractional Brownian motion and thereafter induced series of papers concerning the perpetual American options with fractional Brownian motion (see [17–20]).According to the evident statistical properties of the volatility such as self-similarity and long-range dependence, it is reasonable to apply the fractional Brownian motion in asset pricing. It is clear that traditional Black–Scholes models are not sufficient for capturing many empirical features including non-normality, dependence, and nonlinearity. In this paper, we propose a fractional Black–Scholes model and its hedging strategy for fractional European option based on the fractional risk neutral measure.In the literature, there exist several numerical methods for pricing European and American option. For instance, the power penalty method for pricing American stock options was adopted in [21]. The fitted finite volume method was used for pricing standard European stock options in [22], and it was further expanded to evaluate other types of options (cf., see [23–28]).While models have been developed for pricing carbon emission allowance derivatives, but on how to utilize mixed fractional Brownian motion for pricing carbon financial derivatives still remains an open question. Based on [29], we employ a novel fitted finite volume approach to American carbon emission option pricing under the mixed fractional Brownian motion. We illustrate how to price carbon emission allowance spot options using the fractional Black–Scholes model under framework of mixed fractional Brownian motion.Motivated by the work of [16, 29], our paper presents the trend of carbon emission allowance price and regards the long-term memory process as the underlying characteristics. Specifically, we use the mixed fractional Brownian motion (mfBm) to capture special characteristics of underlying asset. Following Ref. [15], we assume H (Hurst index) > 1/2 because it is the most circumstance that frequently encountered in real financial data. Moreover, we purpose a numerical method for pricing American options under assumption that carbon emission allowance price is driven by a geometric FBM. While obtaining a valuation model of carbon derivatives which duplicate the assets portfolio, we further investigate the availability and feasibility of the numerical method.The structure of this paper is as follows. Section2 brieﬂy outlines the definition and properties of mfBm and then derives the pricing model in the framework of mfBm. In Section 3, we present the numerical algorithm that is combined the fitted finite volume method with the power penalty function. In Section 4, numerical experiments are conducted to demonstrate the usefulness of the method and the convergence of numerical solution is discussed. Section 5 presents conclusions.


## 2. A fractional Brownian Market Model

In this section, we brieﬂy introduce basic definitions and preliminaries. For details, see the relevant papers on the mixed fractional Brownian motion (mfBm) and fractional Ito integral (see [15, 16, 30–33]). For simplicity, we only analyze American call derivatives. Essentially, the methodology can be used for other American put option, and the similar complementarity problems can be also defined analogously.


### 2.1. Formal Definition

As a background, we first brieﬂy describe the definition of fBm.Definition 1.
LetBtH (t > 0) be a fractional Brownian motion (fBm) with Hurst index H∈0,1 given on a filtered probability space Ω,F,P for any, and the mathematical expectation and the variance are, respectively, as follows:(1)EBtH=0,EBtH,BsH=12t2H+s2H−t−s2H.Definition 2.
LetBt,t≥0 be a standard Brownian motion and BtHt>0 be a fractional Brownian motion (fBm) with Hurst index H∈1/2,1 given on a filtered probability space Ω,F,P for any t∈R+ and a mixed fractional Brownian motion of two real constant α and β, and H is a linearly combination of these two different Brownian motions defined as follows:(2)MtH=αBt+βBtH.
Properties that follow the mfBmMtH can be seen in References [7, 33].


### 2.2. The Fractional Black–Scholes PDE

As mentioned above, the fractional Brownian motion can characterize self-similarity and long-range dependence, so it is more suitable to model financial asset price processes. It is also well known that these features are critical for asset pricing [16–19]. Carbon emission has no exception.It is essential to consider long-range dependence in models. The so-called fractional long memory stochastic model takes fractional Brownian motion (fBm) as a random source. For better modeling of underlying asset, we adapt the work of Reference [16] to describe the long memory. As a result, it may capture some important aspects of underlying assets.The underlying asset price is modeled as a mixed fractional Brownian motion. The linear SDE is given as follows:(3)dSt=rStdt+σStdBt+σStdBtH,where St=S0∈R,0≤t≤T,T>0,μ,σ1,σ2,σ3<0, the drift term µ is constant, and Bt and BtH are standard Brownian motion and fractional Brownian motion (fBm), respectively.It is clear that the close solution to the SDE is given as follows:(4)St=S0exprt+σWt+WtH−12σ2t−12σ2t2H.In Figure1, it is noted that the mixed fractional Brownian motion fluctuates more strongly than the single Brownian motion and the price process is nonstationary, which corresponds to the long memory properties in spot price process of carbon emission allowance.Figure 1
The solution of (4) under mixed fractional Brownian motion and Brownian motion.Meanwhile, it poses severe challenges because of the semimartingale property of fractional Brownian motion. However, it is relatively easy to overcome such obstacles by applying the Wick–Ito–Skorohod integral. LetVt,St be the valuation of the carbon emission options at time t, and a type of Black–Merton–Scholes equation can be obtained under the process given in (3).In order to derive the pricing formula in a mixed fractional market, the following assumptions are needed:(1)
The market is frictionless(2)
The carbon emission allowance trading is continuous(3)
The interest rate is taken as given and is constant(4)
There exists no arbitrage opportunityLike the Black–Scholes case, the PDE can be obtained via Wick calculus, which given as follows:(5)∂V∂t+12σ2St2∂2V∂St2+Hσ2t2H−1St2∂2V∂St2+rSt∂V∂t−rV=0,with terminal condition V∗St,T=St−K+ for an European call option and value matching conditionVSt,t=St−K and smooth pasting conditions ∂V/∂St=1 for American option, respectively. It should be noted that the classical Black–Scholes equation has no third terms in (5).LetVS,t denote the value of an American call option contract on a carbon emission spot with striking price K, and V∗S,t is the payoff at the expiry date T. It is well known that American option pricing can be transformed into a problem of linear complementarity. Clearly, the value of option V meets the similar strong linear complementarity as follows:(6)LV≥0,V−S−K+≥0,LV⋅V−V∗=0,where St,t∈0,Smax×0,T with terminal and boundary conditions, S denotes the underlying asset price, and L is a degenerate parabolic partial differential operator defined as follows:(7)LV=−∂V∂t−12σ2St2∂2V∂St2−Hσ2t2H−1St2∂2V∂St2−rSt∂V∂t+rV.The initial condition is(8)VS,0=maxS−K,0,and the boundary conditions are(9)V0,t=0,S⟶0,VS,t=S−K,S⟶0.For purpose of computation, it needs to constrictS on I=0,Smax, where Smax is a big number for accuracy of the solution (see [24]).We also apply the following penalty method for the complementarity problem of (6)–(9):(10)LV−λV∗−V+1/k=0,V0,t=0,VSmax,t=V∗Smax,VS,T=V∗S,where λ is the penalty parameter satisfying λ>1 and a+=max0,a for every a. For notational simplicity, we will subsume the subscript λ in Vλ. In order to apply the finite volume method, we reformulate (10) via self-adjoint transformation:(11)−∂V∂t−∂∂StaSt2∂V∂St+bStV+cv=λV∗−V+k,where(12)a=12σ2+Hσ2t2H−1,b=r−2a=r−σ2−2Hσ2t2H−1,c=r+b=2r−σ2−2Hσ2t2H−1.In particular, we focus on the case that Hurst exponents are in3/4,1 and the corresponding boundary conditions. Except for permanent American options, it is impossible to get an explicit solution and numerical methods for American option pricing are inevitable. In this regard, efficient and accurate numerical methods are critical.


## 2.1. Formal Definition

As a background, we first brieﬂy describe the definition of fBm.Definition 1.
LetBtH (t > 0) be a fractional Brownian motion (fBm) with Hurst index H∈0,1 given on a filtered probability space Ω,F,P for any, and the mathematical expectation and the variance are, respectively, as follows:(1)EBtH=0,EBtH,BsH=12t2H+s2H−t−s2H.Definition 2.
LetBt,t≥0 be a standard Brownian motion and BtHt>0 be a fractional Brownian motion (fBm) with Hurst index H∈1/2,1 given on a filtered probability space Ω,F,P for any t∈R+ and a mixed fractional Brownian motion of two real constant α and β, and H is a linearly combination of these two different Brownian motions defined as follows:(2)MtH=αBt+βBtH.
Properties that follow the mfBmMtH can be seen in References [7, 33].


## 2.2. The Fractional Black–Scholes PDE

As mentioned above, the fractional Brownian motion can characterize self-similarity and long-range dependence, so it is more suitable to model financial asset price processes. It is also well known that these features are critical for asset pricing [16–19]. Carbon emission has no exception.It is essential to consider long-range dependence in models. The so-called fractional long memory stochastic model takes fractional Brownian motion (fBm) as a random source. For better modeling of underlying asset, we adapt the work of Reference [16] to describe the long memory. As a result, it may capture some important aspects of underlying assets.The underlying asset price is modeled as a mixed fractional Brownian motion. The linear SDE is given as follows:(3)dSt=rStdt+σStdBt+σStdBtH,where St=S0∈R,0≤t≤T,T>0,μ,σ1,σ2,σ3<0, the drift term µ is constant, and Bt and BtH are standard Brownian motion and fractional Brownian motion (fBm), respectively.It is clear that the close solution to the SDE is given as follows:(4)St=S0exprt+σWt+WtH−12σ2t−12σ2t2H.In Figure1, it is noted that the mixed fractional Brownian motion fluctuates more strongly than the single Brownian motion and the price process is nonstationary, which corresponds to the long memory properties in spot price process of carbon emission allowance.Figure 1
The solution of (4) under mixed fractional Brownian motion and Brownian motion.Meanwhile, it poses severe challenges because of the semimartingale property of fractional Brownian motion. However, it is relatively easy to overcome such obstacles by applying the Wick–Ito–Skorohod integral. LetVt,St be the valuation of the carbon emission options at time t, and a type of Black–Merton–Scholes equation can be obtained under the process given in (3).In order to derive the pricing formula in a mixed fractional market, the following assumptions are needed:(1)
The market is frictionless(2)
The carbon emission allowance trading is continuous(3)
The interest rate is taken as given and is constant(4)
There exists no arbitrage opportunityLike the Black–Scholes case, the PDE can be obtained via Wick calculus, which given as follows:(5)∂V∂t+12σ2St2∂2V∂St2+Hσ2t2H−1St2∂2V∂St2+rSt∂V∂t−rV=0,with terminal condition V∗St,T=St−K+ for an European call option and value matching conditionVSt,t=St−K and smooth pasting conditions ∂V/∂St=1 for American option, respectively. It should be noted that the classical Black–Scholes equation has no third terms in (5).LetVS,t denote the value of an American call option contract on a carbon emission spot with striking price K, and V∗S,t is the payoff at the expiry date T. It is well known that American option pricing can be transformed into a problem of linear complementarity. Clearly, the value of option V meets the similar strong linear complementarity as follows:(6)LV≥0,V−S−K+≥0,LV⋅V−V∗=0,where St,t∈0,Smax×0,T with terminal and boundary conditions, S denotes the underlying asset price, and L is a degenerate parabolic partial differential operator defined as follows:(7)LV=−∂V∂t−12σ2St2∂2V∂St2−Hσ2t2H−1St2∂2V∂St2−rSt∂V∂t+rV.The initial condition is(8)VS,0=maxS−K,0,and the boundary conditions are(9)V0,t=0,S⟶0,VS,t=S−K,S⟶0.For purpose of computation, it needs to constrictS on I=0,Smax, where Smax is a big number for accuracy of the solution (see [24]).We also apply the following penalty method for the complementarity problem of (6)–(9):(10)LV−λV∗−V+1/k=0,V0,t=0,VSmax,t=V∗Smax,VS,T=V∗S,where λ is the penalty parameter satisfying λ>1 and a+=max0,a for every a. For notational simplicity, we will subsume the subscript λ in Vλ. In order to apply the finite volume method, we reformulate (10) via self-adjoint transformation:(11)−∂V∂t−∂∂StaSt2∂V∂St+bStV+cv=λV∗−V+k,where(12)a=12σ2+Hσ2t2H−1,b=r−2a=r−σ2−2Hσ2t2H−1,c=r+b=2r−σ2−2Hσ2t2H−1.In particular, we focus on the case that Hurst exponents are in3/4,1 and the corresponding boundary conditions. Except for permanent American options, it is impossible to get an explicit solution and numerical methods for American option pricing are inevitable. In this regard, efficient and accurate numerical methods are critical.


## 3. The Penalty Approach and Its Convergence

The penalized nonlinear mfBS equation (10) can only be solved via numerical approximations in practice. Fortunately, there are different numerical methods for PDEs in the literature (see, for example, [21, 22]).


### 3.1. Spatial Grid

Due to the kink in the payoff function, the space discretization of PDE is based on a nonuniform Cartesian mesh. The advantage of a nonuniform grid defined in this section can overcome this difficulty and improve the efficiency [34]. In the S-direction, a nonuniform mesh 0=S0<S1<⋯<Sm1=Smax is chosen with relatively many mesh points in a given interval Sleft,Sright⊂0,Smax containing the strike K. Because the payoff function is discontinuous at the strike price K, the application of mesh generation is employed to overcome numerical difficulties. Let integer m1≥1 and parameter d1≥0. The grid S is partitioned via the transformation:(13)Si=φxi,0≤i≤m1,where the equidistant points xmin=x0<x1<⋯<xm1=xmax with(14)xmin=sinh−1Smin−Sleftd1,xint=sinh−1Smin−Sleftd1,xmax=xint+sinh−1Smax−Srightd1,and(15)φx=Sleft+d1sinhx,xmin≤x≤0,Sleft+d1x,0≤x≤xint,Sright+d1sinhx−xint,xint≤xmax<0.This mesh for S is uniform inside the intervalSleft,Sright instead of outside. The parameter d1 dominates the segmentation of points S. For simplicity, we subsume the subscript S in the notation. The interval I=0,Smax is partitioned into N subintervals:(16)Ii≔Si,Si+1,i=0,1,…,N−1,with 0=S0<S1<⋯<SN=Smax. For each i=0,1,…,N−1, we use the notation hi=Si+1−Si and h=max0≤i≤N−1hi. We also let Si−1/2∈Si−1,Si and Si+1/2∈Si,Si+1 be two midpoints for eachi=1,2,…,N−1. We further divide 0,Smax into new partition by letting S−1/2=S0 and SN+1/2=SN. Integrating both sides of (10) over Si−1/2,Si+1/2, we have(17)−∫Si−1/2Si+1/2∂V∂tdS−SaS∂V∂S+bVSi−1/2Si+1/2+∫Si−1/2Si+1/2cVdS=λ∫Si−1/2Si+1/2V∗−V+1/kdS,for i=1,2,…,N−1, where we use the notation li=Si+1/2−Si−1/2 representing the length of interval Si+1/2,Si−1/2.Using the one-point quadrature, we have(18)−∂Vi∂tli−Si+1/2ρV|Si+1/2−Si−1/2ρV|Si−1/2+ciVili−λV∗−V+1/kli=0.Now, we discuss how to approximate the continuous ﬂuxρV at the midpoint Si+1/2. The approach is to locally approximate the flux of a function via a constant, i.e., a locally nonlinear approximation. Further details of the discussion can be found in [21, 22]. The flux has the following characterization.Case 1.
(approximation ofρ at Si+1/2 for i≥1). The two-point boundary value problem is as follows:(19a)aSV′+bi+1/2V′=0,S∈Ii,(19b)VSi=Vi,VSi+1=Vi+1,where bi+1/2=bxi+1/2,t. Integration of (19a) results in the first-order linear equation:(20)ρiV=aSV′+bi+1/2V=C1,where C1 denotes an additive constant. Solving this linear system gives(21)ρiV=C1=bi+1/2Si+1ηiVi+1−SiηiViSi+1ηi−Siηi,where ηi=bi+1/2/a.Case 2.
(approximation ofρ at S1/2). The approximation of the flux at S1/2 is similar. We need to analyze again the approximation of the flux because (19a) is degenerated. To overcome such difficulty, we consider following ODE with additional degree of freedom:(22a)aSV′+bi+1/2V=C,S∈I1,(22b)VS0=V0,VS1=V1.
Solving this problem analytically gives(23)ρ1V=aSV′+b1/2Vs1/2=12a+b1/2V1−a−b1/2V0,V=V0+V1−V0SS1,S∈I1=0,S1.
Now, using (17) and (22a)–(22b), we define a globally piecewise constant approximation to ρ1V by ρhVη satisfying ρhV=ρiV if S∈Ii∑ for i=0,1,…,N−1.
Thus, by substituting (19a)–(19b) or (22a)–(22b) into (15) and introducing a time-reverse transformation τ=T−t, the option pricing problem can be solved by semidiscretization [35]:(24)∂Vi∂τ=EitVt+dVi,where(25)li=Si+1/2−Si−1/2,Vi=VSi,t,dVi=λV∗Si−VSi,t+1/k,Vt=V1t,V2t,…,VNtT,E1t=e1,1t,e1,2t,0,…,0,Eit=0,…,0,ei,i−1t,ei,it,ei,i+1t,0,…,0,i=2,3,…,N−2,EN−1t=0,…,0,eN−1,N−2t,eN−1,N−1t,e1,0=S14l1a−b1/2,e1,1=b1+1/2S1+1/2S2η1l1S2η1−S1η1,e1,2=−S14l1a+b1+1/2−b1+1/2S1+1/2S1η1l1S2η1−S1η1−c1,ei,i−1=bi−1/2Si−1/2Si−1ηi−1Siηi−1−Si−1ηi−1li,ei,i+1=bi+1/2Si+1/2Si+1ηiSi+1ηi−Siηili,ei,i=−bi−1/2Si−1/2Siηi−1Siηi−1−Si−1ηi−1li−bi+1/2Si+1/2SiηiSi+1ηi−Siηili−ci,for i=2,…,N.
This is first-order linear ODE ofVt. Now, we consider the discretization of (23). Let tii=0,1,…M be the partition points of 0,T with T=t0>t1>t2⋯>tM=0 whose subinterval length Δtm=tm+1−tm<0. Applying the two-level implicit time-stepping method with a splitting parameter θ∈0,1 to (24), the following full discrete system is generated:(26)θEm+1+GmVm+1+θDVm+1=Gm−1−θEmVm−1−θDVm,for μm=0,1,…,M−1, where Em and Em+1 denote the values of Et at tm and tm+1, respectively, and(27)Gm=diagl1−tm,…,lN−1−tm,DVm=dV1ml1,dV2ml2,…,dVN−1mlN−1T.
Obviously, the scheme is Crank–Nicolson scheme asθ=0.5 and implicit Euler scheme as θ=1; both cases are complete stable and must be of second- and first-order accuracy, respectively.
The standard Newton method is adopted for solution to the nonlinear system (26). The term dVi employs the technique proposed in [16] to approximate smooth.Theorem 1.
The matrix of (26) is an M-matrix.Lemma 1.
(consistency). The discretization of (26) is consistent.Lemma 2.
(stability). The discretization of (26) is stable if it is sufficiently small, where(28)h=maxmaxΔSi,maxΔτi.The proofs of these theorem and lemmas are similar to those of [27], and thus we omit this discussion.


## 3.1. Spatial Grid

Due to the kink in the payoff function, the space discretization of PDE is based on a nonuniform Cartesian mesh. The advantage of a nonuniform grid defined in this section can overcome this difficulty and improve the efficiency [34]. In the S-direction, a nonuniform mesh 0=S0<S1<⋯<Sm1=Smax is chosen with relatively many mesh points in a given interval Sleft,Sright⊂0,Smax containing the strike K. Because the payoff function is discontinuous at the strike price K, the application of mesh generation is employed to overcome numerical difficulties. Let integer m1≥1 and parameter d1≥0. The grid S is partitioned via the transformation:(13)Si=φxi,0≤i≤m1,where the equidistant points xmin=x0<x1<⋯<xm1=xmax with(14)xmin=sinh−1Smin−Sleftd1,xint=sinh−1Smin−Sleftd1,xmax=xint+sinh−1Smax−Srightd1,and(15)φx=Sleft+d1sinhx,xmin≤x≤0,Sleft+d1x,0≤x≤xint,Sright+d1sinhx−xint,xint≤xmax<0.This mesh for S is uniform inside the intervalSleft,Sright instead of outside. The parameter d1 dominates the segmentation of points S. For simplicity, we subsume the subscript S in the notation. The interval I=0,Smax is partitioned into N subintervals:(16)Ii≔Si,Si+1,i=0,1,…,N−1,with 0=S0<S1<⋯<SN=Smax. For each i=0,1,…,N−1, we use the notation hi=Si+1−Si and h=max0≤i≤N−1hi. We also let Si−1/2∈Si−1,Si and Si+1/2∈Si,Si+1 be two midpoints for eachi=1,2,…,N−1. We further divide 0,Smax into new partition by letting S−1/2=S0 and SN+1/2=SN. Integrating both sides of (10) over Si−1/2,Si+1/2, we have(17)−∫Si−1/2Si+1/2∂V∂tdS−SaS∂V∂S+bVSi−1/2Si+1/2+∫Si−1/2Si+1/2cVdS=λ∫Si−1/2Si+1/2V∗−V+1/kdS,for i=1,2,…,N−1, where we use the notation li=Si+1/2−Si−1/2 representing the length of interval Si+1/2,Si−1/2.Using the one-point quadrature, we have(18)−∂Vi∂tli−Si+1/2ρV|Si+1/2−Si−1/2ρV|Si−1/2+ciVili−λV∗−V+1/kli=0.Now, we discuss how to approximate the continuous ﬂuxρV at the midpoint Si+1/2. The approach is to locally approximate the flux of a function via a constant, i.e., a locally nonlinear approximation. Further details of the discussion can be found in [21, 22]. The flux has the following characterization.Case 1.
(approximation ofρ at Si+1/2 for i≥1). The two-point boundary value problem is as follows:(19a)aSV′+bi+1/2V′=0,S∈Ii,(19b)VSi=Vi,VSi+1=Vi+1,where bi+1/2=bxi+1/2,t. Integration of (19a) results in the first-order linear equation:(20)ρiV=aSV′+bi+1/2V=C1,where C1 denotes an additive constant. Solving this linear system gives(21)ρiV=C1=bi+1/2Si+1ηiVi+1−SiηiViSi+1ηi−Siηi,where ηi=bi+1/2/a.Case 2.
(approximation ofρ at S1/2). The approximation of the flux at S1/2 is similar. We need to analyze again the approximation of the flux because (19a) is degenerated. To overcome such difficulty, we consider following ODE with additional degree of freedom:(22a)aSV′+bi+1/2V=C,S∈I1,(22b)VS0=V0,VS1=V1.
Solving this problem analytically gives(23)ρ1V=aSV′+b1/2Vs1/2=12a+b1/2V1−a−b1/2V0,V=V0+V1−V0SS1,S∈I1=0,S1.
Now, using (17) and (22a)–(22b), we define a globally piecewise constant approximation to ρ1V by ρhVη satisfying ρhV=ρiV if S∈Ii∑ for i=0,1,…,N−1.
Thus, by substituting (19a)–(19b) or (22a)–(22b) into (15) and introducing a time-reverse transformation τ=T−t, the option pricing problem can be solved by semidiscretization [35]:(24)∂Vi∂τ=EitVt+dVi,where(25)li=Si+1/2−Si−1/2,Vi=VSi,t,dVi=λV∗Si−VSi,t+1/k,Vt=V1t,V2t,…,VNtT,E1t=e1,1t,e1,2t,0,…,0,Eit=0,…,0,ei,i−1t,ei,it,ei,i+1t,0,…,0,i=2,3,…,N−2,EN−1t=0,…,0,eN−1,N−2t,eN−1,N−1t,e1,0=S14l1a−b1/2,e1,1=b1+1/2S1+1/2S2η1l1S2η1−S1η1,e1,2=−S14l1a+b1+1/2−b1+1/2S1+1/2S1η1l1S2η1−S1η1−c1,ei,i−1=bi−1/2Si−1/2Si−1ηi−1Siηi−1−Si−1ηi−1li,ei,i+1=bi+1/2Si+1/2Si+1ηiSi+1ηi−Siηili,ei,i=−bi−1/2Si−1/2Siηi−1Siηi−1−Si−1ηi−1li−bi+1/2Si+1/2SiηiSi+1ηi−Siηili−ci,for i=2,…,N.
This is first-order linear ODE ofVt. Now, we consider the discretization of (23). Let tii=0,1,…M be the partition points of 0,T with T=t0>t1>t2⋯>tM=0 whose subinterval length Δtm=tm+1−tm<0. Applying the two-level implicit time-stepping method with a splitting parameter θ∈0,1 to (24), the following full discrete system is generated:(26)θEm+1+GmVm+1+θDVm+1=Gm−1−θEmVm−1−θDVm,for μm=0,1,…,M−1, where Em and Em+1 denote the values of Et at tm and tm+1, respectively, and(27)Gm=diagl1−tm,…,lN−1−tm,DVm=dV1ml1,dV2ml2,…,dVN−1mlN−1T.
Obviously, the scheme is Crank–Nicolson scheme asθ=0.5 and implicit Euler scheme as θ=1; both cases are complete stable and must be of second- and first-order accuracy, respectively.
The standard Newton method is adopted for solution to the nonlinear system (26). The term dVi employs the technique proposed in [16] to approximate smooth.Theorem 1.
The matrix of (26) is an M-matrix.Lemma 1.
(consistency). The discretization of (26) is consistent.Lemma 2.
(stability). The discretization of (26) is stable if it is sufficiently small, where(28)h=maxmaxΔSi,maxΔτi.The proofs of these theorem and lemmas are similar to those of [27], and thus we omit this discussion.


## 4. Numerical Results

In this section, we conduct some numerical experiments to show the rate of convergence of our scheme and effectiveness of the penalty.The solution of such a nonlinear equation is obtained numerically. The default parameters are given asS=256,K=100,r=0.0210,σ=0.105,T=0.5,andH=0.76. All the computations are conducted in MATLAB.The solution on the finest mesh is first to approximate the “exact solution”Vexact for each instance because analytical solutions are unknown. Then, we compute the ratios of the numerical solutions of the consecutive meshes as follows:(29)ratio⋅∞=VhΔτ−Vexact∞Vh/2Δτ/2−Vexact∞,where Vαβ denotes the solution on the mesh with spatial mesh α and time mesh size β in the domain VhΔτ−Vexact∞=max1≤i≤N,1≤n≤MVhΔτ−Vexact.We have no exact solution to this problem. Hence, we try to solve (26) on the 1024 × 4098 nonuniform mesh nodes in space and time as the exact solution. For the numerical scheme (26), in order to get the theoretical rates of convergence, we computed errors in ⋅h at t=0 and the discrete maximum norm, ⋅h,∞, which are listed in Table 1. The numbers show the rates of convergence. Numerical tests for the fitted finite volume method combined with the Crank–Nicolson scheme are investigated. The convergence rate of the Crank–Nicolson scheme is clearly observed from Table 1. The results indicate that the method is very efficient.Table 1
Errors inL∞−norm at t = 0.
MeshL∞−normRatio16 × 640.140432 × 1280.04792.930964 × 2560.01742.7504128 × 5120.00712.4470256 × 10240.00282.5584512 × 20480.00093.0928In order to investigate the impact of parameters on the optimal exercise timing, we display the optimal exercise boundary for different Hurst exponents in Figure2.Figure 2
The exercise boundary under different Hurst parametersH = 0.5 and H = 0.76.
(a)(b)(c)We plot the evaluation of the carbon emission call options with different Hurst indexes. The plots show that numerical scheme is stable and robust in these two figures. Obviously, when Hurst exponent isH=1/2, it means (10) that is close to optimal exercise boundary with H=0.76 under the standard Brownian motion. Practically, Greeks are an important measure in risk management and options trading. Because Greek can measure different dimensions of risk in option position, we compute Greeks in our models. In Figure 2, we see that the optimal exercise boundary decreases with Hurst exponent H.We also plot the delta and gamma in Figure3 at the last time step of the American option for a set of given parameters. In Figure 3, the option value, delta, and gamma are qualitatively stable and contain no oscillations. This means that our numerical method is robust. Because Hurst parameter H plays a significant role in the mfBm model, we exhibit its impact on evaluation in Figure 4.Figure 3
(a) American option value, (b) Delta, and (c) Gamma.Figure 4
The influence of Hurst exponentH∈0.5,1 on option value.


## 5. Conclusions

Carbon options are popular financial derivatives that play an essential role in financial market and comprehensive ecological improvement. And above all, we need more efficient valuation and accurate quantification that is very important both in theory and practice analysis of the derivative market. The related carbon option pricing theory has been developed on the basis of log-normal stochastic differential models. However, the existing work [23, 24, 27] neglected the self-similarity and long-range dependence in the price process. To capture the long memory property and to rule out the arbitrage in fractional Brownian motion, this paper uses a mfBm to capture the behavior of the carbon allowance spot price and deduces the carbon allowance option pricing model in the mixed fractional Brownian environment. The numerical method was performed to illustrate the convergence and efficiency of the methods. Moreover, the numerical results are clearly in favor of our model and indicate that the method is stable. In conclusion, since mfBm is a well-developed mathematical model of strongly correlated stochastic processes, this model will be an efficient model for pricing carbon financial derivatives.



---

*Source: 6612284-2021-04-07.xml*