Abstract
In financial markets, there exists long-observed feature of the implied volatility surface such as volatility smile and skew. Stochastic volatility models are commonly used to model this financial phenomenon more accurately compared with the conventional Black-Scholes pricing models. However, one factor stochastic volatility model is not good enough to capture the term structure phenomenon of volatility smirk. In our paper, we extend the Heston model to be a hybrid option pricing model driven by multiscale stochastic volatility and jump diffusion process. In our model the correlation effects have been taken into consideration. For the reason that the combination of multiscale volatility processes and jump diffusion process results in a high dimensional differential equation (PIDE), an efficient finite element method is proposed and the integral term arising from the jump term is absorbed to simplify the problem. The numerical results show an efficient explanation for volatility smirks when we incorporate jumps into both the stock process and the volatility process.
1. Introduction
Due to the well-known phenomenon of volatility ‘smile’ and ‘smirk’ exhibited in option pricing processes, many attempts have been made to solve the problem by extending the classical Black-Scholes models [1] and relaxing the assumptions. One of the famous approaches is to replace the constant volatility by a process described by the volatility model such as a local volatility model or a stochastic volatility model, which is widely studied to capture the phenomenon of volatility skew. The local volatility model assumes that the local volatility of the stock is a function of stock price and time . For example, in Dupire’s Model [2], the classical Black-Scholes model is modified to include a time dependent local volatility rather than a constant volatility. The constant elasticity of variance model (CEV model) attempts to capture the stochastic nature of volatility and the leverage effect by assuming , which is firstly developed in [3] and then applied to calibrate and estimate the energy commodity market in [4]. Different from the local volatility model, the stochastic volatility model assumes that the volatility process is related to another stochastic process instead of the stock price process itself. The commonly used model is the Heston model [5], from which Heston generalises the Black-Scholes model to a two-dimensional stochastic model by allowing the volatility to follow a Cox Ingersoll Ross model (CIR) process and derives a semiclosed solution by applying the method of characteristic function. Stein and Stein [6] also promote a stochastic volatility model driven by the OU process. Other stochastic volatility models are proposed for different formation of stochastic volatility. The traditional Heston model [7] assumes that the underlying volatility process is a CIR process with the power of ; the model assumes that the diffusion of volatility process is a flipped CIR process, raising the power of . The process is the combination of the CIR process and the flipped CIR process (see [8]).
However, recent empirical study shows that the single-factor Heston model is overly too restrictive and multifactor stochastic models are required in order to obtain a more accurate result. The multifactor model is based on the modification of the term structure. The term structure of volatility is much more complicated than that in single factor models. The idea that the mean reversion rate of low frequency data is different to that for the high frequency data has been noticed by [9–11]. Thus, multifactor stochastic processes are introduced considering the different frequencies in the observed data (see [12]). Heston decomposes the stochastic volatility into multifactors by the principle analysis regarding the frequency of the volatility. The author suggests that the use of multifactor stochastic volatility may enhance the option pricing model by a large extent, and at least two factors should be taken into consideration in the study of path-independent and path-dependent option pricing problems (see [13]). The concept of time-scale is firstly proposed by Fouque to model the volatility process as a combination of fast-scale and slow-scale process (see [14–16]). In 2008, Fouque proposed a numerical algorithm based on asymptotic approximation and asymptotic homogenization to study the effect of the fast and the slow scale of the volatility OU process on option pricing (see [16]). The definition of time-scale is distinguished by the fluctuation frequency of the volatility process. The fast-scale volatility relates to the highly frequent short period fluctuation, while the slow-scale volatility relates to the less frequent and long term variation. The phenomenon of time-scale can easily be observed by stock prices generated by using the 27 years daily SPX data downloaded from the Chicago Board Options Exchange (CBOE) website, as shown in Figure 1. Slow scale volatility can be tracked from the long period variation, and it does not have to be mean reversion, while the fast scale volatility is the smaller but drastic oscillations between the peak and the bottom.
An alternative approach to capture the leptokurtic features and the implied volatility smile is the jump diffusion model. The jump process can be used to sketch the unexpected abrupt change of stock price within a short period. The pioneering work of Merton in [17] assumes that the asset return process follows a Brownian motion plus a jump process, and the jump process is a compound Poisson process with constant jump intensity and normally distributed jump-size distribution. Different from Merton’s Model, the work of Kou in [18] assumes that the distribution of the jump-size is a double exponential distribution instead of a normal distribution for the simplicity of computation. More recent work proves that combination of the stochastic volatility model and the traditional jump diffusion model leads to more accurate models. Bates introduces the SVJ (stochastic volatility with jumps) model by allowing both jump diffusion and stochastic volatility in the return process. The SVJ model is then extended in [19] to incorporate the jump term not only in the return process, but also in the stochastic process (see [19]). The SVJ model is also studied in [20], in which the authors assume an affine structure of characteristic function and apply fast Fourier transformation to solve the SVJ problem to obtain a semianalytic solution.
Many numerical algorithms have been proposed to study the option pricing. The traditional Black-scholes equation is a convection-diffusion parabolic equation, and the finite difference scheme for solving the equation is studied in detail by [21]. In [22], Hull and White suggest a modification to the explicit finite difference method (FDM) for valuing derivatives, which ensures a more accurate approximation with small time steps, and the approach has been applied to calculate bond options under two different interest rate processes. A FDM scheme to price the PIDE arising from a jump diffusion model is presented in [23], and an explicit-implicit FDM scheme was proposed for solving the PIDE to price the European and Barrier options with Levy process. Convergence and stability are also considered in [23]. Another important work on the application of the FDM method is studied in [24], in which an alternating direction implicit method (ADI) is applied to solve the PIDE arising from the Possion jump. The ADI approach is shown to be unconditionally stable and efficient when it is combined with the fast Fourier transform (FFT) methods (see [24]). The finite difference method (FDM) is applied in [25] to solve the variance swaps problem using the assumption of constant volatility (see [1]), in which the two-dimensional (2D) problem is tranformed to a system of one-dimensional partial differential equations, and the price of variance swap is calculated as the average of all the solutions. The FEM ensures more flexibility and adaptivity of mesh compared to the FDM. The FEM is suitable for pricing almost all option types. Three simple applications of the FEM approach in option pricing are given in [26], including the standard Black-scholes equation, the stochastic volatility model, and the path-dependent Asian option. The FEM is also applied to study the multiasset American type option in [27]. By adding the penalty term with continuous Jacobian and solving the final ordinary differential equation (ODE) with an adaptive variable order and variable step size solver SUNDIALS, the authors prove that their approach is efficient even for multidimensional PDEs.
In this paper, our contribution includes two aspects. Firstly, we extend the multiscale volatility model in [28] by incorporating both multiscale volatility processes and the jump diffusion process to price European options and the expectation of the realised volatility. The jump term is included both in the stock process and in the volatility process, and the correlation effect is also taken into consideration. Secondly, we develop an efficient FEM method to solve the problem numerically. Inclusion of both of the two factors and the jump results in a high dimensional partial integral differential equation (PIDE). Our chosen element is eight-nodal hexahedron, which can be seen as a tensor product of three one-dimensional non-parametric elements. This largely simplifies the problem by absorbing the integral part in only one tensor (one dimensional problem).
The paper is organised into five sections. In Section 1, we briefly introduce the background of the work and our main contribution. Section 2 describes the model of the underlying asset price, with the volatility following a multiscale stochastic process, which incorporates the jump diffusion term. Section 3 presents the numerical algorithm we use to solve the problem. Numerical results are given in Section 4, followed by a conclusion in Section 5.
2. Model Setup
The price of stock is assumed to follow the following stochastic process:where is a function of two factors and which denotes fast/slow scale volatility. If , the volatility process is formed by a CIR process; if , the volatility is a process which can be viewed as a combination of the CIR process and the process, and the assumption is in line with the idea that the volatility should not remain too close to zero (see [8]). It is assumed that and follow the stochastic processes
The concept of fast-scale and slow-scale is distinguished by the frequencies of the observed volatility data, and it is suggested to consider them simultaneously in [12]. Additionally, we assume that the Brownian motion are correlated with the following correlation: , , and for simplicity.
In this paper, we consider both the European option and the variance swap. For the case of European put option, the pay-off function at the maturity time is
Variance and volatility swap are well-known financial derivatives which allow investors to trade the realized volatility against the current implied volatility. Different from European options, variance swap and volatility swap are time-dependent. The pay-off function of a variance swap, as shown in (5), is convex in volatility. This phenomenon indicates that the variance swap will boost the gains and discount the losses, which explains why the variance swap is more attractive than the volatility swap. The difference between the realized volatility and the implied volatility is that the realized volatility is calculated by applying the historical data of option prices, while the implied one is derived from the prices of options.
The realized volatility is commonly approximated by the following two formulas:where denotes the underlying stock price at the th time step. AF is the annualized factor and if the sampling frequency is every month. In this paper, we let as a simplification. The pay-off of the variance swap iswhich is equal to zero under the assumption of zero entry cost. Therefore, the fair strike price can be defined as . As a result, the variance swap pricing problem becomes calculating the expected value of the realized variance in the risk neutral world.
We apply the dimensional reduction technique due to [25] by introducing a new variable driven by the underlying processwhere is the Dirac delta function, which means if and if . The terminal condition becomesFor the reason that we are more interested in the relationship between the maturity time and the strike price, we construct a new variable and then obtainAccording to the Ito formula and (1), we obtain a new processIf the problem in question is an European put option,and the pay-off function isIf the investigated problem is a variance swap, we have two different situations,where if the jump rate follows the double exponential distribution as in Kou’s model with the density ofIn contrast to the model (1) which absorbs the jump in the stock process only, the multidimensional jump process is more interesting. With this motivation, we include the jump process in both the stock price process and the multiscale volatility process, namely,However, incorporating more factors makes the model harder to tackle with, and thus development of an efficient numerical method for high dimensional PIDE is of great importance.
3. Algorithm of FEM
According to the Feyman-Kac theorem, we obtain the following partial differential equation:with the infinitesimal generator of the three-dimensional Markov process . Letting , we obtainwithwhich can be rewritten in vector form bywhere
In order to obtain option price, we have to solve the differential equation (22). However, different from , is a dynamic process which is related to time. Letting , then (22) is divided into different partial differential equations. We can then solve them one by one and then substitute the solutions back into (6) to obtain the .
The weak form of (22) can be written as
Thus, by applying the Green Theorem, we obtainwhich is derived by using the divergence theorem we assume that the test function vanishes on the boundary, and denotes inner product.
Letting , , then we obtain the following ODE system:where the mass matrix , the matrix of the diffusion part , the matrix of the convection part , , and denotes the matrix of the integral part:
The 8-node hexahedral elements can be seen as the tensor product of three one-dimensional linear elements,with , , and denoting the one-dimensional (1D) shape function in each dimension. For example, assume the natural shape functions in one-dimension as
If it is in x-dimension ; if it is in y-dimension ; else, . Thus, the shape functions of the 8-node hexahedral element arewith , , and denoting the natural coordinates of the ith nodes. To be more specific, the 8-node hexahedral element can be expanded into form shown in (see Table 1).
Let denote the integral term,where is a double exponential density function and according to (32) and is determined by the relationship between integers and . Substituting (32) into (28), the integral term can be rewritten aswith function approximating by the finite element interpolation .
The detail proof is in the Appendix. By simple calculation, we obtain
Therefore, can be seen as a Kronecker product of inner products in three dimensions:Moreover,
Let , (27) can be written as
To solve the ODE system (37), we simply apply the backward Euler method, considering its unconditional stability property
4. Numerical Results and Discussion
In this section, we present our numerical results for European options and variance swap by allowing both multiscale volatility and jump properties. Firstly, we start simulating both the stock price process and the multiscale volatility processes to show the motivation of our study. Then we apply the FEM algorithm to solve the three-dimensional PIDE. The validity of our algorithm is verified by comparing our results with the results of the two factors Heston Model in [13]. Pricing of variance swap is also studied in our paper as an application.
4.1. Validity and Motivation of Our Model
To show the motivation of our model, we firstly apply Monte Carlo simulation to generate a sample path of the stock price. Figure 2 is the stock price generated by the models (1) and (2) by the classic Euler-Maruyama Method [29]. As we can see from the figure that the asset process is a martingale process and upward sloping.
In terms of the algorithm validity, we apply our FEM method to solve the model and compare the result with the semi-analytical result shown in [13]. It is seen from Figure 3(a) that our result is well fitted. Figure 3(b) shows the underlying trajectory of the fast scale volatility process, which is highly oscillated due to the small fast-scale rate . The slow scale volatility is simulated in Figure 3(c) with the slow scale rate . The incorporation of jump process in both stock price processes has practical significance, as shown in Figures 4(a) and 4(b).
(a) Comparison of FEM with semianalytic solution
(b) Simulation of fast-scale volatility process
(c) Simulation of slow-scale volatility process
(a) Simulation of jump
(b) Simulation of stock price process with jump
However, analytic solution only exists for some special cases if we can find the characteristic function. For other models, it is not possible to obtain.
4.2. The Effects of Multi-Scale Volatility and Jump Term
Our method is applied to determine the option price of the classical European option model (12) as well as the strike price of variance swap with the payoff function shown in Figure 8(a). To be specific, let , , , , and . Both the fast-scale process and the slow scale process are assumed to be mean reverted process. The parameters we selected are from the calibrated results of [30]. , our model reduces to the original multiscale volatility model by [31]. Parameters in (1) and (2) are shown in Table 2.
Figures 5(a) and 5(c) are the surface plot of the option price and strike price of variance swap when the value of slow-scale stochastic volatility is fixed and equivalent to . If both stock price and volatilities are all variables, we obtain the three-dimensional plot shown in Figures 5(b) and 5(d).
(a) 2D view of European put option
(b) 3D view of European put option
(c) 2D view of strike price of variance swap
(d) 3D view of strike price
To show the validity of our approach, another set of data shown in Table 3 has been applied from [13] with the numerical result shown in Figures 6(a) and 6(b). When , the jump process here is assumed to be a double exponential process with , , and .
(a) 1D view of European put option
(b) 3D view of European put option
(c) 1D view of jump effect
(d) 2D view of jump effect
It can be seen from Figures 6(c) and 6(d) that the jump intensity has significant effect on option price. Hence, our model is more general compared to the multi-factor Heston model. The option price increases with the jump intensity , mainly because the growth of jump intensity leads to large uncertainty and risk exposure rate, which offers investors more possibilities to be in the money.
Different from jumps, the effect of stochastic volatility is a combination result of fast-scale and slow-scale volatility correction. The effects of fast scale rates and slow scale rates are displayed in Figures 7(a) and 7(b), respectively. As we can see from Figure 7(a), the option price increases with the fast-scale rate, while in Figure 7(b), the option price decreases with the slow scale rate, and the effects of fast-scale rate outweigh the effects of the slow-scale rate in a short period.
(a) Fast-scale rate effect
(b) Slow-scale rate effect
(a) Option price comparison with jump including in different process
(b) Relationship between T and strike price
Also, the jump terms can also be incorporated into both the fast scale volatility and the slow scale volatility process. The change of the option price, though small, can be seen from Figure 8(a). In Figure 8(a), MSJ denotes the multiscale stochastic volatility model with jumps in the stock price process, MS1J denotes the MSJ model with one jump term in the fast scale volatility, and MSV2J denotes the MSJ model incorporating jump terms in each of the three processes. For the reason that we expand the dimension by Kronecker method, we can avoid the loop and consequently save a lot of time. To make it more understandable, we compare the time (9.3490) that we assemble 3D matrix by using the nested loop with the time of the Kronecker method , which is much faster.
We also study the fair strike price of variance swap. Figure 8(b) shows the relationship between the strike price and the maturity time of variance swap, which is anticorrelated due to the introduction of fast and slow scale volatility. The fast scale and slow scale rate we choose in this analysis are and separately. The result verifies that volatility provides a measure of risk exposure. The longer the investors hold the contract, the higher risk they are exposed.
5. Conclusions
In this paper, the finite element method and the dimension reduction technique are applied to obtain the approximate solution for the classical European option price and the fair strike price of variance swaps under both multiscale stochastic volatility and jump diffusion process. The time scale rate of stochastic volatility is used to model the short term and the long term perturbation of volatility process. Our numerical results are compared with Monte Carlo simulation and are a good fit. We find that the option price increases with the jump rate and volatility value, which is in line with the reality. In terms of the effects of multiscale volatility, it is a combination result. As assumed in our model, the volatility of the stock process is driven by both the fast scale volatility and the slow scale volatility. The fast scale volatility always relates to the short term volatility with high frequency, while the slow scale volatility relates to the long term volatility and is more smooth. The option price increases with the fast-scale rate and decreases with the slow scale rate, and the effect of slow scale volatility outweighs the effect of fast scale volatility in a long run. Also, the strike price of variance swap is found to be anti-correlated with the maturity time. Volatility is a measure of risk, and the strike price decreases as the maturity time increases.
The significance of this work is in two aspects. Firstly, the exact solution can be obtained only for some specified models. For most partial differential equations, especially the high dimensional ones, closed form solution is impossible to obtain, and hence it becomes necessary to use the numerical approach in order to solve the problem. Secondly, even though stochastic volatility has already been considered in some work, multifactors in volatility are not to be tackled due to the high dimensional difficulties, and we combine both multiscale rate and jump process to make the result more reliable. In addition, the numerical method and dimensional reduction technique proposed in our paper can be applied to solve some other similar three-dimensional pricing problems.
For the future research, we are more interested in the calibration of this model with high-frequency data and the study of the skew effect. The application of the multiscale model in other financial derivatives can also be a future study area.
Appendix
In this appendix, we show the detail of approximating the integration term by interpolation method. We have the following approximation:withwith . To be consistent with the natural element, we map it to the interval of by letting
Thus, we can rewrite (A.2) in the form of being only related to the difference of .
Data Availability
The data used to support the findings of this study are included within the supplementary information files.
Conflicts of Interest
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors acknowledge the support of the Humanities and Social Science Foundation of Ministry of Education of China, Grant no. 17YJC630236. The first author also acknowledges the financial support from Curtin International Postgraduate Research Scholarship (CIPRS) and Chinese Scholarship Council (CSC).