A JUMP-DIFFUSION WITH STOCHASTIC VOLATILITY AND INTEREST RATE

In this study, we present the application of Time C hanged Levy method to model a jump-diffusion process with stochastic volatility and stochastic i nterest rate. We apply the Lewis Fourier transform method as well as the risk neutral expectation pric ing method to derive a formula for a European option pricing. These combining methods give quite a short route to derive the formula and make it efficient to compute option prices. We also show th e calibration of our model to the real market with global and local optimization algorithms.


INTRODUCTION
The success of Black and Scholes (BS) model rests on the ease of computation and traceability: it has a closed form solution and allows for dynamic hedging. However, the BS model fails to explain many aspects of the real distribution of an asset return. Since then, there have been continued efforts to make a better model to describe a model in the real world starting from Merton's jump-diffusion, Heston's stochastic volatility and Bates' stochastic volatility with a jump-diffusion. While the search for better models for dynamics of asset prices continues, the practitioner's community usually focuses more on the utilizability of the models to be just as important as how the models can describe the dynamic of an asset in the real world. The issues are on the computation, calibration and traceability of the model.
In this study, we will apply the Time Changed Levy method to model a jump-diffusion with stochastic volatility and stochastic interest rate which is quite different from the typical method that describes the dynamic of the model by separated components of the stochastic factors. The Time Changed Levy method, however creates the stochastic volatility by changing the time of a pure Levy process to a random stochastic time.
In the computation part, the Lewis' Fourier Transform method was used to calculate the option prices on the complex plane which works seamlessly with the Time Changed Levy method where the measure changed is defined in the complex domain. The option formula by the Lewis (2001) method comes out in a single Fourier integral form which helps to reduce the computation time compared with the other methods which normally generate two integrals.
We also calibrate the model to the real market prices using both global and local search methods to find the minimum of the discrepancies between the market prices and model prices to obtain the optimal parameters of the model. The result is reported in the last section.

Typical Model
The typical risk neutral model for jump-diffusion with stochastic volatility and stochastic interest rate can be described by Equation 1: Here W t , v t W and r t W are the Brownian motions associated to the underlying asset process, the variance process and the interest rate process respectively. The process S t is the underlying asset price process and r t is the instantaneous risk free rate process. The process N i is a Poisson process with jump frequency λ and independent from the other processes. The jump size Y t -1 is distributed as described above. The jump is included in the model to make short term implied volatility curve steep as indicated by the empirical studies. The process v t is the variance process, k is the speed of mean reversion, θ is the mean of long term variance and σ is the volatility of the variance process. The variance process is the square root process known as CIR process. To explain the leverage effect, the negative correlation is usually introduced between the underlying asset and the variance as shown above. The interest rate process is also a CIR process but with different parameters and independent from the other processes.

The Time Changed Levy Approach
We derive the Fourier option pricing using the Lewis method. There is a variety of Fourier Transform Pricing methods but we choose this method as its integration domain is on the complex plane. This complex domain will correspond to the domain for the time changed Levy process. Another nice feature of this method is that it produces a formula in a single integration form compared with the typical approaches which produce two integrations such as the approach in Sattayatham and Intarasit (2011). This single integral reduces computation time of option prices in the calibration process. During our calculation, we also apply the Modular approach, pioneered by Zhu (2010), which employs the rule of independence of characteristic functions to write the characteristic function as product of each characteristic function of an independent stochastic factor. This approach will help us to handle each stochastic factor independently which results in the reduction in the dimensions of the calculation.
Our dynamic of stock price will be an exponential Levy process which is driven as Equation 2: Here S t is the stock price at time t, S 0 = exp(X 0 ) is the price of stock at time t = 0, r t is a risk free rate process and L t is a Levy process with exp(L t ) being a martingale. Let us assume the characteristic function of L t , where α and β are real numbers and z is a complex number. The Fourier transform of the payoff for a European call option strike at K with a payoff function at the maturity, max(S T -K) or (S T -K) + can be computed as: Here exp (X T ) = S T or , Ĥ (z) is defined in the region where the imaginary part of Fourier transform variable z, is greater than 1. The corresponding generalized inverse Fourier transform H(x) for the payoff function is defined below: In Equation 3 and Equation 4, we extend the transform variable z to take a value in the complex domain that is defined in the generalized Fourier transform sense. Given, Ĥ (z) is well defined on the plane where the imaginary part of z is greater than 1, the integration in Equation 4 is just the line integration on a complex plane paralleled to the real axis with any Im(z)>1. From the Fundamental Theorem of Asset Pricing, the no arbitrage condition is equivalent to the existence of a risk neutral measure where a discounted asset price is a martingale. Based on this Fundamental Theorem, we can write the value of European call option at t = 0 as risk neutral expectation of the discounted payoff as:

