2. Limiting implied volatilities

2.1. Wishart multifactor stochastic volatility process

Da Fonseca et al. [Reference Da Fonseca, Grasselli and Tebaldi8] suggested that in an arbitrage-free frictionless financial market, the WMSV can represent the dynamics of a return of a risky asset price $S_t$ ,

(1) \begin{align}\nonumber Y_t &= - \frac12 \textrm{Tr}\left[\Sigma_t\right] \rd t+\textrm{Tr}\left[\sqrt{\Sigma_{t}}\left(\rd W_{t}R^{\top}+\rd B_{t}\sqrt{\mathbb{I} - RR^{\top}}\right)\right],\\\rd\Sigma_{t}&=\big(\Omega\Omega^{\top}+M\Sigma_{t}+\Sigma_{t}M^{\top}\big) \rd t+ \sqrt{\Sigma_{t}} \rd W_{t} Q + Q^{\top}\left(\rd W_{t}\right)^{\top}\sqrt{\Sigma_{t}},\end{align}

where $Y_t = \ln S_t$ , $\mathbb{I}$ is the n-dimensional identity matrix, $\Omega, R, M, Q\in \mathcal{M}_{n}$ (the set of square matrices), and $W_{t}, B_{t}\in \mathcal{M}_{n}$ are composed of $n^{2}$ independent Brownian motions under the risk-neutral measure. Here, the risk-free rate r is assumed to be zero without loss of generality. As noted in [Reference Da Fonseca, Grasselli and Tebaldi8], (1) represents the matrix analogue of the square-root diffusion of the Heston model [Reference Heston22]. The random shocks to the stock return and the volatility process are correlated by the matrix R.

There are additional assumptions on the model parameters. The matrix M is negative semi-definite, which is related to the mean-reverting feature in typical volatility modeling practice. It is also assumed that Q is nonsingular and $\Omega\Omega^{\top}=\beta{Q}^{\top}Q$ with real parameter $\beta>{n-1}$ . These standard conditions ensure that the resulting symmetric matrix $\Sigma_t$ is positive semi-definite. For the rest of this paper, such conditions are treated as part of the WMSV model.

The usefulness of the WMSV model lies in the transform techniques widely studied in the literature, including [Reference Duffie, Pan and Singleton10], still being applicable. The transform formula is well documented in [Reference Da Fonseca, Gnoatto and Grasselli6]. We record here a version that fits our purpose: for a scalar p,

(2) \begin{equation}\ln \mathbb{E}\left[ \textrm{e}^{p Y_t}\right] = p Y_0 + \textrm{Tr}\left[ A(p,t) \Sigma_0\right] + b(p,t).\end{equation}

Here, A(p, t) is a solution to the matrix Riccati differential equation given as

(3) \begin{equation}\frac{\rd A}{\rd t} = A\left( M + p Q^{\top} R^{\top}\right) + \left( M^{\top} + p RQ\right) A + 2 A Q^{\top} Q A + \frac{p(p-1)}{2} \mathbb{I}\end{equation}

with the initial condition $A(p,0) = \textbf{0}$ . Lastly, $b(p,t) = \beta \int_0^t \textrm{Tr}\left[ Q^{\top} Q A(p,s)\right] \rd s$ . There are many results in the literature regarding solutions to the above differential equation. For completeness, we provide a statement that is used in the next section.

Proposition 1. Consider a matrix differential equation for a fixed p:

\begin{equation*} \frac{\rd}{\rd s} \left(\begin{array}{c@{\quad}c} G & F \end{array}\right) = \left(\begin{array}{c@{\quad}c} G & F \end{array}\right) \left( \begin{array}{c@{\quad}c} M + p Q^{\top} R^{\top} & - 2 Q^{\top} Q \\ \\[-7pt] \frac{p(p-1)}{2}\mathbb{I} & - \left(M^{\top} + p RQ\right) \end{array}\right)\!,\end{equation*}

with $G(p,0) = 0$ and $F(p,0) = \mathbb{I}$ . As long as F(p, t) is invertible, (3) has a solution $A(p, s) = F(p,s)^{-1}G(p,s)$ on [0, t]. Conversely, if (3) has a solution A(p,s) on [0, t], then F(p,s) is invertible on [0, t] so that $A(p,s) = F(p,s)^{-1}G(p,s)$ . In such a case, b(p, t) is simplified to

\begin{equation*} b(p,t) = - \frac{\beta}{2} \textrm{Tr}\left[\ln F(p,t) + t\left( M + p Q^{\top} R^{\top}\right) \right]. \end{equation*}

The first statement is already noted in [Reference Da Fonseca, Grasselli and Tebaldi8, Reference Da Fonseca5], where the full details are referred to [Reference Grasselli and Tebaldi20]. The second statement is a slight modification of a result reported in [Reference Choi3, Reference Reid26]. Hence, we have the equivalence between the existence of A and the non-singularity of F(p, t). In addition, the validity of the affine transform formula (2) for general affine processes is proved in [Reference Keller-Ressel and Mayerhofer23]; that is, the region for p in which exponential moments of $Y_t$ exist coincides with the region for p such that the solution to (3) exists at a given time t. Consequently, we see that the blow-up condition of exponential moments at time t is specified as $\det F(p,t) = 0$ .

2.2. Large deviations principle

A large deviations principle for the Heston stochastic volatility model is nicely derived in [Reference Forde and Jacquier11]. The authors compute the limiting cumulant-generating function (CGF)

\begin{equation*}\Lambda(p) = \lim_{t \downarrow 0} t \ln \mathbb{E} \big[ \textrm{e}^{p (Y_t - Y_0)/t}\big]\end{equation*}

on some open interval $(p_-, p_+)$ including 0. Subsequently, its Fenchel–Legendre transform $\Lambda^*(x) = \sup_{p \in (p_-, p_+)} \{px - \Lambda(p)\}$ is shown to determine the small-time limit of implied volatilities under the Heston model.

Lemma 1. Suppose that RQ is symmetric in the WMSV model. We further assume that $\mathbb{I} - R^{\top} R$ is invertible. Then, the limiting cumulant-generating function $\Lambda(p)$ is given by

\begin{equation*} \Lambda(p) = \textrm{Tr}\left[ \left( \cos(p\Gamma) - \Gamma^{-1} \sin (p\Gamma)RQ \right)^{-1} \left( \frac{p}{2}\Gamma^{-1}\sin(p\Gamma)\right) \Sigma_0\right] \end{equation*}

on some open interval $(p_-, p_+)$ which includes 0. Here, $\Gamma$ is a square root of $Q^{\top} (\mathbb{I} - R^{\top} R) Q$ . The scalars $p_\pm$ are the largest negative and smallest positive numbers such that

\begin{equation*} \det\!\big(\!\cos(p\Gamma) - \Gamma^{-1} \sin (p\Gamma) RQ \big) = 0. \end{equation*}

Proof. By definition, the target function is given as follows:

\begin{align*} \Lambda(p) &= \lim_{t\rightarrow 0} t \ln \mathbb{E}\big[ \textrm{e}^{p(Y_t - Y_0)/t}\big] \\[3pt] &= \lim_{t \rightarrow 0} \Big\{ \textrm{Tr}\Big[ t A\Big(\frac{p}{t},t\Big)\Sigma_0\Big] + t b\left(\frac{p}{t},t\right)\Big\} \\[3pt] &= \textrm{Tr}\Big[ \lim_{t\rightarrow 0} t A\Big(\frac{p}{t},t\Big)\Sigma_0\Big] + \lim_{t\rightarrow 0} t b\Big(\frac{p}{t},t\Big), \end{align*}

provided the limits exist. We compute those limits by looking at the solution to the matrix differential equation in Proposition 1. It is easy to see that the solution is given as $ \left( \begin{array}{c@{\quad}c} G\left(\frac{p}{t},t\right) & F\left(\frac{p}{t},t\right) \end{array}\right) = \left( \begin{array}{c@{\quad}c} \textbf{0} & \mathbb{I} \end{array}\right) \textrm{e}^{t\textsf{M} + \textsf{N}_t}$ , with

