ACCURATE PARAMETER ESTIMATION FOR COUPLED STOCHASTIC DYNAMICS

A BSTRACT . We develop and implement an efﬁcient algorithm to estimate the 5 parameters of Heston’s model from arbitrary given series of joint observations for the stock price and volatility. We consider the time interval T separating two observations to be unknown and estimate it from the data, thereby estimating 6 parameters with a clear gain in ﬁt accuracy. We compare the maximum likelihood parameter estimates based on a Euler discretization scheme to analogous estimates derived from the more accurate Milstein discretization scheme; we derive explicit conditions under which the two set of estimates are asymptotically equivalent, and we compute the asymptotic distribution of the difference of the two set of estimates. We show that parameter estimates derived from the Euler scheme by constrained optimization of the approximate maximum likelihood are consistent, and we compute their asymptotic variances. Numerically, our estimation algorithms are easy to implement,and require only very moderate amounts of CPU. We have performed extensive simulations which show that for standard range of the process parameters, the empirical variances of our parameter estimates are correctly approximated by their theoretical asymptotic variances.


1.
Introduction.Derivatives such as futures and options were introduced to be used as hedging tools, and are now heavily traded in the market.Trading prices for these derivatives are indicators of market expectations for the near future.Due to high volatilities, efficient hedging tools require accurate pricing of these derivatives.To evaluate option prices through robust model based inference from asset dynamics data, it is crucial that the underlying joint stochastic dynamics of the asset price and volatility be calibrated correctly.Our objective here is to propose a fast and robust parameter estimation algorithm for such pairs of coupled SDEs, with good quantitative control for the accuracy of these estimates.
Stock price volatilities are notoriously not constant in time [22].Local volatility models consider the volatility to be dependent on time and the underlying asset [11] [12] [24].Stochastic volatility models drive volatility and price dynamics by correlated Brownian processes [17][25] [26][27] [28].Various methods have been proposed for estimation of stochastic volatility models [1][4] [8][14] [29].[1] employs Maximum Likelihood using a Hermite approximation of the likelihood function.[4] develops a weighted non parametric approach to determine the risk-neutral measure of the future states of the market.Bayesian estimation is proposed in [29].We focus our analysis on parameter estimation for the often used pair of coupled SDEs introduced by Heston to price options [17].Feller's "square root" SDE [15] for volatility was first applied by Cox, Ingersoll and Ross to model short term interest rates [9].We do not study here the pragmatic suitability of such SDEs to model the joint dynamics of stock price and volatility, but focus on generating easily computable parameter estimates from finite sets of observed data, and on the simultaneous evaluation of the estimates' accuracy.In this context, there is no closed form expression for the log-likelihood or the joint density of price and volatility.
Heston's pair of coupled SDEs for price and volatility involve 5 unknown parameters.We derive robust parameter estimates from an approximate maximum likelihood principle applied to the Euler discretization scheme for SDEs.The time interval T between consecutive observations, is not fixed a priori but also simultaneously estimated from the data.
We evaluate asymptotic variances for these consistent parameters estimates, and compare their accuracy to a second group of estimates based on the more elaborate Milstein discretization scheme for SDEs.The Milstein scheme approximates underlying diffusion processes with L 2 norm accuracy equivalent to T, while the Euler scheme accuracy is of the order of √ T [21].We derive explicit conditions to force the Euler scheme parameter estimates to remain very close to the Milstein scheme estimators, and we compute the asymptotic distribution of the differences between these two types of estimators.We show that the estimators obtained are consistent.Numerically, our estimation algorithms are easy to implement,and require only very moderate amounts of CPU.Our extensive simulations show that for standard ranges of the process parameters, the empirical variances of our parameter estimates are correctly approximated by their theoretical asymptotic variances.
2. Heston's coupled SDEs for volatility and asset price.Consider Heston's classical coupled SDE system (under market measure) [17] At time t, Y t = exp(S t ) is the asset price and √ X t is its volatility.The squared volatility X t is a mean reverting square root process [9] with mean reversion speed κ > 0 and volatility γ.Processes W 1 (t), W 2 (t) are correlated Brownian motions with E[dW 1 (t)dW 2 (t)] = ρ dt.Parameters µ and θ are respectively the mean rate of return of the asset under market measure, and its long run average variance.Existence and uniqueness of the solution for these SDEs is well known [16].Denote the vector of model parameters by PAR = (µ, κ, θ, γ, ρ) and assume the following classical constraints to be satisfied The last constraint 2κθ > γ 2 is Feller's square root condition [15] ensuring that X t stays a.s.away from 0, so that X 0 > 0, then forces X t to remain positive.Let The conditional density of X t given X s is a non central χ 2 with 2m degrees of freedom.Its steady state density p(x) is given by ν m Γ(m) [9].
By integrating (2) and taking expectations, one easily obtains explicit expressions for E[X t ] and E[X 2  t ], which show that as t → ∞ the mean and variance of X t converge at exponential speed towards θ and γ 2 θ 2κ respectively.
The joint density of {S t , X t }, has no explicit closed form.We want to estimate the 5 model parameters given a finite series of joint observations for the stock price {Y t } and square of volatility {X t }.Note that in practice the volatility X t is not directly observed.Good practical estimates of X t have been proposed and used To focus on parameter estimation, we will assume here that the volatility is directly observed.
Let (U 0 , U 1 , . . .U N+1 ) and (V 0 , V 1 , . . .V N+1 ) be N + 2 joint observations for the log of asset price and the square of volatility at time points We assume that the time interval T = t n+1 − t n is fixed but unknown, to inject an adjustable time scale in model fitting.We want to select optimal T and PAR to achieve the best fit of the data by Heston's SDE system.
3. Approximate log-likelihood based on Euler Discretization.Euler discretization scheme for the SDE system provides the following difference equations, linking the one step differences ∆U n , ∆V n , ∆W 1 (n) and ∆W 2 (n), Due to the presence of a square root term, the coefficients of our SDEs do not satisfy a global Lipschitz condition and therefore standard results on L 2 convergence for the Euler discretization do not apply here.However one can show that, under the parametric restrictions stated above, if the process and its Euler discretized approximation are both killed at the first (random) time the volatility becomes smaller than a fixed ǫ > 0, then weak convergence of the Euler approximation scheme will hold in the associated space of continuous paths, endowed with an adequate natural metric.In particular, in order to ensure that V n remains in (0, ∞), sufficiently small discretization step sizes have to be imposed.The recent paper [19] develops efficient modified Ito-Taylor schemes to overcome both these difficulties, but we have not yet implemented their interesting approximation scheme.From (3) we get The independent random vectors Z n+1 = (∆W 1 (n), γ∆W 2 (n)) are Gaussian with zero mean, variances T and γ 2 T, and covariance ργT.Denote the coordinates of Z n+1 by (z 1 (n + 1), z 2 (n + 1)).The joint density of {Z 1 }, ..., {Z N } is: The maximum likelihood principle leads us to select the 6 unknown parameters T and PAR to maximize the log-likelihood LL given by The bivariate normal density of Z n is completely determined by the covariance matrix of Z n , and an easy computation then yields where A necessary condition to maximize LL is to let all its first derivatives be equal to zero : where z 1 (n), z 2 (n) are as in equation (4).
For the derivative ∂L ∂µ , we use the simpler log-likelihood based only on the SDE verified by S t , which does not make a significant difference in the estimates: Solving the 5 equations (5) − (9) is achieved numerically by gradient descent, in order to compute the maximum likelihood parameter estimates based on Euler discretization.We come back to the detailed study of these estimates below.In principle a more accurate discretization scheme with faster speed of convergence to the true SDE solutions could also be used to generate an approximation of the log-likelihood, and then could provide parameter estimates by maximizing the log-likelihood.We first show why, for parameter estimation, there is no real advantage in using more precise but more complex discretizations.We focus below on the Milstein discretization scheme, which approximates in L 2 the true SDE solution with accuracy of order T, instead of the √ T accuracy provided by the Euler scheme [21] [19].4. Milstein discretization scheme.Equations (1) and (2) in integral form are: Expand the integrals by Stochastic Taylor expansion [6] to obtain the following discretization: where are martingales with mean zero [23], variance As seen above for the Euler scheme, an approximate log-likelihood can be formally computed from the Milstein difference equation, and then one can compute the first derivatives of this approximate log-likelihood to set them equal to zero.The derivatives of the additional terms in the Milstein scheme are zero with respect to all parameters except γ.The functional form of the first order equations for κ, θ, γ, ρ is hence essentially the same for both discretizations.
Denote this system of first order equations by F(PAR, Z) where PAR = (κ, θ, γ, ρ), Z = (z 1 (1), z 1 (2), . . ., z 1 (N), z 2 (1), z 2 (2), . . ., z 2 (N)).We do not need to consider the equation corresponding to µ here since conditions corresponding to µ can be easily analyzed separately.Denote by (PAR E , Z E ) and (PAR M , Z M ) the parameter estimates and data vectors respectively corresponding to the Euler and the Milstein schemes.Define where ∆PAR is a 4 × 1 vector, and all partial derivatives are matrices of adequate dimensions, evaluated at (PAR E , Z E ).By definition of PAR M , the left hand side of the above equation is zero : The difference ∆PAR between our two vector estimates can be computed by inverting the preceding linear system, and of course depends on T. Using this expression for ∆PAR, and for any fixed small tolerance level η, we can compute a corresponding explicit upper bound for T, compelling each coordinate of ∆PAR to have variance less than η.
We also derive the asymptotic distribution of ∆PAR for large N.As N → ∞ the coefficient matrix of ∆PAR divided by N in the preceding linear system converges to the constant 4 × 4 symmetric matrix f , where all coefficients are zero except the following ones Denote the right hand side of the linear system (12) above by R = ∂F ∂Z .∆Z.We show that as N → ∞ the random vector R √ N is asymptotically Gaussian with mean zero and a limit covariance matrix Σ which we compute explicitly.
The two last results just mentioned show that as N → ∞, √ N times the difference √ N (PAR M − PAR E ) between the two groups of parameters estimates deduced from the Milstein and Euler schemes becomes asymptotically Gaussian with mean zero and an explicitly computed covariance matrix.Using these results, we have derived, for fixed large N, an explicit bound τ on T which will force the easily computed Euler scheme parameter estimates to remain reasonably close to the more accurate (but far harder to compute) estimators based on Milstein's discretization.

