Abstract

The coupled additive and multiplicative (CAM) noises model is a stochastic volatility model for derivative pricing. Unlike the other stochastic volatility models in the literature, the CAM model uses two Brownian motions, one multiplicative and one additive, to model the volatility process. We provide empirical evidence that suggests a nontrivial relationship between the kurtosis and skewness of asset prices and that the CAM model is able to capture this relationship, whereas the traditional stochastic volatility models cannot. We introduce a control variate method and Monte Carlo estimators for some of the sensitivities (Greeks) of the model. We also derive an approximation for the characteristic function of the model.

1. Introduction

The classical Black-Scholes-Merton (BSM) model [1] for derivative pricing assumes that the stock price, , has lognormal distribution, with the following governing stochastic differential equation (SDE):where the constants and are the drift and volatility and is a standard Brownian motion. There are several generalizations of the BSM model that relax various model assumptions. One of these assumptions is that the volatility is constant. It has been long known that this is not supported empirically. Two main generalizations of the constant volatility assumption are given by the local volatility and stochastic volatility models. The local volatility models write volatility as a deterministic function of time and stock price, whereas the stochastic volatility models describe the behavior of volatility by another stochastic differential equation. A comprehensive treatment of stochastic volatility models can be found in [2].

In most stochastic volatility models, (1) is replaced bywhere the volatility, , is now a function of another stochastic process , which follows the SDE:Here, the standard Brownian motion is typically assumed to be correlated with the Brownian motion of the asset price equation (2). The specific forms of the functions , vary between different models. For example, Table 1 lists some stochastic volatility models where each model is determined by the choices for , , ,   ( is the instantaneous correlation coefficient of and ).

The Hull-White [3] model uses a lognormal process for the volatility, the Scott [4] and Stein-Stein [5] models use a mean-reverting Ornstein-Uhlenbeck process, and the Heston model [6] uses the Feller (Cox-Ingersoll-Ross) process. In the table, the constant is the rate of the mean reversion, and is the long term mean of . The constant is sometimes called the volatility of the volatility process.

A generalization of these stochastic models, called the multiscale stochastic volatility model, was introduced by Fouque et al. [7] and further studied by Fouque and Han [8]. These models take the following form:The volatility is driven by two stochastic factors: a fast mean-reverting diffusion process and a slowly varying diffusion process . The different scale parameters are given by and . The standard Brownian motions , , have a general correlation structure.

In this paper, we consider a new class of stochastic volatility models, called the “coupled additive and multiplicative” (CAM) noises model. The CAM model writes using two Brownian terms, one additive and one multiplicative. Together with the SDE for the stock price, the model takes the following form:In the rest of the paper, except for the last section, we let ; that is, the logarithm of the volatility follows the process in (7).

The CAM model was first used by Anteneodo and Riera [9] in derivative pricing. The authors computed the density function for the steady state solution of (7) and showed that the density function matches the empirical density of the log-volatility of their stock data very well. The model (7) was used by Sura and Sardeshmukh [10] to model the daily sea surface temperature variations. The authors gave empirical evidence suggesting that the relationship between the kurtosis and skewness of the variations was best explained by the CAM model. Amusan [11] introduced and compared various methods to estimate the parameters of the CAM model from data, compared the empirical and theoretical density functions for stock price data, and compared European option prices when the stock price was modeled using the CAM model and the Ornstein-Uhlenbeck process.

Let us investigate the empirical analysis given by Sura and Sardeshmukh [10] further. Let and denote the “excess” kurtosis (kurtosis of the distribution minus 3; the Gaussian kurtosis) and skewness of . The authors proved that where , are constants. Then they derived bounds on the constants, arguing for weak multiplicative noise, , and in general and . Does empirical financial data support such a relationship between the excess kurtosis and skewness of log-volatility? To answer this, we compute the logarithm of the historical volatility for the CBOE S&P 500 Buy Write Index, for June 30, 2002, to July 31, 2012. Figure 1 gives a scatter plot of the excess kurtosis of the log-volatility data, against its skewness, and plots the parabola . Each point in the scatter plot corresponds to of the log-volatility of daily price of 6 months’ data. Interestingly, we observe an excess kurtosis-skewness relationship , similar to what Sura and Sardeshmukh [10] observed in their ocean temperature data.

