A non-parametric inference for implied volatility governed by a Lévy-driven Ornstein–Uhlenbeck process

1.1 Background

Understanding volatility is of fundamental importance for effective portfolio choice, derivative pricing, and risk management, among other issues. The unobservable nature of volatility leads us to two substantive aspects of the non-parametric identification procedure which we propose: 1) We assume that the stochastic implied volatility follows an OU process driven by a hidden Lévy process, which is observed at discrete time points. 2) We provide a non-parametric estimation of the Lévy measure based on a preliminary estimator of the characteristic function of the stationary distribution. These claims are supported by a conclusive body of empirical studies on the behaviour of implied volatilities of exchange traded options. In particular, using a Karhunen-Loéve decomposition of time series of volatility surfaces, Cont and da Fonseca (2002) performed a joint study of dynamics of implied volatilities for all strikes and maturities of index options. They concluded, among many other interesting findings, that

Relative movements of implied volatilities have little correlation with the underlying asset;

The use of OU processes in volatility modelling to capture the mean reverting behaviour of data is so firmly established in the literature that it hardly requires further motivations (Hull and White (1987), Heston (1993) and Barndorff-Nielsen and Shephard (2001)). Much of the recent research has stemmed from the latter work where the volatility follows a non-Gaussian OU process. Notable references include, among others:

Roberts et al. (2004) who fit SV models with the changes in the log price modelled via a drift-less diffusion and volatilities through a superposition of OU processes with Gamma marginals; Griffin and Steel (2006) who, in addition, investigate models with jumps in the returns, an incorporation of a leverage effect and different risk premium coefficients for different component volatilities; Gander and Stephens (2007) look at such OU models when the marginal of the OU process is generalised inverse Gaussian or inverse gamma and introduce an approximation to fractional Brownian motion in the returns equation. The finance literature is slightly more varied, Li et al. (2008) consider the impact of models which include, in the jumps on the returns, infinite activity Lévy processes (e.g. Applebaum (2009)). In more recent papers, Zhang and Xu (2014) introduce stochastic state variables into volatility dynamics and analyse the influence of state-variable volatile characters on investment stopping boundaries; and Ma et al. (2015) use a discretised version of mean-reverting square root model for valuation of volatility options.

The Lévy density of the asset return process summarises the information about the jumps and has been the main object of study of a large body of statistical work (cf. Todorov (2011)). Todorov and Qin (2017) reasons that the Lévy density in the finite activity jump case can be viewed as the conditional probability of jump arrival of given size. Therefore, the Lévy density of time series of implied volatility data, would contain information both for future expected jump risk as well as for the so-called risk-neutral probability measure, which is of interest both from a theoretical and applied point of view. All of the aforementioned research (except for Todorov and Qin (2017)) exploit distinct parametric specifications for the Lévy measure in the information set. In contrast, the non-parametric volatility models are void of any specific functional form assumptions about the stochastic processes governing the dynamics. Additionally, they differ from the parametric models in their focus on providing measures of the notional volatility, rather than the expected volatility.

In contrast to the existing literature on non-parametric models of stochastic processes in financial mathematics, instead of estimating the drift function and the coefficient function of the Lévy measure (e.g. Soulier (1998), Ait-Sahalia (1996) and Bandi and Renò (2011)), we directly estimate the hidden Lévy measure of the OU process governing the stochastic volatility. For example, Bandi and Renò (2011) assume that the stochastic volatility is governed by a SDE consisting of a Brownian motion correlated with the underlying asset and a jump competent independent of the Brownian motion. Then they identify the relevant coefficient functions (through estimates of the system’s infinitesimal moments) using non-parametric kernel methods for jump-diffusion processes. Despite its efficiency, their method is still limiting because of the assumption of normality in the Brownian motion, as well as the Poisson (or Gaussian) jump distribution. We take one step further and make no distributional assumptions for the finite variation terms.

It is important to realise that traditional estimation techniques such as maximum-likelihood and Bayesian estimation are cumbersome. This is caused by the intractability of the likelihood. By the Lévy-Khintchine theorem, a Lévy process has its natural parametrisation in the Fourier domain and, in general, there exists no closed form for the marginal distribution of a Lévy process. For discrete time observations from a stationary OU-process, the estimation is further complicated by the fact that the observations are dependent, so that the transition density of the process is needed. The latter is even harder to get hold of than the marginal density. For this reason an alternative estimation method is proposed, which we will now explain briefly. As a starting point, we take a sequence of preliminary estimators such that for each point in time, we have an unbiased estimator of the characteristic functions, either almost surely or in probability. The canonical example of such an estimator is the empirical characteristic function, though other estimators are possible. To prove that the empirical characteristic function is an appropriate preliminary estimator, we can use the β-mixing result of Davydov (1973). Our estimator for the canonical function of the Lévy process is a Cumulant M-estimator (CME) that minimises the random criterion function. Our approach is closely related to the seminal work of Jongbloed et al. (2005), where the robust statistical properties of the M-estimator is developed. Since a stationary OU-process is β–mixing, the convergence is guaranteed by an application of Birkhoff’s ergodic theorem (Krengel and Brunel (1985)).

