Accepted Manuscript
Calibration of a path-dependent volatility model: Empirical tests Paolo Foschi, Andrea Pascucci PII: DOI: Reference: To appear in: S0167-9473(08)00506-9 10.1016/j.csda.2008.10.042 COMSTA 4234 Computational Statistics and Data Analysis
Please cite this article as: Foschi, P., Pascucci, A., Calibration of a path-dependent volatility model: Empirical tests. Computational Statistics and Data Analysis (2008), doi:10.1016/j.csda.2008.10.042 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
a Faculty
of Economics, Forl` Campus, University of Bologna ı Piazzale della Vittoria, 15 - 47100 Forl` - Italy. ı
of Mathematics, University of Bologna Piazza di Porta S. Donato, 5 - 40126 Bologna - Italy.
b Department
Abstract
Key words: Option pricing, Hobson & Rogers model, Stochastic volatility, Calibration, Inverse problem
1
Introduction
AC C
Among non-constant volatility models in complete markets, the Hobson and Rogers (HR) model [24] seems to be one of the more appealing. In this model the volatility is a function of the trend of the underlying asset, defined as the difference between the spot price and a weighted average of past prices of the underlying asset. Since no exogenous source of risk is added, market completeness is preserved and the standard arbitrage pricing theory applies. Moreover, preliminary tests performed in [24] show that this model is potentially capable to reproduce the observed smiles and volatility term structure patterns.
∗ Corresponding author. Address: Faculty of Economics, Forl` Campus, University ı of Bologna - P.le della Vittoria, 15 - 47100 Forl` - Italy. Tel. +39 328 7655 498 ı Email address: foschip@gmail.com (Paolo Foschi).
Preprint submitted to Elsevier Science
EP
TE
DM
AN
The Hobson and Rogers model for option pricing is considered. This stochastic volatility model preserves the completeness of the market and can potentially reproduce the observed smile and term structure patterns of implied volatility. A calibration procedure based on ad-hoc numerical schemes for hypoelliptic PDEs is proposed and used to quantitatively investigate the pricing performance of the model. Numerical results based on S&P500 option prices are discussed.
US
CR
October 29, 2008
Paolo Foschi a,∗ and Andrea Pascucci b
IPT
Calibration of a path-dependent volatility model: empirical tests.
ACCEPTED MANUSCRIPT
Dτ = Z τ −
+∞ 0
λe−λs Zτ −s ds,
λ > 0,
dSτ = µ(Dτ )Sτ dτ + σ(Dτ )Sτ dWτ .
CR
where Zτ = log(e−rτ Sτ ) is the discounted log-price and the parameter λ amounts to the rate at which past prices are weighted. We assume that S is the solution to the stochastic differential equation
Uτ − rU + rSUS − λD +
σ2 2 σ2 UD + S USS + 2SUSD + UDD = 0, (3) 2 2
AC C
It should be recalled that the main advantages of the HR model against other stochastic volatility models consists on the fact that only one source of uncertainty drives both the price and the latent process. For this reason the model is able to incorporate path-dependency without loosing the market completeness property [24]. Recent works addressing the applications of stochastic volatility 2
EP
Compared with the assumption σ = σ(τ, Sτ ) of the widespread local volatility models, the path dependence property seems to be very realistic and natural. For instance it is well known that volatility increases after market reversals: this is difficultly captured by a model which only takes into account of the present price of the underlying and cannot relate past to present prices. Moreover a path-dependent model incorporates information on the preceding behavior of investors. Then, in some sense, the model “knows” how investors behave in different market circumstances and can also keep into account of the (positive or negative) trend of the asset. For this reason it seems that the HR model does not need to be frequently recalibrated, which is a known deficiency of local volatility models.
TE
DM
for (S, D, τ ) ∈ R+ × R×]0, T [. As in the Black&Scholes framework, the drift term in (2) does not enter in the valuation PDE while a key role is played by the volatility function σ which is an input of the model and has to be estimated in order to fit market observations.
AN
We denote by U = U (S, D, τ ) the price at time τ of an European option with maturity T , expressed as a function of the state variables: time τ , asset price S and deviation D. Then usual no arbitrage arguments lead to the pricing PDE (cf. formula (4.4) in [24])
US
The coefficients µ and σ in (2) are positive functions satisfying usual hypotheses in order to guarantee that the system of SDEs (1)-(2) is uniquely solvable and the solution (Sτ , Dτ ) is a Markovian process.
IPT
(1) (2)
A simplified version of the HR model is defined as follows: in a Wiener space with one-dimensional Brownian motion W , we denote by S the stock price and by D the deviation of prices from the trend, defined by
ACCEPTED MANUSCRIPT
In Section 3 the inverse problem arising in the calibration is stated as a simple nonlinear least squares problem. Aiming to motivate the model, Hobson and Rogers consider in [24] a volatility function of the form
EP
TE
As a first goal of this paper, in Section 2 we provide accurate numerical schemes ad hoc designed for the degenerate parabolic equation (3). These schemes exploit the intrinsic regularity properties of the PDE derived by the H¨rmander o condition and allow for an efficient implementation.
AC C
√ σ(D) = min η 1 + εD2 , N ,
DM
For these reasons, the HR model raised some interest among researchers. We recall that parameters calibration (λ in the average weight and the volatility function σ) was studied by Platania and Rogers [30], Fig`-Talamanca and a Guerra [17], Hahn, Putsch¨gl and Sass [22]. Sekine [31] obtained a closed-form o expression of the Laplace transform E e−z log Sτ for a special parametrization of σ. More general path dependent volatility models were proposed by Hubalek, Teichmann and Tompkins [25]; Foschi and Pascucci [19] proposed the use of a more general weight function, possibly with compact support, in (1) and analyzed the out-of-sample performance of the path dependent volatility in comparison with standard stochastic volatility models. An extension to the framework of term-structure modeling was given by Chiarella and Kwon [7]. Hahn, Putsch¨gl and Sass [23] considered the HR dynamics in a portfolio opo timization problem. The robustness of the HR model with respect to the data and parameters was studied by Hallulli and Vargiolu [5]. Trifi [34] showed that the HR model is the continuous time limit of an ARCH-type model. Recently the free boundary and optimal stopping problems for American options in a path dependent volatility model were studied in [29] and [14].
AN
US
CR
Another fine feature of the HR model is that, due to some invariance property of (3), a simple change of variables allows to evaluate all plain vanilla option prices corresponding to different strikes and different time-to-maturities in a single run (cf. Remark 1). This considerably speeds up the calibration procedure by PDEs’ techniques. Actually that approach also has the natural advantage of allowing to compute the derivatives with respect to the parameters of the solution which will be useful in the procedure.
for some large constant N and positive parameters ε, η: then they show that the model can indeed exhibit smiles and skews of different directions. In this note we aim to select σ without imposing a priori assumptions on its shape but simply calibrating it to market prices of plain vanilla options. More precisely, 3
IPT
(4)
models to the pricing, hedging and calibration of options includes [3,27,37]. Additional references on estimation, selection and forecasting using stochastic volatility model can be found in [2].
ACCEPTED MANUSCRIPT
in order to maintain the approach as much flexible as possible, we assume that σ is approximated in a space of piecewise cubic Hermite polynomials. In Section 3, we also introduce what we call the Hobson-Rogers recalibrated model which is a refinement of the standard HR model: in this alternative, as suggested by Hobson and Rogers, the level of the offset Dt used to evaluate new derivatives is inferred from the current cross section of option prices, instead of using the observed trajectory of St for its estimation. This new approach seemingly overcomes possible observation or model errors and gives a better fit to market prices.
Acknowledgement. We thank two anonymous referees for their valuable remarks that improved the paper.
2
Kolgorov schemes
For an European call option with strike K, equation (3) is coupled with the initial condition U (S, D, T ) = (S − K)+ .
DM
AN
US
(x, y) ∈ R2 .
Finally, in the last part of the paper, the results of the calibration are tested on a set of S&P500 index options prices and experimental results regarding the fitting of the model to observed prices are presented.
CR
We simplify the pricing PDE by a simple change of variables. We rewrite equation (3) as Lu := a(∂xx u − ∂x u) + (x − y)∂y u − ∂t u = 0, where u = u(x, y, t) is determined by the transformation U (S, D, τ ) = Ke−r(T −τ ) u(r(T −τ )+log(S/K), r(T −τ )+log(S/K)−D, λ(T −τ )), (7) and σ 2 (x − y) . (8) a(x, y) = 2λ By this change of variables, problem (3), (5) is equivalent to the Cauchy problem for (6) in the strip R2 × [0, λT ] with initial condition (6)
EP
TE
AC C
u(x, y, 0) = (ex − 1)+ ,
Remark 1 Problem (6)-(9) is independent of K. Therefore formula (7) gives option prices corresponding to different strikes by solving a unique PDE. 4
IPT
(5) (9)
ACCEPTED MANUSCRIPT
Xu := ∂x u,
Y u := (x − y)∂y u − ∂t u
AN
Since the natural framework for the study of the properties of equation (6) is the theory of subelliptic PDEs, also for the numerical approximation the best results are obtained in a non-Euclidean setting. Indeed it is known that the differential operators (10)
a X 2 − X u + Y u = 0. Therefore, in the numerical solution of the pricing equation by finite difference methods, it is natural and more efficient to approximate the directional derivatives X and Y rather than the usual Euclidean ones. Some “intrinsic” schemes for the HR model were first introduced in [11]. Here we refine those results and introduce some schemes, hereafter called Kolmogorov schemes, whose main features are the following: i) in the discretization of the PDE in a finite region, no boundary conditions on the y-variable are required. This is a significant advantage since there are no obvious financial motivations for imposing conditions on the option price U (S, D, t) for some fixed D; ii) solving the scheme only involves the inversion of a tri-diagonal matrix which leads to a fast and easy implementation. To be more specific, we consider the uniform grid G = {(i∆x , j∆y , n∆t ) | i, j, n ∈ Z, n ≥ 0}, 5 (11)
AC C
EP
TE
DM
are the main (in some intrinsic sense) directional derivatives and we can rewrite equation (6) in the form
US
2.1 Finite difference schemes of Kolmogorov type for the pricing PDE
CR
Due to the additional state variable D on which the option price depends, equation (6) is of degenerate type since the quadratic form associated to the second order part of L is singular. However (6) belongs to the noteworthy subclass of H¨rmander PDEs today called of Kolmogorov or Ornstein-Uhlenbeck o type. For this class a very satisfactory theory has been developed and many sharp analytical results are available even under weak regularity assumptions: see, for instance, [26] for an exhaustive survey on this topic. In particular, in [12] conditions are given for the existence and uniqueness of a classical solution to the Cauchy problem (6)-(9); moreover sharp estimates for u and its derivatives are provided.
IPT
ACCEPTED MANUSCRIPT
and approximate as usual the derivatives ∂x u and ∂xx u by the centered differences and the three-point schemes, respectively: ∂x u(x, y, t) ∼ D∆x u(x, y, t) = and
2 ∂xx u(x, y, t) ∼ D∆x u(x, y, t) =
u(x + ∆x , y, t) − u(x − ∆x , y, t) , 2∆x
Thus, the approximation
US
2 ∂xx u(x, y, t)−∂x u(x, y, t) ∼ D∆x u(x, y, t) − D∆x u(x, y, t) = d1 u(x − ∆x , y, t) + d2 u(x, y, t) + d3 u(x + ∆x , y, t),
CR
u(x + ∆x , y, t) − 2u(x, y, t) + u(x − ∆x , y, t) , ∆2 x (13)
with d1 = 1/∆2 + 1/(2∆x ), d2 = −2/∆2 and d3 = 1/∆2 − 1/(2∆x ), is of order x x x ∆2 . x The second main derivative Y is approximated either by
+ Y u(x, y, t) ∼ Y∆t u(x, y, t) =
u(x, y, t) − u(x, y − (x − y)∆t , t + ∆t ) , (15) ∆t
or by
− Y u(x, y, t) ∼ Y∆t u(x, y, t) =
u(x, y, t) = (1 − γ)u(x, y, t) + γu(x, y + ∆y , t),
TE
where u(x, y, t) denotes the linear interpolation of u at the point (x, y, t) based on the two nearest grid points. Specifically, (17)
AC C
where γ = (y − y)/∆y and y = [y/∆y ] ∆y denoting by [·] the integer part. Since u(x, y, t) approximates u(x, y, t) with an error of the order of ∆y , then the approximations (15) and (16) are of the order of ∆t + ∆y . We remark that interpolation (17) is necessary because (x, y, t) and (x, y − (x − y)∆t , t − ∆t ) cannot both belong to the same uniform grid. In [10] a different change of variables has been proposed in place of (7). That approach allowed for both the points to belong to the grid, but at the cost of imposing the grid size condition ∆y = ∆x ∆t . The discrete operators L+ and L− are defined by G G
2 ± L± u = a(D∆x u − D∆x u) + Y∆x u, G
EP
DM
6
u(x, y + (x − y)∆t , t − ∆t ) − u(x, y, t) , (16) ∆t
AN
IPT
(12) (14) (18)
ACCEPTED MANUSCRIPT
and approximate L in the sense that
L∞
for some positive constant C depending on the L∞ -norms of a, ∂xxx u, ∂y u, 4 ∂x u, Y 2 u, ∂xx Y u, and ∂xxy u.
DM
7
In the calibration procedure we consider the volatility σ as a smooth function of a parameter vector α ∈ Rp for some p ∈ N: more precisely, we assume that + σ = σ(d; α) ∈ C 1 (R × Rp ) and that there exist two positive constants C1 , C2 + such that C1 ≤ σ(d, α) ≤ C2 , |∂αk σ(d, α)| ≤ C2 (d, α) ∈ R × Rp , + k = 1, . . . , p.
Thus we rewrite the dynamics (2) of the price as
TE
dSt = µ(Dt )St dt + σ(Dt ; α)St dWt ,
2
AN
2.2 Calibration and continuous dependence results
US
In the Appendix we formulate the discretization of the PDE (6) by means of (18), as the block bi-diagonal linear system (A.4). It is remarkable that the solution of such a system only requires the inversion of a tri-diagonal matrix which can be computed very efficiently.
CR
Hereafter, we refer to L+ and L− respectively as explicit and implicit schemes G G for the discretization of L. The implicit scheme is unconditionally stable, while ∆2 x the stability condition for the explicit method is given by ∆t ≤ 2 sup a and ∆x < 2.
Theorem 2 Let L(α) be the operator defined by (6) with diffusion coefficient a = a(·; α) and consider the solution u to the Cauchy problem
AC C
EP
and denote by u(·; α) the solution to (6)-(9) with a(x, y; α) = σ (x−y;α) . For 2λ what follows, it will be useful to compute the derivatives of u with respect to the parameters α. A linear system, similar to (A.4), derives from the discretization of the corresponding PDE. Indeed the following result about continuous dependence with respect to the parameters holds.
L(α) u(·; α)
where the initial datum u0 is H¨lder continuous and such that |u0 (x, y)| ≤ o C4 (x2 +y 2 ) C3 e with C4 suitably small. Then, u(z; ·) ∈ C 1 (Rp ) for every z = +
= 0, u(x, y, 0; α) = u0 (x, y),
in R2 ×]0, λT [, in R2 ,
IPT
(20) (21)
Lu − L± u G
≤ C ∆2 + ∆t + x
∆2 y , ∆t
(19)
ACCEPTED MANUSCRIPT
(x, y, t) ∈ R2 ×[0, λT ]. Moreover for k = 1, . . . , p, the derivative v = ∂αk u(·; α) satisfies the PDE L(α) v = −(∂αk a(·; α))(∂xx u(·; α) − ∂x u(·; α)), in R2 ×]0, λT [, with initial condition v(x, y, 0) = 0, for (x, y) ∈ R2 .
The proof of the theorem is based on the following
Proof of Theorem 2. For simplicity, we only consider the case p = 1. We set z = (x, y, t) ∈ R2 ×]0, λT [ and denote by Γ(α) the fundamental solution of L(α) . Since v(z) := u(z; α) − u(z; α0 ) is solution to the problem
we have
EP
u(z; α) − u(z; α0 ) = α − α0
AC C
Therefore, as α goes to α0 , by the dominated convergence theorem combined with Lemma 3, we infer that u is differentiable w.r.t. α and it holds ∂α u(z; α0 ) =
t 0 R2
This concludes the proof. 8
TE
L(α0 ) v
DM
Γ(α0 ) (z; ξ, η, s)·
A detailed proof of Lemma 3 is given in [13]. The proof is rather delicate since it involves the study of some singular integrals related to degenerate parabolic PDEs. We only mention that the continuity of ∂x u and ∂xx u can be proved by using the techniques in [12] (see, in particular, Chap. 5 regarding some potential estimates); estimate (24) can be proved as Theorem 8.2, Chap. V in [15], by using the pointwise estimates for the fundamental solution and its derivatives provided in [12], Theorem 1.4.
AN
= −(a(·; α) − a(·; α0 ))(∂xx u(·; α) − ∂x u(·; α)), v(x, y, 0) = 0,
t R2
0
a(ξ, η; α) − a(ξ, η; α0 ) (∂ξξ − ∂ξ )u(ξ, η, s; α)dξdηds. α − α0
Γ(α0 ) (z; ξ, η, s)∂α a(ξ, η; α0 )(∂ξξ − ∂ξ )u(ξ, η, s; α0 )dξdηds. 2
US
Lemma 3 Under the above assumptions, ∂x u, ∂xx u are continuous functions w.r.t. the variables (x, y, t, α) and there exist some positive constants C5 , C6 and δ > 0, only dependent on C1 , . . . , C4 and the H¨lder constant of u0 , such o that 2 2 eC6 (x +y ) |∂x u(x, y, t; α)| + |∂xx u(x, y, t; α)| ≤ C5 , (24) t1−δ for every (x, y, t) ∈ R2 ×]0, λT [ and α ∈ Rp . +
CR
in R2 ×]0, λT [, in R2 ,
IPT
(22) (23)
ACCEPTED MANUSCRIPT
2.3
Boundary conditions
The numerical solution of (6) by finite-difference methods requires the discretization of the equation in a bounded region and the specification of some initial-boundary conditions. More precisely, we approximate the Cauchy problem (6)-(9) in the cylinder
for some suitably large µ, ν. By transformation (7), this corresponds to the initial-boundary value problem for (3) in the domain of the points (S, D, t) such that Ke−µ−r(T −t) < S < Keµ−r(T −t) , |D − log er(T −t) S/K | < ν and 0 < t < T. The conditions on the parabolic boundary of Q, defined by ∂P Q = ∂Q ∩ {(x, y, t) | t < λT }, are set as follows: u(x, y, 0) = (ex − 1)+ , moreover, we set
AN
i, j ∈ Z,
for x ∈ [−µ, µ], y ∈ [−ν, ν];
US
|i| ≤ i0 ,
CR
|j| ≤ j0,n .
Q = {(x, y, t) | |x| < µ, |y| < ν and 0 < t < λT },
(∂xx u − ∂x u)(±µ, y, t) = 0,
DM
9
for y ∈ ] − ν, ν[, t ∈ ]0, λT [.
We note explicitly that (27) corresponds to condition ∂SS U = 0 in the original variables, that imposes linearity on U at the border. This condition is standard in the Black&Scholes framework [33,35] and it has been analyzed in [28,36]. It is remarkable that the Kolmogorov schemes avoid imposing conditions on the lateral boundary {y = ±ν}, provided that ν is suitably large. To be more specific, let us first introduce some notation. Fixed i0 , j0,n ∈ N for n ∈ N∪{0}, we denote un = u(i∆x , j∆y , n∆t ), i,j
EP
TE
Applying the discrete operator in (14) to un for |i| ≤ i0 − 1 gives i,j
AC C
n n D∆2 un − D∆x un = (d1 un i,j i−1,j + d2 ui,j + d3 ui+1,j ). x i,j
Consider now the discretization (16) and assume that (x, y, t) = (i∆x , j∆y , n∆t ) belongs to the grid. Then we have u(x, y + (x − y)∆t , t − ∆t ) = (1 − γ) un−1 + γ un−1 i,j+k i,j+k+1
IPT
(25) (26) (27) (28) (29)
ACCEPTED MANUSCRIPT
and
where k and γ are, respectively, the integer and fractional part of (x−y)∆t /∆y , that is
Applying the discrete operator L− to un we have G i,j
− aij D∆2 un − D∆x un + Y∆t un = 0, i,j i,j x i,j
where aij = a(i∆x , j∆y ). On the other hand, assuming µ = i0 ∆x , condition (27) is equivalent to
− Y∆t un = 0, i,j
Next we fix j0,N ∈ N and examine the domain of dependence of the set of values
N UN = uN | |i| ≤ i0 , |j| ≤ j0 ; i,j
more precisely, for 0 ≤ n ≤ N − 1, we specify j0,n as the maximum of the set of the indexes j’s such that UN depends on un through conditions (32)-(33). i,j Moreover we set νn = j0,n ∆y . Since, by (30) and (31), it holds j0,n−1 = j0,n + j0,n ∆t + i0 it follows that ∆x ∆t ∆x ∆t + 1 ≤ j0,n (1 + ∆t ) + i0 + 1, ∆y ∆y
and thus νN −n ≤ zn , where zn is defined by the difference equation zn+1 = (1 + ∆t )zn + µ∆t + ∆y , z 0 = νN ,
AC C
which has solution zn = (1 + ∆t )n (y0 + µ + ∆y /∆t ) − µ − ∆y /∆t : indeed, recall that the solution of the difference equation zn+1 = αzn + β, with initial value z0 and α = 1, is given by zn = αn (z0 − z∗ ) + z∗ , where z∗ = β/(1 − α) is the equilibrium value. Finally, we deduce νn = (1 + ∆t )(N −n) (νN + µ + ∆y /∆t ) − µ − ∆y /∆t 10
EP
TE
νn−1 ≤ νn (1 + ∆t ) + µ∆t + ∆y ,
DM
AN
i = ±i0 , |j| ≤ j0,n .
US
|i| ≤ i0 − 1, |j| ≤ j0,n ,
CR
k=
(x − y)∆t = ∆y
i
∆x − j ∆t ∆y
and
γ=
(x − y)∆t −k . ∆y
IPT
(31) (32) (33)
− Y∆t u n = i,j
1 un − (1 − γ)un−1 − γun−1 i,j+k i,j+k+1 , ∆t i,j
(30)
ACCEPTED MANUSCRIPT
and
Thus we have proved the following result. Theorem 4 Assume that
Remark 5 Notice that, under condition (34), the approximation error of L± G in (19) reduces to an order of ∆2 + ∆t . x Moreover, condition (34) ensures that the width of the initial region can be chosen independently of the refinement of the grid. Alternatively, one can solve (6) in the prism {(x, y, t) | |x| < µ, |y| < eλT −t νN + (eλT −t − 1)(µ + C0 ), 0 < t < λT }, rather than in the whole cylinder Q.
3
Calibration
EP
In this section we consider the calibration of the HR model, that is the estimation of the volatility function σ from observed market prices of European options. As u is observed only at a finite number of points, specific restrictions should be imposed on σ in order to obtain a well posed problem. For this reason, as in Subsection 2.2, we assume that σ = σ(·; α) depends on a vector α = (α1 , . . . , αp ) of real positive parameters and denote by u(x, y, t; α) the solution to the Cauchy problem (6)-(9) corresponding to
TE
a=
AC C
This defines a mapping from Rp to C ∞ (R2 ×]0, λT [). The aim of this section + is to develop and test a numerical procedure to “invert” that function. ˆ Let fi be the observed option value at the point zi ≡ (xi , yi , ti ), for i = 1, 2, . . . , M , and let fi (α) be the price given by (7) in terms of the solution u(xi , yi , t; α) of the PDE (6) for a given α at the observation point zi . Since the 11
DM
σ 2 (x − y; α) . 2λ
AN
US
∆y = C0 , (34) ∆t for some constant C0 . Then, in order to approximate the solution u(x, y, λT ) for |x| ≤ µ and |y| ≤ νN , conditions on the lateral boundary {y = ±ν}, where ν = eλT νN + (eλT − 1)(µ + C0 ), are superfluous.
CR
IPT
ν0 = (1 + ∆t )λT /∆t (νN + µ + ∆y /∆t ) − µ − ∆y /∆t ≤ eλT νN + (eλT − 1)(µ + ∆y /∆t ).
ACCEPTED MANUSCRIPT
ˆ εi (α) = fi (α) − fi .
M α∈R+ i=1
min p
ε2 (α)/wi + ρ(α), i
CR
We aim to find α that best fits the data by solving the nonlinear least squares (NLLS) problem (35)
3.1 The dataset
AC C
Following Ait-Sahalia and Lo [1] we do not use quotations on the underlying St and realized dividends δt , because the first ones are affected by synchronization errors and the second ones are not necessarily equal to the expected rates at the time the option is priced. In order to avoid pricing anomalies that can arise from these problems the following procedure is used: at each time t and for each maturity T and strike K, in order to avoid arbitrage opportunities, 12
EP
The calibration procedure here described is applied to a set of European options quotations on the S&P 500 index from the Chicago Board Options Exchange (CBOE). Only calls with time-to-maturity from two weeks to six months are considered, moreover the average of bid and ask prices is used as reference. Observations have been taken each 15 minutes from 11:00 to 14:00 of each trading day in the period from Nov-15-2002 to May-23-2003. The empirical distribution of this dataset, which contains 190397 observations, w.r.t. absolute time, deviation from the trend, time-to-maturity and moneyness is shown in Figure 1.
TE
DM
where wi is some weight given to the i-th observation which will be specified later in Subsection 3.2 and ρ is positive definite regularization function introduced to penalize large estimates of σ and to render the objective function free of non-optimal stationary points. The NLLS problem (35) is solved using the Matlab routine lsqnonlin which is a trust-region method based on the interior-point method described in [8]. The algorithm needs the first order derivatives ∂αk u at the points zi , for i = 1, . . . , M , which are computed by solving a set of p + 1 PDEs (6) and (22): thus the computational cost linearly increases with the number p of parameters involved. In the experiments performed, the non-convexity of the objective function affected only the speed of convergence of the method, but not its robustness.
AN
US
IPT
point zi may not belong to the grid G, the value of fi (α) is approximated by using a linear interpolation of the nearest points of the grid. The error made in fitting the i-th observed value for a given α is denoted by
ACCEPTED MANUSCRIPT
10000
5000
0
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 Week
20000
0
!0.24 !0.22
!0.2
!0.18 !0.16 !0.14 !0.12 !0.1 Deviation from the mean
40000
Freq.
0
0.5
1
1.5
2
2.5 3 3.5 Months to expiration
AN
4
20000
US
4.5 5
!0.08 !0.06 !0.04
0
!0.35 !0.3 !0.25 !0.2 !0.15 !0.1 !0.05 !0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 Moneyness log(S/K)
Figure 1. Distribution of observations. The four panels show the number of observations plotted against time (measured in weeks), deviation from the mean, time to maturity (measured in months) and moneyness.
EP
the following parity relations hold: FT −t = St e(rt,T −δt,T )(T −t) ,
and
AC C
put call UT −t (St , K, rt,T , δt,T ) + Ke−rt,T (T −t) = UT −t (St , K, rt,T , δt,T ) + FT −t e−rt,T (T −t) ,
put call where FT −t , UT −t and UT −t denote, respectively, the future, call and put prices with expiration at T , rt,T and δt,T the corresponding interest and dividend rates. Thus, future prices can be inferred from synchronous quotations of put
TE
DM
13
50000
Freq.
CR
5.5
40000
Freq.
IPT
Freq.
ACCEPTED MANUSCRIPT
to replace option prices U (St , K, rt,T , δt,T ) with ex-rates prices U (FT −t , K, 0, 0). Notice that many authors modify each cross section before the calibration. For example in [6] implied volatilities curves are computed by smoothing splines from market observations and in [4] constrained cubic splines are used to smooth market prices and to enforce non-arbitrage conditions. In both cases the reported results refer to the modified values. Here, we only infer the underlying future price and use raw data for the fitting. Furthermore, differently from other investigations, where transactions or daily closing prices are used, we use a dataset built from intra-day observation of bid-ask prices which do not necessarily correspond to some transaction. To compute the exponential trend Mt and the corresponding deviation Dt we have used closing day prices from October 1982 to September 2002 and then intra-day prices until May 2003. Index and trend computed with λ = 1 are shown in Figure 2. We have chosen λ = 1 as it gives the best reproduction of the term structure of implied volatilities. The period considered in the dataset has an initial decreasing phase followed by sharp rise of the index level.
1400 1200 1000 800 600
TE
Jul97
DM
Jan00
AN
US
Jul02
Figure 2. Price of the S&P 500 index St and the corresponding exponential trend computed with λ = 1 for the years 1997-2004.
AC C
The dependency of the option prices on Dt is evident from Figure 3. The four panels√ plot the implied volatilities against the adjusted log-moneyness log(F/K)/ T − t: as it was noted by Foque et al. [18] the implied volatility √ smile of a cross section of index options is a function of log(F/K)/ T − t only. This relation captures the term structure of option implied volatilities. Each panel refer to a different range of Dt . The base of the smile shifts from 14
EP
CR
Jan05
ert,T (T −t) UT −t (St , K, rt,T , δt,T ) = UT −t (FT −t , K, 0, 0) = f (FT −t , Mt , T − t),
IPT
and call options with the same strike and maturity. In order to obtain reliable values, FT −t is computed as the weighted average of the implied futures over all the available strikes. Log of volumes plus one are used as weights, even if a simple average produces similar results. Next, in order to reduce the number of input parameters of the calibration procedure, we use the homogeneity relation
ACCEPTED MANUSCRIPT
Dev. from trend 0.5 0.45 Implied Volatility 0.4 0.35 0.3 0.25 0.2 0.15
Dev. from trend
Dev. from trend
Dev. from trend 0.5
!0.25 !0.2 !0.15 !0.1 !0.05 !0.25 !0.2 !0.15 !0.1 !0.05 !0.25 !0.2 !0.15 !0.1 !0.05 !0.25 !0.2 !0.15 !0.1 !0.05
!0.5 0 0.5 Moneyness
!0.5 0 0.5 Moneyness
!0.5 0 0.5 Moneyness
!0.5 0 0.5 Moneyness
3.2
Calibration results
AC C
The HR model has been calibrated to the data for two different choices of the volatility function σ 2 (D). Following [20], the first choice consists in a piecewise 2 2 cubic Hermite polynomial σSpline (D) interpolating the abscissae σSpline (Di ) = αi , for 1 ≤ i ≤ 7 at the knots Di ∈ {−1, −2/3, −1/3, 0, 1/3, 2/3, 1}. In Figure 4 the interpolating polynomial and the knots are represented by a dashed line 2 and by circles, respectively. The function σSpline (D) is chosen linear outside the interval [−1, 1] continuous up to the first derivative in D = 1 and D = −1. Further, its positiveness is ensured by the constrains αi ≥ 0, 1 ≤ i ≤ 7, in the NLLS problem (35). Hereafter, we refer to this function as the Spline volatility function and to the corresponding model as the Spline model. 15
EP
TE
Since the observations are not homogeneously distributed with respect to timeto-maturity, moneyness and deviation from the trend, we adopt the following strategy to choose the weights wi in (35). We divide option transactions in 18 groups based on maturity ([0, 3] or ]3, 6] months), log-moneyness (] − ∞, −.1], ] − .1, .1] and ].1, ∞[) and deviation from trend (] − ∞, −.2], ] − .2, −.15], ] − .15, −.1], ] − .1, ∞]). The weight wi of the ith observation has been chosen equal to the number of elements in the corresponding group.
DM
AN
an implied volatility of 0.22 to one of 0.15 and from a log-moneyness of -0.5 to one of -0.25. Notice that, the observations of Dt , like St , are affected by errors, due to synchronization or to market inefficiencies.
US
Figure 3. Effects of the deviation from the trend on marked implied volatilities. The √ implied volatilities are plotted against adjusted log-moneyness log(F/K)/ T − t and grouped by different ranges of D as shown by the bar in the top of each panel.
CR
0.25 0.2 0.15
IPT
0.45 0.4 0.35 0.3
ACCEPTED MANUSCRIPT
10 8 Volatility function 6 4 2 0 !1.5 !1 !0.5 0 0.5 Deviation from trend 1 1.5 HR Spline Spline knots Volatility function 1 0.8 0.6 0.4 0.2 0 !0.4 HR Spline Spline knots
!0.2 0 0.2 Deviation from trend
CR
2 2 Figure 4. Calibrated Volatility functions σHR and σSpline .
√ 2 σHR (D) = min α1 + α2 (D − α3 )2 , 5 ,
To analyze the quality of the fit, four kinds of errors are considered: ˆ • absolute or dollar (valuation) error: εi = fi − fi ; ˆ; • percentage or relative error: εi = εi /fi ˆ • error outside the bid-ask spread, or simply outside error: ˆ • percentage or relative outside error: εOE = εOE /fi ; ˆi i ˆ ˆ where fiask and fibid are the bid and ask prices of the ith observation (cf. [16]). In order to partially eliminate the bias in percentage errors, only calls with price larger than 10$ are considered when computing statistics for relative errors (note that the value of the underlying is 887$ in mean). 16
AC C
EP
ˆ ˆ εOE = sign(εi ) max(fi − fiask , fibid − fi , 0); i
TE
In order to overcome possible observation or model errors in Dt , we have re-calibrated the offset parameter α3 of the HR model for each day of the dataset. More precisely, the parameters α1 and α2 have been kept fixed to the previously estimated values, that is to the values reported in Table 2 and the parameter α3 is calibrated using the observations of each given day: this, (3) results in a series of estimates αt which allows to recover the option implied (3) deviations defined by Dt = Dt + α3 − αt . The two time-series Dt and Dt are shown in Figure 5, where it can be noticed that the difference is often substantial. We refer to this calibration as the HR re-calibrated model.
DM
AN
with the constrains αi ≥ 0, for i = 1, 2 and α3 unconstrained. We refer to this function as the HR volatility function and to the corresponding model as the HR model. The calibrated parameters of the two functions and their standard deviations are reported in Table 2. Figure 4 shows the graphs of the two volatility functions.
US
The second choice consists in
IPT
0.4 0.6
ACCEPTED MANUSCRIPT
!0.05 Esti4ated 74plied !0.1 !0.15 /
t
!0.2 !0.25 !0.3
0 1 2 3 4 5 6 ( 8 * 101112131415161(181*202122232425262( ,ee.
Table 1 contains the following resuming statistics for the three models and the four types of errors: mean error, standard deviation, extreme values, Root Mean Square Error (RMSE) RMSE(εi ) = M and Mean Absolute Error (MAE),
AN
−1
M
i=1
MAE(εi ) = M −1
DM
17
M
i=1
|εi |.
The four panels in Figures 6 and 7, respectively, plot the RMSEs of absolute and relative outside errors for different ranges of absolute time, deviation from the trend, time-to-expiration and moneyness. We first compare the fits of the HR and Spline specifications. The overall performances of both models are good, also considering the length of the period and that the volatility function, and thus the pricing kernel, is kept constant. The good results are confirmed by the histograms of absolute and percentage errors for the HR and Spline calibrations shown in Figure 8. From the third panel of Figures 6 and 7 it should be noticed that both the models well explain the term-to-maturity structure of option prices. This has been pointed as a deficit of stochastic volatility models in [18,21,32].
AC C
The Spline model obtains slightly better results than those of the HR model, but the difference is so small that the HR model should be preferred for its parsimony. It can be seen from Figure 4 that the two volatility functions differ only on their tails and especially on the right ones. However, as indicated by 2 the standard deviations reported in Table 2, σSpline is completely imprecise 2 in that region. In fact, the last two parameters of σSpline are not identifiable
EP
TE
US
1/2
Figure 5. Comparison of two estimators for Dt . The thin line represents the time series computed from historical data by discretizing (1). The thick line corresponds to the value for Dt implied from each daily cross sections of options.
ε2 i
,
CR
IPT
ACCEPTED MANUSCRIPT
Absolute Errors Model Spline HR HR Recal. RMSE 1.853 1.857 1.532 MAE 1.458 1.463 1.229 Mean -0.757 -0.789 -0.613 Std. Dev. 1.692 1.681 1.403 Minimum -5.944 -6.113 -4.948 Maximum 7.049
7.131 4.81
Model Spline HR HR Recal
RMSE 1.153 1.154 0.852
MAE 0.740 0.741 0.533
Mean -0.366 -0.387 -0.271
Std. Dev. 1.094 1.087 0.808
Minimum -5.532 -5.846
US
-4.131 Minimum -34.40% -36.56% -26.00% Minimum -29.57% -31.72%
Percentage Errors Model Spline HR HR Recal. RMSE 5.16% 5.19% 4.33% MAE 3.29% 3.34% 2.73% Mean -0.71% -1.05%
Std. Dev. 5.11%
AN
5.09% 4.30% Std. Dev. 3.35% 3.33%
-0.47%
Percentage Outside Errors Model Spline HR RMSE 3.37% 3.37% MAE Mean Maximum 29.69% 28.66%
1.70% 1.73%
AC C
Next we discuss the performance of the HR re-calibrated model. This approach allows for a further improvement in the fit: both the absolute and relative outside errors are reduced of one third (cf. Table 1). Furthermore, implying Dt from option prices allows to better capture the term structure of option prices (see the third panel in Figures 6 and 7). While the distribution of the residuals of the HR and Spline models, shown in Figure 8, is clearly not normal, that of the HR re-calibrated model seems to be a normal one. This should point out a misspecification of Dt in the first two models that almost disappears in 18
EP
because they do not significantly affect the solution u(x, y, t) of (6) at the observation points.
TE
HR Recal. 2.62% 1.19% -0.26% 2.60% -23.01% 31.62% Table 1 Statistics for absolute (relative) errors and outside errors for the calibrations with the Spline, HR and HR re-calibrated models.
DM
-0.34% -0.56%
CR
6.043 6.125 3.976
Absolute Outside errors
Maximum
IPT
Maximum 33.08% 32.05% 34.06%
ACCEPTED MANUSCRIPT
Outside Errors 2 Spline HR HR Recalibrated
1 0.5 0
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 Week
1.4 1.2 RMSE 1 0.8 0.6 0.4 !0.24 !0.22 !0.2 !0.18
!0.16 !0.14 !0.12 Deviation from trend
US
!0.1 !0.08 4 4.5 5
2 1.5 RMSE 1 0.5 0
0.5
1
1.5
2
2.5 3 3.5 Months to expiration
AN
2 1.5 RMSE 1 0.5
AC C
the HR re-calibrated model. This also suggests that further investigation in the model with respect to the specification of the state variable Dt is needed. With this respect, we mention the extensions considered in [19] and [25] and the use of an additional deviation Dλ2 ,t with a faster decay (a shorter memory) or a second order offset function Dt =
(2) ∞ 0
EP
Figure 6. Plot of RMSEs of absolute outside errors against observation time (measured in weeks), deviation from the mean, time to maturity (measured in months) and moneyness. Results for the Spline, HR and HR re-calibrated models are represented by circles, squares and diamond, respectively.
as already considered in [24], to better explain the market behavior. 19
TE
0 !0.35 !0.3 !0.25 !0.2 !0.15 !0.1 !0.05 !0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 Moneyness log(S/K)
DM
λe−λτ (Zt − Zt−τ )2 dτ,
CR
!0.06 !0.04 5.5
IPT
1.5 RMSE
ACCEPTED MANUSCRIPT
Relative Outside Errors 0.08 Spline HR HR Recalibrated
0.04 0.02 0
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 Week
0.06 0.05 RMSE 0.04 0.03 0.02 0.01 !0.24 !0.22 !0.2 !0.18
!0.16 !0.14 !0.12 Deviation from trend
US
!0.1 !0.08 4 4.5 5 0.25 0.3
0.05 0.04 RMSE 0.03 0.02 0.01 0 0.5 1 1.5 2
2.5 3 3.5 Months to expiration
AN
0.08 0.06 RMSE 0.04 0.02 0
AC C
Figure 9 shows implied volatilities for the HR model using the daily calibrated deviations Dt . As can be seen by comparing Figures 3 and 9, the model is not able to fully reproduce the smiles for deep in- or out-of-the money options. Furthermore, smiles implied by the market are more significant especially for shorter maturities. 20
EP
Figure 7. Plot of RMSEs of relative outside errors against observation time (measured in weeks), deviation from the mean, time to maturity (measured in months) and moneyness. Results for the Spline, HR and HR re-calibrated models are represented by circles, squares and diamond, respectively.
TE
!0.15 !0.1
!0.05
DM
!0
0.05 0.1 0.15 0.2 Moneyness log(S/K)
CR
!0.06 !0.04 5.5 0.35 0.4
IPT
0.06 RMSE
ACCEPTED MANUSCRIPT
4 3.5 3 2.5 Freq. 2 1.5 1 0.5
x 10
4
7 Spline HR HR Recalibrated 6 5 4 3 2 1
x 10
4
Spline HR HR Recalibrated
0 !8
!6
!4
!2
0 Errors
2
4
6
8
0 !0.4
!0.3
!0.2
CR
!0.1 0 0.1 Relative Errors 0.2
Figure 8. Histograms of the errors in the calibration with the HR and Spline volatility models.
αi α1 α2 α3 α4 α5 α6 α7 0.7791 0.4195 0.1308 0.0289 0.0050 9.9862 9.9927
Std. dev. (3.59e-2) (2.01e-3) (1.24e-4) (4.17e-5) (3.54e-4) (2.46e+5) (6.58e+5)
αi
0.0272 0.7114 0.0616
Table 2 Calibrated parameters for the volatility functions σspline and σHR .
3.3 Out-of-sample tests
In order to test out-of-sample performances, the HR model has been calibrated on the observations ranging from Nov-15-2002 to Jan-14-2003 and tested on the successive week, month and three months periods. The calibration produces the volatility function
AC C
and examples of the in- and out-of-sample fits are shown in Figure 10, where the calibrated payoffs for different maturities are represented by the curves ˆ and market prices fi by the crosses.
Statistics for in-sample and out-of-sample tests are given in Table 3 and RMSEs are plotted against time, deviation, maturity and log-moneyness in Figure 11. As expected, out-of-sample performances get worse as the horizon considered increases. Moreover, since the calibration is performed mainly on short 21
EP
σHR (D)2 = 0.0297 + 0.9360(D − 0.0352)2 ,
TE
DM
AN
US
Std. dev. (8.51e-5) (2.11e-3) (6.51e-4)
Spline
HR
IPT
0.3 0.4
Freq.
ACCEPTED MANUSCRIPT
Dev. from trend 0.5 0.45 Implied Volatility 0.4 0.35 0.3 0.25 0.2 0.15
Dev. from trend
Dev. from trend
Dev. from trend 0.5
!0.25 !0.2 !0.15 !0.1 !0.05 !0.25 !0.2 !0.15 !0.1 !0.05 !0.25 !0.2 !0.15 !0.1 !0.05 !0.25 !0.2 !0.15 !0.1 !0.05
!0.5 0 0.5 Moneyness
!0.5 0 0.5 Moneyness
!0.5 0 0.5 Moneyness
!0.5 0 0.5 Moneyness
02!Dec!2002 13:00:00 ! Dt = !0.13 250 Market Prices Calibrated Payoffs 200 250
AN
200 Call Prices 150 100 50 0
US
Figure 9. Effects of the deviation from the trend on marked implied volatilities. The √ implied volatilities are plotted against adjusted log-moneyness log(F/K)/ T − t and grouped by different ranges of the deviation as shown by the bar at the top of each panel.
27!Mar!2003 13:00:00 ! Dt = !0.16 Market Prices Calibrated Payoffs
Call Prices
150
100
50
0
!0.25 !0.2 !0.15 !0.1 !0.05 0 0.05 0.1 Moneyness: log(K/F)
DM
0.15 0.2 0.25
!0.25 !0.2 !0.15 !0.1 !0.05 0 0.05 0.1 Moneyness: log(K/F)
CR
0.25 0.2 0.15
0.15
AC C
EP
maturity options it does not give a very good prediction for longer maturities. Nevertheless the RMSEs in the first panel follow the same pattern of the RMSEs in Figure 6, with an exception at week 17. Table 3 also confirms the validity of the model which, contrary to local volatility models, seems to have a good performance even if it is not continuously re-calibrated. In particular, the quality of the fit seems to be quite acceptable even one month after the calibration period, confirming the theoretical arguments.
TE
Figure 10. Market prices and call price curves on December 2nd (in sample) and March 27th (out of sample). The calibration has been performed over the period November 15th 2002 - January 14th 2003.
22
IPT
0.45 0.4 0.35 0.3
0.2 0.25
ACCEPTED MANUSCRIPT
RMSE in-sample out-of-sample 1 week 1 month 3 months 1.58 2.23 2.69 1.71
MAE 1.34 1.18 1.49 1.75
Mean -0.46 -0.02 0.18 0.60
Std. Dev. 1.65 1.58 2.22 2.63
Minimum -4.28 -3.17 -3.97 -5.16
Maximum
Absolute outside errors RMSE in-sample out-of-sample 1 week 1 month 3 months 0.97 1.67 2.13 0.54 0.84 1.07 0.08 0.35 0.65 0.99 MAE 0.58 Mean -0.07 Std. Dev. 0.99
Minimum -3.29
US
-2.18 -5.06 -3.03 Minimum -20.17% - 5.33% -17.36% -33.17% Minimum -15.27% - 2.16% -14.75%
0.97 1.63 2.02
AN
Std. Dev. 4.06% 3.85% 5.72% 7.88% Std. Dev. 2.50% 2.20% 3.95%
Percentage errors RMSE in-sample out-of-sample 1 week 1 month 3 months 4.04% 5.84% 8.19% 4.06% MAE 2.56% 2.64% 3.78% Mean
DM
1.23% 1.15% 2.23% Mean 0.13% 0.71% 1.16%
-0.17%
4.86%
RMSE in-sample out-of-sample 1 week 1 month
TE
Percentage outside errors Maximum 30.35% 14.44% 36.77%
MAE
2.50% 2.31%
1.14%
AC C
3 months 6.45% 3.16% 2.01% 6.13% -28.16% 45.26% Table 3 In-sample and out-of-sample statistics for percentage errors εi and absolute errors εi of the calibration over the period from Jan-17 to Feb-17-2003. Out-of-sample ˆ statistics refers to the one week, one month and three months periods starting the 18-Feb-2003.
EP
1.16%
4.12%
2.08%
23
CR
Maximum 6.52 4.45 10.01 12.28
IPT
7.58 5.46 10.51 13.29 Maximum 33.14% 19.20% 39.71% 48.38%
Absolute errors
ACCEPTED MANUSCRIPT
Outside Errors 4 3 RMSE 2 1 0 in!sample out!of!sample
0
1
2
3
4
5
6
7
8
9 10 11 12 13 14 15 16 17 18 19 20 21 Week
2.5 RMSE 2 1.5 1 0.5 !0.24 !0.22 !0.2
!0.18 !0.16 Deviation from trend
AN
!0.14 4
US
!0.12 !0.1 0.5 1 1.5
3
5 4 RMSE 3 2 1 0
DM
2
2.5 3 3.5 Months to expiration
4 3 RMSE 2 1
AC C
Figure 11. Plot of RMSEs versus observation time (measured in weeks), deviation from the mean, time to maturity (measured in months) and moneyness. In-sample (from January 17th to February 16th) and out-of-sample (from February 17th to May 23th) RMSEs are represented by circles and squares, respectively.
EP
0 !0.35 !0.3 !0.25 !0.2 !0.15 !0.1 !0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 Moneyness log(S/K)
TE
24
CR
4.5 5 5.5
IPT
ACCEPTED MANUSCRIPT
4
Conclusions
A
Appendix: discretization and linear systems
EP
− Let us now consider the application of the discrete operator Y∆t in (30) to − the vector un . The generic element Y∆t un is the linear combination of the i,j n−1 n n−1 corresponding element in u and two elements, ui,j+k and un−1 . i,j+k+1 of u − Thus, applying Y∆t to un is equivalent to the difference of two linear operators, ∆−1 In and ∆−1 Zn , applied respectively to un and un−1 , where In denotes the t t − identity operator in RIJn . Specifically, the vector with elements Y∆t un is given i,j by
AC C
where Zn ∈ RIJn ×IJn−1 is the matrix such that the entry corresponding to the index i, j of Zun−1 is given by
n−1 (1 − γ)ui,j+k + γun−1 . i,j+k+1
TE
DM
25
Throughout this Appendix we use the notations of Subsection 2.3. The aim is to formulate the discretization of the partial differential equation (6) as a block bidiagonal linear system. We define I = 2i0 + 1, Jn = 2j0,n + 1 and denote by un ∈ RIJn the vector containing the values un for |i| ≤ i0 and i,j |j| ≤ j0,n : those values are sorted by the couple of indices (j, i) in lexicografic order.
1 n (u − Zn un−1 ), ∆t
AN
US
Aiming to test the effectiveness of the model we study its calibration as a nonlinear least squares problem. In this approach we estimate the volatility function both in a parametric and a semi-parametric framework. We perform an empirical experiment on a large dataset of intraday S&P 500 option quotations ranging from November 2002 to May 2003. The results show good performances of the model even when not continuously recalibrated. For this reason the model seems to be a valid alternative to the widely used local volatility models.
CR
In this paper we consider the path dependent volatility model proposed by Hobson and Rogers. In this model option prices can be expressed in terms of solution to a hypoelliptic parabolic PDE of Kolmogorov type. In order to exploit the intrinsic regularity properties of the degenerate PDE, we introduce a non-Euclidean finite differences scheme which allows for an accurate approximation of option prices.
IPT
ACCEPTED MANUSCRIPT
Then it turns out that the linear system (32) can be rewritten in matrix form
for 1 ≤ n ≤ N , where A ∈ RIJn ×IJn is the diagonal matrix with elements aij and
D 0 . . .
Dn =
and
ˇ D=
d2 d3 · · · .. .. .. . . .
CR
d1 d2 d3
ˇ 0 0 ··· D
ˇ 0 ··· 0 ˇ D ··· 0 . .. . . . . . .
0 0 d1 0 0
0 ··· 0
···
0 . . , .
are tridiagonal matrices of order IJn and I, respectively.
Similarly, combining the forward and backward schemes allows to derive the θ-method: θ∆t An Dn un + (1 − θ)∆t Zn An−1 Dn−1 un−1 − Zn un−1 + un = 0, or ¯1 ¯2 An un = An un−1 ,
AN
0 0 . . .
EP
Finally, by setting N = λT /∆t it turns out that considering (A.3) for n = 1, . . . , N and imposing the initial conditions u0 = v0 is equivalent to the linear system
I ¯1 −A2 . . . 0
TE
¯1 ¯2 with An = (In + θ∆t An Dn ) and An = Zn (In − (1 − θ)∆t An−1 Dn−1 ). As usual, the θ-method reduces to the explicit, implicit or Crank-Nicholson schemes when θ = 0, 1 or 0.5, respectively. Notice that the θ-method is unconditionally ¯ stable for 0.5 ≤ θ ≤ 1. The matrices An and Dn have an identical structure, 1 specifically they are block diagonal with tridiagonal blocks. Thus, the computational cost required to solve (A.3) is of the order of IJn . Furthermore, the structure of the matrices can be exploited to design computationally efficient and/or parallel algorithms for the solution of the PDE (6).
DM
0 0 .. . ··· ··· 0 ¯2 ¯1 −AN AN 26
1 ≤ n ≤ N,
US
u0 u1 . . .
0
0 ··· 0
0 0
0
AC C
¯ A1 1 .. .
¯N ¯1 · · · −A2 −1 AN −1 0 ···
0
N −1 u
=
uN
v 0 . . . 0
0
IPT
(In + ∆t An Dn )un − Zn un−1 = 0,
(A.1)
(A.2)
(A.3)
ACCEPTED MANUSCRIPT
or, with the appropriate substitution,
¯ where A is block bidiagonal and u ∈ Rq where ¯
N
q=
n=0
IJn .
References
[3] G. Barone-Adesi, H. Rasmussen, and C. Ravanelli, An option pricing formula for the garch diffusion model, Computational Statistics and Data Analysis, 49 (2005), pp. 287–310. [4] D. S. Bates, Post-’87 crash fears in the S&P 500 futures option market, J. Econometrics, 94 (2000), pp. 181–238. [5] V. Blaka Hallulli and T. Vargiolu, Financial models with dependence on the past: a survey, Applied and Industrial Mathematics in Italy, M. Primicerio, R. Spigler, V. Valente, editors, Series on Advances in Mathematics for Applied Sciences, World Scientific 2005, 69 (2005). [6] M. Broadie, M. Chernov, and M. Johannes, Model specification and risk premiums: The evidence fromthe futures options, tech. rep., Columbia University, 2005.
AC C
[7] C. Chiarella and K. Kwon, A complete Markovian stochastic volatility model in the HJM framework, Asia-Pacific Financial Markets, 7 (2000), pp. 293– 304. [8] T. Coleman and Y. Li, An interior, trust region approach for nonlinear minimization subject to bounds, SIAM Journal on Optimization, 6 (1996), pp. 418–445.
EP
TE
DM
27
[2] A. Amendola and G. Storti, A gmm procedure for combining volatility forecasts, Computational Statistics and Data Analysis, 52 (2008), pp. 3047– 3060.
AN
[1] Y. A¨ ıt-Sahalia and A. W. Lo, Nonparametric estimation of state-price densities implicit in financial asset prices, The Journal of Finance, LIII (1998), pp. 499–547.
US
Existence and uniqueness of the solution are related to the non-singularity of ∆ ¯ An and, in view of the expressions of d1 , d2 , d3 , are clearly guaranteed if ∆ t2 is 1 x suitably small. Stability and numerical stability are driven by the properties ¯2 ¯1 of the matrices An and An .
CR
IPT
(A.5)
¯u ¯ A¯ = v ,
(A.4)
ACCEPTED MANUSCRIPT
[9] F. Corielli, P. Foschi, and A. Pascucci, Analytical approximation for financial derivatives, tech. rep., University of Bologna, 2008. [10] M. Di Francesco, P. Foschi, and A. Pascucci, Analysis of an uncertan volatility model, Journal of Applied Mathematics and Decision Sciences, 2006 (2006), pp. 1–17.
[12] [13]
, On a class of degenerate parabolic equations of Kolmogorov type, AMRX Appl. Math. Res. Express, 3 (2005), pp. 77–116. , A continuous dependence result for ultra-parabolic equations, J. Math. Anal. Appl., 336 (2007), pp. 1026–1041.
[14] M. Di Francesco, A. Pascucci, and S. Polidoro, The obstacle problem for a class of hypoelliptic ultraparabolic equations, Proc. R. Soc. Lond. A, 464 (2008), pp. 155–176.
[16] B. Dumas, J. Fleming, and R. E. Whaley, Implied volatility functions: empirical tests, J. Finance, 53 (1998), pp. 2059–2106.
[18] J.-P. Foque, G. Papanicolau, K. Sircar, and K. Solna, Maturity cycles in S&P 500 volatility, Journal of Computational Finance, 6 (2003), pp. 1648– 1665.
[20] F. Fritsch and R. Carlson, Monotone piecewise cubic interpolation, SIAM Journal on Numerical Analysis, 17 (1980), pp. 238–246.
AC C
¨ [22] M. Hahn, W. Putschogl, and J. Sass, Parameter estimation for stock models with non-constant volatility using markov chain monte carlo methods, Operations Research Proceedings, Springer Berlin Heidelberg, (2007), pp. 227– 232. [23] , Portfolio optimization with non-constant volatility and partial information, Brazilian Journal of Probability and Statistics, 21 (2007), pp. 27– 61.
EP
[21] R. Garcia, E. Ghysels, and E. Renault, The econometrics of option pricing, in Handbook of Financial Econometrics, Y. Ait-Sahalia and L. Hansen, eds., 2005.
TE
[19] P. Foschi and A. Pascucci, Path dependent volatility, Decis. Econ. Finance, 31–1 (2007), pp. 1–20.
DM
28
` [17] G. Figa-Talamanca and M. L. Guerra, Fitting prices with a complete model, Journal of Banking and Finance, 20 (2006).
AN
[15] E. DiBenedetto, Partial differential equations, Birkh¨user Boston Inc., a Boston, MA, 1995.
US
CR
[11] M. Di Francesco and A. Pascucci, On the complete model with stochastic volatility by Hobson and Rogers, Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci., 460 (2004), pp. 3327–3338.
IPT
ACCEPTED MANUSCRIPT
[24] D. G. Hobson and L. C. G. Rogers, Complete models with stochastic volatility, Math. Finance, 8 (1998), pp. 27–48. [25] F. Hubalek, J. Teichmann, and R. Tompkins, Flexible complete models with stochastic volatility generalising Hobson-Rogers, working paper, (2005).
¨ ¨ ´ [27] E. Lindstrom, J. Strojby, M. Broden, M. Wiktorsson, and J. Holst, Sequential calibration of options, Computational Statistics and Data Analysis, 52 (2008), pp. 2877–2891.
[29] A. Pascucci, Free boundary and optimal stopping problems for American Asian options, Finance and Stochastics, XII-1 (2008), pp. 21–41. [30] A. Platania and L. C. G. Rogers, Putting the Hobson&Rogers model to the test, working paper, (2003). [31] J. Sekine, Marginal distribution of some path-dependent stochastic volatility model, preprint, (2008). [32] S. Sundaresan, Continuous-time methods in finance: A review and an assessment, Journal of Finance, 55 (2000), pp. 1569–1622. [33] D. Tavella and C. Randall, Pricing Financial Instruments: The Finite Difference Method, John Wiley & Sons, Inc., 2000. [34] A. Trifi, Issues of aggregation over time of conditional heteroscedastic volatility models: What kind of diffusion do we recover?, Stud. Nonlinear Dyn. Econom., 10 (4), Article 5 (2006). [35] P. Wilmott, S. Howison, and J. Dewynne, Option pricing, Oxford Financial Press, Oxford, 1993.
AC C
[37] J. Yu, Z. Yang, and X. Zhang, A class of nonlinear stochastic volatility models and its implications for pricing currency options, Computational Statistics and Data Analysis, 51 (2006), pp. 2218–2231.
EP
[36] H. Windcliff, P. Forsyth, and K. Vetzal, Analysis of the stability of the linear boundary condition for the black-scholes equation, Journal of Computational Finance, 8 (2004), pp. 65–92.
TE
DM
29
AN
US
[28] M. D. Marcozzi, On the approximation of optimal stopping problems with application to financial mathematics, SIAM Journal on Scientific Computing, 22 (2001), pp. 1865–1884.
CR
[26] E. Lanconelli, A. Pascucci, and S. Polidoro, Linear and nonlinear ultraparabolic equations of Kolmogorov type arising in diffusion theory and in finance, in Nonlinear problems in mathematical physics and related topics, II, vol. 2 of Int. Math. Ser. (N. Y.), Kluwer/Plenum, New York, 2002, pp. 243–265.
IPT