Stochastic volatility models with one noise term, as summarized in Table 1, cannot explain the relationship between kurtosis and skewness that was observed in Figure 1. The logarithm of the volatility in the Scott and Stein-Stein models is Gaussian and thus has constant kurtosis and skewness of three and zero. In the Hull-White and Heston models, the kurtosis and skewness of the stationary distribution does not exist. To the best of our knowledge, there are no analytical results on the kurtosis and skewness of the multiscale stochastic volatility model of Fouque et al. [7], given by (4). We use Monte Carlo simulation to compute 95% confidence intervals for the kurtosis and skewness of the distribution of the volatility and obtain (−0.008, −0.0076) and (−0.0462, −0.046), respectively. Clearly these results cannot explain the empirical financial evidence observed in Figure 1.

To get a qualitative understanding of how the density function of the stock price looks like under the CAM model, we plot in Figure 2 the density function for under the CAM model (solid curve) and the Black-Scholes model (dotted curve). The CAM density plot is computed by Monte Carlo simulation. The constant volatility for the Black-Scholes model is taken as . The parameters used in the CAM model are , , , , , , , and . We observe that the tails of the CAM density are fatter and asymmetric.

The empirical evidence we discussed, as well as the evidence given in [9, 11], encourages the further analysis of the CAM model. In this paper, we develop efficient Monte Carlo estimators for the CAM model and derive estimators for its Greeks. In addition, we derive an approximation to the characteristic function of the CAM model, which enables the efficient calculation of its density function and option prices.

The rest of the paper is organized as follows. In Section 2, we develop a control variate method for European call option pricing under the CAM model. This method is based on a control variate method developed by Fouque and Han [8] for the multiscale stochastic volatility model. We compare randomized QMC methods with MC and the control variate method with the crude MC implementation. In Section 3, we introduce Monte Carlo estimators for estimating the sensitivities (Greeks) of option prices under the CAM model, which are of fundamental importance for hedging strategies. In Section 4, we use fast Fourier transform to compute option prices and present numerical results.

2. Martingale Control Variate Method

Fouque and Han [8] introduced a control variate method (called the “martingale control variate method”) to price options under the multiscale stochastic volatility model given by (4). In this section, we develop a similar control variate method for the CAM model and present numerical results for European options with MC and RQMC implementations.

We start with the Euler discretization of the CAM model. The model has a general correlation structure for the Brownian motions , , , with , , , and for . We rewrite these processes in terms of three independent standard Brownian motions , , :

Then, the CAM model can be rewritten in the following form:

We first derive the unique probability measure , equivalent to , such that the discounted stock price is a martingale under the new measure . Then, one can compute European style option prices with expiry and payoff , using the formula .

To this end, define From the multidimensional Girsanov theorem (see [12]), under certain integrability conditions, and are independent standard Brownian motions under the (risk-neutral) measure defined by

We choose as . The processes , , can be any adapted processes, under the aforementioned integrability condition. Then, under this new risk-neutral measure, the SDEs (9) becomewhere The functions , are given as The three Brownian motions , , and remain independent under the new measure. In our numerical results, we set for simplicity.

We now rescale the various parameters in terms of a parameter : , , , , and . Our equations become where

Under the CAM model, the price of a European type option with the integrable payoff function and expiry is given bywhere denotes the expectation with respect to , conditioned on the current states , . Crude Monte Carlo simulation estimates the option price at time 0 by where is the total number of independent sample paths and denotes the th simulated stock price at time .

Similar to the derivation given by Fouque and Han [8], we apply Itô’s lemma to the discounted option price (assuming the necessary smoothness), then integrate the process from time to the maturity , and obtain the following martingale representation:where the martingales are defined by Following [8], we replace in the above equations, by an approximate option price. A good choice is the Black-Scholes option price , where is the volatility at the mean-reverting value ; that is, . The martingale control variate estimator is defined as

In the simple case when the instant correlation coefficients , , and are zero, the variance of the control variate estimator can be computed as follows: In particular, it is possible to show that the variance of the control variate estimator is bounded. The following result was proved by Fouque and Han [8] for the multiscale stochastic volatility model and Huang [13] for the CAM model.

Theorem 1. Assume the payoff function is continuous and piecewise smooth, and the discounted European style option price is as a function of the parameters. Then, for any fixed initial state, there exists a constant such that, for ,

The martingale control variate method can offer significant error reduction, especially when combined with randomized quasi-Monte Carlo. Figure 3 plots European call option estimates using crude Monte Carlo and the martingale control method, against the sample size . The plot on the left uses a pseudorandom sequence (Mersenne twister, [14]), and the plot on the right uses a random-start randomly permuted Halton sequence (Rasrap; [15, 16]). Table 2 displays the sample standard deviation of 50 option estimates when pseudorandom and rasrap sequences are used with the crude Monte Carlo (labeled as MC and RQMC) and with the martingale control (labeled as CV-MC and CV-RQMC) implementations. Three different sample sizes are considered: 20K, 60K, and 100K. The last column gives the ratio of the sample variance of CV-MC to CV-RQMC. In these numerical results, the following parameters are used: , , , , , , , , , , , and . The discretization uses 4 time steps. We choose the function in (5) through (7) as in the numerical results. The numerical results indicate significant advantages of the control variate method over the crude implementation and RQMC over MC.