The CME is calculated by a high dimensional maximisation problem, for which there is no known analytical solutions. Apart from the classic expectation maximisation of Dempster et al. (1977) maximum, there are a number of recent innovations in optimisation in the field of statistics which can be used. Here, we employ the support reduction algorithm introduced in Groeneboom and Wellner (1992) and further developed in Groeneboom and Wellner (1992), Groeneboom et al. (2001) and Groeneboom et al. (2008).

After obtaining the non-parametric estimation of the distribution, we will examine the model’s explanatory power using Kolmogorov-Smirnov and Anderson-Darling goodness-of-fit tests. A fundamental issue in statistical analysis is testing the fit of a particular probability model to a set of observed data. A simulation approximation to the null distribution of the test provides a convenient and powerful means of testing model fit. In this paper we consider the acceptance-rejection method for the simulation part. Finally, we compare the forecasting performance of the estimation procedure, using out-of-sample observations, with two very well known stochastic volatility models; namely the Heston (1993) model and the Barndorff-Nielsen and Shephard (2001) (BNS) model.

The remainder of this paper is structured as follows. Section 2 explains the theoretical framework of our model and method for estimating the Lévy measure. Section 3 provides the details of the estimation via the support reduction algorithm, as well as the out-of-sample test to analyse the forecasting power of our method. Section 4 presents our conclusion.

We purpose to derive a nonparametric estimator of the process X (t), which is observed at equal time intervals between [0, T]. Each observation time is defined by t n m = m Δ n ( m = 0 , . . . , n - 1 ) . The invariant probability distribution of X, denoted by π, has characteristic function ψ which can be expressed explicitly in terms of k in the following way

ψ k ( t ) = ∫ e thity π ( dy ) = exp ( ∫ 0 ∞ ( e ity - 1 ) k ( y ) y dy ) , t ∈ T , y > 0 .(2.5)

Here, k (.) is a right-continuous, decreasing canonical function that describes the process Z by the relationship k (x) = ν (x, ∞) (x > 0). The special scaling in the defining equation for X, (2.1), allows for a separate estimation of λ and k. The former estimation problem is deferred to Section 2.2.

In this section we outline a new estimation method for the canonical function k. As pointed out in the introductory section, maximum likelihood techniques are hampered by the lack of a closed form expression for the density of π in terms of the “natural” parameter k (the density of a nondegenerate self-decomposable distribution always exists, see Sato (1999), Theorem 27.13). Even if such an expression would exist, we would in fact need more, namely, the transition densities of X, in terms of k. This motivates our choice for another type ofM-estimator.

This estimator is constructed by firstly defining a preliminary estimator ψ ˜ n for ψ₀, which denotes the “true” underlying characteristic function. Any reference to the true distribution will be denoted by a subscript 0. For example, F₀ denotes the true underlying distribution function of X (1) and k₀ denotes the true canonical function. In what follows, we choose ψ ˜ n such that for each n, ψ ˜ n is a ch.f and ∀ t ∈ ℝ , ψ ˜ n ( t ) → ψ 0 ( t ) , as n→ ∞. Here, the convergence is either almost surely, or in probability. Hence, ψ ˜ n estimates ψ₀ consistently pointwisely.

Jongbloed et al. (2005) discuss that the convergence and consistency of such estimator relies on three important properties of the OU process, namely, the Markov property, self-decomposability and β - mixing. We have proved these properties of (2.1), at the beginning of this section.

Let ℒ¹ (η) be the space of integrable functions on ℛ⁺, for the measure η defined by η ( dy ) = 1 ∧ y y dy , y ∈ R + .

The definition of the measure η precisely suits the integrability condition on k, which can now be formulated by ||k|| _η < ∞. Then, let K : = {k ∈ ℒ¹ (η) : k (y) ≥0} for the measure η. Thus, we suppose K ⊆ ℒ¹ (η) is a convex cone, containing the functions of non-degenerate self-decomposable distributions.

Suppose ψ ˜ n is non-vanishing for sufficiently large n and thus admits a distinguished logarithm. Then a CME is defined as

