PAPER:
The calibration of the Heston stochastic volatility model using
filtering and maximum likelihood methods
Lorella Fatone, Francesca Mariani, Maria Cristina Recchioni, Francesco Zirilli
- Abstract
- Introduction
- Estimation
problem using synthetic data Estimation Movie
1, Estimation
Movie 2
- Estimation
problem and forecasting experiments using real data Forecasting Movie 1, Forecasting Movie 2
- References
Warning: If your browser does not allow you to see the content of this website click
here to download a
word file that reproduces the material contained here.
In [1] we study the
problem of obtaining accurate estimates of the parameters, of the initial
stochastic variance and of the risk premium parameter of the risk neutral measure
of the Heston stochastic volatility model from the observation at discrete
times of the stock log-returns and of the prices of a European call option on
the stock. This problem is an inverse problem known in the literature as
calibration problem. As a byproduct of the solution of the calibration problem
we develop a tracking procedure that can be used to forecast the stochastic
variance and the stock price. From a mathematical point of view the problem
considered is formulated as a constrained optimization problem where the
objective function is the logarithm of the likelihood associated to the
parameter, the initial stochastic variance and the risk premium parameter
values given the observed stock log-returns, option prices and observation
times, this function is called (log-)likelihood function. The evaluation
of the (log-)likelihood function associated to a given choice of the
parameter, of the initial stochastic variance and of the risk premium parameter
values requires the solution of a filtering problem for the Heston model. An
accurate and computationally efficient solution of this filtering problem is
necessary for a satisfactory solution of the calibration problem for the Heston
model. A similar problem has been considered in [2] and [3]. The aim of this paper is to extend and improve the
results obtained in [2] with respect to the formulation of
the problem, the accuracy of the solution obtained and the computational efficiency
of the solution method. In particular in comparison with [2]
we reformulate the problem adding to the quantities that must be estimated the
risk premium parameter and to the observations the option price at the initial time.
The addition of the risk premium parameter to the quantities that must be
estimated makes the problem considered more realistic and makes interesting the
analysis of time series of real financial data using the method proposed. The
addition of the option price at the initial time to the data used in the
solution of the problem is natural and improves significantly the estimate of
the initial stochastic variance. Moreover we simplify the expression of some
formulae used in [2] in the computation of the
(log-)likelihood function and we improve the optimization method employed to
solve the maximum likelihood problem introducing some ad hoc preliminary
optimization steps. Finally a new, easy to compute, formula that gives the ``Heston"
option price is derived. The use of this new formula reduces substantially the
computational cost of evaluating the (log-)likelihood function when
compared to the cost of the (log-)likelihood function evaluation in [2]. Some numerical examples of the calibration problem using
synthetic and real data are presented. As real data we use the 2005 historical
data (precisely the data of the period Jan. 3, 2005, May 11, 2005) of the US
S&P500 index and of the corresponding option prices. Very good results are
obtained forecasting future values of the S&P500 index and of the
corresponding option price using the solution of the calibration problem and
the forecasting and tracking procedure mentioned above. In this website some
auxiliary material useful in the understanding of [1]
including some animations and some numerical experiments is shown. A more
general reference to the work of the authors and of their coauthors in
mathematical finance is the website: http://www.econ.univpm.it/recchioni/finance.
The
problem that we are considering is the
calibration of the Heston stochastic volatility model including the estimation
of the initial stochastic variance and of
the risk premium of volatility implicit in option prices. The main differences
between the work contained in [1] and the work presented in [2] are: i)
the addition to the quantities that must be estimated of the risk premium
parameter contained in the risk neutral measure of the Heston model and the
addition to the observations used in the estimation problem of the observation
at the initial time of the option price, ii) the reduction of the formula used
to compute the ``Heston" option price from being a three dimensional
integral that must be evaluated numerically to being a one dimensional integral
that must be evaluated numerically, iii) an improved optimization method that
includes some ad hoc preliminary steps that is used to maximize the likelihood
function. These differences lead to a substantial improvement of the accuracy
of the estimates obtained of the model parameters, of the initial stochastic
variance and of the risk premium parameter and to great savings in the
computational cost of solving the calibration problem. Note that the addition
of the risk neutral parameter to the quantities that must be estimated makes
the problem considered more relevant from the practical point of view and
together with the improved computational efficiency in the solution of the
problem makes possible the use of the method proposed to analyze time series of
real financial data. The work presented in [2] corresponds
to choose the risk premium parameter equal to zero. In particular we show that
the addition of the initial option price to the data and the improved
optimization method determine a substantial improvement of the quality of the
estimate of the initial stochastic variance obtained and, consequently, an
improvement of the quality of the estimates of all the remaining quantities. In
fact the numerical experience shows that the accuracy of the estimate of the
initial stochastic variance is crucial to obtain a good solution of the
calibration problem and satisfactory forecasted values of the stock and option
prices. On the other hand the numerical evaluation of a one dimensional
integral in order to compute the option prices instead than the three dimensional
integral evaluated numerically in [2] is responsible for
the substantial reduction of the computational time needed to solve the
calibration problem. The improved numerical efficiency makes practical the use
of observation samples bigger than those used in [2] and
this gives an extra possibility of improving the quality of the results
obtained.
Let
us give some details. Assuming the stock prices and their stochastic
volatilities or variances as state variables and the stock prices and
option prices observed at discrete times as observations, we study the
problem of forecasting (tracking) the unobserved and time-varying stochastic
volatility or variance when the model parameters, initial condition and risk
premium parameter are known using the observations (filtering problem)
and the problem of estimating the model parameters, the unknown initial
stochastic volatility or variance and the risk premium parameter using the
observations of the stock price and of the price of a European call option on
the stock (estimation problem).
Let
St and vt be the state variables at time t>0, that is,
respectively, the stock price and the stochastic variance at time t>0, the Heston model prescribes the following dynamics for the
state variables (St,vt), t>0 (see [4]):
where Wt1, Wt2,
t>0, are standard Wiener processes such that W01=W02=0, dWt1, dWt2
are their stochastic differentials and <dWt1dWt2
>= rdt where < × > denotes the mean of × and rÎ[-1,1] is a constant known as correlation parameter; m is a real constant related to the behaviour of
the mean of the stock log returns, associated to the stochastic process St
, t>0. Note that the stochastic
variance vt, t>0, is
assumed to be described by a mean-reverting process with speed g³0 and parameters q³0 and e³0, that is by equation (2). The equations (1),
(2) must be equipped with an initial condition, that is:
Here 
and 
denote random variables concentrated
with probability one in a point; for simplicity, we continue to denote these
points respectively with and .
We note that , the stock price at t=0, is
observable while , the stochastic variance at t=0,
cannot be observed. The quantities m, q, e, g
and r are real constants and are the
model parameters that we want to estimate. Hence, the parameters and unknown
initial condition component of the Heston model that we want to estimate
are m, q, e, g, r,, moreover we want to estimate the risk premium parameter of the risk neutral
measure l. The meaning of the parameter l is explained below. We assume that 2gq/e2 ³1, this assumption guarantees that
if is positive with probability one, vt, solution of
(1), (2), remains positive with probability one for t>0.
We suppose to have the observations of the
prices of the stock and the prices of a European call option of known maturity
and strike price on the underlying whose
price St, t>0, satisfies (1), (2), (3), (4) at prescribed
discrete times t0= 0<t1<t2<…<tn<tn+1
=+¥, let Ft= {
: ti £t}, t>0, be the set of the observations up
to time t>0 of the stock prices St and of the option prices Ct
that is for i=1,2,...,n let 
be the observation of
the stock price and 
be the observation of
the option price at time ti and let 
be the vector made of
the model parameters, of the initial stochastic variance and of the risk premium
parameter. We note that since we consider option prices as data and we evaluate
the option prices under the risk neutral measure we add to the parameters and
to the initial condition component of
the Heston model (1), (2), (3), (4) to be estimated the risk premium parameter l [4] that appears in the risk
neutral measure used to price the options. Note that as shown in [4]
the risk neutral measure is obtained
from the statistical measure associated to (1), (2), (3), (4) with a simple
explicit transformation. We use the notation t0 =0 and F0={
}, this notation simplifies some of the formulae that follow.
The option prices considered refer to the prices of an European call option on
the underlying whose prices are described by St, t>0, with known
strike price E and maturity time T such that T>tn, we denote with
C(S,v,t;E,T,Q) the price under the risk neutral
measure of this option in the Heston
stochastic volatility model when 0£t<T, St=S and vt=v. We assume that the option prices observed are
given by the theoretical option prices in the Heston model affected by a
Gaussian error, that is: for i=1,2,...,n, =C(, 
,ti; E,T,Q)+ui, where ui is an additive Gaussian error
term; that is we assume that ui is sampled from a Gaussian variable
with mean zero and known variance 
, and C(,,ti; E,T,Q) is the option price in the Heston model
assuming known the (unobserved) variance and the (observed)
stock price at time ti.
Note that we assume and known without error.
Let us briefly describe how we get the
forecasted values of the state variables 
, 
, t>0 (see [1], [2] for further details). Let us introduce the stock log returns:
(5) without loss of generality, we can assume S0==1. From (1), (2), (3), (4) and (5) it is easy to see that (xt,vt),
t>0, satisfy the following equations:
with
initial condition:
(8)
(9)
Given , we can compute the
joint probability density function 
(x,v)Î(-¥,+¥)
(0, +¥) , t>0, of the state variables 
, that is of the stock
log return and of the stochastic variance, conditioned to the observations Ft,
t>0, and to the initial condition (8), (9)
as solution of the equation presented in [1], [2] , that is: for
i=0,1,...,n:
(10)
with
initial condition: for i=0:
(11)
and for
i=1,2,...,n:
(x,v)Î(-¥,+¥)( 0,+¥) .
(12)
where for i=1,2,...,n: 
=log(Si/) and 
(x,v,ti;E,T, Q) is the value in the Heston model
at time t=ti (see [1] Section 3) when 
, 
of a call option on
the stock whose price is St=exp(xt), t>0, with strike price E and maturity time
T.
Once we know the joint probability density 
(x,v,t|Ft,Q), (x,v) Î (-¥,+¥)(0, +¥), t>0, solution of (10), (11), (12) we can
forecast the values of the random variables (xt,vt),
t>0, using the formulae:
The
quantities 
are used to forecast
the stock prices i.e. we forecast with 
, t>0, and the stochastic variance i.e. we forecast with 
, t>
In order to reconstruct the vector from a set of
observations made at time t=ti,
i=0,1,…,n, we use the maximum likelihood method that consists in searching the
vector QÎM that maximizes the following (log-)likelihood function:
(15)
where 
is the set of the admissible vectors Q and 
, i=0,1,…,n-1, is the following function:
(16)
We note that (15) represents the sum of the
logarithms of the integrals of the product between the transition probability
densities 
and the weight
functions 
, i=0,1,...,n-1, and of the logarithm of 
. The use of these weight functions in the expression
of the likelihood function (15) has the
purpose of exploiting the fact that also the option prices , i=0,1,...,n, are
observed, assigning bigger weights to the values of the v variable that make more
likely the observation of at time t=ti,
i=0,1,...,n. The function F(Q) represents the (log-)likelihood of the parameter
vector Q
given the observations Ft, t>0, and 
, 
at time t=t0.
A solution of the estimation problem (i.e. the
problem of calibrating the parameters of Heston model (1), (2), (3), (4), the
stochastic initial variance and the risk premium parameter l) can be obtained solving the following problem:
(17)
using an iterative method (such as, for
example, a variable metric steepest ascent method) starting from an initial
guess Q0ÎM. Problem (17) is the maximum (log)-likelihood
problem that has been announced.
Further
details can be found in the website http://www.econ.univpm.it/pacelli/mariani/finance/w1
3. Estimation problem using
synthetic data
In this Section
we describe some numerical experiments obtained applying the numerical method
proposed in [1] and summarized in Section 2
to reconstruct the vector Q on some test problems. We remind that once known Q, we can forecast the stock prices, the stochastic
variance and the option prices solving the filtering problem, that is using
formulae (13), (14) and the appropriate formula for the option price.
The first numerical experiment consists in the calibration
of the Heston model through the solution of
problem (17) . We use synthetic data. The time unit is a ``year"
made of 252.5 trading days. We consider
a temporal interval of 36 consecutive trading days and n=9 observation times in
this time period. Note that we choose t0=0, and the observation
times ti equally spaced, that is ti=ti-1+D,
i=1,2,...,9, where D=4/252.5. The choice of the vector Q that generates the trajectory of
the initial value problem (6), (7), (8), (9) used to generate the data is: =
. Moreover we choose =0 and =10-4, i=0,1,…,9. We remind that the vector 
is expressed in annualized
unit. Similarly all the remaining choices of the vector 
made in
the numerical experiments are expressed in annualized unit. Table 1 shows the
synthetic data generated integrating numerically one trajectory of (6), (7),
(8), (9) and using the appropriate change of variable and option price formula.
Starting from 
(0.5,8.22,0.15,0.0067,-0.634,0.3,0.01)
, we solve the optimization
problem (17) associated to the data of Table 1. Note that we have: 
|
i |
ti |
|
|
|
0 |
0 |
0 |
0.5099 |
|
1 |
4/252.5 |
1.2355e-1 |
0.6410 |
|
2 |
8/252.5 |
6.8935e-2 |
0.5804 |
|
3 |
12/252.5 |
-2.9809e-2 |
0.4791 |
|
4 |
16/252.5 |
1.2049e-2 |
0.5201 |
|
5 |
20/252.5 |
4.4983e-2 |
0.5538 |
|
6 |
24/252.5 |
-7.7662e-4 |
0.5065 |
|
7 |
28/252.5 |
-1.3411e-2 |
0.4936 |
|
8 |
32/252.5 |
-1.2461e-1 |
0.3894 |
|
9 |
36/252.5 |
-7.1041e-2 |
0.4377 |
Table 1. The synthetic data of the first
experiment
The numerical procedure evaluates the objective
function (15) in a point 
solving the filtering problem (10), (11), (12) and evaluates
the gradient of the objective function (15) in a point using elementary finite
differences formulae so that to compute the function (15) and its seven
dimensional gradient in a point it is necessary to
solve eight filtering problems.
The iterative procedure used to solve the
optimization problem (17) is a variable metric steepest ascent algorithm that
generates a sequence of vectors 
k=1,2,..,itmax, and in
our numerical experience we fix the maximum number of iterations itmax=200. We
use the following stopping criterion:
where
is a given tolerance. In this experiment we choose =10-6 and the maximization procedure that starts
at 
stops at the k=92-th iteration. In Table 2 we
show the values of F(
), and of 
, when k=0, 92
and we show the value of F(). Note that the maximum likelihood procedure finds as estimated value of the vector
the vector 
.
|
k |
F( |
|
|
0 |
-929.14 |
0.392 |
|
92 |
43.17 |
0.041 |
Table 2. Value of the
(log)-likelihood function and relative distance between and 
when 
, k=0, 92
Figure 1 shows the tracking of the stochastic variance
obtained with formula (14) (dashed dotted line) when =
and when =
, the (blue) solid line represents the “true trajectory”
obtained integrating numerically one trajectory of (6), (7), (8), (9) with 
. The substantial improvement obtained in the tracking of the
stochastic variance passing from 
to 
shows that the maximum
likelihood problem (17) is a valid formulation of the calibration problem.
Figure 1. Tracking of the
stochastic variance
Estimation
Movie 1, click here to download a digital movie showing the tracking
of the stochastic variance trajectory when
Q = Qk ÎM , k=0,1,…,92, remind that these Q are vectors obtained during the iterative optimization
procedure used to solve problem (17) with the data of Table 1. Moreover the
movie shows the true trajectory of the stochastic variance corresponding to the parameter vector Q=
obtained integrating numerically one trajectory of (6), (7), (8), (9), this trajectory is the one
used to generate the data of Table 1.
The second experiment consists in analyzing a series of daily observations during a
period of two years, that is a series made of 505 observation times. Remind
that we have considered a “year” made of 252.5 trading days. We consider 505
observation times corresponding to 1010 observations 
,, t=ti=i/252.5, i=0,1,…,504, where we have chosen
t0=0, x0=
=0, v0=
=0.5, T=2.5, E=0.5, 
, i=0,1,2,…,504. The data are generated integrating
numerically with Euler method one trajectory of (6), (7), (8), (9) for the two
years period. We choose the vector 
associated to the
stochastic differential system that generates the data to be equal to the
vector 
in the first year of
data (i.e. when t=ti , i=0,1,…,252) and to the vector 
in the second year of
data (i.e. when t=ti, i=253,254,…,504) where 
and 
are given by: 
=(0.026,5.94,0.306,0.01159,-0.576,0.5,0),and
=(0.4,2,0.01,0.01,0.5,0.012,0.7)
.
Click
on the names to download the file
containing the observations (samplex.txt)
, (sampleC.txt), made at time ti=i/252.2,
i=0,1,…,504 used in the numerical experiment.
We solve problem (17) using 8 consecutive
observations that is the observations ,, made at time t=ti,j i=0,1,…,7 when j=0,1,…,62, that
is for j=0,1,…,62 the time values ti,j=t8j+i, i=0,1,2,…,7
belong to the same time window j. We solve 63 problem (17) corresponding to the
data associated to the data to the 63 time windows described above starting
from different initial guesses. In Figure 2 the diamond represents the vector and the square represents the vector . The unit
length in Figure 2 is 
. The origin is the mid point of the segment
that joins the points and and the length of the segment that joins the
diamond and the square is the Euclidean
norm of the vector - . The vector
determined solving problem (17) using the optimization procedure are
represented by stars. A point P=Q is plotted around when the quantity 
is smaller than the
quantity 
, otherwise it is plotted around . The point
P belonging to the cluster associated to is plotted on the left or on the right of the
point that represents according to the sign of the first component of 
. A similar statement
holds for 
and for a point P=
belonging
to the cluster associated to .
Figure 2. Reconstruction of the parameter
vectors and
Estimation
Movie 2, click here to download a
digital movie showing the points determined solving the 63 problem (17) described above as a function
of the time window of the observations.
Figure
2 and the Estimation Movie 2 show, as expected, that
the points obtained solving problem
(17) are concentrated on two disjoint segments one on the left and one
on the right of the origin and that they form two disjoint clusters, that is,
the solution of the 63 problems (17) corresponding to the time windows
described previously shows that two sets
of parameters seem to generate the data studied. Moreover as expected the
points obtained using data contained in the first year of observation cluster
around and the points obtained using data contained in the second year of
observation cluster around 
.
4. Estimation problem and
forecasting experiments using real data
In this section we show some numerical results
obtained applying the calibration procedure to real data.We analyze the daily
closing values of the U.S. S&P 500 index and the corresponding bid prices of a European call option on the
U.S. S&P 500 index with maturity date December 16, 2005 and strike price E=1200 during the period of about four months going from
January 3, 2005 to May 11, 2005. Due to
the number of data actually available in 2005 the time unit is a “year” made of 253 trading
days. Note that we have selected an option that during the period considered is
at the money and that this option has a significant volume traded. We have chosen
as initial time t=t0=0 the day January 3, 2005 and we have
considered 81 windows of 10 daily consecutive observation times. In particular
the observation times are chosen as ti,j=(i+j)/253, i=0,1,…,9, j=0,1,…,80
that is in each window we consider daily observations and when we move from a
window to the next one we discard the observations made in the first day of the
window and we add the next daily observations
as last observations to build the new
window. In the observation times t=ti,j , i=0,1,…,9, j=0,1,…,80 the
index j labels the window and the index i labels the observation dates in the window
j. In each window we solve problem (17) choosing =5.10-4 and we choose the starting
point taking into account the implied volatility associated to the option
prices, the standard deviation and the mean value of the historical log-returns
data. The calibration procedure, applied to the eighty-one windows considered
provides eighty-one parameter vectors. We use the parameter vector obtained using
the data contained in the j-th window (i.e. observation times ti,j =(i+j)/253, i=0,1,…,9) to forecast
the value of the S&P500 index and the bid price of the option the day next to the last observation time of the j-th window (see
Figures 3 and 4), that is we forecast St, C at time t=(j+10)/253. This
is done for the j-th window, j=0,1,…,80. The forecasts are made using formula
(13) and a similar formula for the option price (see [1]
Section 3 and Section 4 ). We note that in the eighty-one calibration and
forecasting problems solved the mean relative error made in the forecasts of
the “next day” S&P500 index value
and option bid price are 0.0063 and 0.0785 respectively. We note that in order
to forecast the option price we use the value of the stochastic variance
obtained with the tracking procedure, that is obtained with formula (14).
Figure 3. Tracking of the U.S.
S&P 500 index
Figure 4. Tracking of the bid
option price on the U.S. S&P 500 index
Figures 3 and 4 show the forecasted and the observed values of the index and of the bid price of the option versus time
(expressed in days) respectively. The solid line and the square represent the
“true” values provided by the historical
data while the dotted line and the stars are the forecasted values obtained
using the procedure described previously.
The
forecasts shown in Figure 3 and Figure 4 are “next day” forecasts. We consider
also forecasts two, three, four and five days in the future obtained with the
vectors estimated previously.
In Table 2 we show the mean relative errors made in the forecast of the
S&P500 index value eindex and of the bid price of the option eoption
obtained forecasting up to five days in the
future. The mean relative errors are computed averaging the errors made in the
solution of the eighty-one minus the number
of days in the future of the forecasts problems considered. For example when we make “next
day” forecasts the average is made on eighty values, when we made “two days” in
the future forecasts the average is made on seventy-nine values and so on. The
errors are computed comparing the forecasts with the historical data.
|
Days in the future |
eindex |
eoption |
|
1 |
0.0063 |
0.0785 |
|
2 |
0.0083 |
0.1010 |
|
3 |
0.0097 |
0.1181 |
|
4 |
0.0109 |
0.1296 |
|
5 |
0.0118 |
0.1435 |
Table 3. Mean relative errors on forecasted values
of the SP&500 and of the corresponding option price
Forecasting Movie 1,
click here to download a digital movie showing the
one and four days in the future forecasted values of the SP&500 index as function of time. The blue boxes are the “true”
(historical) values, the red stars are the one day in the future forecasted
values and the green stars are the four days in the future forecasted values.
Figure 5
shows the true bid option prices on the US S&P500 index provided by the
historical data (blue solid line and squares) and the bid prices of the option
forecasted as explained below.
Figure 5. Daily estimates
of bid option prices on the US S&P
500 index
Forecasting Movie 2, click here to download a digital movie showing the bid option prices forecasted as explained
below as a function of time. The blue
boxes are the “true” (historical) values and the red stars are the forecasted
values.
Finally
Figure 5 and the Forecasting Movie 2 show the bid prices of the option computed
using the value of the stock price actually observed and the value of the
stochastic variance obtained as one day in the future forecast derived from the
parameters obtained solving the
calibration problem with the data observed in the previous ten days. In this case mean relative error made is e1option=0.0298.
[1] L. Fatone, F. Mariani, M.C. Recchioni, F.
Zirilli: “The calibration of the Heston stochastic volatility model using
filtering and maximum likelihood methods”, in Proceedings of Dynamic Systems and Applications, G.S.Ladde, N.G.Medhin,
Chuang Peng, M.Sambandham Editors, Dynamic Publishers, Atlanta, USA, 5, 170-181 (2008).
[2] F. Mariani, G. Pacelli, F. Zirilli:
"Maximum likelihood estimation of the Heston stochastic volatility model
using asset and option prices: an
application of nonlinear filtering theory", Optimization Letters, 2, 177-222 (2008). (http://www.econ.univpm.it/pacelli/mariani/finance/w1 contains downloadable software).
[3] L.
Fatone, F. Mariani, M.C. Recchioni, F. Zirilli: " Maximum likelihood
estimation of the parameters of a system of stochastic differential equations
that models the returns of the index of some classes of hedge funds",
Journal of Inverse and Ill-Posed Problems, 15, 329-362 (2007). (http://www.econ.univpm.it/recchioni/finance/w5
contains downloadable software).
[4] S.L. Heston, “A closed-form solution for
options with stochastic volatility with applications to bond and currency
options”, Review of Financial Studies, 6,
327-343 (1993).