Euler-Discretized Hull-White Stochastic Volatility Model

# Asymptotics for the Euler-Discretized Hull-White Stochastic Volatility Model

Dan Pirjol  and  Lingjiong Zhu Department of Mathematics
Florida State University
Tallahassee, FL-32306
United States of America
27 April 2016. Final: 4 July 2017
###### Abstract.

We consider the stochastic volatility model , with uncorrelated standard Brownian motions. This is a special case of the Hull-White and the (log-normal) SABR model, which are widely used in financial practice. We study the properties of this model, discretized in time under several applications of the Euler-Maruyama scheme, and point out that the resulting model has certain properties which are different from those of the continuous time model. We study the asymptotics of the time-discretized model in the limit of a very large number of time steps of size , at fixed and , and derive three results: i) almost sure limits, ii) fluctuation results, and iii) explicit expressions for growth rates (Lyapunov exponents) of the positive integer moments of . Under the Euler-Maruyama discretization for , the Lyapunov exponents have a phase transition, which appears in numerical simulations of the model as a numerical explosion of the asset price moments. We derive criteria for the appearance of these explosions.

###### Key words and phrases:
linear stochastic recursion, Lyapunov exponent, phase transitions, critical exponent, large deviations, central limit theorems.
###### 2010 Mathematics Subject Classification:
60G99,60K99,82B26,60F10,60F05

## 1. Introduction

Stochastic volatility models are widely used in financial practice for modeling the dynamics of the volatility surface. Some of the most popular models are affine models with stochastic volatility such as the Heston model, and models where the volatility is the exponential of a Gaussian process such as a Brownian motion or Ornstein-Uhlenbeck process.

In this paper we will study the stochastic volatility model defined by the process

 (1) dSt=σtStdWt (2) dσt=ωσtdZt

with initial condition , where and denote independent standard Brownian motions. The model parameter is a positive real constant.

This model is a particular case of the Hull-White model [17]

 (3) dSt=rStdt+√YtSt(ϱdWt+√1−ϱ2dBt) (4) dYt=μYtdt+ζYtdWt

where are uncorrelated standard Brownian motions. This reduces to the model (1),(2) by identifying and taking and .

The case of zero correlation has received special attention in the literature because of its analytical tractability [14, 32]. This is also a particular realization of the SABR model [16], corresponding to the so-called log-normal SABR model

 (5) dSt=σtSβtdWt (6) dσt=ωσtdZt,E[dWt,dZt]=ϱdt.

The model is also a limiting case (zero mean reversion) of the Scott model [28, 6], which corresponds to assuming that is the exponential of an Ornstein-Uhlenbeck process.

In practice the original SABR model as formulated in [16] is used mostly in parametric form for interpolating swaption or caplet volatilities, due to the unrealistic assumption of zero mean reversion for the volatility process. Nevertheless, due to its simplicity the model was studied extensively and many analytical results are available. The properties of this model were studied in continuous time in [19, 1, 23]. We briefly summarize a few results.

The asset price is a strict martingale only if the correlation is non-positive [19, 23], see also [29, 4]. This is a necessary condition if the diffusion (1) is to be used to model the price of a tradeable asset. The model (1), (2) has also moment explosions. Moment explosions in stochastic volatility models have been studied in a wide class of models, see [1, 12, 10]. Define the explosion time of the -th moment of the asset price with as

 (7) T∗(q)=sup{t≥0:E[Sqt]<∞}