\begin{equation*} \textsf{M} = \left( \begin{array}{c@{\quad}c} M & \textbf{0} \\ \textbf{0} & - M^{\top}\end{array}\right)\!,\qquad \textsf{N}_t = \left( \begin{array}{c@{\quad}c} p \textsf{A} & t \textsf{B} \\ \textsf{c}(t)\mathbb{I} & -p\textsf{A} \end{array}\right)\!, \end{equation*}

where $\textsf{A} = Q^{\top} R^{\top} = RQ$ , $\textsf{B} = - 2 Q^{\top} Q$ , and $\textsf{c}(t) = \frac{p(p-t)}{2t}$ . For notational simplicity, we simply write $\textsf{c}$ if there is no risk of confusion.

We will later show that $\textsf{M}$ does not affect the limit as t decreases to zero. Provided this is valid, it is enough to compute $\lim_{t \rightarrow 0} \textrm{e}^{\textsf{N}_t}$ . Straightforward computations yield that

\begin{equation*} \textsf{N}_t^{2k} = \left(\begin{array}{c@{\quad}c} \textsf{Z}_t^k & \Pi_k \\ \textbf{0} & \textsf{Z}_t^k\end{array}\right)\!, \qquad \textsf{N}_t^{2k+1} = \left(\begin{array}{c@{\quad}c} * & * \\ \textsf{c} \textsf{Z}_t^k & -p \textsf{Z}_t^k\textsf{A} \end{array}\right)\!, \end{equation*}

where $\textsf{Z}_t = p^2 \textsf{A}^2 + \textsf{c} t \textsf{B}$ and $\Pi_k$ is given by $\textsf{c} \Pi_k = p\big( \textsf{A} \textsf{Z}_t^{k} - \textsf{Z}_t^{k}\textsf{A}\big)$ , with $\Pi_1 = pt(\textsf{A}\textsf{B} - \textsf{B}\textsf{A})$ . Since $\lim_t \textsf{Z}_t = p^2 \textsf{A}^2 + \frac{p^2}{2} \textsf{B} =\!:\,\textsf{Z}_0$ , the matrix norm $|\textsf{Z}_t|$ is bounded by some constant for all small t. Consequently, by looking at the lower block matrices of $\textrm{e}^{\textsf{N}_t}$ , we obtain

\begin{align*} G\left(\frac{p}{t},t\right) &= \textsf{c} \sum_{k=0}^{\infty} \frac{1}{(2k+1)!} \textsf{Z}_t^k, \\ F\left( \frac{p}{t}, t\right) &= \sum_{k=0}^{\infty} \frac{1}{(2k)!} \textsf{Z}_t^k - p \sum_{k=0}^{\infty} \frac{1}{(2k+1)!} \textsf{Z}_t^k \textsf{A}. \end{align*}

Since the convergence of infinite sums is uniform on a finite interval, we get the following limits:

\begin{align*} G^*(p) \,:\!= \lim_{t \rightarrow 0} t G\left(\frac{p}{t},t\right) &= \frac{p^2}{2} \sum_{k=0}^{\infty} \frac{1}{(2k+1)!} \textsf{Z}_0^k, \\ F^*(p) \,:\!= \lim_{t \rightarrow 0} F\left( \frac{p}{t}, t\right) &= \sum_{k=0}^{\infty} \frac{1}{(2k)!} \textsf{Z}_0^k - p \sum_{k=0}^{\infty} \frac{1}{(2k+1)!} \textsf{Z}_0^k \textsf{A}. \end{align*}

On the other hand, we notice that $\textsf{Z}_0 = -p^2 Q^{\top}\left(\mathbb{I} - R^{\top} R\right)Q$ . Hence, it is represented by $-p^2 \Gamma^2$ . The final expression is easily obtained by the definitions of matrix sine and cosine.

Let us turn our attention to $\textrm{e}^{t\textsf{M} + \textsf{N}_t}$ . We note that

\begin{align*} &\textrm{e}^{t \textsf{M} + \textsf{N}_t} - \textrm{e}^{\textsf{N}_t} \\ &\quad= \sum_{l=0}^{\infty} \frac{1}{(2l)!}\left\{ \left( t\textsf{M} + \textsf{N}_t\right)^{2l} - \textsf{N}_t^{2l} \right\} + \sum_{l=0}^{\infty} \frac{1}{(2l+1)!}\left\{ \left( t\textsf{M} + \textsf{N}_t\right)^{2l+1} - \textsf{N}_t^{2l+1} \right\}\\ &\quad= \sum_{l=0}^{\infty} \frac{1}{(2l)!}\left\{ \left( t\textsf{M} + \textsf{N}_t\right)^{2l} - \textsf{N}_t^{2l} \right\} + \sum_{l=0}^{\infty} \frac{1}{(2l+1)!} \left(t \textsf{M} + \textsf{N}_t\right)^{2l} t\textsf{M} \\ &\qquad + \sum_{l=0}^{\infty} \frac{1}{(2l+1)!} \left\{ \left( t\textsf{M} + \textsf{N}_t\right)^{2l} - \textsf{N}_t^{2l} \right\} \textsf{N}_t. \end{align*}

It is easy to verify that the matrix norm $|(t\textsf{M} + \textsf{N}_t)^2|$ is bounded by a constant for all small t values. This implies that the second term converges to zero as t decreases to zero. Indeed, $(t\textsf{M} + \textsf{N}_t)^2 = t^2 \textsf{M}^2 + \textsf{M} (t\textsf{N}_t) + t\textsf{N}_t \textsf{M} + \textsf{N}_t^2$ , and each term is at most $\mathcal{O}(1)$ in t.

Regarding the first and third terms on the right-hand side, we are concerned with their lower left, say (i), and lower right, (ii), block matrices of size $n \times n$ as they are relevant to matrices G and F. Simple calculations show that

\begin{align*} \includegraphics{inline1.png}\\[-30pt]\end{align*}

As a result, it can be verified that (i) and (ii) are of $\mathcal{O}(1)$ and $\mathcal{O}(t)$ , respectively. Therefore, as t decreases to zero, (i) multiplied by t does not affect the limit $G^*(p)$ . Likewise, $F^*(p)$ is independent of (ii). Similar arguments can be made for the third summation. Hence, the limit in the statement is valid for any M.

Now we notice that $F^*(0) = \mathbb{I}$ , and thus there is a maximal open interval $(p_-, p_+)$ including zero such that $\det F^*(p) \neq 0$ . The equivalence of this and the finiteness of $\Lambda(p)$ comes from the limit

\begin{equation*} \Lambda(p) = \textrm{Tr}\left[ F^*(p)^{-1} G^*(p) \Sigma_0\right] = \frac{1}{\det F^*(p)} \times \text{some non-blow-up function of \textit{p}}, \end{equation*}

and $\det F^*(p) = \lim_t \det F\left( \frac{p}{t},t\right)$ , in addition to that the matrix A exists if and only if $\det F$ is nonzero. Hence, for $\Lambda(p)$ to be finite, each $\det F\left(\frac{p}{t},t\right)$ must be nonzero for all small t values with a nonzero limit. Alternatively, one may apply arguments similar to [Reference Forde and Jacquier11, Lemma B.1]. The last statement then follows.

So far we have not discussed the limit $\lim_{t \rightarrow 0} t b\left(\frac{p}{t},t\right)$ . However, from Proposition 1, we notice that

\begin{align*} \lim_{t \rightarrow 0} b\left(\frac{p}{t},t\right) &= -\frac{\beta}{2} \lim_t \textrm{Tr}\left[ \ln F\left( \frac{p}{t},t\right) + tM + p Q^{\top} R^{\top} \right] \\ &= -\frac{\beta}{2} \textrm{Tr}\left[ \ln F^*(p) + pQ^{\top} R^{\top}\right] \end{align*}

as long as $p \in (p_-,p_+)$ . As a result, this term does not appear in $\Lambda(p)$ . □

Remark 1. The above lemma imposes the assumption that RQ is symmetric. This condition is arguably mild. In [Reference Da Fonseca, Grasselli and Tebaldi8] it was shown that the stochastic correlation $\rho_t$ between the stock noise and the volatility noise depends on the product RQ according to