k ˆ n = arg min k ∈ K ∫ | log ψ k - log ψ ˜ k | 2 w ( t ) dt ,(2.6) where w is an integrable weight function with compact support S _w . We call this estimator a cumulant-M-estimator (CME).

Note that in (2.6) we are using the cumulant function (i.e. the distinguished logarithm of a characteristic function), instead of the characteristic function, directly. This is because, apart from the issue whether this estimator is well defined, one disadvantage of using characteristic function is that the objective function is non-convex in k (convexity being desirable from a computational point of view).

Now let G denote the set of cumulants corresponding to K. Then the cumulant corresponding to a particular canonical function k₁, or characteristic function ψ₁ is denoted by g₁. Define L : K → G by [ L ( k ) ] ( t ) = ∫ 0 ∞ ( e ity - 1 ) k ( y ) y dy .

Jongbloed et al. (2005) proved that k is determined by the values of g = L (k) restricted to the compact support, S, under the assumption w is symmetric and strictly positive around the origin. Now let L² (w) denote the space of square integrable functions with respect to w (t) dt, such that for f, a measurable function of time, ∫|f (t) |²w (t) dt< ∞. Since elements of G are continuous and w is compactly supported, G ⊆ L² (w). Additionally, the mapping L : K → G is continuous, onto and one-to-one.

Next, define the mapping Γ _n  : G → [0, ∞) by Γ n ( g ) : = | | g - log ψ ˜ n | | w 2 = ∫ | g ( t ) - log ψ ˜ n | 2 w ( t ) dt .

Then the cumulant M-estimator, k ˆ n is the minimiser of Γ _n  ∘ L, such that [ Γ n ∘ L ] ( k ) : = | | L ( k ) - log ψ ˜ n | | w 2 = ∫ | [ Lk ] ( t ) - log ψ ˜ n | 2 w ( t ) dt , k ∈ K , over an appropriate subset of K, such that the estimator is well-defined. To find precise conditions, we will go along the following two steps:

Use the fact from Hilbert-space theory that every non-empty, closed, convex set in L² (w) contains a unique element of smallest norm. Since Γ _n is a squared norm, it suffices to specify a closed, convex subset G′ of G to obtain a unique minimiser of γ _n over G′.

Let R > 0, and (i) K _R is a compact, convex subset of ℒ¹ (η), as well as (ii) G _R is a compact, convex subset of L² (w). Since Γ _n has a unique minimiser over G _R and to each G _R belongs a unique member of K _R , there exists a unique minimiser of Γ _n L over K _R . More precisely, Let g ˆ n = arg min g ∈ G R Γ n ( g ) and k ˆ n = arg min k ∈ K R [ Γ n L ] ( k ) . The arguments for the existence and uniqueness of k ˆ n can be established based on the fact that L : K _R  → G _R is onto and one-to-one, to each g ∈ G _r . Thus, there corresponds a unique k ∈ K _R such that L (k) = g.

For numerical purposes we will approximate the convex cone K by a finite-dimensional subset. For 1 ≤ j ≤ N, set Θ = {θ₁, . . . , θ _N } such that θ _j  = jh corresponds to grid points of the mesh h. Then, set 𝒰 _Θ  : = {u _θ , θ ∈ Θ} and define a basis function by u θ ( y ) : = 1 [ 0 , θ ) ( y ) , y ≥ 0 , z θ ( t ) = [ Lu θ ] ( t ) = ∫ 0 θ e itu - 1 u du , t ∈ R .

Now define a convex cone K _Θ by K Θ = { k ∈ K | k = ∑ i = 1 N α i u θ i , α i ∈ [ 0 , ∞ ) , 1 ≤ i ≤ N } .

Then define a unique sieved estimator by

k ˇ n = arg min k ∈ K Θ Γ n ∘ L ( k ) = arg min α i ∈ [ 0 , ∞ ) | | ∑ i = 1 N α i z θ i - log ψ ˜ n | | w 2 .(2.7)

The consistency of k ˇ follows from the consistency of ψ ˜ n , as proven by Jongbloed et al. (2005). Further, Δ may be fixed or Δ ↓ 0. Since many procedures in the statistical theory for stochastic processes are tailor-made either for fixed Δ or for Δ ↓ 0, we must note that the presented algorithm is robust with respect to this sampling frequency. An elaborate discussion on the consistency of this CME is provided in Jongbloed et al. (2005), to which interested readers may refer.

Suppose X_1Δ, X_2Δ, …, X _nΔ are discrete-time observations from the stationary OU process, where Δ is a constant. The latter condition is very important, otherwise, nΔ→ ∞ should be a necessary condition to identify the drift, since for a fixed time horizon even with continuous observations Girsanov’s theorem shows that the corresponding laws are equivalent and thus their distance is fixed as long as nΔ is bounded.