3. Estimating Option Sensitivities

The sensitivities of an option (Greeks) are important for practitioners for hedging. This information is used to manage the risk of an asset. There are two main “direct” Monte Carlo methods to estimate sensitivities: pathwise and likelihood ratio method. In general, these methods give unbiased estimators for the sensitivities. An extensive treatment of the Monte Carlo estimation of sensitivities is given by Glasserman [17].

Consider the price of a European call option: where is the expiry, is the risk-free interest rate, and is the exercise price. The expectation is taken under the risk-neutral measure. The stock price process, , follows the CAM model. Discretize the time interval , uniformly, as , and . In this section, we derive pathwise estimators for the Greeks delta and rho, which are defined as the derivatives and , respectively. Since we are working with a discretization of the SDE, these estimators are not unbiased. The pathwise estimator of delta is given by the following theorem.

Theorem 2. The pathwise estimator of delta for the European call option is given by

Proof. Sufficient conditions for the existence of the pathwise derivative and the interchangeability of the differentiation and expectation operations are given in Glasserman (page 394, [17]). We first verify these four conditions. Consider the Euler discretization of the CAM model:(i) From (25), we have for each . The solution of this recursion implies that, at each , exists with probability 1. (ii) We have ; that is, points of nondifferentiability occur with probability zero. (iii) The payoff function is Lipschitz continuous with respect to . (iv) The last condition holds with .
Having verified the existence and interchangeability conditions, we derive the pathwise estimator of delta for European call options. From the chain rule, Differentiating both sides of (25) with respect to , we getContinuing recursively, we obtain (), . By induction, we conclude .

The next theorem gives the pathwise estimator for rho, the derivative of the option price with respect to the risk-free interest rate.

Theorem 3. The pathwise estimator of rho for the European call option is given by where can be computed recursively as and .

Proof. The existence of pathwise derivative and the justification of interchange of operations can be shown similar to the proof of Theorem 2. From the chain ruleWe differentiate both sides of (25) with respect to :One more step in the recursion givesand the last step of the recursion gives .

4. Characteristic Function and the Discrete Fourier Transform

In this section, we derive an approximation to the characteristic functionfor the CAM model. Once we have the characteristic function, we can compute the probability density function and option prices efficiently using the fast Fourier transform approach introduced by Carr and Madan [18]. We will work in the risk-neutral measure and thus take , where is the risk-free interest rate, in (5) of the CAM model.

We apply a log-transformation for the asset, that is, put , and write the SDE for the new variable:and let follow the process in (7); that is,Assume that the discounted characteristic function has the following form: where represents the time to maturity of a European style option and   and   are functions to be determined.

Consider a general conditional expectation function, at least twice differentiable, of the formFrom Itô’s lemma

By construction, is a martingale, and thus . Then, the deterministic part of (39) must be zero, which yields the Fokker-Planck forward equation:

We make the following approximations to linearize the coefficients of this equation so that an analytical solution is possible. This approach is based on [19]. To this end, we let and replace by its linear approximation at the mean-reverting volatility ; that is, . We also choose to linearize some of the terms (although it is a priori possible that -values could be negative in the original SDE, these values are far from the focus of our linearization, and ignoring this issue introduces no worse errors than the linearization itself). The partial differential equation (40) is simplified as

We next substitute the function (37) for in (41), which reduces the partial differential equation into two ordinary differential equations: with initial conditions , The solution of the system of ordinary differential equations is (details can be found in Huang [13]) where

The characteristic function can be used, together with a fast Fourier transform approach introduced by Carr and Madan [18], to compute the density function and option prices efficiently. In the following numerical results, we use this fast Fourier transform approach; details can be found in Huang [13]. To get a qualitative understanding of how the density function of the stock price looks like under the CAM model, we plot in Figure 4(a) the density, together with the lognormal density for the Black-Scholes model (the dotted curve), where the constant volatility is taken as . The parameters used in the CAM model are , , , , , , and . We observe that the tails of the CAM density are fatter and asymmetric. Figure 4(b) shows a close-up of the density function at its peak, while the correlation coefficient of the multiplicative and additive noises takes various values. We observe that the density function concentrates more at its mean as increases.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.