\begin{equation*}\rho_t = \frac{\textrm{Tr}[RQ\Sigma_t]}{\sqrt{\textrm{Tr}[\Sigma_t] \textrm{Tr}\left[Q^{\top} Q\Sigma_t\right]}\,}.\end{equation*}

Symmetric RQ does not affect the flexibility of this model to reflect the stochastic skew effect. The same condition appears in other related works. For instance, in [Reference Gnoatto and Grasselli18] or [Reference Da Fonseca5], the authors derive the explicit Laplace transform for the WMSV model under technical conditions including symmetric RQ.

The CGF $\ln \mathbb{E} \left[ \textrm{e}^{p(Y_t - Y_0)/t}\right]$ is convex in p for each t, by Hölder’s inequality. Therefore, its limit $\Lambda(p)$ is convex as well. On the other hand, its derivative $\Lambda'(p)$ can be computed as

(4) \begin{equation} \Lambda'(p) = \textrm{Tr}\left[ F^*(p)^{-1}\left( - {F^*}'(p) F^*(p)^{-1} G^*(p) + {G^*}'(p) \right) \Sigma_0\right] = \frac{f(p)}{\det F^*(p)^2}\end{equation}

for some non-blow-up function f. Here, the derivative of a matrix with respect to a scalar p is computed by componentwise differentiation. Recall that $G^*$ and $F^*$ are the limits of $tG\left(\frac{p}{t},t\right)$ and $F\left(\frac{p}{t},t\right)$ defined in Lemma 1. Consequently, $|\Lambda'(p)| \rightarrow \infty$ at the boundary of $(p_-,p_+)$ as long as f(p) does not vanish at the boundary points. Based on this essential smoothness of $\Lambda$ and its convexity, we obtain the next result.

Proposition 2. Assume that all of the conditions in Lemma 1 hold, and that f(p) of (4) does not vanish at $p_\pm$ . Then, $(Y_t - Y_0)$ satisfies the large deviation principle as t approaches zero, with the rate function $\Lambda^*(x) = \sup_{p \in (p_-,p_+)} \{px - \Lambda(p)\}$ for all $x \in \mathbb{R}$ . As a result, if $\ \textrm{Tr}[\Sigma_0] > 0$ , then