With a slight abuse of notations set X _i  : = X _iΔ . The chain X _n satisfies the first order auto-regressive relation X n = e - λ Δ X n - 1 + W n ( λ ) , where (W _n  (λ))  _n is an i.i.d. sequence of random variables distributed as W λ : = ∫ 0 Δ e λ ( u - Δ ) dZ ( λ u ) , with common infinitely divisible distribution.

Let θ = e^-λΔ and denote the true parameter by θ₀. As in Jongbloed et al. (2005), define the estimator θ ˆ n = min 1 ≤ k ≤ n X k X k - 1 .

Then, θ ˆ n ≥ θ 0 , for each w. Hence, θ ˆ n is biased upwards. However,

Theorem 4. θ ˆ n is consistent: θ ˆ n → P θ 0 as n→ ∞.

Proof. Let ɛ > 0. Since { | θ ˆ n - θ 0 | > ɛ } = { θ ˆ n > θ 0 + ɛ } , p ( n , ɛ ) : = P ( | θ ˆ - θ 0 | > ɛ ) = P ( X k / X k - 1 > θ 0 + ɛ , ∀ k ∈ { 1 , ⋯ , n } ) = P ( W k ( λ ) > ɛ X k - 1 , ∀ k ∈ { 1 , ⋯ , n } )

Define N n : = ∑ k = 1 n 1 { X k - 1 > 1 } , then p ( n , ɛ ) = ∑ j = 0 n P ( W k ( λ ) > ɛ X k - 1 , ∀ k ∈ { 1 , ⋯ , n } | N n = j ) P ( N n = j ) ≤ ∑ j = 0 n ( P ( W 1 ( λ ) > ɛ ) ) j P ( N n = j ) , where the inequality holds since {W _k  (λ)} _k≥1 is an i.i.d. sequence. Since W₁ (λ) has support [0, ∞) (Sato (1999), Corollary 24.8), α _ɛ  : = P (W₁ (λ) > ɛ) ∈ [0, 1). This gives p ( n , ɛ ) ≤ ∑ j = 0 ∞ α ɛ j P ( N n = j ) .

By dominated convergence, lim n → ∞ p ( n , ɛ ) ≤ ∑ j = 0 ∞ α ɛ j [ lim n → ∞ P ( N n = j ) ] . Finally, note that lim n → ∞ P ( N n = j ) = 0 , which concludes the proof. □

Then λ ˆ n = - log θ ˆ n , then λ ˆ n → P λ 0 , as n→ ∞, where λ₀ denotes the true value of λ. Nielsen and Shephard (2003) provided the detailed asymptotic analysis, which showed that θ ˆ n equals the maximum likelihood estimator for the model.

This approach is critically interesting, from the application point of view, compared to the parametric cases. The advantage is the flexibility and the ability to ‘let the data speak for themselves’. In other words, our data adaptive approach allows us to estimate the density function of the process for stochastic volatility, which is closer to the ‘true’ distribution. However, problems that can usually be solved using standard techniques in the parametric case, can be much more difficult in the non-parametric situation. For example, a M-estimator is defined as the minimiser of a random criterion function over an appropriate parameter set. In parametric models, estimates can be computed explicitly or approximated using some numerical technique; whereas, in non-parametric models, the computational issues often boil down to high dimensional constrained optimisation problems. From the practical point of view, this means higher computational expense, which may be a major deterrent to adopt to implementation. To overcome this problem, we use the support-reduction algorithm based on Groeneboom et al. (2008) for the computation of the estimator. See Section 3.1 further details.

The M-estimator is defined as the minimiser of a random criterion function over the parameter set. In parametric models, estimates can often be computed explicitly or approximated using some numerical technique for solving (low dimensional) convex unconstrained optimisation problems. Nevertheless, in our non-parametric case, we face a high dimensional constrained optimisation problem. Apart from the classic expectation maximisation of Dempster et al. (1977), there are a number of recent innovations in optimisation in the field of statistics. Here, we use the support-reduction (SR) algorithm, introduced in Groeneboom and Wellner (1992) and further developed in Groeneboom and Wellner (1992), Groeneboom et al. (2001) and Groeneboom et al. (2008). A version of this algorithm, known as the iterative convex minorant algorithm, was also used in Jongbloed et al. (2005) for the numerical approximation of the CME introduced in Section 2.1. However, the method of Groeneboom et al. (2008) (which is used in this paper) provides faster convergence, when computing non-parametric M-estimators in mixture models. Within the optimisation theory, the SR algorithm can be classified as a specific instance of an active set method. Within the field of statistical computing, the algorithm fits in the class of vertex direction algorithms. An algorithm, related to our support reduction algorithm, can be found in Meyer (1997), which led to the idea of the iterative spline algorithm for convex regression, as described in Groeneboom et al. (2001).

