MPhil Thesis August 15, 2013
Option Pricing and Time-varying Jump Risk Premia
Andrei Lalu
MPhil Candidate Tinbergen InstituteSupervisors:
Prof. dr. Peter Boswijk
Prof. dr. Roger Laeven
Dept. of Quantitative Economics Dept. of Quantitative Economics University of Amsterdam University of Amsterdam
Abstract
This paper develops a semi-closed form option pricing approach in the context of a fully parametric model for asset returns with jumps. The model used is more flexible in accommo-dating asset price jump patterns than the classic time homogeneous Poisson jump diffusion. The stochastic jump intensity in the model self-excites as a result of jumps occurring, so the model can accommodate clustered jumps. The paper presents a way to estimate model parameters using joint time series of spot and option prices. We find evidence in favor of a time varying jump risk premium in a S&P 500 index data-set.
1
Introduction
Extant research in financial econometrics has documented the need to use both stochastic volatility and jumps when modeling asset returns. This has consequences not only for the way we conduct inference on time series of equity returns but also for option pricing models. In the turmoil of financial markets during 2008 asset price crashes occurred more frequently than a Poisson jump component driven model would suggest. Assuming investors attach a risk premium for return volatility as well as for stock price jumps we propose a parametric framework and an identification strategy aimed at disentangling the time-varying jump risk premium from the diffusive risk premium using time series of asset prices and option prices.
Large drops in asset prices occur more often than sensibly parametrized standard Brownian-driven statistical models would suggest. To account for these discontinuities in asset return series - a point process is incorporated in asset return models, and, following the seminal work of Merton (1976), the typical point process choice is the homogeneous Poisson process. Compared to a standard diffusion, a Poisson jump augmented diffusion process is better suited at modeling sparse asset return discontinuities. However it is often the case, especially during financial turmoil periods, that asset price crashes tend to cluster over short time spans of days or even hours. The standard Poisson diffusion could not justify such a repeated clustering of return jumps over time. Aït-Sahalia et al. (2013) introduce a novel way of modeling this complex interplay between jumps over time and across different asset markets by introducing self-exciting jump processes, known as Hawkes processes after Hawkes (1971). By modeling the point process component as a Hawkes process we aim to better accommodate asset return jump patterns in our parametric study of jump risk premia.
The study of risk premia constitutes another central topic in finance. Classical portfolio theory (Markowitz (1952)) is centered around explaining the equity risk premium, i.e. the compensation that rational investors require for return variation risk - variation in prices. Still, volatility patterns implied from option prices using the Black-Scholes model suggest that market participants are more ’afraid’ of the extreme negative moves than this model would predict. The distribution implied from options prices is negatively skewed whereas in the Black-Scholes model it should essentially be symmetric (Bates (2000)). Empirical financial econometrics studies provide a lot of evidence that asset return variance changes over time (e.g. Bollerslev et al. (1994)), and this variation naturally introduces an additional risk source from holding assets. Similarly, provided that a jump process is necessary to justify return distributions, the question of jump risk and its market pricing arises naturally. Furthermore, given that empirical events suggest the jump point process is better modeled as a time inhomogeneous point process, questions regarding the possible dynamics of this jump risk also arise.
The thesis intends to contribute to a vast literature aimed at explaining the observed dif-ferences between the time series behavior implied from option prices and the actual price of the underlying asset through the introduction of a jump risk component (e.g. Bates (2000), Pan (2002), Eraker (2004), Bollerslev and Todorov (2011)). These studies, out of which some
are based on parametric and some on non-parametric techniques, find evidence of a jump risk premium embedded in equity premiums. Most of the existing work is geared towards estab-lishing the existence and gauging the (relative) size of jump risk premiums. Some of these studies also employ state-dependent intensities for the jump process, however the intensity of the jump process is linked to other contemporaneous variables while our approach links the jump process to its own history. Therefore, this thesis aims to model the dynamics of the jump risk premium through parametric assumptions which allow for a more flexible dynamic of the premium compared to existing studies to date on the topic.
In order to be able to identify a rich parametric model and to make risk premium-related assumptions, we have to resort to using as much data on the underlying process as is available. Using the time series of the underlying together with the option prices allows us to devise an identification strategy for the model. However, this approach gives rise to another set of challenges given that option pricing is only feasible under a risk neutral measure. Unlike the Black-Scholes set-up, the presence of random jump sizes in each asset makes our set-up an incomplete market with respect to the universe composed of the underlying stock, the finite number of option contracts and the money market account, i.e. the state price density in our set-up is not unique. Under standard arbitrage free assumptions we choose a suitable candidate pricing kernel which prices jump risks and diffusive risks1. Under this assumption, we can use option prices to retrieve the parameter-dependent instantaneous volatility and jump intensity , i.e. the unobservable state variables, in a similar way in which volatility is implied from option prices in the Black-Scholes set-up. The implied intensity and volatility series can then be used together with the asset return series to estimate parameters in a straightforward econometric exercise (using GMM).
Using a data-set consisting of spot and option prices for the S&P 500 index covering the turbulent period around Lehman’s bankruptcy, we find statistically significant evidence in favour of accepting the hypothesis that investors demand a time varying jump risk premia as it can be reflected from option price data.
The remainder of the thesis develops as follows: Section 2 details the properties of the Hawkes self-exciting jump process. Section 3 elaborates on the model and the candidate pricing kernel and Section 4 continues with details about the approach we took for option pricing. Section 5 then outlines the parameter estimation procedure. Sections 6 is devoted to an empirical exercise on the S& P500. Lastly, Section 7 briefly concludes the thesis. Some technical details are omitted from the text, but are presented in appendices.
1Other risk factors such as interest rate uncertainty, dividend risk or liquidity concerns could potentially
have a very different impact on option prices. We only focus on the jump risks and diffuse risks to be able to investigate the differences and the dynamic behavior of the two. An alternative approach to this problem could explore preference-based equilibrium pricing such as in the seminal paper of Lucas (1978).
2
Self-exciting Jump Processes
In order to model the time-dynamics of the jump component, we consider using a Hawkes point process. The Hawkes process2 (N, λ) is a Markov process, where N is a counting process with intensity driven by λ. Similar to a Poisson process, this intensity describes the Ft conditional
mean jump rate per unit of time, i.e
P [Nt+∆− Nt= 0|Ft] = 1 − λt∆ + o(∆) P [Nt+∆− Nt= 1|Ft] = λt∆ + o(∆) P [Nt+∆− Nt> 1|Ft] = o(∆) (2.1)
The special feature of the Hawkes process however is its self-exciting property, which determines an interdependency between N and λ. This aspect is encompassed in the following stochastic integral equation for the intensity process:
λt= ¯λ + t Z −∞ δe−kλ(t−s) dN s (2.2)
The intensity of N is determined by λ and, in turn, λ also changes in response to a jump in the counting process N . If we interpret the model in the context of an asset jump, ¯λ > 0 represents
the intensity until a first jump occurs, whereas δe−kλ(t−s), with k
λ > δ ≥ 0 measures the impact
of a jump on the intensity process. The impact of the jump decays over time at an exponential rate with speed kλ. The compensated process Nt−
t
R −∞
λsds is a local martingale. As shown in
Appendix A.1, the intensity process (2.2) has the equivalent model formulation:
dλt= kλ(¯λ − λt) dt + δdNt (2.3)
Jump clustering can be captured by this model by choosing an appropriate size for the δ parameter. This intensity model prescribes a parsimonious way to incorporate the jump history in the time t intensity variable. Sample paths from a simulated Hawkes process illustrate the dynamics of the pair (N, λ) in Appendix A.2.
In the study in which they introduce the use of univariate and multivariate Hawkes processes3 for the modeling of the jump component in a continuous time model for asset returns, Aït-Sahalia et al. (2013) find econometric evidence in favour of such a specification with statistically significant self- and cross excitation jump effects for a panel of global equity indices. Errais et al. (2010) employ a multivariate specification of this point process to model clustered portfolio defaults. Hawkes processes are also used for the (high-frequency) modeling of trade book order arrivals and microstructure noise (e.g. Bacry et al. (2013)).
2
All throughout the thesis we refer to probability spaces of the type (Ω, F , P) where Ft is a complete infor-mation filtration satisfying all the usual conditions, see e.g. Shreve (2004).
3This type of point process has been employed in the academic literature to model earthquakes and their
3
The Model
The model describes univariate equity return dynamics and is an adaptation of the classic Bates (2000) model. While the model intrinsically features three risk sources, in our parametrization volatility risk is not priced4, we only focus on parametrizing the risk premia prices for the diffusive shocks to the asset return series and the jump risks.
3.1 The Model under P (DGP)
In a suitably defined probability space, assume that the stock price process, its volatility and the intensity of the point process have the following dynamics:
dSt= rt+ ηVt+ (µ − µQ)λt Stdt + p VtStdWt1+ ZtdNt− µλtStdt (3.1) dVt= kv( ¯V − Vt)dt + σv p Vt ρdWt1− q 1 − ρ2dW2 t (3.2) dλt= kλ(¯λ − λt)dt + δdNt (3.3)
This specification uses the familiar Heston (1993) model for stochastic volatility and the Hawkes point process to model the jumps. We assume here that the (continuos) risk free interest rate5
rt is a deterministic processes with respect to the filtration and that W1, W2 are standard
adapted Brownian motions.
This model captures two salient features of stock price dynamics: stochastic volatility and price jumps. Stochastic volatility is modeled through the process V defined in (3.2), which is also known as a CIR process from Cox et al. (1985). This volatility process is mean reverting towards a long run mean ¯V with a mean reversion rate kv. The specification renders the total
Brownian component in V to be correlated with the Brownian component in (3.1), via the correlation parameter ρ, capturing the leverage effect coined by Black (1976).
Some studies (e.g. Eraker et al. (2003), Eraker (2004)) justify the inclusion of a jump component in the volatility process, to better accommodate volatility dynamics, possibly even synchronising the jumps in volatility with the ones in the stock price process. We do not include this variation in the stochastic volatility equation out of tractability concerns in the estimation approach.
Pertaining to the jump component in our model, Z here denotes a serially independent random variable governing the jump size. Conditional on a jump arriving, the stock price jump is St= St−exp(Zt), where Ztis normally distributed with mean µj and standard deviation σj.
This implies that the mean relative jump size is µ = E [exp(Zt) − 1] = exp(µj+12σj2) − 1. The
jump times are determined by the counting process Nt with intensity λt, the pair forming the 4
Adding a volatility risk premium would be possible. As volatility is unobserved and is itself volatile, including such a premium could be justified. As we are mainly interested in the time dynamic of the jump risk premium, for parsimony and identification purposes, we do not pursue the addition of a volatility risk premium.
5
Possibly also considering a (continuous) dividend rate process e.g. qtwhich is also deterministic with respect to the filtration Ft
Hawkes process detailed in the previous section. The details about the terms in the drift of the stock price process (3.1) follow in the next sub-section after the candidate pricing kernel and the corresponding model under the risk neutral probability measure are detailed.
3.2 Candidate Pricing Kernel & The Model under Q
Given the model design (3.1)-(3.3), the market setting is rendered incomplete with respect to the underlying, the finite number of option contracts available and the money market account, so any pricing kernel in this set-up would not be unique. We restrict our attention to a candidate kernel which prices jump risks and diffusive price shocks and maintains the dynamics of the stochastic volatility and the intensity process in the affine model class as it is defined by Duffie et al. (2000). Complete details about the choice of measure change can be found in Appendix B. Letting Q denote the equivalent martingale measure associated with the chosen candidate pricing kernel, the stock price process is determined by the following dynamics:
dSt= rtStdt + p VtStdWt1,Q+ ZtQdNt− µQλtStdt (3.4) dVt= kv( ¯V − Vt)dt + σv p Vt ρdWt1,Q−q1 − ρ2dW2,Q t (3.5) dλt= kλ(¯λ − λt)dt + δdNt (3.6)
Here, (W1,Q, W2,Q) are standard Brownian motions under Q. The intensity process is the same as under P, and hence the counting process Nt is not affected by the measure change.
Conditional on a jump event taking place, the risk neutral mean relative jump size is µQ = EQ
h
exp(ZQ) − 1i= exp(µQj +12σ2j) − 1, with a different mean than under the physical measure. The last term in (3.4) is the corresponding compensator for the jump process under Q.
Market price of risks
Comparing the stock price drift in the specifications under each measure, the pricing of the risk factors becomes intuitive. The market price of diffusive risks (Brownian shocks) is pinned down by the parameter η in (3.1). It is akin to the risk-return trade-off in the CAPM framework, the larger η is, the higher the premium for diffusive risks is. It is time-varying, as the volatility of the process changes over time according to (3.2) and hence the premium ηVtchanges over time.
The jump risk premium is captured by the difference between the mean relative jump size under P, i.e. µ and the mean relative jump size under Q, i.e. µQ. This implies that the only premium investors attach to jumps is due to jump size uncertainty, while the timing of the jumps plays no role. This is a simplification we adopt for the sake of parsimony, as the sample identification of both a jump timing and a jump size risk premium would be challenging. Therefore the expected excess return which investors demand in exchange for bearing stock price jump risks is (µ − µQ)λt, which is time varying as λt changes over time.
4
Option Pricing
In order to conduct inference on the model parameters using option price data together with data on the underlying, we have to be able to price derivative contracts based on the model. To achieve this, we derive the discounted transform (characteristic function) of (log St, Vt, λt)T
-the state vector and invert it to find -the price of -the derivative contract at t = 0 assuming that the initial state vector (log S0, V0, λ0)T and the model parameters are known. The computation
of the transform can be carried out in closed form up to the solution of a system of ordinary differential equations. Its derivation in the following subsection is based on the results obtained in the general setting of affine jump diffusions developed in Duffie et al. (2000). The general approach for an affine jump diffusion process is summarized in Appendix C.1.
4.1 Affine Set-up and Transform Based Pricing
Denoting log St ≡ Yt, the system of stochastic differential equations (under Q) which defines
the state vector of our model, i.e. Xt≡ (log St, Vt, λt)T can be written in the following matrix
form: d YT Vt λt | {z } ≡Xt = rt− √ Vt 2 − µ Q λt kv( ¯V − Vt) kλ(¯λ − λt) | {z } ≡µ(Xt) dt + √ Vt 0 0 σvρ √ Vt σv p 1 − ρ2√V t 0 0 0 0 | {z } ≡σ(Xt) dWt1 dWt1 0 | {z } ≡dW3×1 t + ZQ 0 δ | {z } ≡Z3×1 t dNt (4.1) In the notation of Duffie et al. (2000), this system fits into the class of an affine jump diffusion and can be written as:
dXt= µ(Xt)dt + σ(Xt)dWt3×1+ Zt3×1· dNt (4.2)
Given the initial state vector X0 and the model parameters denoted by χ, a European option
contract with maturity T is priced under the risk neutral probability measure Q at time t as:6
C(d, c, T, χ) = Eχ exp − T Z t rsds ed·Xt− c+ (4.3)
This can be rewritten making use of the Fourier inversion of the discounted transform as:
C(d, c, T, χ) = Gd,−d(− log c : X0, T, χ) − cG0,−d(− log c : X0, T, χ) (4.4) where Ga,b(y : X0, T, χ) = ψχ(a, X0, 0, T ) 2 − 1 π ∞ Z 0 Im ψχ(a + ivb, X0, 0, T )e−ivy v dv (4.5) 6
and where the discounted transform is: ψχ(u, Xt, t, T ) = Eχ exp − T Z t rsds eu·Xt|Ft (4.6)
Proposition: The transform is equivalent to:
ψχ(u, Xt, t, T ) = eα(T −t)+β1(T −t)Yt+β2(T −t)Vt+β3(T −t)λt where α(T − t) and βi(T − t) are solutions to the system of ODE:
˙ β1 = 0 ˙ β2 = − β1 2 (1 − β1) − kvβ2+ σvρβ1β2+ 1 2σ 2 vβ22 ˙ β3 = −(eµ Q j+ σ2 j 2 − 1)β1− k λβ3+ eβ1µ Q j+ β2 1σ2j 2 +δβ3− 1 ˙ α = −r + rβ1+ kvV β¯ 2+ kλλβ¯ 3
subject to the initial conditions β = (u1, 0, 0)T and α = 0.
The proof of the proposition follows from the application of Proposition 1 in Appendix B of Duffie et al. (2000), using the system matrices detailed here in the thesis in Appendix C.3.
This system of ordinary differential equations can be solved using accessible numerical meth-ods (e.g. Runge-Kutta), a full analytical solution7 is not possible due to the non-linearities in
the differential equation involving ˙β3.
4.2 Benchmarking the Pricing Approach in a Monte Carlo Exercise
The option pricing approach outlined in the previous section involves finding the numerical solution of a system of ordinary differential equations as well as the evaluation of an indefinite integral, both steps potentially inducing pricing errors. In order to benchmark the results of the transform pricing approach we conducted a small Monte Carlo investigation by simulating paths of the stock process under the risk neutral measure, and computing the sample average of a discounted call contract payoff at maturity.
To simulate each stock price path we employed a (Euler) discretization scheme, so, based on a set of starting values (initial state vector) St0, Vt0, λt0 the discretized version of the model at
7The first equation is trivially solved by β
1(T − t) = u1, while a full analytical solution for the second one
a time point tj = tj−1+ ∆ based on its values at tj−1 is: Stj = Stj−1 + (rt− µ Q j λt)Stj−1∆t+ p VtStj−1∆W 1 t + Stj−1(e µQj+σj∗N (0,1)− 1)∆N tj Vtj = Vtj−1+ kv( ¯V − Vtj−1)∆t+ σv q Vtj ρ∆Wt1+ q 1 − ρ2∆W2 t λtj = λtj−1+ kλ(¯λ − λt)∆t+ δ∆Nt (4.7) One immediate problem with this scheme above is that Vt could become negative at some time
point with non-zero probability. To alleviate this problem we opted to truncate Vtj−1 by taking max(Vtj−1, 0) at each time point, a procedure referenced as ’absorption’ or ’full truncation’ in
the simulation literature (e.g. Lord et al. (2010) for more details).
The simulation approach also has its drawbacks compared to transform based pricing method, as it introduces discretization error into call prices. However, to keep this error to a minimum and the simulations parsimonious we opted for a fairly small a time8 increment ∆t = 10−3 and simulated N = 10000 distinct asset price paths to price each call option. Derivative prices (for European style options) were simply obtained by averaging the payoffs at maturity over all simulated paths and discounting this average. The parameter values under which both the simulation and the transform based pricings were done can be found in Appendix C.4.
Table 1 shows the average pricing difference between the Monte Carlo based price and the transform based price over 10 simulations for each different maturity and strike-price ratio.
T = 0.1 T = 0.25 T = 0.5 T = 1 K/S = 0.8 0.0469 -0.0343 0.1046 0.0109 K/S = 0.9 -0.0205 -0.0172 0.0837 0.0178 K/S = 1.0 -0.0255 0.0016 0.0644 0.0184 K/S = 1.1 0.0071 0.0207 0.0503 0.0164 K/S = 1.2 -0.0279 0.0229 0.0261 0.0294
Table 1: Average Pricing Differences
While there is some difference between the two pricing approaches, there seems to be no consistently large difference or a specific bias towards under or overpricing. The differences are likely due to the discretization error in the Monte Carlos which is probably too rough of an approximation for the jump process, with bigger consequences for contracts written with a longer maturity, where the differences seem slightly more pronounced than for the shorter maturity ones.
5
Model Estimation
In this section we outline an estimation strategy for the model parameter set which we denote as χ =µQj, σj, kv, ¯v, σv, ρ, kλ, ¯λ, δ, η, µj
. Note that the number of parameters is nχ= 11. The
8
estimation uses time series data on the underlying together with a panel of option prices. Besides the relatively large number of parameters, the biggest difficulty in estimation arises because two out of the three variables in the state vector are not observable, i.e. we only observe the stock price, while the volatility and intensity series are never observable.
One possible approach to identify the parameters is by minimizing the price differences between the actual observed call prices and the transform based theoretical ones by varying the parameters. This approach would identify parameters under the risk neutral measure, but wouldn’t give any insights about the risk premia in the model. Using the stock price dynamics together with option prices would instead allow for that. Assuming that we could somehow observe the full state vector at a series of (∆-spaced) time points t∆, t2∆, t3∆, . . . tn∆
we could use the equations of the model under P to write down a set of Fti∆-conditional moment conditions which would also identify the risk premium parameters η and µj.
Following Black-Scholes, the concept of implied volatility allows us to gauge the unobservable constant volatility parameter from market prices of traded derivative contracts. The concept can be extended to the more complicated model set-up we are in, so we can back out the implied-volatility state Vnχand the implied-intensity state λχn. As opposed to the Black-Scholes case where the implied volatility does not depend on any other parameters in the Black-Scholes model, in our case, the model-implied state variables depend on the full parameter vector
χ. When implied from option prices at the true parameter vector, say χ0, the implied-states
will be different from when they are implied from a different parameter vector χ. Given this identification we can devise an estimation strategy for the model parameters by choosing an initial set of parameters χ and then, by using the option pricing framework to imply Vnχand λχn, we evaluate moment conditions derived from the model under P based on the full state vector (log St, Vtχ, λχt)T and the choice of χ. We change the choice of χ and go through the process again until a moment based criterion function is minimized. This type of estimation approach was coined in Pan (2002) as implied-state GMM. In Pan (2002) the jump intensity is allowed to time-vary, however it is modeled as a multiple of the volatility, i.e. a λt= λ∗Vt, so the state
vector only contains the observable log-stock price and the unobservable volatility, being more parsimonious at the expense of a stronger assumption about the intensity process dynamics, which is essentially linked to a Brownian component.
5.1 Estimation Procedure
Given a time series for the underlying and a panel of call prices {Sn, Cn} observed at some
discrete time intervals {0, ∆, 2∆, 3∆, . . . , N ∆}, we use the transform based pricing to invert9 the pricing relationship for a parameter set χ and obtain the implied states at time point n∆:
Vnχ = f (Cn(K, Sn, T ), χ)
λχn = f (Cn(K, Sn, T ), χ) (5.1) 9The function f (·) generically denotes this inversion
The implied states are then used in the sample moment condition derived on the basis of the model (3.1)-(3.3): g(χ) = 1 N N X n=1 h (Sn, Vnχ, λχn) (5.2)
The GMM estimator ˆχ is:
ˆ
χ = arg min
χ g(χ)
TW
ng(χ) (5.3)
Where {Wn} is an adapted sequence of nh× nh positive definite weight matrices converging in the limit as n → ∞ to a positive definite matrix W .
5.2 GMM-based Inference
This section details the selection of moment conditions used and discusses the optimal instru-ment choice in this context. Given the known dynamics of the state variables10 we can derive approximations for the one time-period ahead expectation of raw moments for different state variable changes. We used the following raw moments:
RM1= E [d log St] = rt+ Vt(η −12) − µQλt ∆ RM2= Ed log St2 = rt+ Vt(η −12) − µQλt 2 ∆2+ Vt∆ + λtJ(2)∆ + 2 rt+ Vt(η −12) − µQλtλtJ ∆2 RM3= Ed log St3 = 3 rt+ Vt(η −12) − µQλtVt∆2+ 3 rt+ Vt(η −12) − µQλtJ(2)λt∆2+ + 3VtλtJ(1)∆2+ λtJ(3)∆ + o(∆3) RM4= Ed log St4 = 4 rt+ Vt(η −12) − µQλtλtJ(3)∆2+ 6VtλtJ(2)∆2+ λtJ(4)∆ + o(∆3) RM5= E [d log Vt] = kv V − V¯ t∆ RM6= E d log V2 t = k2 v V − V¯ t 2 ∆2+ σ2 vVt∆ RM7= E [d log λt] = kλ λ − λ¯ t ∆ + δλt∆ RM8= Ed log λ2t = k2λ λ − λ¯ t 2 ∆2+ 2kλ ¯λ − λtδλt∆2+ δ2λt∆ RM9= Ed log λ3t = 3kλ λ − λ¯ tδ2λt∆2+ δ3λt∆ + o(∆3) RM10= E [d log StdVt] = rt+ Vt(η −12) − µQλtkv V − V¯ t∆2+ σvVtρ∆ + kv V − V¯ tλtJ(1)∆2 RM11= E [d log Stdλt] = kλ λ − λ¯ t rt+ Vt(η −12) − µQλt∆2+ kλ λ − λ¯ tλtJ(1)∆2+ + rt+ Vt(η −12) − µQλtλtδ∆2+ δJ(1)λt∆ RM12= E [d log Stdλt] = kv V − V¯ tkλ λ − λ¯ t∆2+ kv V − V¯ tδλt∆2 RM13= Ed log StdVt2 = rt+ Vt(η −12) − µQλtσ2vVt∆2+2kv V − V¯ tσvVtρ∆2+J(1)λtVtσv2∆2+o(δ3) RM14= Ed log St2dVt= kv V − V¯ tVt∆2+kv V − V¯ tJ(2)λt∆2+2 rt+ Vt(η −21) − µQλtσvVtρ∆2+ . + 2J(1)λtσvVt∆2 RM15= Edλ2td log St= rt+ Vt(η −12) − µQλtδ2λt∆2+ δ2J(1)λt∆2+ o(∆3) RM16= Edλtd log S2t = kλ λ − λ¯ tVt∆2+ J(2)λtkλ λ − λ¯ t∆2+ Vtδλt∆2+ J(2)δ∆2
where J(i)is ith moment of the log-normal jump size distribution with mean µj and standard deviation σj
By subtracting the empirical counterpart ( ˜RM ) from the theoretical raw moment, we get a valid
moment condition with zero Fn-conditional expectation when evaluated at the true parameter set χ0. Moment condition i at time point n is:
hni(Sn, Vnχ, λnχ) = RMin−RM˜ in (5.4) 10
In order to have stationary moment conditions we use the log-stock price as a state variable instead of the stock price level, see Appendix C.2 for the state vector SDEs
The choice of the moment conditions is somewhat intuitive, the ’individual’ moment conditions based on a single state variable contribute to the identification of parameters that characterize the respective state variable, whereas ’cross’ moment conditions based on two state variables help identify the variables such as the Brownian correlation ρ and to better pin down the intensity process.
As we are using more moment conditions than parameters in the model, we should try to find an optimal instrument set which would lead us to the more efficient moment conditions. We follow11 Hansen (1985) and introduce the time point n instrument matrix An ∈ Fn of size
nχ× 16 defined as :
Aχn= DTn × Covχn(hn(Sn, Vnχ, λχn))
−1
(5.5) Where Dn is the conditional score of the moment condition with respect to the parameter set:
Dn= −
∂hn(Sn, Vnχ, λχn)
∂χ (5.6)
The instrument is hence the conditional score multiplied by the inverse of the conditional co-variance of the realized moment condition. Intuitively, the conditional coco-variance corrects for conditional keteroskedasticity and adjusts the scaling of the original moment conditions, while the Jacobian measures the ’sensitivity’ of each moment condition to the model parameters, weighing it in the optimal moment conditions. Hansen (1985) shows that this choice of in-struments leads to the asymptotically efficient estimator, however this result does not hold in our set-up. This is because the conditional score only picks up the sensitivity of the moment conditions to the parameters by treating the implied states as given, when these are in fact a function of the parameters as well as a function of the data. Given that the derivatives prices are obtained numerically, adjusting the Jacobian for the sensitivity of the implied states with re-spect to the parameter set is infeasible in our case. Hence at most we get more efficient moment conditions by using this instrument set-up. Hence the ’pseudo’-optimal moment conditions are:
˜
hn(Sn, Vnχ, λχn) = Anχhn(Sn, Vnχ, λχn) (5.7)
The GMM estimator becomes: ˆ χ = arg min χ g(χ)˜ TW ˜g(χ) (5.8) where ˜ g(χ) = 1 N N X n=1 ˜ h (Sn, Vnχ, λχn) (5.9)
The optimal weight matrix W is a consistent estimate of
S = Eh˜h (Sn, Vnχ, λχn) ˜h (Sn, Vnχ, λχn)T
i−1
(5.10)
and the asymptotic covariance matrix of the estimator is as in standard GMM:
Ω = ∆−1D˜TS−1D˜ (5.11) where ˜D denotes the score of the pseudo-optimal moment conditions with respect to the
pa-rameter vector.
5.3 Practical Issues
A first practical issue that arises when implementing the estimation approach is choosing the number of derivative contracts to imply states from. While in theory two call contract prices would suffice, after running the estimation procedure on a simulated data-set we reached the conclusion that using four or more call price contracts allows the backing out of implied states with a minimum of four decimal precision. In the empirical exercise we used four near at the money call contracts at each time point to back out the implied states to compute the criterion function.
Then there is the issue of computing time cost. Optimizing a criterion function which depends on the implied states that change with every parameter set change requires that for each GMM criterion function evaluation all the implied states for each time point get re-calculated. This in practice takes a lot of computing time and makes it almost impossible to use very long series of time points, albeit for the asymptotic properties of GMM estimator to kick in this would be desirable. To speed up the estimation we used at every time point four call contracts of the same maturity (as near to 30 calendar days as possible), hence reducing the number of transform calculations to the minimum (for each of the four call prices at every time point the transform values can be reused as long as the maturity is identical).
Last but not least, imposing restrictions on the parameter set domain is important in order to achieve sample identification. Imposing the Feller condition (2kvv ≥ σ¯ 2v) which ensures
that volatility does not reach the zero lower boundary and the stationarity conditions for the intensity parameters (i.e. kλ > δ) turn out to be very useful restrictions in the criterion function
minimization. Depending on the time-length of the data-set used, it can be useful to force the mean jump sizes to be negative µQj < µj < 0 as well as the correlation coefficient: ρ < 0. In
the empirical exercise that follows suit the convergence of the GMM criterion function to the same minimum point occurred irrespective of the starting values chosen, albeit selecting ’bad’ ones in the sense that they are away from the minimum parameter vector significantly increased computing time.
6
Empirical Application
6.1 Data description
In this applied estimation procedure we used joint spot and options data on the S&P 500 index from the Option Metrics database (accessible on our campus via WRDS). The data had daily frequency and provided for every trading day the best bid and ask close prices for all Chicago Board Options Exchange call contracts traded on the S&P 500 index. The data-set also provided the interest rate levels and expected dividend yields for each day. Taking the latter into account is important as they mainly impact the drift rate of the stock process which also contains the risk premiums we seek to investigate.
In order to keep computing time costs bearable we selected a time span of 95 trading days starting August 1, 2008 and ending December 22, 2008. The period covers the bankruptcy of Lehman Brothers which occurred on September 15, 2008 and its aftermath. This particular time span was chosen deliberately as it includes many significant drops in the S&P 500 index level which happened in rapid succession and could easily be reconciled with the self-exciting jump framework we use. It is worth mentioning that the applicability of the set-up is general and should fit through non-turbulent periods as well as turbulent ones, but given the computing time costs involved we choose a period which comprises many possible jumps to boost the identification of the jump dynamics and risk premium while still using a small length dataset.
For every trading day in the period four option contracts with strike/price ratios ranging between 0.9 − 1.1 and with the same time to maturity were selected. The strike price ratio is justified by the need to have as less microstructure noise in the call prices as possible while still varying the strike range. The maturity of the contracts was chosen every day to be as close as possible to 30 calendar days and was kept the same for all four contracts to save computing time.
6.2 Results
The GMM estimates of the model for the daily S&P 500 index, Aug 1, 2008 to Dec 22, 2008:
η kv V¯ σv ρ kλ λ¯ δ µj µQj σj
6.697 6.929 0.039 0.244 -0.391 21.31 0.30 13.18 -6.39% -11.10% 2.79% (2.9) (3.1) (0.02) (0.09) (0.02) (5.3) (0.23) (4.8) (3.82) (2.5) (0.83) The obtained parameter estimates are somewhat in line with other estimations conducted in the literature. For example, Andersen et al. (2002) for the S&P 500 find the estimates kv =
3.4, ¯V = 0.016, σv = 0.17, ρ = −0.33, while Pan (2002) finds for a model with jumps the
estimates kv = 7.1, ¯V = 0.0134, σv = 0.28, ρ = −0.52, η = 3.1.
Most parameter estimates are statistically significant at conventional confidence levels, a notable exception being the long run mean of the intensity process ¯λ and the long run mean of
4.71%. This would imply that investors are risk averse towards jumps and view jump as being ’more’ negative than they are. The diffusion risk premium is large, more than six times the daily standard deviation. This could be an artefact of the particularly volatile sample we use for the estimation.
The plots hereafter show the state variables (observed stock price, implied daily standard deviation and implied jump intensity) as implied from option prices when using parameter estimates.
A peculiar behavior can be noticed in the volatility plot. Volatility swings are frequent around big jump times, which can likely indicate a mis-estimated volatility process. Around the jump dates, the jump component is so important in size relative to the Brownian component to the extent to which over a small sample size of roughly one hundred trading days as is ours, the fitted model cannot distinguish between the two sources of variation. Using a sample with a longer time span is however difficult due to the high computing time costs.
6.3 Model-generated Volatility Surfaces
Equipped with the estimation results from the previous section we investigate the implications of the model for the (out of sample) cross section of option prices. Given that we used only 4 option prices cross-sectionally at each time point we can use some of the remaining options in our data-sample to benchmark how the model implied volatility structure12 for different moneyness levels compares to the (interpolated) volatility structure implied from option prices. We used the interpolated implied volatilities provided by Option Metrics and contained in our data-set. The plots below show the model implied volatility term structure fit for two days in our sample: September 16, 2008 the day after the turbulent default of Lehman Brothers and December 22, 2008 the last day covered in our sample. The aim was to see how the model fits during high jump intensity times compared to low jump intensity times. The plots confirm the
12
We imply volatilities under the Black-Scholes set-up, this approach is often used in the literature to benchmark option model mis-pricings relative to observed market prices.
intuition from the previous subsection, i.e. that the model fits better during turbulent times, e.g. soon after jumps occur, however during calm periods, when volatility plays a more impor-tant role than the jump component, pricing imprecisions occur. The source of the pricing errors could be attributed to the small sample sized used which could lead to imprecise parameter esti-mates, but it could also be due to the parametrization of the model and/or the parametrization of the probability measure change. Further investigation is necessary, however this extension is not addressed at this stage in this paper, being a topic reserved for future research. Appendix D presents a couple of model generated volatility surfaces based on the estimates in this section.
7
Conclusions
We proposed using a fully parametric model for asset returns which is more flexible in accom-modating price jump patterns than are the classic time homogeneous Poisson jump diffusions and jump diffusions in which the jump intensity is completely determined by another stochastic state variable such as the volatility process. After assuming a parametric form for the state price density we derived a semi-closed form approach to option pricing based on this model and devised a way to estimate its parameters using the joint time series of the spot and option prices. After implementing this procedure for the S&P 500 index we find evidence in favour of a time varying jump risk premium. Provided the sample span could be increased with more manageable computing time costs, the approach developed yields filtered series of volatility and jump intensities which could serve for practical applications.
Appendix A.1
Intensity Process Parametrization
Consider the intensity process parametrization using the notation in Hawkes (1971), such that the pair (N, λ) is a Markov process:
λt= ¯λ + t
Z −∞
βe−α(t−s)dNs (?)
To make sure the intensity process is stationary, consider taking its expectation:
E(λt) = ¯λ + E t Z −∞ βe−α(t−s)dNs = ¯λ + λ t Z −∞ βe−α(t−s)ds = ¯λ + λ ∞ Z 0 βe−αudu = ¯λ − λβ1 α e −αu|∞ 0 = λ∞− β αλ ⇒ λ = α α − βλ∞
Ensuring stationarity and non-negativity of the intensity process implies that α > β .
To obtain from the above form the stochastic differential equation governing the dynamics of λ take the differential of (?) with respect to time13
dλt dt = λ∞+ t Z ∞ βe−α(t−s)dNs 0
Taking derivative of the integral14 results in the following expression:
dλt dt = t Z −∞ ∂ ∂tβe −α(t−s) dNs+ ∂t ∂tβe −α(t−s) dNs|s=t= −α t Z −∞ βe−α(t−s)dNs+ β dNt dt From (?): λt= ¯λ + t Z −∞ βe−α(t−s)dNs⇒ t Z −∞ βe−α(t−s)dNs= λt− ¯λ
Substituting back in the differential expression for the integral we get the mean-reversion type expression for the instantaneous dynamics of the intensity of process:
dλt= −α(λt− ¯λ)dt + βdNt⇒ dλt= α(¯λ − λt)dt + βdNt
13Another way to prove this is by applying integration by parts for Stieltjes integrals in the expression for λ 14
Using Leibniz’s rule: ∂
∂α b(α) R a(α) f (α, x)dx = b(α) R a(α) ∂ ∂αf (α, x)dx − a(α) ∂α f (α, a) + ∂b(α) ∂α f (α, b)
Appendix A.2
Simulated paths from a Hawkes jump process
(a)
(b)
(c)
Sample simulated series with parameters: (a) T = 1; ¯λ = 5; α = 10; β = 5;
(b) T = 1; ¯λ = 5; α = 2; β = 1;
Appendix B
Candidate Pricing Kernel
In order to disentangle the diffusive price risk from the jump risk, consider a state price density that connects the model under the physical probability measure outlined in (3.1)-(3.3) to its version under the risk neutral probability measure (3.4) - (3.6) of the following form:
ξt= exp − Z t 0 rsds exp − Z t 0 γ1,sdWs1+ Z t 0 γ1,sdWs2 −1 2 Z t 0 γ1,s2 ds − Z t 0 γ1,s2 ds × × exp X s∗ i≤t Zs∗i (B.1)
γ1 and γ2 denote the market price of risk for the two Brownian motions which are defined as:
γ1,s= η p Vs, γ2,s= − ρ p 1 − ρ2η p Vs (B.2)
The jump part of the kernel is designed to keep model estimation tractable under both prob-ability measures, albeit the case for a richer specification can easily be made. The state price density jumps at the same time with the underlying asset, hence the last part of the right hand side term in (B.2), i.e.: exp P
s∗
i≤t
Zs∗i
!
, only ’kicks in’ when the underlying jumps as well, hence the s∗i time points denote the times at which the stock price experiences a jump. The price density jump sizes denoted by Zs∗
i are assumed to be i.i.d. normally distributed with mean
µ∗j and variance σj∗. Also we assume that the kernel jump sizes Zs∗i are independent from all other stochastic variables in the model such as the Brownian motions, the jump times, and all previous jump sizes in the kernel and in the log stock price, and only correlated ’contempora-neously’ with the jump size variable Zsi in the stock price. This allows us to keep the timing of the jumps (hence the intensity process dynamics) the same under both measures, while only varying the mean of the jump distribution between the risk neutral specification and physical probability specification of the model. More precisely, assume that the contemporaneous cor-relation between the stock price jump and the kernel jump is denoted by ρ∗. This means that we can express the mean relative jump size of the kernel as µ∗ = expµj+ σjσ∗jρ∗+12σj2
. By constraining the mean relative price jump size in the kernel to zero, i.e. forcing µ∗j+12(σ∗j)2= 0 we restrict the jump risk premia to be jump-size related and not jump-timing related.
To prove that (B.1) is a state-price density we show that the deflated processes of the stock (d ˜S) and the money market account (d ˜B) are local martingales (more on why this rules out
rule ( and assuming the deterministic dividend process is known w.r.t. the filtration): d ˜St= d ξStexp t Z 0 qsds = p Vt− γ1,t ˜ StdWt1− γ2,tS˜tdWt2+ + (exp(Z∗+ Z) − 1) ˜St−dNt− λtµ∗S˜tdt d ˜Bt= d ξ exp t Z 0 rsds = −γ1,tB˜tdWt1− γt,2B˜tdWt2+ (exp(Z ∗ ) − 1) ˜Bt−dNt (B.3)
Given the assumption previously made, that µ∗ = expµj+ σjσ∗jρ∗+12σj2
, it becomes clear that the first equation in B.3 is a local martingale (as the jump part is locally compensated) and using the other assumption µ∗j+12(σ∗j)2= 0, the second equation is also a local martingale. Finally, note that if we define the density process to be ξtexp(
t
R 0
rsds) (which, trivially, is also a
local martingale) it defines an equivalent martingale measure Q, under which:
Wt1,Q= Wt1+ t Z 0 γ1,sds Wt2,Q= Wt2+ t Z 0 γ2,sds (B.4)
Using (B.4) and taking into account the change in the mean expected jump size, the dynamics of the state variables under Q are those defined in (3.4)-(3.6).
Appendix C.1
Affine Jump Diffusion - The General Case
Given a Markov process - state vector X of size n × 1 defined on a domain D with values in Rn which follows a general affine jump diffusion of the form:
dXt= µ(Xt)dt + σ(Xt)dWt+ dZt
where W is a n × 1 vector of standard Brownian motions and Z is a jump process with non-negative intensity process {λ(Xt) : t ≥ 0}, so the intensity process could potentially depend on the state vector. The jumps that characterize the jump process Z have a fixed probability distribution ν which is independent of the state vector. The affine form of the other components is as follows:
• µ(x) = K0+ K1x, for K0∈ Rn and K1 ∈ Rn×n
• (σ(x)σ(x)T)
ij = (H0)ij+ (H1)ij· x for H0 ∈ Rn×n and H1 ∈ Rn×n×n, hence H1 is a tensor
of order n
• λ(x) = l0+ l1· x for l0∈ R and l1∈ Rn
• The jump size distribution ν is such that it is determined by the Fourier transform well defined for c ∈ Cn and such that θ(c) = R
Rn
exp(c · x)dν(z)
• The discount rate can be modeled as an affine function of the state R(x) = ρ0+ ρ1· x for
ρ0 ∈ R and ρ1 ∈ Rn, or as a constant by setting ρ1 = 0
The transform ψχ: Cn× D × R+× R+→ C of Xt conditional on the filtration Ftdefined by:
ψχ(u, Xt, t, T ) = Eχ exp − T Z t R(Xs)ds eu·Xt|Ft
where χ denotes the collection of parameters under which the expectation is computed: χ = (K, H, l, θ, ρ) Duffie et al. (2000) prove that under regularity conditions this transform is
ψχ(u, Xt, t, T ) = eα(t)+β(t)·x
where the vector functions α(t) and β(t) are the unique solutions to the ODEs: ˙
β(t) = ρ1− K1Tβ(t) − 1/2β(t)TH1β(t) − l1(θ(β(t)) − 1)
˙
α(t) = ρ0− K0β(t) − 1/2β(t)TH1β(t) − l0(θ(β(t)) − 1)
The application of this to option pricing arises by manipulating the expected value of a contin-gent claim on the state vector:
C(d, c, T, χ) = Eχ exp − T Z t R(Xs)ds ed·Xt − c+ = Eχ exp − T Z t R(Xs)ds ed·XT1 −d·XT≤− log c − − cEχ exp − T Z t R(Xs)ds e0·XT1 −d·XT≤− log c
Consider defining the general function
Ga,b(y : X0, T, χ) = Eχ exp − T Z t R(Xs)ds ea·XT1b·XT≤y
So the contingent claim price can be rewritten as:
C(d, c, T, χ) = Gd,−d(− log c : X0, T, χ) − cG0,−d(− log c : X0, T, χ)
Because Ga,b(y) is an increasing function one can apply a Fourier transform to it, essentially treating it like a distribution measure.
The Fourier transform of Ga,b is
Ga,b(z) =
+∞ Z −∞
eizydGa,b(y)
This integral can be computed explicitly and then using Fourier inversion the option price can be computed. Duffie et al. (2000) show that
Ga,b(v : X0, T, χ) = ψχ(a + ivb, X0, 0, T )
Hence to price the derivative, one first determines the transform ψχ(a + ivb, X0, 0, T ) by solving
the ODEs and then using the inverse Fourier transform back out the ’measure functions’ using
Ga,b(y : X0, T, χ) = ψχ(a, X0, 0, T ) 2 − 1 π ∞ Z 0 Im ψχ(a + ivb, X0, 0, T )e−ivy v dv
Appendix C.2
State Vector SDEs
Applying It¯o’s lemma and rewriting the process using the log-stock price Yt≡ log St as one of
the state variables in the system we get the following system of stochastic differential equations characterizing the full state vector (log St, Vt, λt)
dYt= (rt− √ Vt 2 − µ Qλ t)dt+ p VtStdWt1+ ZdNt dVt= kv( ¯V − Vt)dt + σv p Vt ρdWt1+ q 1 − ρ2dW2 t dλt= kλ(¯λ − λt)dt + δdNt
Appendix C.3
Our Model in Affine Jump Diffusion Form
We identify the coefficient matrices in our model in the form of an AJD (affine jump diffusion) using the notation in Duffie et al. (2000) as follows:
• µ(Xt) = K0+ K1Xt ⇒ K0= rt kvV¯ kλλ¯ K1= 0 −12 −µQ 0 −kv 0 0 0 −kλ • σ(Xt)σ(Xt)T ij = (H0)ij+ (H1)ij · Xt ⇒ σ(Xt)σ(Xt)T = Vt 1 σvρ 0 σvρ σ2v 0 0 0 0 ⇒ H0 = 0 0 0 0 0 0 0 0 0 H1 = k = 1; 0 0 0 0 0 0 0 0 0 k = 2; 1 σvρ 0 σvρ σv2 0 0 0 0 k = 3; 0 0 0 0 0 0 0 0 0 • λ(Xt) = l0+ l1· Xt ⇒ l0= 0 l1= 0 0 1
Appendix C.4
Monte Carlo Pricing Simulation Details
When implementing the Monte Carlo simulations and the transform based pricing, the process parameters were set at:
Drift Stoch. Vol. Intensity log St rel. jump Intensity jump
r = 0.0319 kν = 6.21 kλ= 10 µQj = −0.05 δ = 5
¯
V = 0.1 λ = 0.5¯ σj = 0.05
σv = 0.61
Besides initial state values, the ranges of maturities and strike prices of the call option contracts used in this simulation study were:15
Option Type Maturities K X0 initial state vector
Call T=0.10 80 log(S0) = log(100) T=0.25 90 V0= 0.1
T=0.50 100 λ0 = 2
T=1.00 110 120
Some sample simulated paths of the stock price, volatility and intensity process under the risk neutral measure:
(a) (b)
(c) (d)
(e) (f)
Sample simulated series: (a)&(b) Stock price paths; (c)&(d) Volatilities; (e)&(f) Intensities
15
Appendix D
Model Generated Volatility Surfaces
The following plots show the model generated volatility surfaces on September 16, 2008 and December 22, 2008, using the estimated parameters from Section 6.2 and the implied volatility and intensity series on the respective dates. They are generated under the assumption of a flat interest rate structure and a flat dividend yield structure, and thus not immediately suitable for comparison with traded-options’ implied volatility surface on the respective days.
(a) Sep 16, 2008 - Model-generated Volatility surface
References
Aït-Sahalia, Y., Cacho-Diaz, J., and Laeven, R. J. (2013). Modeling financial contagion using mutually exciting jump processes.
Andersen, T. G., Benzoni, L., and Lund, J. (2002). An empirical investigation of continuous-time equity return models. The Journal of Finance, 57(3):1239–1284.
Bacry, E., Delattre, S., Hoffmann, M., and Muzy, J.-F. (2013). Modelling microstructure noise with mutually exciting point processes. Quantitative Finance, 13(1):65–77.
Bates, D. S. (2000). Post-’87 crash fears in the s&p 500 futures option market. Journal of
Econometrics, 94(1-2):181–238.
Black, F. (1976). Studies of stock price volatility changes.
Black, F. and Scholes, M. (1973). The pricing of options and corporate liabilities. The Journal
of Political Economy, pages 637–654.
Bollerslev, T., Engle, R. F., and Nelson, D. B. (1994). Arch models.
Bollerslev, T. and Todorov, V. (2011). Tails, fears, and risk premia. The Journal of Finance, 66(6):2165–2211.
Cox, J. C., Ingersoll Jr, J. E., and Ross, S. A. (1985). A theory of the term structure of interest rates. Econometrica: Journal of the Econometric Society, pages 385–407.
Duffie, D. (2010). Dynamic asset pricing theory. Princeton University Press.
Duffie, D., Pan, J., and Singleton, K. (2000). Transform analysis and asset pricing for affine jump-diffusions. Econometrica, 68(6):1343–1376.
Eraker, B. (2004). Do stock prices and volatility jump? reconciling evidence from spot and option prices. The Journal of Finance, 59(3):1367–1404.
Eraker, B., Johannes, M., and Polson, N. (2003). The impact of jumps in volatility and returns.
The Journal of Finance, 58(3):1269–1300.
Errais, E., Giesecke, K., and Goldberg, L. R. (2010). Affine point processes and portfolio credit risk. SIAM Journal on Financial Mathematics, 1(1):642–665.
Ghysels, E. and Goldberger, S. (2012). Conditional asset pricing models with instruments closer to agentŠs information set. Technical report, Technical report, University of North Carolina. Hansen, L. P. (1985). A method for calculating bounds on the asymptotic covariance matrices
of generalized method of moments estimators. Journal of Econometrics, 30(1):203–238. Hawkes, A. G. (1971). Point spectra of some mutually exciting point processes. Journal of the
Heston, S. L. (1993). A closed-form solution for options with stochastic volatility with applica-tions to bond and currency opapplica-tions. Review of Financial Studies, 6(2):327–343.
Lord, R., Koekkoek, R., and Dijk, D. V. (2010). A comparison of biased simulation schemes for stochastic volatility models. Quantitative Finance, 10(2):177–194.
Lucas, R. E. (1978). Asset prices in an exchange economy. Econometrica: Journal of the
Econometric Society, pages 1429–1445.
Markowitz, H. (1952). Portfolio selection. The Journal of Finance, 7(1):77–91.
Merton, R. C. (1976). Option pricing when underlying stock returns are discontinuous. Journal
of Financial Economics, 3(1):125–144.
Pan, J. (2002). The jump-risk premia implicit in options: Evidence from an integrated time-series study. Journal of Financial Economics, 63(1):3–50.
