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

+is 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-developedmathematical model of strongly correlated stochastic processes, this model will be an efficient model for pricing carbon financial derivative.


Introduction
e 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. e 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][4][5][6][7][8][9]). In the literature, Guy Jumarie [10] derived two fractional Black-Scholes equations and endeavored to obtain the solutions of these pricing equations [11][12][13][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][18][19][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]. e 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][24][25][26][27][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. e structure of this paper is as follows. Section 2 briefly 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.

A fractional Brownian Market Model
In this section, we briefly 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][31][32][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.

Formal Definition.
As a background, we first briefly describe the definition of fBm. Definition 1. Let B H t (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: Definition 2. Let B t , t ≥ 0 be a standard Brownian motion and B H t (t > 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: Properties that follow the mfBm M H t can be seen in References [7,33].

e 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][17][18][19]. Carbon emission has no exception.
It is essential to consider long-range dependence in models.
e 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. e underlying asset price is modeled as a mixed fractional Brownian motion. e linear SDE is given as follows: It is clear that the close solution to the SDE is given as follows: In Figure 1, 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.
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. Let V(t, S t ) be the valuation of the carbon emission options at time t, and a 2 Discrete Dynamics in Nature and Society 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) e market is frictionless (2) e carbon emission allowance trading is continuous (3) e interest rate is taken as given and is constant (4) ere exists no arbitrage opportunity Like the Black-Scholes case, the PDE can be obtained via Wick calculus, which given as follows: with terminal condition V * (S t , T) � (S t − K) + for an European call option and value matching con-ditionV(S t , t) � S t − K and smooth pasting conditions (zV/zS t ) � 1 for American option, respectively. It should be noted that the classical Black-Scholes equation has no third terms in (5). Let V(S, 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: where (S t , t) ∈ (0, S max ) × (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: e initial condition is and the boundary conditions are For purpose of computation, it needs to constrict S on I � [0, S max ], where S max 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): where λ is the penalty parameter satisfying λ > 1 and [a] + � max 0, 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: where a � 1 2 In particular, we focus on the case that Hurst exponents are in ((3/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.

The Penalty Approach and Its Convergence
e 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]).

Spatial Grid.
Due to the kink in the payoff function, the space discretization of PDE is based on a nonuniform Cartesian mesh. e 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 � S 0 < S 1 < · · · < S m 1 � S max is chosen with relatively many mesh points in a given interval [S left , S right ] ⊂ [0, S max ] 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 m 1 ≥ 1 and parameter d 1 ≥ 0. e grid S is partitioned via the transformation: where the equidistant points and is mesh for S is uniform inside the interval [S left , S right ] instead of outside. e parameter d 1 dominates the segmentation of points S. For simplicity, we subsume the subscript S in the notation. e interval I � [0, S max ] is partitioned into N subintervals: with 0 � S 0 < S 1 < · · · < S N � S max . For each i � 0, 1, . . . , N − 1, we use the notation h i � S i+1 − S i and h � max for i � 1, 2, . . . , N − 1, where we use the notation l i � S i+(1/2) − S i− (1/2) representing the length of interval (S i+(1/2) , S i− (1/2) ).
Using the one-point quadrature, we have Now, we discuss how to approximate the continuous flux ρ(V) at the midpoint S i+ (1/2) . e 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]. e flux has the following characterization.
Case 1 (approximation of ρ at S i+(1/2) for i ≥ 1). e twopoint boundary value problem is as follows: Integration of (19a) results in the first-order linear equation: where C 1 denotes an additive constant. Solving this linear system gives where η i � (b i+(1/2) /a).
Case 2 (approximation of ρ at S (1/2) ). e approximation of the flux at S (1/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: Solving this problem analytically gives Now, using (17) and (22a)-(22b), we define a globally piecewise constant approximation to us, 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]: where for i � 2, . . . , N. is is first-order linear ODE of V(t). Now, we consider the discretization of (23). Let t i (i � 0, 1, . . . M) be the partition points of [0, T] with T � t 0 > t 1 > t 2 · · · > t M � 0 whose subinterval length Δt m � t m+1 − t m < 0. Applying the two-level implicit time-stepping method with a splitting parameter θ ∈ [0, 1] to (24), the following full discrete system is generated: for μ m � 0, 1, . . . , M − 1, where E m and E m+1 denote the values of E(t) at t m and t m+1 , respectively, and 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. e standard Newton method is adopted for solution to the nonlinear system (26). e term d(V i ) employs the technique proposed in [16] to approximate smooth.

Theorem 1.
e matrix of (26) is an M-matrix.

Lemma 2 (stability). e discretization of (26) is stable if it is sufficiently small, where
e proofs of these theorem and lemmas are similar to those of [27], and thus we omit this discussion.

Numerical Results
In this section, we conduct some numerical experiments to show the rate of convergence of our scheme and effectiveness of the penalty.
e solution of such a nonlinear equation is obtained numerically. e default parameters are given as S � 256, K � 100, r � 0.0210, σ � 0.105, T � 0.5, and H � 0.76. All the computations are conducted in MATLAB. e solution on the finest mesh is first to approximate the "exact solution" V exact for each instance because analytical solutions are unknown. en, we compute the ratios of the numerical solutions of the consecutive meshes as follows: where V β α denotes the solution on the mesh with spatial mesh α and time mesh size β in the domain 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. e numbers show the rates of convergence. Numerical tests for the fitted finite volume method combined with the Crank-Nicolson scheme are investigated. e convergence rate of the Crank-Nicolson scheme is clearly observed from Table 1. e results indicate that the method is very efficient.
In order to investigate the impact of parameters on the optimal exercise timing, we display the optimal exercise boundary for different Hurst exponents in Figure 2.
We plot the evaluation of the carbon emission call options with different Hurst indexes. e plots show that numerical scheme is stable and robust in these two figures. Obviously, when Hurst exponent is H � (1/2), it means (10) Discrete Dynamics in Nature and Society 5 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 Figure 3 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. is 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.

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. e 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 longrange 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. e 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.
Data Availability e numerical simulation data used to support the findings of the study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.
Acknowledgments e paper was supported by "the Fundamental Research Funds for the Central Universities", South-Central University for Nationalities (CSY18013).