For the model (1), (2) the explosion time is given by [19, 23]

 (8) T∗(q)={+∞,ϱ<ϱ∗(q)0,ϱ>ϱ∗(q)

where the critical correlation is

 (9) ϱ∗(q)=−√1−1q.

For all positive moments the critical correlation is negative , such that at zero correlation all these moments explode in zero time.

In practical implementation using Monte Carlo approaches, stochastic volatility models are simulated in discrete time. For this purpose, the stochastic differential equation (1), (2) is discretized in time, using one of several available time discretization schemes [21, 20]. The simplest scheme is the forward Euler time discretization, or Euler-Maruyama discretization [21]. This can be applied either directly to the stochastic differential equation for , or to that for . We will call these schemes the Euler and log-Euler schemes, respectively. The same treatment can be applied to the stochastic differential equation for . Although the latter can be solved exactly for this case, we consider also its discretization as an illustration for more complicated volatility processes, where an exact treatment is not available.

The Euler-Maruyama time discretizations of the stochastic volatility model (1), (2) have distinctive properties which can be different from those of the continuous-time model. For example, the asset price is a true martingale for any correlation (Proposition 28), in contrast to the continuous time model where this property holds only for . Also, under Euler-Maruyama discretization of , the moments of the asset price are finite for any parameter values . On the other hand, in continuous time, as noted above, all moments with explode in zero time [19, 23].

One surprising feature of the Euler-discretized model is that the positive integer moments with have a sudden rapid increase for sufficiently large or simulation time step . This phenomenon is well known to practitioners, and is known to appear in simulations of stochastic volatility models with log-normally distributed volatility. See [18] for an informal discussion. For a discussion in the context of Monte Carlo simulations of the Hull-White model with arbitrary correlation, see Sec. 5.2 in [11], where the explosion of the variance of the asset price is controlled by imposing an upper bound on the values of the stochastic volatility process .

This phenomenon introduces difficulties in the estimation of the the error of Monte Carlo pricing of payoffs , since a very large variance of the asset price may lead to a very large variance of the payoff [13]. Large values of the higher moments can have also direct relevance for pricing certain instruments. For example, in fixed income markets, the second moment of forward Libor rates is relevant for pricing certain instruments such as Libor payments in arrears [1].

The explosion of moments observed for the time discretization of the model (1), (2) appears in a wider class of models. It was observed in a discrete time stochastic compounding process with multipliers proportional to a geometric Brownian motion [24]. The positive integer moments of the compounding process were observed to explode for sufficiently large values of the volatility or time step . We emphasize that the explosion is to very large (but finite) values, and the moments remain strictly finite, as expected in a discrete time setting. This phenomenon was studied in [25] using large deviations theory, where it was shown that the moment explosion is due to a discontinuous behavior of the Lyapunov exponents . The existence of this limit requires that the model parameters are rescaled with such that a certain combination is kept finite in the limit. We will use in this paper a similar approach to study the phenomenon of moment explosion in the Euler discretized version of the stochastic volatility model (1), (2).

We study in this paper the large asymptotics of as in the model (1), (2) discretized in time under various applications of the Euler-Maruyama scheme. In usual applications of the Euler discretization one is interested in the limit of a very large number of time steps at fixed maturity . The limit considered here is different, as we take at fixed and . As explained in the next section, this covers the fixed-maturity limit fixed, with small and large (denoted as Regime 3 below). In addition, this scaling includes other regimes, corresponding to large maturity , small (Regime 1), and small maturity , large (Regime 2).

As mentioned, the Euler-Maruyama time discretized versions of the stochastic volatility model (1) can have different properties from those of the continuous time model. The different limit considered here introduces also differences in the asymptotics of the discrete time model compared to those of the usual Euler-Maruyama discretized model. For example the limits may be different from the corresponding large time limit of the continuous time model.

The motivation for adopting the specific large scaling of this paper is that the growth rates (Lyapunov exponents) of the positive integer moments exist and are finite under this scaling. We obtain explicit results for the Lyapunov exponents under this special scaling of the model parameters using large deviations theory, and study their functional dependence on the model parameters. We find that, under the application of the Euler-Maruyama scheme to , the Lyapunov exponents have non-analyticity in the model parameters which is similar to the phase transition studied in [25]. This phenomenon is responsible for the numerical explosions of the moments of the asset price observed in numerical simulations of the model, which have implications for the Monte Carlo simulation of the model as discussed above.

Section 2 introduces the different Euler-Maruyama schemes, and the scaling of the model parameters under the large limit considered in this paper. The asymptotic properties of the different schemes are discussed separately: the Euler-Log Euler scheme (Sec. 3), the Log Euler-Log Euler scheme (Sec. 4), Log Euler-Euler scheme (Sec. 5) and Euler-Euler scheme (Sec. 6). Section 7 presents a detailed comparison of the asymptotics of these schemes with the known results of the continuous time model, and we demonstrate good agreement between the properties of the phase transition for Euler-Log Euler scheme obtained from the asymptotic analysis of Section 4 with exact numerical simulations of the model. This agreement demonstrates the practical usefulness of our asymptotic results, as they give thresholds for the numerical explosion of the moments observed in numerical simulations of the model.

## 2. Euler-Maruyama Time Discretizations

Stochastic volatility models are usually simulated in practice in discrete time. Finite grid simulations discretize the time, asset, and volatility, while Monte Carlo simulation discretize only in time.

Consider the simulation of the one-dimensional Itô stochastic differential equation for the -dimensional stochastic vector for

 (10) dXt=σ(Xt,t)dWt+b(Xt,t)dt

with initial condition and coefficients which are globally Lipschitz functions. This ensures the existence of strong solutions [21, 22]. is a standard Brownian motion.

We divide the interval into time steps with uniform length . We would like to simulate the SDE (10) on the sequence of discrete time steps

 (11) 0=t0

The simplest time discretization is the explicit Euler scheme, or Euler-Maruyama discretization, which is defined by the recursion

 (12) Xi+1=Xi+σ(Xi,ti)(Wi+1−Wi)+b(Xi,ti)(ti+1−ti)

with initial condition . The convergence properties of the Euler-Maruyama scheme were studied in [30, 3, 15].

We can apply the Euler-Maruyama discretization either to or to the log-price , and same for . Upon taking all possible combinations, we have altogether four possible discretizations: Euler-log-Euler scheme, Log-Euler-log-Euler scheme, Euler-Euler scheme and Log-Euler-Euler scheme. For instance in Euler-log-Euler scheme, the first “Euler” refers to discretization of the asset price and the second “log-Euler” refers to the discretization of the volatility and so on and so forth.

###### Definition 1 (Euler-log-Euler scheme).

The Euler-Maruyama scheme (or simply Euler-log-Euler scheme) for the discretization of the SDE (1) is defined by

 (13) Si+1=Si(1+σi√τεi) σi+1=σiexp(ω(Zi+1−Zi)−12ω2τ).

We denote for simplicity , and and , with i.i.d. Gaussian variables with mean and variance independent of i.i.d. Gaussian variables with mean and variance .

The log-Euler discretization for coincides with the exact solution for the volatility at the start of the time interval since is a constant

 (14) σi=σ0eωZi−12ω2ti.

This scheme has the disadvantage that the asset price can become negative. An alternative scheme which preserves the positivity of the asset price is defined by applying the Euler scheme to . This is given by the following recursion:

###### Definition 2 (Log-Euler-log-Euler scheme).

The Euler-Maruyama scheme for (or the log-Euler-log-Euler scheme) for the discretization of the SDE (1) is defined by

 (15) Si+1=Siexp(σi(Wi+1−Wi)−12σ2iτ) σi+1=σiexp(ω(Zi+1−Zi)−12ω2τ).
###### Definition 3 (Euler-Euler scheme).
 (16) Si+1=Si(1+σi√τεi) σi+1=σi(1+ω√τVi).
###### Definition 4 (Log-Euler-Euler scheme).
 (17) Si+1=Siexp(σi(Wi+1−Wi)−12σ2iτ) σi+1=σi(1+ω√τVi).

Throughout the paper, we assume the following:

###### Assumption 5.

We assume that

 (18) β:=12ω2n2τ, (19) ρ:=σ0√τ,

are fixed positive constants111Note the change of notation for from the equation (5). In the remainder of the paper will be defined as in (18)..

The model parameters may depend on . Our results should hold also under the less restrictive conditions , , but for the sake of simplicity we will use the definitions in Assumption 5.

In usual applications of the Euler-Maruyama discretization the model parameters are kept fixed as . However, as shown above, in the continuous time case, the moments will be infinite for any . On the other hand, under the scalings (18), (19) the Lyapunov exponents of these moments will be seen to be finite. (The moments themselves will approach infinity as since .) Thus the scalings (18), (19) ensure that a well-defined limit exists for the growth rates of the moments (Lyapunov exponents).

Our assumptions (18), (19) include the following asymptotic regimes:

• Regime 1. When and are fixed constants, the volatility of volatility is of order and the maturity is of order . So this is the small volatility of volatility and large maturity regime.

• Regime 2. When is of the order , then such that is of the order and . This corresponds to the large initial volatility and small maturity regime.

• Regime 3. When is of the order , then is of the order , and is of the order . This corresponds to the fixed maturity regime.

The accuracy of the asymptotic results for the Lyapunov exponent increases with the number of steps . We illustrate the application of the asymptotic results with numerical examples in Section 7, where we show that the asymptotic limit gives a reasonably good approximation for the finite result of the Lyapunov exponents for values of as low as 40.

## 3. Euler-Log-Euler Scheme

### 3.1. Lyapunov Exponents of the Moments

We would like to compute the moments of the asset price in the time discretizations introduced in the previous section. We consider in this section the scheme (13). Let us recall that this scheme is defined by the recursion

 (20) Si+1=Si(1+σi√τεi) σi+1=σiexp(ω(Zi+1−Zi)−12ω2τ).

with initial condition .

The th moments of , are given by

 (21) E[(Sn)q] =Sq0E[n−1∏k=0(1+σ0√τeωZk−12ω2tkεk)q] =Sq0E⎡⎣n−1∏k=0[q/2]∑j=0q!(2j)!(q−2j)!(2j−1)!!(σ0√τ)2je2j(ωZk−12ω2tk)⎤⎦,

where we took the expectations over and used the identity for the th moments of Gaussian random variables

 (22) E[(1+σε1)q]=[q/2]∑j=0q!(2j)!(q−2j)!(2j−1)!!σ2j.

Since and , we can express the th moments by

 E[(Sn)q]=Sq0m(q)nE[n−1∏k=0elogρYk+ωZkYk−12ω2tkYk],

where are i.i.d. random variables such that

 (23) P(Yk=2j)=1m(q)q!(2j)!(q−2j)!(2j−1)!!,j=0,1,2,…,[q/2].

is the normalizer defined as

 (24) m(q)=[q/2]∑j=0q!(2j)!(q−2j)!(2j−1)!!.

Now, observe that and . Therefore,

 (25) E[n−1∏k=0elogρYk+ωZkYk]e−β ≤E[n−1∏k=0elogρYk+ωZkYk−12ω2tkYk] ≤E[n−1∏k=0elogρYk+ωZkYk].

Hence

 (26) limn→∞1nlogE[(Sn)q]=logm(q)+limn→∞1nlogE[n−1∏k=0elogρYk+ωZkYk],

if the limit exists. Indeed, using Large Deviations theory we will show that this limit exists, and is given by the following result.

###### Theorem 6.

For any , exists and it can be expressed in terms of a variational formula

 (27) λ(ρ,β;q)=supg∈Gq{g(1)logρ+β∫10(g(1)−g(x))2dx−∫10Iq(g′(x))dx},

where , where

 (28) fq(θ):=logE[(1+ε1eθ)q]=log[q/2]∑j=0q!(2j)!(q−2j)!(2j−1)!!e2jθ,

where is a normal random variable with mean and variance and

 (29) Gq :={g:[0,1]→[0,2[q/2]],g(0)=0, g is % absolutely continuous and 0≤g′≤2[q/2]}.
###### Proof.

The proof will be given in Section 8. ∎

###### Remark 7.

We have studied the limit for any non-negative integer . Observe that can take negative values, therefore is not well-defined for non-integer . If is a negative integer, is not well defined for any Gaussian random variable since is not Lebesgue integrable for any open interval including . Therefore, we restrict ourselves to the study of non-negative integer moments.

###### Remark 8.

Recall that

 (30) E[(1+ε1eθ)q]=[q/2]∑j=0q!(2j)!(q−2j)!(2j−1)!!e2jθ,

which equals to the th moment of a normal random variable with mean and variance . It is known that

 (31) E[N(μ,σ2)q]=(−2σ2)q2U(−q2,12,−μ22σ2),

where is a confluent hypergeometric function and the function’s second branch cut can be chosen by multiplying with . For general , it seems to be difficult to get a more explicit expression. But it is possible at least for .

(i) For ,

 (32) I2(x)=supθ∈R{θx−log(1+e2θ)}=12xlogx+12(2−x)log(2−x)−log2.

(ii) For ,

 (33)

(iii) For ,

 (34) I4(x) =supθ∈R{θx−log(1+(42)⋅1⋅e2θ+(44)⋅3⋅1⋅e4θ)} =supθ∈R{θx−log(1+6e2θ+3e4θ)} =x2logη4(x)−log(1+6η4(x)+3(η4(x))2),

where .

(iv) For ,

 (35) I5(x) =supθ∈R{θx−log(1+(52)⋅1⋅e2θ+(54)⋅3⋅1⋅e4θ)} =supθ∈R{θx−log(1+10e2θ+15e4θ)} =x2logη5(x)−log(1+10η5(x)+15(η5(x))2),

where .

The formula for is complicated but the limits for large , small , and large are more tractable. We have the following result.

###### Proposition 9.

(i) The limit for the Lyapunov exponent is

 (36) λ(ρ,0;q)=sup0≤x≤2[q/2]{xlogρ−Iq(x)}=fq(logρ).

(ii) The Lyapunov exponent is bounded from above and below as

 (37) β43[q/2]2+2[q/2]logρ−Iq(2[q/2]) ≤λ(ρ,β;q)≤43β[q/2]2+fq(logρ).

(iii) These bounds give the asymptotic behavior in the large limit

 (38) limβ→∞λ(ρ,β;q)β=43[q/2]2,

and in the large limit

 (39) limρ→∞∣∣∣λ(ρ,β;q)−β43[q/2]2−2[q/2]logρ+Iq(2[q/2])∣∣∣=0.
###### Proof.

The proof will be given in Section 8. ∎

### 3.2. Variational Problem

Consider the variational problem appearing in Theorem 6. This can be expressed equivalently in terms of the functional

 (40) Λ[f]≡logρ∫10f(x)dx+β∫10(∫1xf(y)dy)2dx−∫10Iq(f(x))dx

defined in terms of a function subject to the constraints

 (41) 0≤f(x)≤2[q/2].

The function is related to the function appearing in Theorem 6 as . The rate function is given by

 (42) Iq(x)≡supθ∈R{θx−fq(θ)},

where is given by Equation (28).

The variational problem for gives an integral equation

 (43) δΛ[f]δf=logρ+2β∫10K(y,z)f(z)dz−I′q(f(y))=0,

where the kernel is . This can be transformed into an ordinary differential equation by taking successive derivatives with respect to . The equation (43) is written as

 (44) logρ+2β∫y0zf(z)dz+2βy∫1yf(z)dz−I′q(f(y))=0.

Take one derivative with respect to

 (45) 2β∫1yf(z)dz−ddyI′q(f(y))=0.

Taking another derivative gives

 (46) 2βf(y)+d2dy2I′q(f(y))=0.

The function is given by the solution of this ordinary differential equation with boundary conditions

 (47) logρ=I′q(f(0)),f′(1)=0.

The first condition follows by taking in (44), and the second condition is obtained by taking in (45).

We will prove next that the equation (46) can be formulated in such a way that it only requires the cumulant function but not its Legendre-Fenchel transform . Introduce a new unknown function defined by

 (48) h(y)=I′q(f(y)).

This is inverted as

 (49) f(y)=fq(h(y)).

The proof of this relation follows from the observation that where is the value of which achieves the supremum in the definition of the rate function .

In conclusion, the equation satisfied by the function is given by the following result.

###### Proposition 10.

The function satisfies the second order differential equation

 (50) h′′(y)=−2βfq(h(y))=−V′(h(y))

where we defined the potential function as

 (51) V(h)=2βfq(h).

The boundary conditions for the equation (50) are

 (52) h(0)=logρ,h′(1)=0.
###### Proof.

This follows directly from combining Equations (46) and (49). ∎

This reduces the variational problem to that of finding the solution of an ordinary differential equation. The equation (50) with boundary conditions (52) is solved straightforwardly. We start by noting that the quantity

 (53) E=12[h′(y)]2+V(h(y))=V(h(1))

is a constant of motion. This can be used to express the derivative in terms of as

 (54) dydh=1√2(V(h(1))−V(h(y))=12√β1√fq(h(1))−fq(h)

so integrating from to we obtain

 (55) y=12√β∫h(y)h(0)=logρdx√fq(h(1))−fq(x).

Taking in this relation gives an equation for

 (56) Fq(h(1);ρ)=2√β

where

 (57) Fq(a;ρ)=∫alogρdx√fq(a)−fq(x).

Once is known, the function can be found by substituting its value into (55) and solving for .

Another important simplification is that can be expressed only in terms of . This reduces the variational problem to that of finding the extremum of a real function of one variable. This is given by the following result

###### Proposition 11.

The functional can be expressed as

 (58) Λ[f]=fq(h(1))−1√β∫h(1)logρ√fq(h(1))−fq(x)dx.
###### Proof.

The proof is given in Section 8. ∎

###### Remark 12.

The variational problem for the functional defined in (40) generalizes the variational problem considered in [25] to a wider class of rate functions . The problem considered in [25] corresponds to . The simpler result (58) agrees with the result for the functional in Proposition 7 of [25], upon substituting into (58).

### 3.3. Numerical Study and Analytical Approximations

The Lyapunov exponent given by the solution of the variational problem in Theorem 6 displays the phenomenon of phase transition. This is manifested as a discontinuity of the partial derivatives and at points along a curve in the plane.

This phenomenon was studied in [25] for the case of the moment generating function . It can appear when the equation (56) has multiple solutions for , and the optimal value of switches between two of these solutions. For sufficiently small values of below a critical value the equation (56) has multiple solutions, while for the solution is unique.

The properties of the phase transition have been studied in [25] using both numerical and analytical methods. The existence of a phase transition for sufficiently small and its absence for sufficiently large have been proved rigorously in Section 7 of [25]. These proofs can be extended immediately to the case of the moments, for which the corresponding moment generating functions and have a similar functional dependence to the function considered in [25]. For the higher moments the properties of the phase transition can be studied numerically.

We present in this section a numerical study of the Lyapunov exponent and its phase transition for the first few positive integer moments .

#### 3.3.1. The second moment (q=2)

The numerical solution of the variational problem can be simplified by taking as independent variable instead of . This is related to as shown in (139). Expressed in terms of this variable, the Lyapunov exponent of the moment is obtained by taking the supremum of the following functional over

 (59) Λ2(d) =βd2+log(1+ρ2) −βd3(1+ρ2)∫10y2dy1+ρ2−eβd2(y2−1).

We computed the Lyapunov exponent of the second moment using this relation. We show in Figure 1 plots of vs for several values of . These results show that for sufficiently small , below a critical value , the Lyapunov exponent has discontinuous partial derivatives and at a point . This corresponds to a phase transition of the Lyapunov exponent.

The phase transition appears when the supremum over of switches between two different values of . We show in Figure 3 the phase transition curve . This curve ends at a critical point with coordinates .

An analytical approximation for the function can be obtained in the limit. For the integral is approximated by Lemma 28 in [25] as

 (60) ∫10y2dy1+ρ2−eβd2(y2−1)=11+ρ2{13−12βd2log(ρ21+ρ2)+O(β−2)}.

This gives the approximation valid in the region .

 (61) Λ2(d)=−13βd3+βd2+12dlogρ21+ρ2+log(1+ρ2).
###### Lemma 13.

The cubic polynomial in in (61) reaches its supremum at or , according to the following conditions

 (62) supdΛ2(d)=⎧⎪ ⎪⎨⎪ ⎪⎩Λ2(d∗2) for β>−23logρ21+ρ2,Λ2(d∗1) for β<−23logρ21+ρ2.
###### Proof.

The derivative is a quadratic polynomial in

 (63) Λ′2(d)=−βd2+2βd+12logρ21+ρ2.

We are interested in the region , so we have . For we have in fact the stronger result that for all . This implies that is decreasing on and thus its supremum is reached at .

For , the derivative becomes positive in a region centered on , where are the roots of quadratic equation . has a local maximum at the largest root . This defines . (The smallest root corresponds to a minimum of .)

Under what conditions is a global supremum for ? This condition can be written as

 (64) Λ2(d∗2) =23β(1+12βlogρ21+ρ2)d∗2+16logρ21+ρ2+log(1+ρ2) ≥Λ2(d∗1)=log(1+ρ2).

We used the fact that is a solution of to eliminate cubic and quadratic terms in on the left side. This inequality can be further written as

 (65) 2R3+3R2−1=(R+1)2(2R−1)≥0

where . This is satisfied for . Thus we conclude that is a global supremum of provided that , and otherwise the global supremum is realized at . ∎