The algorithm discussed here is based on isotonic regression techniques as can be found in Härdle (1989). It is mostly used to compute shape-restricted estimators of distribution functions in semiparametric models. We apply the algorithm to compute M-estimators in mixture models, using unconstrained optimisations iteratively. It is to be noted that in each step the algorithm adds one support knot among several possible support knots to the existing iterate, resulting in a sparse next iteration. These reduced knots causes a low-dimensional unlimited optimisation problem during the next iteration. A fundamental advantage of this method is that in problems where the solution is a sparse mixture, the model can be handled well, due to the speed at which the low-dimensional optimisations are performed.

We seek to find an approximate solution for k ˇ n , numerically and characterised by variational inequalities. The existence of a unique maximiser over the cone K _Θ is proven in Groeneboom et al. (2008), Lemma 1. Most importantly, they showed that k ˇ n minimises Γ _n  ∘ L (k) in (2.7) iff

[ D k ˇ n ( Γ n ∘ L ) ] ( u θ j ) { ≥ 0 ∀ j ∈ { 1 , . . . , N } , = 0 ∀ j ∈ J ,(3.8) where J : = {j ∈ {1, . . . , N} |α _j  > 0} refers to the set of positive α-weights, and [D _a Φ] (b) is a directional derivative function from the right D _a Φ in the direction of a defined by [ D a ( Φ ) ] ( b ) : = lim ɛ → 0 [ Φ ] ( b + ɛ a ) - [ Φ ] ( b ) ɛ .

Note that whenever (Φ) (b)< ∞, convexity of Φ guarantees the existence of [D _a  (Φ)] (b). Using straight forward, but tedious, calculations we can show [ D k ˇ n ( Γ n ∘ L ) ] ( u θ j ) = 2 〈 L u θ - g ˜ n , L k ˇ n 〉 w .

With a slight abuse of notations, let k ^J from the cone K _Θ be the current iterate from each iteration. Then suppose k ^J  = ∑_j∈Jα _j u _{θ j} .

The optimisation algorithm in Groeneboom et al. (2008) prescribes to evaluate the optimality of k ^J (3.8). If k ^J is not optimal, then there exist an i ∈ {1, . . . , N} ∖ J with a negative directional derivative, which means u _{θ i} provides a direction of descent for the objective function. Therefore, our initial aim is to determine the set of grids with negative directional derivatives, all of which are non-optimal and provide a direction of descent. Then using a search algorithm, we pick a direction of descent and find a new iterate. In this stage, we must recalculate the optimal weights, using the standard least square technique. For doing so, we differentiate the objective functions with respect to α _j  (j ∈ {1, . . . , m}) and set them to zero. This gives the following system

〈 z θ i , z θ j 〉 w α _ = 〈 z θ i , g ˜ n 〉 w , i , j = 1 , . . . , m ,(3.9)

Equation (3.9) is easily proved by showing that the left-hand-side is symmetric and non-singular. This is briefly as follows.

Let A : = 〈z _{θ i} , z _{θ j} 〉 _w and α _ j denote the j-th the column of A. We aim to show that if ∑ i = 1 m = h i α _ i = 0 _ then all h j ∈ ℝ , j ∈ { 1 , ⋯ m } are zero. Now, ∑ k = 1 m h k w θ k ( t ) = 0 for all t ∈ S, where w θ ( t ) = t d dt z θ ( t ) = cos ( θ t ) - 1 + i sin ( θ t ) .

Then also, d p dt p ∑ k = 1 m h k w θ k ( t ) | t = 0 = ∑ k = 1 m h k ( i θ k ) p = 0 , p = 1 , ⋯ , m .

Rewriting this in a linear system we find that ( i θ 1 i θ 2 ⋯ i θ m - θ 1 2 - θ 2 2 ⋯ - θ m 2 ⋮ ⋮ ⋱ ⋮ ( i θ 1 ) m ( i θ 2 ) m ⋯ ( i θ m ) m ) ( h 1 h 2 ⋮ h m ) = ( 0 0 ⋮ 0 ) .

Since the matrix in this display is a generalised Vandermonde matrix, its determinant is unequal to zero. Thus, h₁ = h₂ = ⋯ = h _m  = 0.

Finally, we evaluate the new iterate by (3.8) for optimality, and repeat the process until the optimal solution is found.

Theorem 1 in Groeneboom et al. (2008) gives conditions to guarantee that the sequence of iterates {k⁽ⁱ⁾}  _n (generated by the support-reduction algorithm) indeed converges to the solution of our minimisation problem. Since these conditions are met in our case, we have ( Γ n L ) ( k ( i ) ) ↓ ( Γ n L ) ( k ˇ n ) , i → ∞

The support-reduction algorithm is summarised as follows

Step 1. In the current iterate, add a basis function u _θ for Γ _n L. Due to the linearity condition of L,

[ Γ n L ] ( k + ɛ u θ ) - [ Γ n L ] ( k ) = ɛ c 1 ( θ , k ) + 1 2 ɛ 2 c 2 ( θ ) ,(3.10)

where c 1 ( θ , k ) = 2 〈 Lk - g ˜ n , Lu θ 〉 w , and c 2 ( θ ) = 2 | | Lu θ | | w 2 .

We implement the Divide and Conquer (D&C) algorithm to search for the solution. In this algorithm, a successful join of points is determined when the new solution decreases the objective function and is not outside of K _Θ . Throughout the optimisation process, we may find solutions that contain support points in too close vicinity of each other. This is because the grids are not continuous. To prevent such situation we locally optimise the support by scanning for closely located points. If we found such support points, we determine a convex combination of them that minimises Γ _n  ∘ L (k). This may have a smoothing effect that does not interfere with the convergence condition specified in Groeneboom et al. (2008). Additionally, the extra step can reduce the convergence speed dramatically; therefore, we will only implement the replacement when we find two points with two new directions of descend, that cause the algorithm zigzag between them.

Another issue that might arise is that the solution α _i may generate a measure ∑α _i u _{θ i} , not necessarily belonging to ℛ⁺ as α _i might be negative for some i. If this function happens to belong to ℛ⁺, it is a new iterate. Otherwise, one travels as far as possible along the line segment connecting the current iterate with this infeasible measure. A remedy for this problem is offered in Groeneboom et al. (2008), which we implement in our algorithm as part of Step 2. For this purpose, we use a sequence of unrestricted minimisations and support-reduction, to obtain a new subset of support-points that minimises Γ _n  ∘ L (k).

Finally, in relation to the computational efficiency the support-reduction estimation has a very low run time. Even though the calculation of matrix c₂ (θ) needs a long time (in fact it depends on the total number of mesh grids), the algorithm works highly efficiently, once it is initialised, the determination of the vector is fast.

In this section, we approximate the canonical functions for VDAX and VIX, using the cumulant M-estimator. Further, we approximate the intensity parameter λ for the implied volatility of both indices. Table 1 presents the results for the support-reduction algorithm applied to the volatility data set. Additionally, Fig. 2 plot ψ ˜ ( t ) for each of the volatility indices.

OPEN IN VIEWER

Fig.2

These figures exhibit the plots for the estimated canonical functions and characteristic functions for both VIX and VDAX. The shape of the characteristic functions remind us of log-normal distribution, therefore, we test the hypothesis that whether the data sets were drawn from log-normal distributions. The results are provided in Table 1.

Next, we want to obtain the stationary distribution function f _n using fast Fourier transform. Note that ψ ˜ ( t ) = exp ( ∫ 0 ∞ ( e ity - 1 ) k ˜ ( y ) y dy ) , = exp ( ∑ j = 1 M α j z θ j ( t ) ) , thus, f ˜ ( y ) = 1 2 π ∫ - ∞ ∞ e - ity + ∑ α j z θ j ( t ) dt .

The characteristic function of a self-decomposable distribution does not need to be integrable. Such condition means that there is no guarantee that we can calculate the Fourier transform analytically. Here, we use an alternative method based on Schorr (1975), where we numerically approximate the Fourier transform. We only present the results for the numerical approximation, without the proof. Interested readers may refer to Schorr (1975) for the detailed calculations of similar examples. f ( y ) = ∑ κ = - ∞ + ∞ a κ exp ( i κ wy ) , where the parameters are defined as follow: a k = w 2 π ∫ - ∞ + ∞ f ( u ) exp ( - i κ wu ) du = w 2 π ψ ˜ ( - κ w ) . Here, w = 2 π T , T = T₊ - T_-, where T + = ( - 2 ln ( ɛ + ) ) 1 2 σ , T - = ( - 2 ln ( ɛ - ) ) 1 2 σ , where σ is the standard deviation of the data, so Fourier approximation follows

f ˜ ( y ) = ∑ κ = - ∞ ∞ T - 1 exp ( ∑ j = 1 M α j ∫ 0 jhw κ e - iz - 1 z dz + i κ wy ) .(3.12) Schorr (1975) shows that f ˜ ( y ) is well defined.

We consider 60 grid points for the numerical approximation with lags of 0.05; therefore, we will have 60 basis functions. Moreover, for estimation of density functions based on the Fourier transform, T is considered to be 12 and κ takes values from -150 to 150. This interval is sufficiently large enough to obtain a reasonably accurate approximation.

Next, using the MLE method discussed in the previous section, we estimated the intensity parameters, λ, which are also reported in Table 1, then, in Fig. 3 we plot the non-parametrically estimated distributions for both indices, overlaying the histogram of the datasets.

OPEN IN VIEWER

Fig.3

These figures exhibit the plots for the non-parametrically estimated density functions for VIX and VDAX, overlaying the histogram of the datasets.

OPEN IN VIEWER

Fig.4

These figures present a comparison between different stochastic volatility models (i.e. non-parametric, Heston and BNS models) and their forecasting abilities. Out-of-sample simulated paths are indicated by blue lines and the times series (real) of daily implied volatility for the same time periods are indicated by red lines.

In order to compare the performance of our non-parametric model with some conventional parametric cases; namely, log-normal, Inverse Gaussian, Burr, and Gamma. At the beginning of the paper we argued that the flexibility of the non-parametric framework is a considerable advantage of our approach to model implied volatility, compared to the parametric cases. Here, we are interested to compare the statistical power of our model with the above mentioned parametric models. A fundamental issue in statistical analysis is testing the fit of a particular probability model to a set of observed data. A simulation approximation to the null distribution of the test provides a convenient and powerful means of testing model fit. Tu (2006) proposed a new Monte Carlo-based methodology to construct this type of approximation when the model is semi-structured. They showed that when there are no nuisance parameters to be estimated, the non-parametric Monte Carlo test can exactly maintain the significance level and when nuisance parameters exist, this method can allow the test to asymptotically maintain the level.

Generating random values from the estimated density function (3.12) requires some effort, since the quantile function is not available in closed form and the inversion method is, therefore, not possible. As an alternative, since f ˜ ( y ) is unimodal, we consider the acceptance-rejection method, which is a type of Monte Carlo method based on the observation that to sample a random variable one can sample uniformly from the region under the graph of its density function. Acceptance-rejection method assumes the existence of a density g and the knowledge of a constant c such that f ˜ ( y ) ≤ cg ( y ) , ∀ y . Therefore, in order to apply the algorithm, one needs to be able to determine a suitable constant c. Besides, it is necessary to specify an appropriate random variable with density g, such that cdf G easily invertible and g has to dominate f ˜ .

The basic idea is to search for a density function g, from which we already have an efficient algorithm to generate from, but also such that the function g is “close” to f ˜ . Therefore, the ratio f ˜ ( y ) / g ( y ) is bounded by c, that is

sup { f ˜ ( y ) / g ( y ) } ≤ c .(3.13)

As in Devroye (1986), the algorithm for generating random variables distributed as f ˜ is

Step 1. Determine g and its parameters such that it dominates f ˜ . The parameters of g are set such that the first four moments of g are as close as possible to the first four moments of f ˜ .

For the density functions of both of the datasets, we consider a wide range of candidates for g and their parameters specifications, such that the dominance condition is met. Choosing the most appropriate g is a heuristic process, for which our best guide is the constant c. It can be mathematically presented that the acceptance rate is 1/c, then the expected number of iterations of the algorithm required must be c (cf. Sigman (2007)). If this is not met, we must try a different g.

We found that for simulating (3.12), log-normal distribution is a suitable choice for g. Simulating from this distribution is relatively simple and fast, which is especially important since, in order for the f ˜ ( y ) to be dominated, the factor c needs to grow with the kurtosis. As a result, the number of rejections also increases and, if the algorithm to generate from the dominating density is not sufficiently fast, computation becomes slow. Table 2 reports the parameter specifications of the log-normal distributions used for each data set, as well as the acceptance rate and the boundary constant, c. We can see that in all cases the inverse of the acceptance rate is very close to c, which verifies the appropriateness of choosing the log-normal distribution. We thus have two Kolmogorov-Smirnov (KS) tests and two Anderson-Darling (AD) tests, applied to the simulated distributions. The KS test is based on the maximum distance between the cumulated probability functions of the simulated density and the sample, whereas the AD test gives more weight to the tails than does the K-S test. In all four scenarios we accept the null hypothesis at the significance levelof 0.05.

Next, we compare if our non-parametric estimation of the density better fits the data better than some conventional parametric cases. Table 3 reports this comparison using log-normal, Inverse Gaussian, Burr, and Gamma distributions. At a significant level of 0.05, we can clearly reject the null hypothesis that our samples (both VIX and VDAX) are from Inverse Gaussian, Burr or Gamma distributions. Nevertheless, the log-normal appears to be a contester for our model. For both VIX and VDAX we can accept the null hypothesis that our sample is from a log-normal distribution. In terms of the KS-test, the log-normal distribution has even slightly out performed our model, evident from their p-values. However, the p-values of the AD tests suggests that our non-parametric model might out perform the log-normal distribution in modelling the tail behaviours. The results are too close to conclusively identify the better model, however, due to the general advantageous of non-parametric approaches (discussed in Section 1), our method is still preferred.

Model selection depends not only on the goodness-of-fit of a model to the data, but also on the objective of the analysis. For forecasting, a model that is best in the in-sample fitting does not necessarily provide more accurate forecasts. Thus, many people use the performance of out-of-sample forecasts to aid the selection of a statistical model.

Statistical tests of a model’s forecast performance are commonly conducted by splitting a given data set into an in-sample period, used for initial parameter estimation and model selection, and an out-of-sample period, used to evaluate forecast performance. Empirical evidence based on out-of-sample forecast performance is generally considered more trustworthy than evidence based on in-sample performance which can be more sensitive to outliers and data mining (White (2000)). Out-of-sample forecasts also better reflect the information available to the forecaster in “real time”. This has led many researchers to regard out-of-sample performance as the “ultimate test of a forecasting model” (Stock and Watson (2007)).

We divide the data for both VIX and VDAX into two sub-periods, namely, the estimation sub-sample which is the first period, and the forecasting sub-sample which is the remainder of the data set. Given T data points, say x₁, …, x _T , we divide the data as {x₁, …, x _n } and {x_n+1, …, x _T }, where n is the initial forecast origin.

Suppose, three competing models: A) the non-parametric model in the current work B) the Heston’s stochastic volatility model (Heston (1993)) and C) the BNS pure jump stochastic volatility model (Barndorff-Nielsen et al., (2002))) with Gamma marginal distribution as in Griffin and Steel (2006).

