Abstract
Empirical evidence shows that single-factor stochastic volatility models are not flexible enough to account for the stochastic behavior of the skew, and certain financial assets may exhibit jumps in returns and volatility. This paper introduces a two-factor stochastic volatility jump-diffusion model in which two variance processes with jumps drive the underlying stock price and then considers the valuation on European style option. We derive a semianalytical formula for European vanilla option and develop a fast and accurate numerical algorithm for the computation of the option prices using the fast Fourier transform (FFT) technique. We compare the volatility smile and probability density of the proposed model with those of alternative models, including the normal jump diffusion model and single-factor stochastic volatility model with jumps, respectively. Finally, we provide some sensitivity analysis of the model parameters to the options and several calibration tests using option market data. Numerical examples show that the proposed model has more flexibility to capture the implied volatility term structure and is suitable for empirical work in practice.
1. Introduction
Options’ pricing has played an important issue in the general theory of asset pricing since the celebrated work of Black and Scholes [1] (called by BS model), who assume the asset returns satisfying a Geometric Brownian motion with constant parameters. However, option pricing formulas in the BS model do not perform well empirically in practice because of the stylized facts, such as the asymmetric leptokurtic features and the volatility smile. Due to the limitations of the BS model, there has been a wide range of research carried out to improve upon this model. The BS model has been extended to incorporate either constant volatility plus jumps or stochastic volatility with/without jumps to capture these stylized facts. Merton [2] was the first to introduce a normal jump-diffusion model which can be consistent with the extreme movements sometimes observed in financial time series and lead to the leptokurtic feature, volatility smile, and closed-form solutions for options. Heston [3] employed the square-root process for modeling the variance of the underlying asset upon the BS model and showed how to use the Fourier transform to obtain the analytical formula of vanilla options. Bates [4], Scott [5], Zhang and Wang [6, 7], Deng [8], and Pillay and O’Hara [9] combined Heston’s stochastic volatility model with Merton’s normal jump-diffusion model to model the asset price dynamics and obtained the explicit formulas for European options. Eraker et al. [10] and Chernov [11] investigated Heston’s stochastic volatility model incorporating jumps both in returns and volatility. Duffie et al. [12] suggested a generalized affine jump-diffusion model (AJD), in which the drift vector, covariance matrix, and jump intensities all have affine dependence on the state vector. Consequently, the AJD model has become a main topic for research studies on equity returns and on currencies, among other applications (see Glasserman and Kim [13] and and Kaeck and Alexandex [14], for example). In the case of single-factor stochastic volatility model, one of the parameters of the model is the correlation coefficient between the underlying asset and its instantaneous variance. However, the assumption of a single factor in this model does not provide sufficient flexibility to account for the stochastic behavior of the volatility skew observed in option markets and leads to incorrect pricing formulas for deep in-and-out-of-the money European options.
In a parallel development, two-factor stochastic volatility models (TSV) have generated attention from the option-pricing literature. Christonffersen et al. [15] and Fouque et al. [16] extend Heston’s model to generate a two-factor stochastic volatility model, which is able to account for the stochastic correlation between the asset return and its instantaneous variance. Another extension to Heston’s model is proposed by [17] who consider a matrix Wishart diffusion process to introduce a correlation structure between the asset returns and its volatility factors (denoted as the MSV model). Gourieroux [18] and Gourieroux and Sufana [19] consider the Wishart Autoregressive (WAR) diffusion processes to model volatility process of the underlying asset. Recently, Muhle-Karbe et al. [20] develop a multivariate Ornstein–Uhlenbeck stochastic volatility model based on Lévy process. As demonstrated in [17], the main advantage of the MSV model is that it can calibrate short-term and long-term volatility levels better than a single-factor one. Recently, the MSV model has been applied to value exotic options (see Wong and Chan [21], Romo [22], and Zhong and Deng [23], for example).
Although the MSV model may be useful for option pricing, it is unsuccessful in the modeling of certain financial asset returns. There is empirical evidence suggesting that volatility jumps cannot be recaptured by a diffusion model only (see Jacod and Todorov [24]). Bates [25] and Chernov et al. [26] find that neither one volatility factor model with jumps nor affine two-factor specifications are well equipped to handle the very short-lived but erratic extreme events typically found during financial crisis. For these reasons, this paper aims to extend the MSV model incorporating jump components in asset returns and its volatilities, which can improve the pricing performance of the AJD and MSV models, and also offer a potentially convenient framework for modeling term structure.
This paper contributes to the literature in the following ways. We derive a semianalytic formulas for the prices of European vanilla options in a two-factor stochastic volatility jump-diffusion model by using the martingale method and Fourier transform technique. This enables us to understand how our proposed model can be adjusted to fit volatility smiles and the real data. We also obtain a fast and accurate numerical solution for the European option pricing by the Fast Fourier transform (FFT) technique. This allows us to assess the impacts of the jump components in the two-factor volatility on option price.
The outline of this paper is as follows. Section 2 introduces the model. The characteristic function of the log-stock price and the semianalytical formula for the European call option are derived. An approximation solution for this option pricing by FFT is also provided. Computational results and empirical analysis are reported in Section 3. Section 4 concludes the paper. In this section, we proposed the following model for the asset price dynamics which is driven by two-factor stochastic volatility jump-diffusion models. Under this proposed model, we derive the semianalytic form for option price and then provide a fast and accurate numerical method using FFT technique.
2. Asset Price Dynamics
Consider an arbitrage-free, frictionless financial economy consisting of two assets, namely, a bank account and a risky stock , which are traded continuously over a finite time interval . Let be a complete probability space with filtrated -field , on which are defined a four-dimensional standard Wiener process , three-dimensional Poisson processes , and three-dimensional random jump size sequences . Assume that all uncertainties on this probability space are mutually independent and is generated by these uncertainties.
We assume that the joint dynamics of the stock price and its variance processes are specified under some risk-neutral probability measure (see [12]) by the following system of stochastic differential equation:where and are constant risk-free interest and dividend yields, respectively, and are stochastic variances with the mean reversion speeds and , the long-term variances and , and the volatilities of the variances and , respectively, and and denote the correlation coefficients between the Wiener processes driving the log-stock price and the two variance processes, respectively. The Poisson processes are assumed to have constant intensities for . The jump size distribution has an associated transform of , such that , where the jump transform is defined bywith and
All parameters in this specification are assumed to be constant. To ensure positivity of the parameters must be chosen such that for . We will refer this specification as the TSVIJ model.
The TSVIJ model nests some popular models in option pricing. Without considering jumps in both returns and volatilities, the TSVIJ model can be reduced to the TSV models (called also the double Heston model) considered by Chernov et al. [26] and LeBaron [27]. If we consider a single-factor volatility process (denoted the SVIJ model for short) only in (1), the SVIJ model discussed by Chernov [11] can be regarded as a particular case of the TSVIJ model. When we consider the single-factor volatility process without jumps, the TSVIJ model can be reduced to the SVJ model discussed by Bates [4]. Without jumps in both returns and variance process and with considering a single-factor volatility process only, the TSVIJ model is degenerated to the SV model considered by Heston [3]. We have the normal jump-diffusion model (denoted as Merton model) considered by Merton [2] without considering the stochastic volatility model in (1).
The TSVIJ model specified in (1), unlike the SVIJ and SV model, accounts for a richer variance-covariance structure. In particular, the total conditional variance of the stock return process is the sum of the two variance factors plus jumps:whereas, as shown by Eraker et al. [10], the conditional variance of the variance processes and the correlation between the stock return and the variance process are given byrespectively. Note that the correlation between the stock return and variance is determined by and and depends on the current levels of both the factors and jumps. It is an important scale for modeling the leverage effect which implies a stochastic correlation.
There are two important advantages of the TSVIJ model with respect to the SVIJ/or SV model. First, the TSVIJ model displays not only stochastic variance but also stochastic correlation between stock return and variance, and this feature potentially enables the model to capture fluctuations in option skewness. Second, the TSVIJ model provides more flexibility to modeling volatility term structure; one of the factors determines the correlation between short-term stock return and variance process, while the other factor determines long-term stock return and variance process.
2.1. The Characteristic Function
Given the asset price dynamics specified in (1), it is possible to obtain the characteristic function of the log-stock price. Let . Applying the generalized Itô formula to the first equation in (1), the evolution of the log-stock price is described by
In the following, we use the martingale method and generalized Itô formula to derive the characteristic function whose solution is present in the lemma below.
Lemma 1. Suppose that the dynamics of processes and are given by (1). Define the conditional characteristic function of aswhere is a real-valued argument. Then,where , and are given bywith , , , and for .
Proof. See Appendix A.
With the characteristic function available, European options can be valued using the Fourier transform method (FT). There are increasing research efforts in option pricing using the FT method. For example, Carr and Madan [28] illustrated the fundamental idea of using the Fast Fourier Transform (FFT) for computing European vanilla options based on the characteristic function of the log-asset value. Duffie et al. [12] provide a comprehensive survey on the FT approach to option pricing in a wide range of stochastic processes. Dempster and Hoog [29] developed an approximation method for pricing European options on spreads using the two-dimensional FFT. Pillay and O’Hara [9] and Zhang and Wang [6, 7] also use the FFT technique to pricing European options under the stochastic volatility model with/without jumps. Wong and Zhao [30] use a modified FFT method (FRFT), namely, the fractional FFT, to value currency options. Černý [31] presented a detailed discussion on the implementation of FFT to option pricing. In Section 2.2, we examine how these methods can be used to our case.
2.2. Semianalytical Formula
Now, we consider a European vanilla call option pricing. Let be the strike price and the maturity of a European call option. The price at time 0 of the European call option, denoted by , is computed as the discounted risk-neutral conditional expectation of the terminal payoff as follows:
Let . Then, the call option function (12) can be expressed as a function of in integral form as follows:where is the probability density function of the process . We rewrite (13) as two expectations of indicator functions so that
For the first term in (14), we choose the spot price as numeraire and switch from probability measure to on . The Radon–Nikodym derivative is then given by
Under this new equivalent martingale measure , the call option price (14) can be restated as
The two probabilities in (16) can be calculated using the Fourier inversion transform:where denotes the real part, and and are defined bywhich are the corresponding characteristic functions of under the measures and , respectively. Substituting (17) and (18) into (16), we summarize the result as follows.
Theorem 1. Assuming the dynamics of the asset price specified by (1), the semianalytic formula for the European call option can be given by
From the Fourier transform formula, the probability density function of the log-asset price is given by
2.3. FFT Method
Note that the option price in (20) does not decay to 0 as . Therefore, the FFT method cannot be used directly to compute the integrals in (20). As in Carr and Madan [28] and Dempster and Hoog [29], is multiplied by an exponentially damping factor for so that it is square-integrable in over . Define a modified call option price function . Then, the Fourier transform of is
Given the Fourier transform function , the integral can be calculated by the inverse Fourier transform:
Set , where denotes the integration step width and is the total number of grid points and is assumed to be an arbitrary power of 2. Then, invoking the trapezoid rule, the Fourier integral in (23) can be approximated by the summation:
Recall that a one-dimensional Fast Fourier transform (FFT) computes for any complex (input) array, , and the following (output) array of identical structure is obtained:for all . In order to use the FFT algorithm to evaluate the sum in (14), we define a grid of size along the modified log-strike price by , where for . Choosing gives the following values of the sum on the grid as
Substituting the representations of and in (25) and using the choice of , we obtain
Using Simpson’s rule for numerical integration, define a sequence of weighting factors by
Then, (26) can be rewritten as
This can be computed from the expression of the FFT (25) by taking the input array assuch that
The result is an approximation for the European call option price at different (log) strikes given by
All that remains is to choose the value of . Carr and Madan [28] suggest that is chosen such that , that is, . Therefore, is well behaved.
Note that the FFT algorithm above depends on the intrinsic value of the option. However, the call options with very short maturity or out-of-the money options have no intrinsic value. This causes the integrand in the Fourier inversion to be high oscillate. As shown in Carr and Madan [28] and Zhang and Wang [6, 7], an FFT based on the time value of the option has to be obtained.
Let be the time value of an out-of-the-money option; then,
Similar to (22), the Fourier transform of the modified time value is given bywhere . So, the time value, , can be calculated by the inverse Fourier transform:
Similar to (31), the FFT algorithm for the out-of-the-money options iswhere .
Remark 1. In practice, one may determine an upper bound on the damping parameter from the analytical expression for the characteristic function and the condition . Carr and Madan [28] found that one-fourth of this upper bound serves as a good choice for ( for calls, for puts) and the impact of the damping factor on the option price is not significant.
3. Numerical Experiments
In this section, we present some numerical results. First, we report the numerical experiments that were performed on pricing options using two techniques: the FFT and the Monte Carlo (MC) simulation method. Next, we use the FFT to further investigate the volatility smile and the distributional properties implied by our proposed model (TSVIJ). We analyze further the effects of volatility jump components and correlation coefficients on option prices under the TSVIJ model. Finally, we describe the market data used for empirical study and estimate the risk-neutral model parameters from option prices on a set of days in various models including the TSVIJ, the SVIJ, SV, Merton, and the BS models as benchmarks.
3.1. Performance of FFT
Having presented the FFT algorithm as outlined above, it would be interesting to investigate the performance of FFT implemented to price option under our proposed model. We have also implemented the MC method for comparison purposes. For the FFT method, according to Remark 1 and by means of the given setting of model parameters in Table 1, we used the damping parameter , , and . For the MC simulation, we apply the Euler scheme for the numerical simulation of the log-asset price and variance processes. In order to achieve high accuracy, the simulation is repeated 10,000 times and the time step is kept at 0.02. Other parameter values used in the numerical computation are listed in Table 1. Two methods were implemented in MATLAB R2015b and on the PC with Intel(R) Pentium(R) D CPU 2.80 GHz. Results of the numerical experiment are presented in Table 2.
Table 2 reports the results for 9 call option prices with different strike prices (from deep-in-the-money (ITM) to at-the-money (ATM) and deep-out-of-money (OTM)) obtained using the FFT method and the Monte Carlo simulations. Column one lists the strike price for the options. Columns two to five list the FFT and the MC simulation prices for both and . The numbers in parentheses are difference between the two methods. If we use the Monte Carlo simulation as the benchmark, we can see from Table 2 that the absolute percentage differences of the FFT method are less than for all cases.
The FFT method takes less than 2 seconds to produce 150 option prices corresponding to different strike prices. The Monte Carlo simulation takes more than 90 seconds for each option price. These numerical experiments show that our analytic solution is correct, and the FFT method is accurate and efficient.
3.2. Implied Volatility
The clear advantage of the FFT method enables us to further investigate the volatility smile implied by our proposed model. For the implied volatility calculations, we use two terms to maturity date: 0.5 year (short-term) and one year (long-term). We use the model parameters taken from Table 1 and vary the strike price from 50 to 150, to generate an implied volatility graph. Figure 1 displays the implied volatility curves for three different models, including the Merton, SVIJ, and TSVIJ models.
(a)
(b)
As can be seen from Figure 1, for both maturities, the curvature of the implied volatility curves for the TSVIJ model are significantly decreased, and the difference between the Merton, SVIJ, and TSVIJ implied volatility curves is quite large. This should not be a surprise, as jumps in two variance processes increase the conditional variance of the underlying asset. Furthermore, the TSVIJ model generates very greater the left-tail skewness, especially for short maturity, which should be valid for the real-world option market and the best framework describing asset price dynamics. This indicates the advantage of the two-factor stochastic volatility model with jumps, with respect to the jumps only in returns (the Merton model) and the single-factor stochastic volatility model with/without jumps (the SVIJ/SV model), is that the two-factor stochastic volatility model with jumps in volatility provides more flexibility to capture the volatility smile since one of the factors determines the correlation between short-term asset’s returns and variance, whereas the other factor determines the correlation between long-term asset’s returns and variance.
Next, we examine how the volatility smile varies with the model parameter values of the variance process. We focus on investigating the effects of three parameters, two volatility jump components, and , and correlation coefficient, , on the volatility smile. We compute the prices for call option with year maturity by FFT under the TSVIJ model; then, we take those prices as market prices of option and compute the implied volatilities using the BS model. The results are plotted in Figures 2–4.
First, we can seen from Figure 2 that an increase in , the arrive intensity of the volatility jumps, increases the implied volatility. This implies that the BS model systematically undervalues all options. A similar argument holds for , the mean of the volatility jumps (see Figure 3). Second, the levels of the implied volatilities depend on the sign and magnitudes of correlation coefficient between the volatility and asset (see Figure 4), whereas the levels of the implied volatilities at the lower strikes drop with increasing correlation and ones at higher strikes go up. This feature is consistent with the models with stochastic volatility without jumps in volatility. Thus, the proposed model maintains some properties of the models being of stochastic volatility. We stress, however, that this model enables option prices to fit the forward price term structure exactly.
3.3. Probability Densities
We now go one step further and examine the distributional properties of the TSVIJ model using the FFT method. It is also instructive to compare the short-term and long-term probability density functions of TSVIJ, SVIJ, and Merton’s models. From (21), the FFT algorithm of can be approximated byfor . The results are plotted in Figures 5–8, which examine the values of the probability density function for a range of values of two parameters when all the others are kept constant and taken from Table 1.
(a)
(b)
Figure 5 shows how the positive correlations of variance process with the stock return (which is calculated in equation (3)) create a fat-right tail and thin-left tail distribution of the log-return. Positive correlations lead to a high variance when the stock return increases, which fattens the right tail of the density function. This is consistent with those in Wong and Zhao [30].
Figures 6 and 7 provide the impacts of the volatility jump components, , and in the TSVIJ model on probability densities of the log-return. Figure 6 shows that increasing the arrive intensities of volatility jumps, and , will decrease the kurtosis only, but create a fat-right-tail and thin-left-tail. A similar argument holds for the means of the volatility jumps, and (see Figure 7).
Figure 8 provides comparisons of the short-term ( year) and long-term ( year) probability densities of the TSVIJ, SVIJ, and Merton’s models. Figure 8 shows that the probability density of the TSVIJ model is distinctly different; it has very high peak and skewness features compared to those of the SVIJ and Merton models for both maturities. Moreover, for both maturity density curves in Figure 8 still show that the density of the TSVIJ model attenuates heavy tails on both sides. Therefore, the TSVIJ model can correct the pricing bias in the BS model and has a potential application in practice.
3.4. Sensitivity Analysis
As said previously, the TSVIJ model provides more flexibility than the single-factor specification to model the volatility term structure and probability density of asset returns. This section considers the effects of the model parameters, including the correlation coefficients and , the arrive intensities and of volatility jumps, and the means and of volatility jumps, on the option prices. We analyze the prices of T = 1 year call option only. Model parameter values are taken from Table 1 apart from , and . The numerical results are displayed in Table 3.
From Table 3, we can see that (i) option prices are very sensible to and . Since and indicate the means of the volatility jumps, this sensitivity can be considered as the sensitivity of option prices to their volatilities, which are verified by equation (2). A similar argument holds for the arrive intensities and of the volatility jumps. (ii) Option prices with different strike prices have different sensitivities to the correlations and , whereas in-the-money options () decrease in values with increasing correlations and , at-the-money options (), and out-the-money options () go up. This finding is consistent with those analyzed by Zhang and Wang [6] and among others.
3.5. Calibration of the TSVIJ Model to Market Data
With a successful implementation of the option pricing formula (15) under the TSVIJ model, we can apply the above calculation method to calibrate the TSVIJ model and estimate the model parameters on a set of market data of standard European-style call options by minimizing the difference between market prices and model prices in a least squared error fit.
The dataset used in this study consists of AAPL plain vanilla call option prices traded on the NASDAQ. These option prices are based on the Apple stock without paying dividends, having a value of $ 364.11 on July 5, 2020. We employ the market quotes over the period from July 10, 2020, to July 24, 2020, as the in-sample data to calibrate the model parameters, and use those of July 31, 2020, for the out-of-sample test. The riskless dollar interest rate , in this period. To ensure sufficient liquidity in option trading, we disregard option trading volume lower than 100 in the sample data. Finally, the option sample contains 128 call options, with 4 trading days and maturies from 5 days to 26 days. The description of AAPL options data are reported in Table 4.
The TSVIJ model has 17 parameters that need estimation, viz., . The estimation procedure is to minimize the least squared error over the parameter space , i.e., we evaluatewhere is the number of the sample data and and denote the market price and the model price at time of a call option, respectively, with maturity and strike price . The model price depending on the parameter space is taken as input and calculated by the FFT method above. As the penalty function, we choose the distance to the initial parameter vector (see Mikhailov and Nögel [32]):
We employ the Lsqnonlin function (least-squares nonlinear) in MATLAB to implement (38). To ensure the success of the local optimization scheme (38), we need effective choices of the initial parameters . Based on reasonable results from the literature, we set the initial parameters used the empirical results studied by Duffie et al. [12]. In addition, we also define lower and upper bounds for the optimal parameters which avoid possible solutions that, while mathematically feasible, are not acceptable in finance sense in the calibration. The resulting calibration results are reported in Table 5.
Table 5 gives the estimated risk-neutral parameters and pricing performances for the TSVIJ, the SVIJ, SV, Merton, and the BS models using the in-sample data. For comparative purposes, we provide the root mean square error (RMSE), the average relative percentage error (ARPE), and the CPU time. Considering the RMSE and ARPE values, the TSVIJ model best fits the market prices in the in-sample data, followed by the SV, SVIJ, Merton, and the BS models. Taking a closer look at the estimates, the jump-related parameter in the TSVIJ model is significant, suggesting that modeling the asset dynamics and its variance processes cannot ignore the jump components.
Using these results in Table 5, the model predicted values and its comparison with market prices of the out-of sample data are shown in Table 6. As the table shows, the TSVIJ model provides a good fit for most traded options in a lower RMSE value.
4. Conclusion
This paper proposes a multifactor variance processes model (TSVIJ model) for European option pricing, with jumps in both asset’s returns and variances, which is a reasonable and natural extension of the two-factor stochastic volatility model, Merton’s jump-diffusion model, and the affine jump-diffusion model. Based on the martingale method and Fourier transform technique, we also obtain closed-form solutions for the characteristic function and European call option. Using FFT, we obtain accurate and efficient numerical algorithm for pricing option within this proposed model. Comparisons of both the implied volatility curves and probability densities under alternative models show that the proposed model has better pricing performance. The analysis of the proposed model empirically shows its ability to handle pricing problem much better than the SVIJ, SV, Merton, and the BS models. Once the model parameters have been calibrated to fit market prices, the proposed model can be used to price other products that are not actively traded in the market. An analysis of this model also reveals that the volatility jump components have a strong impact on option prices. Although the independent arriving jumps in volatilities and jumps in returns are considered in the proposed model, we recognize that it could also be extended to jumps occurring simultaneously in both returns and variances by slightly modifying the framework of this model.
Appendix
A. The proof of Lemma
From the conditional characteristic function of in (6), then for an arbitrary date , we obtain
Note that, since , . The law of iterated expectations implies thatwhich shows that is a martingale. The martingale property of implies that
Applying the generalized Itô lemma to , we obtain
Taking conditional expectation with respect to on both sides of (A.4) under measure and using (A.3), it is straightforward to show that the function satisfies the following partial differential-integral equation:with terminal condition
For the dynamics of processes (2), (3), and (6), the coefficients are all affine structure in nature. It can be shown that the characteristic function has the exponential affine form:where , and are deterministic function of , and . From (A.6), it is clear that
Substituting (A.7) into (A.5), we get the following system of ODEs for , and (we suppress the last two arguments):subject to boundary conditions (A.8).
We now consider equation (A.9). Letting yieldswith the initial condition . Integrating both sides over in (A.12), we obtain
We know
Following the formula above and invoking the condition yieldswhich is specified in (7). The same arguments mentioned above can be used to equation (A.10) for the solution of . With and available, from equation (A.11), we obtain
Hence, the proof of the lemma is completed.
Data Availability
All data used to support of the findings for this study are included within this article.
Conflicts of Interest
The author declares that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the NSF of China (Grant no. 11461008) and Guangxi Natural Science Foundation (Grant no. 2018GXNSFAA281016).