### PAPER:

The calibration of the Heston stochastic volatility model using filtering and maximum likelihood methods

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.

1. Abstract

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.

2. Introduction

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>0. Once known the forecasted values of the stock price and of the stochastic variance, of the Heston parameters and of the risk premium parameter we can forecast the value of the option price using a formula derived in [1] Section 3 . This formula  involves only the numerical evaluation of  a one dimensional integral.

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

### 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 downloada 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 downloada 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.

5. References

[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).

Entry n. 8237