\begin{equation*} \Lambda^*(k) = \left\{ \!\!\!\begin{aligned} & - \lim_{t \rightarrow 0} t \ln \mathbb{P}\Big(Y_t - Y_0 > k\Big), & {if}\ k \geq 0, \\ & - \lim_{t \rightarrow 0} t \ln \mathbb{P}\Big(Y_t - Y_0 < k\Big), & {if}\ k \leq 0. \end{aligned} \right. \end{equation*}

Proof. The large deviations principle follows from the well-known Gärtner–Ellis theorem [Reference Dembo and Zeitouni9]. The Fenchel–Legendre transform $\Lambda^*$ is a good rate function. A version in our context for $\{Y_t - Y_0\}$ is

\begin{equation*} - \inf_{x > k} \Lambda^*(x) \leq \liminf_{t \rightarrow 0} t \ln \mathbb{P}( Y_t - Y_0 > k) \leq \limsup_{t \rightarrow 0} t \ln \mathbb{P}( Y_t - Y_0 > k) \leq - \inf_{x \geq k} \Lambda^*(x) \end{equation*}

for $k \geq 0$ . This version is similar to the one summarized in [Reference Forde, Jacquier and Mijatović14]. We refer the reader to an extended version available at the authors’ websites. Since $\Lambda$ is convex, $\Lambda^*$ is also convex. Furthermore, the essential smoothness of $\Lambda$ implies it is a closed proper convex function whose convex dual $\Lambda^*$ is closed and convex with $\Lambda^{**} = \Lambda$ . Theorem 26.3 of [Reference Rockafellar27] then tells us that $\Lambda^*$ is strictly convex.

Straightforward computations yield $F^*(0) = \mathbb{I}$ , $G^*(0) = {G^*}'(0)= \textbf{0}$ , and ${G^*}''(0) = \mathbb{I}$ , and this results in $\Lambda(0) = \Lambda'(0) = 0$ and $\Lambda''(0) = \textrm{Tr}[\Sigma_0]$ . Therefore, if $\textrm{Tr}[\Sigma_0] > 0$ there is a small interval around 0, say $\mathcal{I} = [-\varepsilon, \varepsilon]$ , on which $\Lambda'(p)$ is monotonically increasing and $\Lambda''(p)$ is strictly positive. Outside of $\mathcal{I} $ , the convexity of $\Lambda$ implies $\Lambda'$ is nondecreasing. Now let us consider $x \in [\Lambda'(-\varepsilon), \Lambda'(\varepsilon)]$ . Note that this last interval includes 0 in its interior. For such x, $\Lambda^* (x) = p^* x - \Lambda (p^*)$ , where $p^* = p^*(x)$ is a unique solution to the equation $\Lambda'(p) = x$ . It is clear that $p^*(0) = 0$ .

The nondegeneracy of $\Lambda''$ around 0 then implies the differentiability of $p^*(x)$ thanks to the implicit function theorem. Finally, we observe that

\begin{equation*} {\Lambda^*}'(x) = {p^*}'(x)x + p^*(x) - \Lambda'(p^*){p^*}'(x) = p^*(x), \qquad \Lambda^*(0) = {\Lambda^*}'(0) = 0. \end{equation*}

From the strict convexity of $\Lambda^*$ , we can conclude that ${\Lambda^*}'$ is positive on the positive real axis and negative on the negative real axis. In turn, we see that $\Lambda^*$ is nondecreasing on the positive real axis and nonincreasing on the negative real axis. As a consequence, $-\inf_{x \geq k} \Lambda^*(x) = - \inf_{x > k} \Lambda^*(x) = -\Lambda^*(k)$ for $k \geq 0$ . We can draw a similar conclusion for $k \leq 0$ . The desired statement is immediate. □

In the above proposition, the assumption that f(p) does not vanish at the boundary points is a mild one as it is expected that the two nonlinear functions f and $\det F^*$ do not share zeros in common settings. We note that the large deviations results above are the extension of [Reference Forde and Jacquier11, Theorem 2.1] to the WMSV model. Since their subsequent analyses are general in nature, we adopt those results in order to connect Proposition 2 to implied volatilities in the next subsection.

2.3. First-order expansion

For the reader’s convenience, we record some relevant results.

Proposition 3. (Forde and Jacquier [Reference Forde and Jacquier11].) Assume that all of the assumptions in Proposition 2 hold. Then, we have the following small-time behaviors of vanilla option prices and implied volatilities:

1. • for an out-of-the-money call option with the log-moneyness $x = \ln \frac{K}{S_0} \geq 0$ ,

   \begin{equation*} \Lambda^*(x) = - \lim_{t \rightarrow 0} t \ln C(S_0, K, t); \end{equation*}

2. • for an out-of-the-money put option with the log-moneyness $x \leq 0$ ,

   \begin{equation*} \Lambda^*(x) = - \lim_{t \rightarrow 0} t \ln P(S_0,K, t); \end{equation*}

3. • for the option implied volatility $\sigma_t(x)$ with the log-moneyness $x \neq 0$ ,

   \begin{equation*} I(x) = \lim_{t \rightarrow 0} \sigma_t(x) = \frac{|x|}{\sqrt{2\Lambda^*(x)}\,}. \end{equation*}

Here, $C(S_0,K,t)$ and $P(S_0, K, t)$ are the respective prices of call and put options with initial price $S_0$ , strike K, and maturity t.

Theorem 1. Assume that all of the conditions in Proposition 2 hold. Then, in some neighborhood of zero for the log-moneyness $x = \ln \frac{K}{S_0}$ , the following expansion is valid: with $y = \frac{x}{\textrm{Tr}[\Sigma_0]}$ ,

\begin{align*} I(x) &= \sqrt{\textrm{Tr}[\Sigma_0]} \bigg[ 1 + \frac12 \frac{\textrm{Tr}[RQ\Sigma_0]}{\textrm{Tr}[\Sigma_0]} y + \bigg( \frac16 \frac{\textrm{Tr}[Q^{\top} Q \Sigma_0]}{\textrm{Tr}[\Sigma_0]} \\ &\quad + \frac13 \frac{\textrm{Tr}\left[(RQ)^2 \Sigma_0\right]}{\textrm{Tr}[\Sigma_0]} - \frac34 \frac{\textrm{Tr}[RQ\Sigma_0]^2}{\textrm{Tr}[\Sigma_0]^2}\bigg) y^2 + \mathcal{O}(y^3)\bigg]. \end{align*}

Proof. The formula above is derived by tedious but straightforward computations of Taylor expansions. In the proof of Proposition 2, we already argued that $\Lambda'$ is smooth and strictly monotone in a neighborhood of $p=0$ . Its convex dual $\Lambda^*$ is then given by $\Lambda^*(x) =\break p^*(x)x - \Lambda(p^*(x))$ , where $p^*$ is a smooth solution to $\Lambda'(p) = x$ in a small neighborhood of $x = 0$ . Then, Taylor expansions yield

\begin{align*} \Lambda(p) &= \frac12 \textrm{Tr}[\Sigma_0] p^2 + \frac12 \textrm{Tr}[RQ\Sigma_0] p^3 + \frac16 \textrm{Tr}\left[ \left( Q^{\top} Q + 2(RQ)^2\right)\Sigma_0\right] p^4 + \mathcal{O}(p^5), \\ p^*(x) &= \frac{1}{\textrm{Tr}[\Sigma_0]} x - \frac32 \frac{\textrm{Tr}[RQ\Sigma_0]}{\textrm{Tr}[\Sigma_0]^3} x^2 + \mathcal{O}(x^3). \end{align*}

For notational convenience we denote $\Lambda(p) = p^2 \sum_{i=0}^2 a_i p^{i} + \mathcal{O}(p^5)$ and $p^*(x) = x \sum_{i=0}^2 b_i x^i + \mathcal{O}(x^4)$ . It can be shown that

\begin{equation*} b_2 = \frac92 \frac{\textrm{Tr}\left[ RQ\Sigma_0\right]^2}{\textrm{Tr}\left[\Sigma_0\right]^5} - \frac23 \frac{\textrm{Tr}\left[\left(Q^{\top} Q + 2(RQ)^2\right)\Sigma_0\right]}{\textrm{Tr}\left[\Sigma_0\right]^4}. \end{equation*}

Recall that in Proposition 3, for all sufficiently small nonzero log-moneyness x, we have $I(x)^2 = \frac{x^2}{2\Lambda^*(x)}$ . On the other hand, the above expansions give us

\begin{equation*} \frac{2\Lambda^*(x)}{x^2} = \frac{2}{x^2}\Big( p^*(x)x - \Lambda\left(p^*(x)\right)\!\Big) = c_0 + c_1 x + c_2 x^2 + \mathcal{O}(x^3) , \end{equation*}

for some constants $c_i$ . Those values are readily obtained by collecting relevant terms carefully as functions of the $a_i$ and $b_i$ . To be explicit, they are given by

\begin{equation*}c_0 = \frac{1}{\textrm{Tr}[\Sigma_0]}, \quad c_1 = - \frac{\textrm{Tr}[RQ\Sigma_0]}{\textrm{Tr}[\Sigma_0]^3}, \quad c_2 = \frac94 \frac{\textrm{Tr}[RQ\Sigma_0]^2}{\textrm{Tr}[\Sigma_0]^5} - \frac13 \frac{\textrm{Tr}\left[ \left(Q^{\top} Q + 2(RQ)^2\right)\Sigma_0\right]}{\textrm{Tr}[\Sigma_0]^4}. \end{equation*}

The reciprocal of the above expansion is not difficult to get:

\begin{align*} I(x)^2 &= \frac{1}{c_0} - \frac{c_1}{c_0^2} x + \left( \frac{c_1^2}{c_0^3} - \frac{c_2}{c_0^2}\right) x^2 + \mathcal{O}(x^3), \\ I(x) &= \frac{1}{\sqrt{c_0}\,} - \frac{c_1}{2 c_0^{1.5}} x + \frac{3c_1^2 - 4 c_0c_2}{8 c_0^{2.5}} x^2 + \mathcal{O}(x^3). \end{align*}

The formula in the statement then easily follows. □

The expansion in Theorem 1 reduces to that of the Heston model in [Reference Forde and Jacquier11, Theorem 3.2]. We also note that the above result, in particular $I(x)^2$ , is analogous to one of the main results in [Reference Da Fonseca and Grasselli7, Proposition 3.6] for the WMSV model. Their method is based on the direct expansion of the call option price in terms of the vol-of-vol scale factor $\alpha$ . In fact, the two formulae are the same if $\alpha = 1$ and RQ is symmetric. Nevertheless, we find the large deviations approach is useful because, first, the validity of the expansion is obtained without introducing an extra control parameter $\alpha$ , and second, information about the tail probability behaviors of $Y_t$ is provided via the convex dual of the limiting CGF $\Lambda(p)$ .

3. Small-time smiles

The program in [Reference Forde, Jacquier and Lee13] applies to the WMSV model as well. This extension is the main focus of the current section.

Lemma 2. Assume that all of the conditions in Proposition 2 hold. We further assume that the matrix M is symmetric. Then, for each $p \in (p_-, p_+)\backslash \{0\}$ , we have

\begin{equation*} \mathbb{E}[ \textrm{e}^{\frac{p}{t}(Y_t - Y_0)}] = U(p) \textrm{e}^{ \frac{\Lambda(p)}{t}} (1 + \mathcal{O}(t)) \end{equation*}

as t decreases to zero, where $\Lambda(p)$ is the limiting CGF and U(p) is given by

\begin{align*} &\ln U(p) \\ &\quad= \textrm{Tr}\Big[ {F^*}^{-1} \Big[ -\frac12 \sin p\Gamma \cdot \Gamma^{-1} + \frac12 (p\textsf{D}_1 - \textsf{D}_2) \Big] \Sigma_0\Big] \\ &\qquad +\ \textrm{Tr}\Big[ {F^*}^{-1} \Big[ \sin p\Gamma \cdot \Upsilon_0^{\prime} + \frac1p \sin p\Gamma \cdot \Gamma^{-1} M + \Big(\textsf{D}_1 - \frac1p\textsf{D}_2\Big)RQ \Big] {F^*}^{-1}G^* \Sigma_0 \Big] \\ &\qquad -\ \frac{\beta}{2}\textrm{Tr}\left[ \ln F^* + p RQ\right]. \end{align*}

Here, $F^*(p)$ , $G^*(p)$ , and $\Gamma$ are as in Lemma 1. The new symbols $\textsf{D}_1$ , $\textsf{D}_2$ , and $\Upsilon_0^{\prime}$ satisfy

\begin{align*} \textsf{D}_1 &= \cos p\Gamma \cdot \Upsilon_0^{\prime} \Gamma^{-1}, \\ \textsf{D}_2 &= \sin p\Gamma \cdot \Gamma^{-1}\Upsilon_0^{\prime} \Gamma^{-1}, \\ \textrm{vec} \, \Upsilon_0^{\prime} &= - (\Gamma \oplus \Gamma)^{-1} \textrm{vec}( MRQ + RQ M + Q^{\top} Q ) , \end{align*}

where vec is the vectorization operator and $\oplus$ means the Kronecker sum.

Proof. In the proof of Lemma 1, we showed that the matrix M does not affect the limiting CGF. However, for small-time smiles, we need to calculate the behaviors of the CGF for small t values, which turn out to depend on M. Let us repeat the arguments of the lemma as follows. Recall that the solution $A\left(\frac{p}{t}, t\right)$ is $F^{-1}G$ , where

\begin{equation*} \left( \begin{array}{c@{\quad}c} G\left(\frac{p}{t},t\right) & F\left(\frac{p}{t},t\right) \end{array}\right) = \left( \begin{array}{c@{\quad}c} \textbf{0} & \mathbb{I} \end{array}\right) \textrm{e}^{\textsf{L}_t}, \qquad \textsf{L}_t = \left( \begin{array}{c@{\quad}c} tM+p\textsf{A} & t\textsf{B} \\ \textsf{c}\mathbb{I} & - (tM + p\textsf{A})\end{array}\right)\!,\end{equation*}

where $\textsf{A} = RQ$ , $\textsf{B} = -2Q^{\top} Q$ , and $\textsf{c} = \frac{p^2}{2t} - \frac{p}{2}$ . For notational convenience, we denote $(tM + p\textsf{A})^2 + \textsf{c} t \textsf{B}$ by $\textsf{Z}_t$ . Then, thanks to the symmetry of M, it can be readily shown that

\begin{equation*}\textsf{L}_t^{2k} = \left( \begin{array}{c@{\quad}c} \textsf{Z}_t^k & \Pi_k \\ \\[-6pt] \textbf{0} & \textsf{Z}_t^k \end{array}\right)\end{equation*}

for some matrix $\Pi_k$ . More explicitly, it is given by $\textsf{c}\Pi_k = (tM + p\textsf{A}) \textsf{Z}_t^k - \textsf{Z}_t^k (tM+p\textsf{A})$ . We can rewrite

\begin{equation*}\textsf{Z}_t = \textsf{Z}_0 + t \textsf{C} + t^2 M^2, \qquad \textsf{Z}_0 = p^2 \textsf{A}^2 + \frac{p^2}{2}\textsf{B}, \qquad \textsf{C} = p\Big(M\textsf{A} + \textsf{A} M - \frac12 \textsf{B}\Big).\end{equation*}

From these calculations, we can find the matrices G and F as

\begin{align*} G\left(\frac{p}{t},t\right) &= \textsf{c} \sum_{k=0}^{\infty} \frac{1}{(2k+1)!} \textsf{Z}_t^k, \\ F\left(\frac{p}{t},t\right) &= \sum_{k=0}^{\infty} \frac{1}{(2k)!} \textsf{Z}_t^k - \sum_{k=0}^{\infty} \frac{1}{(2k+1)!} \textsf{Z}_t^k(tM + p\textsf{A}).\end{align*}

Under the assumed conditions, $-\textsf{Z}_t$ converges to a nonsingular $-\textsf{Z}_0 = (p\Gamma)^2$ as t approaches zero. Here, $\Gamma$ is as given in Lemma 1. Since $\textsf{Z}_t$ is symmetric, we can find a nonsingular $\Upsilon_t$ such that $\textsf{Z}_t = - \Upsilon_t^2$ for some symmetric and nonsingular $\Upsilon_t$ , and for all sufficiently small t. Clearly, $\lim_{t \rightarrow 0} \Upsilon_t = p\Gamma$ . Consequently, more compact representations are possible:

\begin{align*} G\left(\frac{p}{t},t\right) &= \textsf{c} \sin \Upsilon_t \cdot \Upsilon_t^{-1}, \\ F\left(\frac{p}{t},t\right) &= \cos \Upsilon_t - \sin \Upsilon_t \cdot \Upsilon_t^{-1}(tM+p\textsf{A}).\end{align*}

Let us now make the next observations: direct differentiations yield

\begin{align*} G\left(\frac{p}{t},t\right) &= \textsf{c} \sin \Upsilon_0 \cdot \Upsilon_0^{-1} + \textsf{c} \textsf{D} t + \textsf{c}\mathcal{O}(t^2)\\ &= \frac1t G^*(p) - \frac{p}{2}\sin \Upsilon_0 \cdot \Upsilon_0^{-1} + \frac{p^2}{2} \textsf{D} + \mathcal{O}(t), \\ F\left(\frac{p}{t},t\right) &= F^*(p) - \sin \Upsilon_0 \cdot \Upsilon_0^{\prime} t - \textsf{D} t(tM+p\textsf{A}) - \sin \Upsilon_0 \cdot \Upsilon_0^{-1} M t + \mathcal{O}(t^2) \\ &= F^*(p) - \left[ \sin \Upsilon_0 \cdot \Upsilon_0^{\prime} + \sin \Upsilon_0 \cdot \Upsilon_0^{-1} M + p \textsf{D}\textsf{A} \right] t + \mathcal{O}(t^2).\end{align*}

Here, $\textsf{D}$ is defined as $\frac{\rd}{\rd t} \sin \Upsilon_t \cdot \Upsilon_t^{-1} \Big|_{t=0}$ and computed as

\begin{equation*}\textsf{D} = \cos\Upsilon_0 \cdot \Upsilon_0^{\prime} \Upsilon_0^{-1} - \sin \Upsilon_0 \cdot \Upsilon_0^{-1} \Upsilon_0^{\prime} \Upsilon_0^{-1}.\end{equation*}

We note that $\Upsilon_0^{\prime}$ satisfies the special case of Sylvester’s equation, $\Upsilon_0^{\prime} \Upsilon_0 + \Upsilon_0 \Upsilon_0^{\prime} = -\textsf{C}$ , by differentiating the defining equation of $\Upsilon_t^2$ . Thanks to the symmetry and the positive definiteness of $\Upsilon_t$ for small t, it is known that there exists a unique solution $\Upsilon_0^{\prime}$ , and the solution is obtained by $\text{vec} \, \Upsilon_0^{\prime} = - \left( \Upsilon_0 \otimes \mathbb{I} + \mathbb{I} \otimes \Upsilon_0 \right)^{-1} \text{vec} \, \textsf{C}$ , where $\otimes$ is the Kronecker product operator and vec is the vectorization operator.

Using the definition of Kronecker sum, we get $\Upsilon_0 \oplus \Upsilon_0 = \Upsilon_0 \otimes \mathbb{I} + \mathbb{I} \otimes \Upsilon_0$ . It is then immediate to see that $\text{vec} \, \Upsilon_0^{\prime} = - (\Gamma \oplus \Gamma)^{-1} \text{vec}( M\textsf{A} + \textsf{A} M - \frac12 \textsf{B})$ . We also write, for notational simplicity,

\begin{eqnarray*}\textsf{D} &\,=\,& \frac{1}{p} \textsf{D}_1 - \frac{1}{p^2} \textsf{D}_2, \qquad \left\{\begin{split} \textsf{D}_1 = &\cos (p\Gamma) \Upsilon_0^{\prime} \Gamma^{-1} , \\ \textsf{D}_2 = &\sin (p\Gamma) \Gamma^{-1} \Upsilon_0^{\prime}\Gamma^{-1}.\end{split}\right.\end{eqnarray*}

Finally, it only takes several simple algebraic manipulations until we arrive at

\begin{align*} A\left(\frac{p}{t},t\right) &= \frac1t {F^*}^{-1}G^* + {F^*}^{-1} \Big[ -\frac12 \sin p\Gamma \cdot \Gamma^{-1} + \frac12 (p\textsf{D}_1 - \textsf{D}_2) \Big] \\ &\quad +\ {F^*}^{-1} \Big[ \sin p\Gamma \cdot \Upsilon_0^{\prime} + \frac1p \sin p\Gamma \cdot \Gamma^{-1} M + \Big(\textsf{D}_1 - \frac1p\textsf{D}_2 \Big)\textsf{A}\Big] {F^*}^{-1}G^* \\&\quad +\ \mathcal{O}(t).\end{align*}

This first-order expansion of A combined with the known formula $b\left(\frac{p}{t},t\right)$ in the proof of Lemma 1 helps us get

(□) \begin{align} \ln \mathbb{E} & \left[ \textrm{e}^{\frac{p}{t}(Y_t- Y_0)}\right] \nonumber\\ & = \frac1t \Lambda(p) + \textrm{Tr}\Big[ {F^*}^{-1} \Big[ -\frac12 \sin p\Gamma \cdot \Gamma^{-1} + \frac12 (p\textsf{D}_1 - \textsf{D}_2) \Big] \Sigma_0\Big] \nonumber\\ & \quad + \textrm{Tr}\left[ {F^*}^{-1} \left[ \sin p\Gamma \cdot \Upsilon_0^{\prime} + \frac1p \sin p\Gamma \cdot \Gamma^{-1} M + \left(\textsf{D}_1 - \frac1p\textsf{D}_2 \right)\textsf{A}\right] {F^*}^{-1}G^* \Sigma_0 \right] \nonumber\\ & \quad - \frac{\beta}{2}\textrm{Tr}\left[ \ln F^* + p\textsf{A}\right] + \mathcal{O}(t). \end{align}

It is possible to obtain an asymptotic expansion of $\ln U(p)$ which becomes handy for approximating small-time smiles. Since the derivation is tedious and long, we briefly sketch some relevant computations.

Corollary 1. Assume that all of the conditions in Lemma 2 hold. Then, in a small neighborhood of zero for p, the following asymptotic expansion is valid:

\begin{align*} \ln U(p) &= - \frac12 \textrm{Tr}[\Sigma_0] p + \left\{ - \frac12 \textrm{Tr}[RQ \Sigma_0] +\frac12 \textrm{Tr}\left[ M\Sigma_0\right] + \frac{\beta}{4}\textrm{Tr}\left[ Q^{\top} Q\right] \right\} p^2 \\ &\quad + \left\{ \textrm{Tr}\left[ \left(- \frac16 \Gamma^2 - \frac16 \Gamma^2 \Upsilon_0^{\prime} \Gamma^{-1} - \frac12 (RQ)^2 \right)\Sigma_0 \right] \right. \\ &\quad \left. +\ \frac12 \textrm{Tr}\left[ \left(\Gamma \Upsilon_0^{\prime} + MRQ + RQM\right) \Sigma_0\right] + \frac{\beta}{6} \textrm{Tr}\left[ Q^{\top} Q RQ \right] \right\} p^3 + \mathcal{O}(p^4) , \end{align*}

where $\Gamma$ and $\Upsilon_0^{\prime}$ are as in Lemma 2.

Proof. Useful expansions are given as follows:

\begin{align*} {F^*}^{-1} &= \mathbb{I} + \textsf{A} p + \left( \textsf{A}^2 + \frac12 \Gamma^2\right) p^2 + \mathcal{O}(p^3), \\ G^* &= \frac{p^2}{2} \mathbb{I} - \frac{1}{12} \Gamma^2 p^4 + \mathcal{O}(p^6), \\ \sin p\Gamma &= \Gamma p - \frac16 \Gamma^3 p^3 + \mathcal{O}(p^5), \\ \cos p\Gamma &= \mathbb{I} - \frac12 \Gamma^2 p^2 + \frac{1}{24} \Gamma^4 p^4 + \mathcal{O}(p^6).\end{align*}

Then, expansions for $\textsf{D}_1$ and $\textsf{D}_2$ can be calculated accordingly. Regarding the last term in $\ln U(p)$ of Lemma 2, we do the formal expansion of $\ln F^* = \alpha_0 + \alpha_1 p + \cdots$ and find matching coefficients via

\begin{equation*}\textrm{e}^{\ln F^*} = F^* = \mathbb{I} - \textsf{A} p - \frac12 \Gamma^2 p^2 + \frac16 \Gamma^2 \textsf{A} p^3 + \mathcal{O}(p^4).\end{equation*}

This leads to $\alpha_0 = \textbf{0}$ , $\alpha_1 = - \textsf{A}$ , $\alpha_2 = - \frac12 \Gamma^2 - \frac12 \textsf{A}^2$ , and $\alpha_3 = - \frac{1}{12} \Gamma^2 \textsf{A} - \frac14 \textsf{A} \Gamma^2 - \frac13 \textsf{A}^3$ . Then, we eventually find that

\begin{align*} - \frac{\beta}{2} \textrm{Tr}\left[ \ln F^* + \textsf{A} p \right] &= - \frac{\beta}{2} \textrm{Tr}\left[ \alpha_2 p^2 + \alpha_3 p^3 \right] + \mathcal{O}(p^4) \\ &= \frac{\beta}{4} \textrm{Tr}\left[ Q^{\top} Q\right] p^2 + \frac{\beta}{6} \textrm{Tr}\left[ Q^{\top} Q \textsf{A}\right] p^3 + \mathcal{O}(p^4).\end{align*}

Carefully collecting relevant terms, we obtain the desired result. □

Proof. Under the given assumptions, it is argued in the proof of Proposition 2 that there is a small interval $x \in [\Lambda'(- \varepsilon), \Lambda'(\varepsilon)]$ for a small $\varepsilon>0$ with the following properties. First, this interval includes 0 in its interior. Second, $\Lambda''$ is strictly positive on $[-\varepsilon, \varepsilon]$ . Third, for such x, there is a unique solution $p^*(x)$ to $\Lambda'(p) = x$ thanks to the strict monotonicity of $\Lambda'$ .

Let us denote $z = p^*+\ri q$ . The CGF of $(Y_t - Y_0)$ is defined as $\Lambda_t\left(\frac{z}{t}\right) = \ln \mathbb{E}[ \textrm{e}^{\frac{z}{t}(Y_t - Y_0)}]$ . Then, observe that

\begin{equation*} \left| \mathbb{E}\left[ \textrm{e}^{\frac{z}{t}(Y_t - Y_0)}\right] \right| = \left| \exp\left(\textrm{Re}\,\Lambda_t\left(\frac{z}{t}\right) + \ri \, \textrm{Im}\,\Lambda_t\left(\frac{z}{t}\right)\right) \right| = \exp\left(\textrm{Re}\,\Lambda_t\left(\frac{z}{t}\right)\right)\!. \end{equation*}

On the other hand, the left-hand side of the above equation is clearly less than or equal to

\begin{equation*} \mathbb{E}\Big[ \Big| \textrm{e}^{\frac{z}{t}(Y_t - Y_0)} \Big|\Big] = \mathbb{E}\Big[ \textrm{e}^{\frac{p^*}{t}(Y_t - Y_0)}\Big] = \exp\Big( \Lambda_t\Big(\frac{p^*}{t}\Big) \Big). \end{equation*}

Since this holds for all t, we obtain the limiting result after multiplying t, $\textrm{Re} \, \Lambda(p^*+\ri q) \leq \Lambda(p^*)$ , making $q = 0$ a maximum.

In order to see that $q = 0$ is a unique maximum, we note that

\begin{equation*} h''(0) = \textrm{Re} \Big( \Lambda''(p^*+\ri q) \ri^2 \Big) \Big|_{q = 0} = - \Lambda''(p^*). \end{equation*}

The strict convexity of $\Lambda$ on $[-\varepsilon, \varepsilon]$ implies that $- \Lambda''(p^*) < 0$ , so that the function h becomes strictly concave at $q = 0$ . This proves the uniqueness of the maximum. □

Based on the two lemmas above, the procedure of [Reference Forde, Jacquier and Lee13] can be applied. For the reader’s convenience, we summarize its outline. Assume that the log-moneyness $x = \ln \frac{K}{S_0}$ is in the small neighborhood of Lemma 3. Then, for sufficiently small t, we have

\begin{equation*} \frac{1}{S_0} \mathbb{E}\left[\left(\textrm{e}^{Y_t} - S_0 \textrm{e}^x\right)^+\right] = (1 - \textrm{e}^x) \textbf{1}_{\{x < 0\}} - \frac{\textrm{e}^x t}{2\pi} \textrm{Re}\left( \int_{\mathscr P} \textrm{e}^{\ri xu/t} \phi_t\left(- \frac{u}{t}\right)\left(\frac{1}{u^2} + \mathcal{O}(t)\right) \rd u\right)\!,\end{equation*}

where $\phi_t(z) = \mathbb{E}\left[ \textrm{e}^{\ri z(Y_t - Y_0)}\right]$ for a complex number z with $- \textrm{Im}(z) \in (p_-, p_+)$ . The integration path ${\mathscr P}$ is from $-\infty + \ri p^*(x)$ to $\infty + \ri p^*(x)$ . Then, Lemma 2 implies that

\begin{align*} & \int_{\mathscr P} \textrm{e}^{\ri xu/t} \phi_t\left(- \frac{u}{t}\right)\left(\frac{1}{u^2} + \mathcal{O}(t)\right) \rd u \\ &\quad= \int_{\mathscr P} \textrm{e}^{\ri xu/t} U(-\ri u) \textrm{e}^{\Lambda(-\ri u)/t} (1 + \mathcal{O}(t))\left(\frac{1}{u^2} + \mathcal{O}(t)\right) \rd u \\ &\quad= \int_{\mathscr P} \textrm{e}^{- H(u)/t} \frac{U(-\ri u)}{u^2} \rd u (1 + \mathcal{O}(t)) ,\end{align*}

where $H(u) = - \ri xu - \Lambda(-\ri u)$ . With nonzero x, Lemma 3 is then utilized to apply the Laplace expansion of [Reference Olver25, Theorem 7.1, Chapter 4]. More specifically, $\textrm{Re}(H(u) - H(u_0))$ is positive on ${\mathscr P}$ for $u_0 = \ri p^*(x)$ . The first-order expansion of the resulting formula is

\begin{equation*}2 \sqrt{\pi t} \textrm{e}^{-H(u_0)/t} \frac{U(-\ri u_0)}{u_0^2 \sqrt{2F''(u_0)}\,}(1 + \mathcal{O}(t)) = - 2\sqrt{\pi t} \textrm{e}^{-(xp^* - \Lambda(p^*))/t} \frac{U(p^*)}{ {p^*}^2 \sqrt{2\Lambda''(p^*)}\,} (1 + \mathcal{O}(t)).\end{equation*}

By definition of the convex dual of $\Lambda$ , we see that $xp^* - \Lambda(p^*) = \Lambda^*(x)$ . The next result is simply the WMSV version of [Reference Forde, Jacquier and Lee13, Theorem 3.1].

Proposition 4. Assume that all of the conditions in Lemma 2 hold. Then, for a nonzero log-moneyness $x = \ln \frac{K}{S_0}$ in a small neighborhood of zero, the asymptotic behavior for European call options is given by

\begin{equation*} \frac{1}{S_0} \mathbb{E}\left[\left(\textrm{e}^{Y_t} - S_0 \textrm{e}^x\right)^+\right] = (1 - \textrm{e}^x) \textbf{1}_{\{x < 0\}} + \textrm{e}^{- \Lambda^*(x)/t} \left( \frac{A(x)}{\sqrt{2\pi}\,} t^{1.5} + \mathcal{O}(t^{2.5}) \right)\!, \end{equation*}

where $A(x) = \frac{\textrm{e}^x U(p^*(x))}{p^*(x)^2 \sqrt{\Lambda''(p^*(x))}\,}$ as t decreases to zero.

Comparison of the above expansion with the call price expansion under the Black–Scholes model leads us to the asymptotic expansion of the small-time smile as in [Reference Forde, Jacquier and Lee13, Theorem 4.2]: for nonzero and small x,

\begin{equation*}\sigma_t^2(x) \approx I(x)^2 + \frac{2 I(x)^4}{x^2} \left[ \ln \frac{A(x) x^2}{I(x)^3} - \frac{x}{2}\right]t.\end{equation*}

Here, I(x) is as in Proposition 3. Since this step is straightforward, we refer the reader to the main reference. This expression still depends on the implicitly defined function $p^*(x)$ . Nevertheless, we can utilize the asymptotic results obtained so far. Specifically, recall that we derived expansions for $\Lambda(p)$ , $p^*(x)$ , and I(x) previously. By plugging in estimates for $\Lambda''$ and U using Lemma 1 and Corollary 1, we get an approximation to A(x). And Theorem 1 provides us with an approximation to I(x). Alternatively, an asymptotic expression can be produced by expanding the necessary quantities to suitable degrees.

Theorem 2. Assume that all of the conditions in Lemma 2 hold. Then, for a nonzero log-moneyness $x = \ln \frac{K}{S_0}$ in a small neighborhood of zero, the small-time smile is approximated by

(5) \begin{equation} \sigma_t^2 \approx I(x)^2 + 2 I(x)^4 (d_0+d_1x) t \end{equation}

for some constants $d_0$ and $d_1$ . In particular, $d_0$ is equal to

\begin{equation*} \frac14 \frac{\textrm{Tr}\left[RQ\Sigma_0\right]}{\textrm{Tr}\left[ \Sigma_0\right]^2} + \frac12 \frac{\textrm{Tr}\left[ M \Sigma_0\right]}{\textrm{Tr}\left[\Sigma_0\right]^2} + \frac{\beta}{4} \frac{\textrm{Tr}\left[ Q^{\top} Q\right]}{\textrm{Tr}\left[\Sigma_0\right]^2} + \frac38 \frac{\textrm{Tr}\left[RQ\Sigma_0\right]^2}{\textrm{Tr}\left[\Sigma_0\right]^4} - \frac16 \frac{\textrm{Tr}\left[ \left(Q^{\top} Q + 2(RQ)^2\right)\Sigma_0\right]}{\textrm{Tr}\left[\Sigma_0\right]^3}. \end{equation*}

Proof. For notational convenience we write $\ln U(p) = u_0 p + u_1 p^2 + u_2 p^3 + \mathcal{O}(p^4)$ . We also retrieve the expressions $\Lambda(p) = \sum_{i=0}^3 a_i p^{i+2} + \mathcal{O}(p^6)$ and $p^*(x) = \sum_{i=0}^3 b_i x^{i+1} + \mathcal{O}(x^5) = x \tilde p(x)$ . Lastly, $I(x)^{-2} = \frac{2\Lambda^*(x)}{x^2} = \sum_{i=0}^3 c_i x^i + \mathcal{O}(x^4)$ . Then, we have

\begin{align*} &\left[ \ln \frac{A(x) x^2}{I(x)^3} - \frac{x}{2}\right] \\ &\quad= \frac{x}{2} + \ln U(p^*) + 2 \ln x - 2 \ln p^* - \frac12 \ln \Lambda''(p^*) - 3 \ln I(x)\\ &\quad= \frac{x}{2} + \ln U(p^*) - 2 \ln \tilde p(x) + \frac12 \ln {p^*}'(x) + \frac32 \ln I(x)^{-2} \\ &\quad= \frac{x}{2} + u_0(b_0 x + b_1 x^2 + b_2 x^3) + u_1 (b_0^2 x^2 + 2b_0b_1 x^3) + u_2 b_0^3 x^3 + \mathcal{O}(x^4) \\ &\qquad -\ 2 \ln \big[ b_0 + b_1 x + b_2 x^2 + b_3 x^3 + \mathcal{O}(x^4)\big] \\ &\qquad +\ \frac12 \ln \big[ b_0 + 2 b_1 x + 3 b_2 x^2 + 4 b_3 x^3 + \mathcal{O}(x^4)\big]\\ &\qquad +\ \frac32 \ln \big[ c_0 + c_1 x + c_2 x^2 + c_3 x^3 + \mathcal{O}(x^4)\big].\end{align*}

Here we utilized the relation $\Lambda'(p^*) = x$ so that $\Lambda''(p^*){p^*}'(x) = 1$ by differentiating with respect to x. Our previous computations show us that $a_0 = \frac12 \textrm{Tr}[\Sigma_0]$ and $b_0 = c_0 = \frac{1}{\textrm{Tr}[\Sigma_0]}$ . We also have $u_0 = - \frac12 \textrm{Tr}[\Sigma_0]$ . This simplifies the formula a bit and leads us to

\begin{align*} \bigg[ \ln \frac{A(x) x^2}{I(x)^3} - \frac{x}{2}\bigg] &= u_0(b_1 x^2 + b_2 x^3) + u_1 (b_0^2 x^2 + 2b_0b_1 x^3) + u_2 b_0^3 x^3 \\ &\quad -\ 2 \ln \left[ 1 + \frac{b_1}{b_0} x + \frac{b_2}{b_0} x^2 + \frac{b_3}{b_0} x^3 + \mathcal{O}(x^4)\right] \\ &\quad +\ \frac12 \ln \left[ 1 + 2 \frac{b_1}{b_0} x + 3 \frac{b_2}{b_0} x^2 + 4 \frac{b_3}{b_0} x^3 + \mathcal{O}(x^4)\right] \\ &\quad +\ \frac32 \ln \left[ 1 + \frac{c_1}{b_0} x + \frac{c_2}{b_0} x^2 + \frac{c_3}{b_0} x^3 + \mathcal{O}(x^4)\right].\end{align*}

It is easy to check that $c_1 = \frac23 b_1$ . Hence, the Taylor expansions of the log functions cancel out x terms because $- 2 b_1 + b_1 + \frac32 c_1 = 0$ . Finally, $\frac{1}{x^2} \left[ \ln \frac{A(x) x^2}{I(x)^3} - \frac{x}{2}\right]$ boils down to $d_0 + d_1 x + \mathcal{O}(x^2)$ for some suitable constants $d_0$ and $d_1$ .

Carefully collecting relevant terms and utilizing the relationships $b_1 = \frac32 c_1$ and $b_2 = 2c_2$ , we can see that $d_0 = u_0 b_1 + u_1 b_0^2 + \frac{b_2}{4 b_0} - \frac{b_1^2}{3 b_0^2}$ , from which we obtain the expression in the statement. □

In the previous theorem it is possible to compute $d_1$ explicitly, in which case $b_3$ and $c_3$ are required. Depending on the timescale for which the user applies implied volatility expansions, however, the expansion $I(x)^2 + 2I(x)^4 d_0 t$ might suffice. In our numerical tests in the next section we use $d_0$ only with the quadratic expansion in Theorem 1 for I(x).

4. Numerical results

In this section we conduct several numerical comparisons to test the effectiveness of the derived formulae in Theorems 1 and 2. We consider a book of European call prices provided by Wharton Research Data Services on the S&P 500 Index (SPX) quoted on June 7, 2019. The time to maturity ranges from two weeks to four weeks. We only use option quotes with log-moneyness between $-0.05$ and $0.05$ , since near-the-money options tend to be more liquid, especially for short-maturity options. Finally, this gives us a sample of 174 option prices.

Let us denote the limiting implied volatility I(x) by $\sigma_0$ , the implied volatility under the WMSV model by $\sigma^{\text{model}}_t$ , the market implied volatility by $\sigma^{\text{market}}_t$ , and the small-time smile approximation in Theorem 2 by $\sigma^{\text{approx}}_t$ . One contribution of Theorem 2 lies in a simple calibration of the WMSV model parameters by minimizing the mean square error (MSE) of $\sigma^{\text{approx}}_t$ , i.e.

(6) \begin{equation}\min \frac{1}{N}\sum_{j=1}^{m}\sum_{i=1}^{n_j}\Big(\sigma^{\text{market}}_{t_j}(S_0,K_i)-\sigma^{\text{approx}}_{t_j}(S_0,K_i)\Big)^2 ,\end{equation}

where $N=n_1+n_2+\cdots+n_m$ , m is the number of maturities, and $n_j$ is the number of options with maturity $t_j$ . The calibration results based on this method are as follows:

\begin{equation*}\Sigma_0 = \Bigg[\begin{array}{c@{\quad}c} 0.0107 & -4.9502\times10^{-4} \\-4.9502\times10^{-4} & 0.0107\end{array}\Bigg],\qquad M = \Bigg[\begin{array}{rr} -14.6632& -0.0295 \\ -0.0295 & -14.6625\end{array}\Bigg] ,\end{equation*}

\begin{equation*}Q = \Bigg[\begin{array}{r@{\quad}r} 0.1559 & 0.2118 \\ 0.2094 & 0.1548\end{array}\Bigg],\qquad R = \Bigg[\begin{array}{r@{\quad}r} -1.5915 & -0.0977 \\ -0.1205 & -1.5919\end{array}\Bigg], \qquad \beta = 5.9386 .\end{equation*}

Based on those calibrated parameters, we find that the MSE in implied volatility is $3.1692\times 10^{-5}$ , and the MSE in price (normalized by $S_0$ ) is $9.2544\times 10^{-8}$ . We further illustrate the difference between model prices for the above parameters and market prices in Figure 1, and compare model implied volatilities $\sigma^{\text{model}}_t$ and market implied volatilities $\sigma^{\text{market}}_t$ in Figure 2. To obtain the implied volatility $\sigma_t^{\text{model}}$ , we first numerically calculate option prices, denoted by $C^{\text{model}}$ , via the widely used Fourier transform method for the WMSV model, and then find $\sigma_t^{\text{model}}$ with which the Black–Scholes price is identical to the price $C^{\text{model}}$ . For Fourier methods, we refer the interested reader to any of the many available references, e.g. [Reference Da Fonseca, Grasselli and Tebaldi8, Reference Carr and Madan2, Reference Duffie, Pan and Singleton10]. Nevertheless, for the reader’s convenience, we document the details of our MATLAB^® implementation, which can be found on the corresponding author’s personal web page or can be requested via email.

Figure 1: Comparison of the WMSV model prices and market prices of the SPX call options with three different maturities quoted on June 7, 2019.

Figure 2: Comparison of implied volatilities derived from the WMSV model prices and market prices of the SPX call options quoted on June 7, 2019.

Since the original paper [Reference Da Fonseca, Grasselli and Tebaldi8] presents detailed sensitivity analyses of parameters on option prices, we focus on the numerical performance of our implied volatility expansions. Figure 3 compares the small-time smile approximations with the implied volatility $\sigma_t^{\text{model}}$ with respect to the log-moneyness $x\in[-0.05, 0.05]$ based on the calibrated parameters listed above. To take a close look at the differences near the money we provide Figure 4, which focuses on the range $x\in[-0.01,0.01]$ . Also, using the same range for x, Figure 5 presents the relative errors of $\sigma_0$ , the relative errors of $\sigma_t$ , and the relative errors of the Black–Scholes prices based on $\sigma_t^{\text{approx}}$ with respect to the model price $C^{\text{model}}$ . These figures clearly confirm our mathematical results in Theorems 1 and 2. Firstly, it turns out that in most cases $\sigma_t^{\text{approx}}$ outperforms $\sigma_0$ in terms of approximating $\sigma_t^{\text{model}}$ . Secondly, as time to maturity decreases, both $\sigma_t^{\text{approx}}$ and $\sigma_0$ tend to be closer to the model implied volatility. Lastly, our approximation works well for near-the-money options.

Figure 3: Comparison of model implied volatility $\sigma^{\text{model}}_t$ , limiting implied volatility $\sigma_0$ , and small-time smile approximation $\sigma^{\text{approx}}_t$ for the three shortest maturities where log-moneyness $x\in[-0.05,0.05]$ .

Figure 4: Comparison of model implied volatility $\sigma^{\text{model}}_t$ , limiting implied volatility $\sigma_0$ , and small-time smile approximation $\sigma^{\text{approx}}_t$ for the three shortest maturities where log-moneyness $x\in[-0.01,0.01]$ .

Figure 5: Relative differences |model_val - approx_val|/model_val for (i) limiting implied volatility $\sigma_0$ , (ii) small-time smile $\sigma^{\text{approx}}_t$ , and (iii) Black–Scholes prices from $\sigma^{\text{approx}}_t$ for the three shortest maturities where log-moneyness $x\in[-0.01,0.01]$ .

One may consider calibrating model parameters using the existing approaches in the literature. The most widely adopted method is to minimize the weighted mean square error of numerically calculated prices $C^{\text{model}}$ , i.e.

(7) \begin{equation}\min \frac{1}{N}\sum_{j=1}^{m}\sum_{i=1}^{n_j}\nu_{ij}^{-2}\Big(C^{\text{market}}(S_0,K_i,t_j)-C^{\text{model}}(S_0,K_i,t_j)\Big)^2 ,\end{equation}

where $\nu_{ij}$ is the vega of the corresponding option, and $C^{\text{market}}$ stands for the option’s market price. As mentioned in [Reference Da Fonseca and Grasselli7], using $\nu_{ij}^{-2}$ for the weight is the market practice and improves the calibration by putting more weight on short-maturity options. However, the calibration errors based on (7) show that our approximation (5) contributes to better calibration of the WMSV model parameters. In particular, compared to our method (6), the Fourier-based calibration (7) turns out to increase both the MSE in implied volatility and MSE in price by 9.06% and 43.56%, respectively. Also, it is worth noting that the optimization problem (6) requires substantially less computational cost than the problem (7). Therefore, at least for near-the-money options with short maturities, our method could be a good alternative in calibrating WMSV model parameters.

Before we end this section, it is worth mentioning the parametric restrictions imposed by the assumptions in Theorems 1 and 2. Even though they are certainly limitations, the flexibility of the multidimensional model makes them less harmful. Indeed, when we use the calibrated parameters using our modeling assumptions as a warm start, the resulting calibration results without any parametric constraints in our theorems deviate only a little.