Constrained parameter estimation based on Euler discretization.
As just shown, the Euler discretization, combined with an adaptive computable upper bound τ on T, generates controllably accurate parameter estimators for the Heston SDEs.The constrained estimation problem is now reduced to the maximization of LL(T, PAR) as a function of T and PAR under the following constraints: where T < 2 κ ensures convergence of the Euler scheme [18].The constraint on T depends on the parameter vector PAR.But PAR itself is unknown and its estimator depends on T. So we implement a sequence of alternating estimations by gradient descent for PAR and T, iterating until the estimates stabilize, according to the following steps : • Start with initial value PAR 0 = (µ θ κ γ ρ) • kth iteration generates estimate PAR k from PAR k−1 by 1-step gradient descent 6. Results of empirical tests.We have performed many numeric simulations of the preceding iterative computation of estimates, and we have always observed that the estimation scheme is actually convergent.We observe that even for moderate values of N the estimates are good.We present below the estimates and empirical variances of these estimators for different values of N. The tests were made on simulated diffusion processes with simulation step δ = 10 −3 , (U 0 , U 1 . . .U N+1 ), (V 0 , V 1 , . . .V N+1 ) observed at t 0 , t 1 , . . ., t N+1 with the true value of T = t n+1 − t n = .005= 5δ.The empirical variances of our parameter estimates were computed over 50 simulated trajectories, and then compared to the theoretical variances obtained above.We observe a good fit between empirical and theoretical variances.

Asymptotic variance of the Euler scheme parameter estimators.
, PAR) • Continue till the sequence of estimates PAR k and T k stabilize 5.3.