JMSS
The expectation E Q is the expectation under a risk neutral measure. The third line is derived from the second line by replacing the payoff function with the corresponding generalized Fourier transform in Equation 4. In Equation 5, here we suppose the interest rate process is independent from the other processes, therefore we can write the expectation out from the other terms.
Here we will apply the Time Changed Levy method to generate a model with stochastic volatility and stochastic interest rate as in Equation 1. That is our model is driven by: where, r t , y t , λ k y and v t are defined as in Equation 1 and t T X is the time changed Levy process which will be defined later. Compared with Equation 2, the Levy part is Equation 7: The term compound Poisson process. This time changed Levy process t T X is constructed by two stochastic processes, a subordinator and an underlying process. The subordinator, as an increasing process, is a function of calendar time t to stochastic time T t . The underlying process is normally a pure Levy process. For our case, the subordinator T t is now defined to be a process Equation 8: where, v s is defined as in Equation 1. The underlying Levy process is the risk neutral Brownian motion with drift rate r t Equation 9: t t t t dS S (r dt dW ) = + (9) whose the log return can be described by: As mentioned previously, the time changed Levy process By assumption of the independence between the time changed Levy process and the jump process, we may write the characteristic function of Levy process T L ( z) φ − from Equation 5 as the product of the characteristic function of the time changed Levy process and the characteristic function of the compensated compound jump process Equation 12: as the characteristic functions of the time changed Levy process and of the compensated Poison process respectively. So we now need to calculate each component of Equation 5 that are The characteristic of the time changed Levy process is derived as: From the second line to the third line above, we apply the measure change defined as the complex value Radon-Nikodym derivative (the detail of this measure can be found in chapter 8 of Zhu (2010) where Equation 14: and X (z) ψ denotes the characteristic exponent of the underlying process t T X . This measure allows us to write the characteristic function of the correlated process (in our case, the underlying process is designed to correlate with the variance) as the Laplace transform under the new measure. That is Equation 15: Here the characteristic exponent of the Levy process ψ x (z) of the process ( t  (2010)). The dynamic of v t is: The variable vM t W is the Brownian motion associated to the variance process under this new measure. Here is the derivation of the above equation. Based on Girsanov Theorem, given a measurable space (Ω,F,Q), an Ito process for the dynamic of v t : Denote M t as a exponential martingale under measure Q defined by: Then M T can be expressed as: By the assumption, the underlying process is correlated with the variance process or v t t dW dW dt = ρ .
Then we have the following results: • M t defines the Radon Nikodym derivative. That is: Then we can solve for the characteristic of the time changed Levy process as: Here is the derivation of The equation for D(t) is called Riccati equation which is a nonlinear ordinary differential equation. To make it simple we will write the equation for D(t) as: The solution for D(t) will be 1 u ' D(t) R u = − where u will satisfy the following auxiliary differential equation: The general solution for the above equation is: The last part is the characteristic of the compensated compound Poison process which is computed below Equation 19:

Calibration
Calibration is the process to obtain a model's parameters that match to the market prices of the options. This model's parameters generally will be different from parameters estimated by the statistical methods. The statistical parameters reflect the past dynamic of the underlying asset and not guarantee that the models built out of this parameters are arbitrage free. Contrary to the statistical method, the parameters from the calibration are arrived with the principles of no arbitrage at least for the traded options that included in the calibration process. The differences of these two methods reflect in the investors risk preferences, hedging costs and views of the participants in the market which is not be able to be captured by a statistical method.
Given we can observe the prices of the option from the market, the equation for the value of the option based on the risk-neutral valuation is Equation 20: Science Publications

JMSS
Here θ p is the set of the model parameters, T and K are the time to maturity and strike of the observed option respectively. The dynamic of S T is described as a parametric model under the risk neutral measure. If we can obtain the prices of options at any time T for all K, we could determine the parameters of the dynamic of S T by the above equation. But this is impossible in the real market where we have limited prices of option in any maturity, one possible way is to minimize the discrepancies between the available market prices and model prices generated from a parametric model. Therefore the problem of calibration has transformed to an optimization problem for the least square of the discrepancies.
The scenario is that the market prices consists of the price of European call options spanning a set of expiration dates T 1 ,…,T N and for each T i , the market quotes for strike K i1 ,…,K iM . The least square method to find the minimum of different of prices and the parameters of the model can be described by Equation 21: The function F(θ) is the objective function with parameter θ p . The function C M (θ p ,T i ,K ij ) and C(T i ,K ij ) are the value of the call option generated by parameters θ p and the observed price at maturity T i and strike K ij respectively. The variable w ij is the weight associated to the confidence of the observed price which varies with the vega of option. The detail of weighing can be found in chapter 4.2 of Poklewski-Koziell (2012). This least square problem will have at least a solution given the domains of parameters are compact.
To find the above minimum problem, we employ the simulated annealing method which is one of the most efficient method to find a global optimum. In this part, we will calibrate the models to the DAX index option prices on July 5, 2002 as shown in Table 1. Based on Equation 17-19, our model will have 12 parameters, that are v 0 , k , θ, σ , ρ , λ , u J , v J , r 0 ,α , ω and β.
We run the calibration algorithm in MATLAB optimization toolbox which provides both a simulated annealing algorithm and gradient based optimization algorithms. We compare our model with The Jump-Diffusion with Stochastic Volatility model (JDSV).

Calibrating of The Model
The simulated annealing method is the global optimization method that replicates the way the metal is heated to a suitable temperature and cooled down slowly to get the optimum structure. The suitable temperature is called the initial temperature T 0 and the way the temperature is reduced is called the annealing schedule. The algorithm is summarized below.
Step 1: Initialization Set the initial solution x 0 and the initial temperature T 0 which is high enough for the acceptance probability of 0.80 to 0.95.

Step 2: Perturbation
Generate the new solution x i according to the designed probability distribution and determine the difference of objective function ∆f = f(x i )-f(x i-1 ).

Step 3: Acceptance Determination
The new point is tested to accept if it falls on the following criterior: where, r is a uniformly distributed random number between [0,1]. The term exp(-K∆f/T) is referred to as the acceptance probability where K is normally equals 1.
Repeat step 2 and 3 until the equilibrium is reached (not much improvement for the objective function in this temperature).

Step 4: Annealing
The temperature is slowly reduced with a specific schedule to zero.
Repeat step 2 to step 4 until the objective function reach a specified goal.
The simulated annealing method differs from the gradient based methods in that it can avoid trapping in the local minimum by allowing the candidate point to be accepted even the objective function is worse than the existing objective function as shown in step 2. The accepting probability depends on the temperature which is high when the temperature is high and low in when the temperature is at low level.

Results
The simulated annealing tool, in the MATLAB's optimization toolbox 2012, provides a lot of options that cater to natures of the problem including the temperature setting, annealing schedule and stopping criterior. We

JMSS
have run the optimization the numbers of times and set the initial temperature at 400. For the random search, the Boltzmann generating function seems to be the best choice to obtain the minimum point to our problem. We compare our model (JDSVSI) with the Jump Diffusion with Stochastic Volatility (JDSV) with has 8 parameters, that are v 0 , k , θ, σ, ρ, λ, u J and v J . The data points that are the in-sample on this calibration are all the data in Table 1 except for the columns of 1 month and 18 months that we use as the out-of-sample points. Due to the nature of simulated annealing method, the search for the candidate points being random can locate only the neighborhood of the minimum. We need the gradient based method or call here "local search" to fine tune the minimum point.
The minimum of objective functions of JDSV and JDSVSI are just approximately 1.20 and 1.10 times higher compared to the BS's one percent error in implied volatility as shown in Table 2. These errors are consistent with the bid/offer spread in the real market which is just about a little higher than 1 percent. The errors are due to the factors that we do not recognize and do not include in the models, the data error and the landscape of the domain of the problem.
The parameters of the calibration of both models are shown in Table 3 and 4.     The Table 5 and 6 show the square discrepancies of both models' implied volatility with the corresponding BS implied volatility. The second last line presents the sum of the square error and the last line is the error adjusted by the total of square error for all columns. As expected, the JDSV is better for the short tenor but the JDSVSI is better in the longer tenor in accordance with the finding from Abudy and Izhakian (2011). The total errors of the out-of -sample data points that shows on the last line of the column of 1 month and l8 months of both tables are within the average 0.125 except for the error for 1 month of JDSVSI model.

CONCLUSION
The combining methods that we implement give a short route and straightforward option pricing formula compared with the existing ones which normally are derived by riskless portfolio partial differential equations or high dimensional risk neutral expectation method. From the calibration results, our model is better than JDSV for the longer tenors but not for the short one.