Let h be the forecast horizon, i.e. we are interested in 1-step to h-step ahead forecasts. The out-of-sample forecasting evaluation works as follows:

Let m = n be the initial forecast origin. Fit models A, B and C. Note that for the non-parametric case, we need to repeat the support-reduction algorithm for the new window. For brevity, we do not discuss the calibration of models B and C. Interested readers may refer to Schorr (1975) and Griffin and Steel (2006).

In this way, we will have (T - n - 1) 1-step ahead forecast errors for each model, and (T - n - 2) 2-step ahead forecast errors for each model, and so on. For l-step ahead forecasts of the model, we then compute the following three functions in the models out-of-sample forecast performance evaluation. MSE ( l ) = 1 T - l ∑ t = 1 T - l ( x t + l 2 - x ˆ t 2 ( l ) ) 2 , R 2 LOG ( l ) = 1 T - l ∑ t = 1 T - l ( log ( x t + l 2 x ˆ t 2 ( l ) ) ) 2 , PSE ( l ) = 1 T - l ∑ t = 1 T - l ( x t + l 2 - x ˆ t 2 ( l ) x ˆ t 2 ( l ) ) 2 , where l = 1, 2, …, h. These three criteria were all suggested by Griliches et al. (1994). The Mean Square Error (MSE) is the average of the squared deviations of the forecast from the implied volatility. Hence, one large deviation is given a much higher weighting than a sum of small deviations even if the sum of the deviations is equal to the one time large deviation. When evaluating volatility forecasting performance this can seem quite illogical since in general one large deviation is not more troublesome than a sum of small deviations that sum up to the size of the large deviation. Moreover, single outliers will have a significant impact on the MSE criteria. The R²LOG still assigns higher weighting to large deviations but they are not as penalised as in the case of MSE. The Percentage Squared Errors (PSE) measures the average of the squared percentage deviation. The deviations are expressed as a percentage of the forecasted volatility. Hence, the PSE takes into account the fact that it is harder to be accurate, in an absolute sense, when estimating high variances and thus measures the relative error as a percentage to account for this.

For all the three loss functions a smaller value is preferred. Table 4 reports the average MSE, R²LOG, and PSE for the two compared parametric implied volatility models, as well as the proposed non-parametric model. The results, unequivocally, suggest that our non-parametric model provides a better out-of-sample fit (in terms of MSE, R²LOG, and PSE) in all cases for both VIX and VDAX.