0.1 | 69.0697 | 69.594 | -0.52425
0.2 | 64.6821 | 62.9238 | 1.7583
0.3 | 62.3638 | 61.0968 | 1.267
0.4 | 61.121 | 60.7674 | 0.35363
0.5 | 60.4721 | 60.6394 | -0.16725
0.6 | 60.1788 | 60.5942 | -0.41539
0.7 | 60.1131 | 60.6083 | -0.49521
0.8 | 60.2009 | 60.6637 | -0.46274
0.9 | 60.3968 | 60.7424 | -0.34556
1 | 60.6711 | 60.8272 | -0.15604
1.1 | 61.004 | 60.9022 | 0.10177

Notes: Time to maturity $0.25$, $\theta=0.2$, $\kappa=4.0$, $\omega=0.1$,
$r=0.0953$, $\rho=-0.5$. Black-Scholes volatility $\sigma_{0}=\sigma=0.2$.
Panel A: Strike Price = 100 and spot volatility $\sigma=0.2$. Panel B: Stock
price = strike price = 100.

Table 11: Approximating $\Gamma$ for the SZ-model

Panel A:
---
Stock price | Fourier Transform, % | KM approx., % | Diff., pp.
80 | 0.57373 | 0.7438 | -0.17007
85 | 2.0189 | 2.0166 | 0.0022704
90 | 3.9028 | 3.4479 | 0.45492
95 | 4.4477 | 4.2612 | 0.18651
100 | 3.581 | 3.8587 | -0.27766
105 | 2.3676 | 2.6194 | -0.25175
110 | 1.3973 | 1.4584 | -0.061154
115 | 0.77065 | 0.7635 | 0.0071461
120 | 0.40811 | 0.3993 | 0.0088082

Panel B:
---
$\sigma(t=0)$ | Fourier Transform, % | KM approx., % | Diff., pp.
0.1 | 0.046477 | -0.077595 | 0.12407
0.2 | 0.03581 | 0.038587 | -0.0027766
0.3 | 0.028296 | 0.028236 | 6.0567e-05
0.4 | 0.023168 | 0.023232 | -6.4091e-05
0.5 | 0.019529 | 0.019976 | -0.00044687
0.6 | 0.016836 | 0.017531 | -0.00069402
0.7 | 0.014772 | 0.015584 | -0.00081176
0.8 | 0.013142 | 0.013988 | -0.00084634
0.9 | 0.011823 | 0.012652 | -0.0008288
1 | 0.010735 | 0.011519 | -0.00078471
1.1 | 0.0098209 | 0.010544 | -0.00072277

Notes: Time to maturity $0.25$, $\theta=0.2$, $\kappa=4.0$, $\omega=0.1$,
$r=0.0953$, $\rho=-0.5$. Black-Scholes volatility $\sigma_{0}=\sigma=0.2$.
Panel A: Strike Price = 100 and spot volatility $\sigma=0.2$. Panel B: Stock
price = strike price = 100.

Table 12: Approximating SZ’s $\mathscr{V}$

Panel A:
---
Stock price | Analytic $\mathscr{V}$ | KM approx. | Diff., abs.
80 | 0.86888 | 0.49785 | 0.37103
85 | 3.5542 | 3.5823 | -0.028132
90 | 8.2878 | 9.6389 | -1.3511
95 | 11.7159 | 12.9976 | -1.2817
100 | 11.6291 | 11.2967 | 0.33233
105 | 9.2815 | 8.3014 | 0.98002
110 | 6.4823 | 6.2093 | 0.27295
115 | 4.163 | 4.3411 | -0.1781
120 | 2.5337 | 2.4634 | 0.070352

Panel B:
---
$\sigma(t=0)$ | Fourier Transform, % | KM approx., % | Diff., abs.
0.1 | 10.438 | 11.2173 | -0.77928
0.2 | 11.6291 | 11.2967 | 0.33233
0.3 | 12.1307 | 87.6963 | -75.5656
0.4 | 12.3724 | -1048.4084 | 1060.7808
0.5 | 12.4988 | -35868.8908 | 35881.3896
0.6 | 12.5661 | -372544.5672 | 372557.1333
0.7 | 12.5994 | -2388099.2386 | 2388111.8381
0.8 | 12.6113 | -11380315.5007 | 11380328.112
0.9 | 12.6084 | -44060486.9252 | 44060499.5336
1 | 12.5946 | -145970936.3363 | 145970948.931
1.1 | 12.5724 | -428003736.4368 | 428003749.0092

Notes: Time to maturity $0.25$, $\theta=0.2$, $\kappa=4.0$, $\omega=0.1$,
$r=0.0953$, $\rho=-0.5$. Black-Scholes volatility $\sigma_{0}=\sigma=0.2$.
Panel A: Strike Price = 100 and spot volatility $\sigma=0.2$. Panel B: Stock
price = strike price = 100.

##  6 Commodity futures with stochastic volatility

All of the models considered so far had in common that the underlying itself
was assumed to follow a geometric Brownian motion. Only stochastic volatility
was assumed to be mean-reverting. While this is a fairly standard assumption
in the modeling of financial derivatives like options or futures e.g. on
stocks, the prices of commodities are usually modeled by mean-reverting
processes. Such mean-reversion of the underlying is well documented for a wide
range of commodities. Lutz [2009, Chapter 2] provides a good discussion on the
empirical evidence and theoretical justifications of the mean-reversion
property of commodity prices.
In the following, I consider a model due to Lutz [2009] in which both, the
price of the underlying as well as the stochastic variance, perform mean-
reversion.

$\displaystyle dX(t)=$
$\displaystyle\left(\eta\left(\alpha-X(t)\right)-\frac{1}{2}v(t)\right)+\sqrt{v(t)}dW_{1}(t)$
(57) $\displaystyle dv(t)=$
$\displaystyle\kappa\left(\theta-v(t)\right)+\omega\sqrt{v(t)}dW_{2}(t)$ (58)
$\displaystyle dW_{1}(t)dW_{2}(t)=$ $\displaystyle\rho dt$ $\displaystyle
S(t)=$ $\displaystyle e^{X(t)}$

Stochastic variance is modeled by a Heston-liked square-root process, which in
the previous experiments had been the most compatible with KM’s approximation.
The commodity price itself follows an Ornstein-Uhlenbeck process, where the
above representation results from setting $X(t)=ln(S(t))$ and then applying
Ito’s lemma. This combination of processes is especially interesting since
KM’s approximation yielded quite good results for the square-root variance
process but was unstable for Ornstein-Uhlenbeck volatility. The model above
combines both processes. While in the previous sections I always considered
plain vanilla call options, I will instead consider Futures contracts in this
section. In the absence of market frictions futures and forward prices are
equal if the interest rate is deterministic.353535See Cox et al. [1981,
Proposition 3], p. 325. As in the previous section I will assume a constant
interest rate, and hence future and forward prices will be considered the same
in the following.

###  6.1 A baseline model for commodity futures

Since I will focus on Futures in this section, the Black-Scholes model can not
serve as the baseline model here. However, there are several choices of
possible models of commodity futures prices available in the literature. A
suitable baseline to approximate the above futures price would be Model 1 from
Schwartz [1997], which I will refer to as Schwartz model in the following. The
model assumes that the log of the commodity price follows an Ornstein-
Uhlenbeck process with constant volatility

$\displaystyle dx(t)$
$\displaystyle=\kappa\left(\alpha-x(t)\right)dt+\sigma_{0}dW(t)$ (59)
$\displaystyle\alpha$ $\displaystyle=\mu-\frac{\sigma_{0}^{2}}{2\kappa}$
$\displaystyle S(t)$ $\displaystyle=e^{x(t)}$

Note that here it is assumed that the $x(t)=ln(S(t))$ follows an Ornstein-
Uhlenbeck process, not $S(t)$ itself. Hence, the processes for the underlying
differs slightly between the true and the baseline model. It is well known
that a variable that follows an Ornstein-Uhlenbeck process has a Gaussian
distribution. Hence, by using standard results from statistics it is easily
verified that $x(t)$ has the following distribution363636See Schwartz [1997],
p. 926.

$\displaystyle x\sim N\left(\underbrace{e^{-\kappa T}x+\left(1-e^{-\kappa
T}\right)\alpha}_{=\mathbb{E}(x)},\underbrace{\frac{\sigma_{0}^{2}}{2\kappa}\left(1-e^{-2\kappa
T}\right)}_{Var(x)}\right)$ (60)

Since the price of a futures contract with maturity $T$ equals the expectation
of the price of the underlying at that maturity date, the price is given
by373737See for this model Schwartz [1997], p. 927, or Cox et al. [1981], p.
339 in general.

$\displaystyle F(S,T)=\mathbb{E}(S(T))=exp\left[e^{-\kappa
T}x+\left(1-e^{-\kappa
T}\right)\alpha+\frac{1}{2}\frac{\sigma_{0}^{2}}{2\kappa}\left(1-e^{-2\kappa
T}\right)\right]$ (61)

###  6.2 KM expansion

KM’s series expansion can again be computed according the principles outlined
in section 2.2. It is straightforward to derive drift vector and the
covariance matrix of the commodity model in (57) and (58)

$\displaystyle\mu(X,t)$
$\displaystyle=\left[\eta\left(\alpha-X(t)\right)-\frac{1}{2}v(t)~{}~{}\kappa\left(\theta-v(t)\right)\right]^{\prime}\text{
,}$ (62) $\displaystyle\sigma^{2}(X,t)$
$\displaystyle=\begin{bmatrix}v(t)&\rho\omega v(t)\\\ \rho\omega
v(t)&\omega^{2}v(t)\end{bmatrix}$ (63)

Whereas the model’s PDE is

$\displaystyle\dfrac{\partial F}{\partial
t}+\left(\eta\left(\alpha-X(t)\right)-\frac{1}{2}v(t)\right)\dfrac{\partial
F}{\partial X}+\kappa\left(\theta-v(t)\right)\dfrac{\partial F}{\partial v}+$
(64) $\displaystyle+\frac{1}{2}\left(v(t)\dfrac{\partial^{2}F}{\partial
X^{2}}+\omega^{2}v(t)\dfrac{\partial^{2}F}{\partial v^{2}}\right)+\rho\omega
v(t)\dfrac{\partial^{2}F}{\partial X\partial v}=0$

Or, in terms of the operator $\mathscr{L}$ given in (5), simply

$\displaystyle\mathscr{L}F(X,t)=0$ (65)

With the boundary condition $F(X,T)=X(T)$. Note that the PDE in (65) lacks the
$rF(X,t)$ term usually appearing in the PDEs of other derivatives.
It is well established in the financial literature that the price of every
derivative should be determined by the fundamental PDE (see also (4))

$\displaystyle\mathscr{L}F(X,t)+c(X,t)=rF(X,t)$ (66)

Where $c(X,t)$ denotes a cash flow from the underlying asset. In all of the
previous sections I assumed $c(X,t)\equiv 0$, as this was also done by KM.
However, doing so in case of Futures might be misleading. Cox et al. [1981,
Propsoition 7] states that the price of a Futures contract equals the value of
an asset that receives continual payments of $c(X,t)=rX(t)$ and additionally
$F(X,T)=X(T)$.383838See Cox et al. [1981], p. 338. Consequently (66) reduces
to (65). Friedman [1975, Theorem 6.1] shows that a direct solution to (65)
would be given by $F(X,t)=\mathbb{E}_{t}\left[X(T)\right]$. By applying
Fourier transform techniques Lutz [2009] provides an analytic (semi-)closed-
form solution to this expectation. However, KM’s closed-form approximation
approach instead is based on the Feynman-Kac representation of a PDE
describing the difference between the true and the baseline model. By denoting
with $\Delta F(X,t)=F(X,t)-F_{0}(X,t)$ the difference between the Futures
price under the true and the Schwartz model, the relevant PDE can be written
as

$\displaystyle\mathscr{L}\Delta F(X,t)+\delta(X,t)$ $\displaystyle=0$ (67)
$\displaystyle\text{s.t. }\quad\quad\quad\quad\quad\quad d(X)$
$\displaystyle=0$

Where $\delta(X,t)$ is the same as in (17), and thus denoting the difference
in the driving forces between the true and the baseline market. The boundary
condition $d(X)$ is equal to zero, since both models, the true and the
baseline, have the same final payoff. In order to find a solution to (67)
similar to (16) recall the outline of KM’s approach in section 2.2, there it
was mentioned that this representation follows from Karatzas and Shreve [1991,
Theorem 7.6]. Further note that (67) differs from the general pricing bias PDE
in (12) only by the term $R(X,t)\Delta F(X,t)$. Such that effectively in (67)
the instantaneous short-term rate is $R(X,t)=0$. Hence, setting $R(X,t)=0$ and
$d(X)=0$ in (15) yields

$\displaystyle
F(X,t)=F_{0}(X,t)+\int_{t}^{T}\mathbb{E}_{t,X}\left[\delta(X,s)\right]dt$ (68)

The case shown here is analogous to the case described in Karatzas and Shreve
[1991, p. 397]. Therein $R(X,t)=\delta(X,t)=0$ whereas $d(X)\neq 0$. Then the
solution to (67) could be found again through Friedman [1975, Theorem 6.1] to
be $\Delta F(X,t)=\mathbb{E}_{t}\left[d(X)\right]$. Analogous to that case
(68) is the solution to (67). Based on this representation of the solution
KM’s series expansion is given as (see section 2.2)

$\displaystyle
F(X,t)=F_{0}(X,t)+\sum_{n=0}^{N}\dfrac{\delta_{n}(z,t)\left(T-t\right)^{n+1}}{(n+1)!}$
(69)

Since $R(X,t)$ was set to zero to derive (68) the corrective terms develop
according to

$\displaystyle\delta_{n}(X,t)=\mathscr{L}\delta_{n-1}(X,t)\text{ for }n>0$
(70)

Whereby the initial pricing error $\delta_{0}$ is still given by (14). Note
that the expectation in (68) alternatively could have been evaluated by Monte
Carlo integration. However, in such a case the obtained solution would not
have a closed-form and hence there would be no advantage compared to Monte
Carlo integration to evaluate the direct solution mentioned before.
Using the covariance matrix in (63) together with the definition in (14) it is
straightforward to derive the initial pricing error as

$\displaystyle\delta_{0}=\frac{1}{2}\left(v(t)-\sigma_{0}^{2}\right)\dfrac{\partial^{2}F_{0}}{\partial
X^{2}}-\frac{1}{2}v(t)\dfrac{\partial F_{0}}{\partial X}$ (71)

While the first element in the initial pricing error is the same convexity
adjustment as in the previous experiments, there is second element which
accounts for the different drift terms of the SDE of the underlying in the
baseline and the true model. Table 13 below shows how the series of corrective
terms developed according to the rule in (70).

Table 13: Iterations of the pricing error for the Commodity Futures model

n |  Pricing error $\delta_{n}(S,t)$
---|---
0 |  $\frac{1}{2}\left(v(t)-\sigma_{0}^{2}\right)\dfrac{\partial^{2}F_{0}}{\partial X^{2}}-\frac{1}{2}v(t)\dfrac{\partial F_{0}}{\partial X}$
1 |  $\dfrac{\partial\delta_{0}}{\partial t}+\left(\eta\left(\alpha-X\right)-\frac{1}{2}v(t)\right)\dfrac{\partial\delta_{0}}{\partial X}+\kappa\left(\theta-v(t)\right)\dfrac{\partial\delta_{0}}{\partial v}+\frac{1}{2}v(t)\dfrac{\partial^{2}\delta_{0}}{\partial X^{2}}+\rho\omega v(t)\dfrac{\partial^{2}\delta_{0}}{\partial X\partial v}$
2 |  $\dfrac{\partial\delta_{1}}{\partial t}+\left(\eta\left(\alpha-X\right)-\frac{1}{2}v(t)\right)\dfrac{\partial\delta_{1}}{\partial X}+\kappa\left(\theta-v(t)\right)\dfrac{\partial\delta_{1}}{\partial v}+\frac{1}{2}\left(v(t)\dfrac{\partial^{2}\delta_{1}}{\partial X^{2}}+\omega^{2}v(t)\dfrac{\partial^{2}\delta_{1}}{\partial v^{2}}\right)+\rho\omega v(t)\dfrac{\partial^{2}\delta_{1}}{\partial X\partial v}$
$\vdots$ |  $\vdots$
N |  $\dfrac{\partial\delta_{N-1}}{\partial t}+\left(\eta\left(\alpha-X\right)-\frac{1}{2}v(t)\right)\dfrac{\partial\delta_{N-1}}{\partial X}+\kappa\left(\theta-v(t)\right)\dfrac{\partial\delta_{N-1}}{\partial v}+\frac{1}{2}\left(v(t)\dfrac{\partial^{2}\delta_{N-1}}{\partial X^{2}}+\omega^{2}v(t)\dfrac{\partial^{2}\delta_{N-1}}{\partial v^{2}}\right)+\rho\omega v(t)\dfrac{\partial^{2}\delta_{N-1}}{\partial X\partial v}$

Notes: $\partial^{2}\delta_{0}/\partial v^{2}=0$

###  6.3 Numerical accuracy

To assess the accuracy of KM’s approximation for this model I use the
parameter values suggested by Lutz [2009]: $S=80$, $\bar{S}=85$ (i.e.
$\alpha=log(85)$), $\eta=1.0$, $\omega=0.2$, $v(t)=0.04$, $\kappa=1.0$,
$\theta=0.05$, and $\rho=-0.5$ together with times to maturity of 12, 6, and 3
months. Following the same argumentation as in the previous sections the
nuisance parameter of the baseline model is set to $\sigma_{0}=\sqrt{v(t)}$.
Figure 16 visualizes the MC result as well as its 95% confidence interval
together with the results of the KM approximation for $N=0$ to $N=4$
corrective terms for these parameter values. The simulated value of
$\mathbb{E}_{t}\left[e^{X(T)}\right]$ is 81.8091 with 95% confidence interval
of $\left[81.7681~{};~{}81.8500\right]$.393939For the same set of parameters
Lutz [2009] reports a MC result of 81.8016 and a 95% confidence interval of
$\left[81.7941~{};~{}81.8090\right]$. Lutz [2009] also reports an analytic
price for theses parameters of 81.8008. This means that my own MC results
overestimates the analytic price by 0.0101%. Hence, I believe that the use of
my MC procedure won’t have any significant influence on the drawn conclusions.
The results of KM’s approximation overestimate the MC for each of the
considered orders of the approximation. Also non of the KM approximations lies
inside the estimated confidence interval. Nevertheless, the approximation
errors are rather small, ranging between 0.2387% and 0.4131% depending on the
order of the approximation. However, the series expansion fails to converge to
the reference value after five corrective terms. An approximation using only
the initial corrective term would lead to the most accurate result
(Approximation error: 0.2387%). The precision is then reduced to it’s lowest
level by including one more element in KM’s expansion (Approximation error:
0.4131%). However, with including more corrective terms there seems be slight
improvements. Overall, the behavior of the series expansion for the first five
elements suggests that there is a significantly larger number of corrective
terms required to achieve clear convergence. While the approximation errors
are such low that the accuracy seems to be sufficient for most practical
applications, this hints at a problematic aspect of KM’s approach in general:
It is not always clear which number of corrective terms yields an optimal
approximation.
Figure 17 shows the result of the approximation when reducing the time
maturity to one month. The overall pattern is the same as before. An
approximation only using the initial corrective term would yield the most
precise approximation. Adding only one additional corrective term yields the
least precise result. Adding more corrective terms then improves the
approximation slightly again, whereas the results for $N=2$ to $N=4$ are
identical. For a maturity of one month the results of KM’s approximation lie
within the 95% confidence interval for each considered order of the
approximation. Note that even though the use of all corrective terms would not
be optimal the accuracy is nevertheless extremely high, with an approximation
error of just 0.0136% for the $N=4$ approximation. However, using only the
initial corrective term yields an approximation error as small as 0.0057%.
Figure 18 depicts the result of KM’s approximation for a time to maturity of
one year. The series expansion behaves as expected in that case. The accuracy
is lower than both previous examples. Again the first order approximation
yields the most and the second order approximation yields the least precise
result. Non of the approximations falls inside the 95% confidence interval of
the MC simulation. However, also with a maturity of one year all approximation
errors remain remarkably low. For the approximation of order $N=4$ the error
is just 0.9521% and for the expansion of order $N=0$ the error is 0.6856%.

Figure 16: MC result vs. KM approximation - maturity 0.5 years

Notes: KM’s approximation compared to MC result for the stochastic variance
commodity futures model. A discription of the MC approach is provided in
Appendix A.5. Parameter values: $S=80$, $\bar{S}=85$ (i.e. $\alpha=log(85)$),
$\eta=1.0$, $T-t=0.5$, $\omega=0.2$, $v(t)=0.04$, $\kappa=1.0$, $\theta=0.05$,
and $\rho=-0.5$. Baseline model is Model I from Schwartz [1997].

Figure 17: MC result vs. KM approximation - maturity 0.25 years

Notes: KM’s approximation compared to MC result for the stochastic variance
commodity futures model. A discription of the MC approach is provided in
Appendix A.5. Parameter values: $S=80$, $\bar{S}=85$ (i.e. $\alpha=log(85)$),
$\eta=1.0$, $T-t=0.25$, $\omega=0.2$, $v(t)=0.04$, $\kappa=1.0$,
$\theta=0.05$, and $\rho=-0.5$. Baseline model is Model I from Schwartz
[1997].

Figure 18: MC result vs. KM approximation - maturity 1.0 years

Notes: KM’s approximation compared to MC result for the stochastic variance
commodity futures model. A discription of the MC approach is provided in
Appendix A.5. Parameter values: $S=80$, $\bar{S}=85$ (i.e. $\alpha=log(85)$),
$\eta=1.0$, $T-t=1.0$, $\omega=0.2$, $v(t)=0.04$, $\kappa=1.0$, $\theta=0.05$,
and $\rho=-0.5$. Baseline model is Model I from Schwartz [1997].

##  7 Conclusion

The main theme of this thesis was the analysis of the behavior of the closed-
form approximation approach developed by Kristensen and Mele [2011]. In terms
of derivatives pricing models I focused on stochastic volatility models as
they yield a convenient set-up to test the accuracy and convergence behavior
of the approximations. I applied KM’s approach to obtain closed-form
approximation of derivative prices as well as Greeks.
KM’s approximation yields remarkably precise approximations for the Heston and
the CEV model over a wide range of parameter choices. The series expansion
converges fast, such that five corrective terms are sufficient to approximate
asset prices in these models. Since the approximation has closed-form it is
straightforward to obtain greeks, where the performance in the approximation
of greeks in the Heston is equally well as the approximation of option prices.
Especially in the case of the CEV model KM’s approximation reveals it strength
in terms of computational efficiency. While computation of reference values
via Monte Carlo simulation requires substantial computation time, KM’s closed-
form approximation yield results virtually immediately.
Both models, CEV and Heston, were approximated for two different sets of
parameters. For the Heston model the approximation mostly performs better for
ITM than for OTM options. Whereby this difference in accuracy becomes more
pronounced the fewer corrective terms in the series expansion are considered.
Additionally the effect of moneyness on the accuracy of the approximation
decreases with the level correlation and with time to maturity, whereas the
latter had a more significant impact on put than on call options. The
approximation of the CEV model appears to be less robust against changes in
the model parameter, especially for high levels of elasticity of variance. An
exception seems to be the sensitivity to correlation, which is almost offset
by a high elasticity of variance.
The approximation of the SZ model performed significantly worse. Specifically,
time to maturity influences the accuracy of the approximation much stronger
than it was the case for the Heston model. While in the Heston model the
approximation remains stable for maturities of up to two years, in the SZ
model instability already occurs for maturities longer than six month. Since
the only difference between the two models is the variance/volatility process,
the instability is most likely related to the usage of the Ornstein-Uhlenbeck
process in the SZ model. The Heston model is affine in the sense of
exponentially linear characteristic functions, while the SZ model as well as
the CEV model are non-affine in this sense. The analyzes of the approximations
for the CEV model indicate that accuracy of KM’s approximation is lower for
non-affine models as it is for affine models. Nevertheless, the convergence
behavior of the series expansion is similar in the CEV and the Heston case,
suggesting the instability might be not related to the non-affine structure of
the SZ model. These results suggest that the convergence of KM’s series
expansion might be not uniform over time. A possible explanation of this could
be found in Kimmel [2008], who links diverging behavior of series expansion
solutions to PDEs to singularities for some maturities. Kimmel [2008] suggests
to extend the area of convergence of such series expansion by replacing the
time variable by a non-linear function of time. KM mention that the approach
developed by Kimmel [2008] could be used to improve convergence of their
approximations, but do not describe how this can be done. Kimmel [2008]
applies his time-change approach only to univariate diffusions and indicates
that an extension to multivariate diffusion might be not straightforward in
many cases.404040See Kimmel [2008], p. 38. Hence, the combination of the
approaches of Kimmel [2008] and KM might be a promising possibility for future
research. If only maturities below six months are considered, the
approximation of the SZ model is overall more precise if correlation is
zero.414141As already noted, in this case the SZ model is identical to the S&S
model. However, while in the zero correlation case an approximation using five
corrective terms yields the most precise results, an approximation using six
corrective terms is superior if correlation is non-zero. Similar to the Heston
and the CEV model cases the approximation performs better for ITM than for OTM
options, whereby again this effect increases with the level of correlation.
The accuracy in the approximation of greeks in the SZ model mimics the
accuracy in the approximation of option prices. However, the approximation
seems to be unstable with respect to spot volatility.
Finally, I applied KM’s approximation to a model of commodity future prices
which models the price of the underlying through an Ornstein-Uhlenbeck process
and stochastic variance through an Heston-like square-root process. The series
expansion was developed for up to five corrective terms, generally showing
neither clear convergence nor clear divergence for maturities of up to one
year. It seems like the approximation is converging to a slightly biased
futures price instead to the true price. However, approximation errors
remained below one percent. Most noteworthy, for the approximation of the
commodity futures model the approximation using only the first corrective
terms always yields the highest accuracy. Considering that also for the SZ
model with zero correlation the approximation using six terms is less precise
than the approximation using five terms and that for the CEV and the Heston
model in some cases also the approximation using four terms was slightly more
accurate than the one using five terms, this indicates that it might be
difficult to determine the optimal order of approximation without solving the
model under investigation also through another procedure to obtain reference
values.
Nevertheless, KM’s approach yields a viable alternative to Monte Carlo
solutions or Runge-Kutta solutions for the futures model due to significantly
higher computational efficiency. The closed-form approximation could also be
considered as an alternative to the analytic solution of the model, which is
only available in terms of Kummer functions of the first and second
kind.424242See Lutz [2009], p. 64. While these functions might be difficult to
implement in some programming languages, KM’s series expansion is
straightforward to implement in most languages.
Overall, KM’s series expansion yields satisfying results for all models
considered in this thesis except for the SZ model. Due to the very abstract
nature of the underlying assumptions it appears to be difficult to point at
specific reason for differences in the convergence behavior or sensitivities
of the approximation towards specific parameter values. Even though the
implementation of the series expansion is straightforward, the derivation is
tedious and prone to error if it is done manually. However, even if symbolic
mathematical software such as Maple is used, the computation is restricted to
only five terms, since the corrective terms are increasing in length too fast
exceeding the computers memory. A major issue in this matter might be the
structure of the chosen baseline model. For most of the examples in this
thesis I used the Black-Scholes model as baseline. As mentioned in one of the
previous sections if only five corrective terms are used, this already
requires the computation of derivatives of the Black-Scholes model up to an
order of ten, which effectively means taking derivatives and cross-derivatives
of a normal distribution up this order. Although if symbolic mathematical
software is used, the series expansion usually can only be developed for a
view terms. Approaches like of finite difference approximation for derivatives
of the baseline model, as it as done by Garcia [2013], ease mathematical
difficulties in developing the approximation but significantly slow down the
computation speed. Additionally, the method then loses its feature of yielding
pricing formulas in closed-form and thus the approach would lose its two main
advantages. Hence, I would suggest three areas of future research regarding
KM’s approach. (i) The combination of Kimmel [2008]’s time-change approach
with KM’s series expansion to improve robustness regarding the time to
maturity. (ii) A closer investigation of the iterative procedure of computing
the corrective terms to derive shorter expressions for the pricing bias and
thus enabling the computation of more terms in the series expansion. (iii) The
development of a tractable test procedure of the underlying assumptions for
specific models. This could be connected with an investigation of criteria to
determine an optimal order of the approximation for specific models. (iv)
Finally, a closer investigation of the interaction between true model and
baseline model to determine optimal choices of the baseline.

##  Appendix A Reference solutions

This appendix provides descriptions on the Fourier transform solutions of the
Heston [1993] model and the Schöbel and Zhu [1999] model as well as on the
Monte Carlo solutions used for the CEV model and the commodity futures model.

###  A.1 Fourier transform solution of the Heston model

The analytic prices of the Heston model are obtained via Fourier transforms.
The Matlab code which I use to compute these prices is provided by Crisostomo
[2014], who uses a formulation of the models characteristic functions that had
been suggested by Gatheral [2006] instead of Heston’s original formulation.
One particular drawback of Heston’s original formulation of the characteristic
functions is the possibility of branch-cuts when evaluating the complex valued
logarithm appearing in the characteristic functions. These branch-cuts might
be problematic in the present context as I will compare the performance of
KM’s approximation under a variety of different parameter values. Hence, if
issue of possible branch-cuts would be ignored, the validity of the analytic
reference values could not be ensured for all parameter choices and thus
potentially spoil the conclusions drawn from the analysis of the approximation
errors. However, by using the characteristic function from Gatheral [2006] the
argument of the logarithm in the characteristic function never crosses the
negative real axis, such that branch cuts do not occur.434343See Gatheral
[2006], p. 20. Nevertheless, Gartherals formulation of the model is equivalent
to Heston’s original formulation.444444See e.g. Zhu [2010], p. 55. Crisostomo
[2014] uses the following pricing framework for the Heston model

$\displaystyle C_{Analytic}(S,t)$ $\displaystyle=S_{0}P_{1}-e^{-r(T-t)}KP_{2}$
(72) with $\displaystyle P_{1}$
$\displaystyle=\frac{1}{2}+\frac{1}{\pi}\int_{\infty}^{0}real\left(\frac{e^{-i\phi
log(K)}\varphi(\phi-i)}{i\phi\varphi(-i)}\right)d\phi$ $\displaystyle P_{2}$
$\displaystyle=\frac{1}{2}+\frac{1}{\pi}\int_{\infty}^{0}real\left(\frac{e^{-i\phi
log(K)}\varphi(\phi)}{i\phi}\right)d\phi$

Where the characteristic function $\varphi(\cdot)$ is given in Crisostomo
[2014]. Note that this approach slightly deviates from the one in Gatheral
[2006]. Both Gatheral [2006] and Heston [1993] derive two distinct
characteristic functions for the probabilities $P_{1}$ and $P_{2}$, while
Crisostomo [2014] uses only one characteristic function for both
probabilities. Also Crisostomo [2014] applies Gatherals method of deriving the
characteristic function to $log(S(t))$, instead of $log(S(t)/K)$ as in
Gatheral [2006]. However, Crisostomo [2014, Appendix B] provides a prove for
the equivalence of the two approaches.
In order to evaluate the real integrals in (72) direct integration via Matlabs
integral() function is used, which employs global adaptive quadrature to
numerically evaluate integrals.454545See Matlab help file for more details.
Zhu [2010, Ch. 4.3 to 4.5] provide a good overview on integration algorithms
in the application of Fourier transform methods.

It is straightforward to obtain a (semi-)closed-form solution via Fourier
transforms for European puts in the Heston framework464646See e.g. Zhu [2000],
p. 35.

$\displaystyle Put(S,t)$ $\displaystyle=e^{-r(T-t)}KP_{2}-S_{0}P_{1}$ (73)
with $\displaystyle P_{1}$
$\displaystyle=\frac{1}{2}-\frac{1}{\pi}\int_{\infty}^{0}real\left(\frac{e^{-i\phi
log(K)}\varphi(\phi-i)}{i\phi\varphi(-i)}\right)d\phi$ $\displaystyle P_{2}$
$\displaystyle=\frac{1}{2}-\frac{1}{\pi}\int_{\infty}^{0}real\left(\frac{e^{-i\phi
log(K)}\varphi(\phi)}{i\phi}\right)d\phi$

Where the characteristic function $\varphi$ is the same in (72). I implemented
the Heston Puts by adapting the Matlab code from Crisostomo [2014].

####  A.1.1 Analytic Greeks in the Heston model

The computation of these so called Greeks for the Heston model is
straightforward, since integration and differentiation of the characteristic
function can be interchanged. Note that I again use the formulation of the
characteristic function developed by Crisostomo [2014]. Recalling the general
structure of call option prices $C(S,t)=SP_{1}-e^{-r(T-t)}KP_{2}$ and the
definitions 1.) to 4.) in section 2.3 the analytic greeks of the Heston model
are defined as 474747See e.g. Zhu [2000], p. 36.

$\displaystyle\Delta_{S}~{}=$ $\displaystyle~{}P_{1}$ (74)
$\displaystyle\Gamma_{S}~{}=$ $\displaystyle~{}\dfrac{\partial P_{1}}{\partial
S}$ (75) $\displaystyle\mathscr{V}~{}=$ $\displaystyle~{}S\dfrac{\partial
P_{1}}{\partial v}-e^{-r(T-t)}K\dfrac{\partial P_{2}}{\partial v}$ (76)
$\displaystyle\text{with }\dfrac{\partial P_{j}}{\partial h}~{}=$
$\displaystyle~{}\frac{1}{\pi}\int_{0}^{\infty}\dfrac{\partial\varPsi_{j}(\phi)}{\partial
h}d\phi$ $\displaystyle\varPsi_{1}~{}$
$\displaystyle=~{}real\left(\frac{e^{-i\phi
log(K)}\varphi(\phi-i)}{i\phi\varphi(-i)}\right)$
$\displaystyle\varPsi_{2}~{}$ $\displaystyle=~{}real\left(\frac{e^{-i\phi
log(K)}\varphi(\phi)}{i\phi}\right)$

Where $h$ equals either $S$ or $v$, and $j=1,2$.

###  A.2 Fourier transform solution of the SZ model

The analytic prices for both, the S&S- and the SZ-model had been obtained via
Fourier transforms, whereas the S&S prices can be computed by simply setting
$\rho=0$ in the characteristic function. The solution has the same structure
as the Fourier solution of the Heston model in (72). However, Schöbel and Zhu
[1999] use two distinct characteristic functions to derive the probabilities
$P_{1}$ and $P_{2}$. In order to implement the Fourier solution I adapted
Matlab code that I used during the Numerical Methods in Finance PC-Lab classes
at University of Tuebingen. The Matlab program simply implements the two
characteristic functions as given in Schöbel and Zhu [1999] and uses global
adaptive quadrature484848Numerical integration via Matlab’s integral()
function. The original PC-lab code used a simple trapezoidal algorithm via
Matlab’s trapz() function. However, while yielding numerically identical
results, I found integral() to be significantly faster. to numerically
evaluate the involved integrals. Since the characteristic functions of the SZ-
model also include a complex valued logarithm the same issue of possible
branch-cuts arises as in the Heston model. While there I could conveniently
use a reformulation of the characteristic functions to ensure stability of the
reference solution, such a reformulation is, to my knowledge, not readily
available for the SZ-model. Hence, I use a simple correction algorithm to
adjust the complex logarithm whenever it’s argument crosses the negative real
axis. Appendix A.3 provides a description of the used algorithm as well as the
corresponding Matlab code.

###  A.3 Digression on complex logarithms

This appendix is an adapted excerpt from two assignments of the 2015
’Numerical Methods in Finance’-class, which I wrote together with F. Slezak
and A. Berg.
The issue of branch-cuts arose in the context of the Fourier transform
solution to the Heston model in section 3 as well as the SZ model in section
5, since in both cases the computation of the option price includes the
evaluation of a complex valued logarithm. It can be shown that the logarithm
of a complex number $z$ has the following form

$\displaystyle w$ $\displaystyle=log(z)=log(|z|)+i(Arg(z)+2k\pi)$ (77)
$\displaystyle\textit{with }k=\pm 0,\pm 1,\pm 2,...$

Where $Arg(z)$ denotes the principal argument of $z$, i.e. $-\pi\leq
Arg(z)<\pi$. Most commercial software packages like Matlab restrict
calculation of complex logarithms to their principal values, i.e. $k$ is set
to zero in equation (77). As any value of $k$ would suffice to recover
$e^{w}=z$, the choice of $k=0$ would not be a problem in an isolated
computation.494949See Zhu [2010], p. 100. However, Schöbel and Zhu [1999] were
among the first to observe that this leads, at least for extreme parameter
values, to discontinuities in the characteristic functions of their and other
stochastic volatility models and thus to a wrong integration.505050See Schöbel
and Zhu [1999], p. 28 and Lord and Kahl [2010], p. 672. The main issue in
implementing equation (77) in order to correctly compute the complex logarithm
is that the appropriate $k$ value can only be found if the path of the complex
number is known.515151See Zhu [2010], p. 100. Every time $z$ crosses the
negative part of the real axis in Figure 19 the value of $k$ needs to be
adjusted.

$Re(z)$$Im(z)$$Z_{0}$$\varphi$IIIIIIIV Figure 19: $z_{0}$ path around the
origin

While Gatheral [2006] provided a reformulation of Heston’s characteristic
function such that the argument of the complex logarithm never crosses the red
line in 19, I am unaware of any such reformulation for the SZ model. In order
to obtain a numerically stable reference solution to the SZ model I used a
correction algorithm which adjusts the number $k$ in (77) every time the
complex argument of the logarithm crosses the red line in Figure 19.525252The
used Matlab code was taken from the same Assignment as the content of this
appendix section and also slightly adapted.

###  A.4 MC simulation for the CEV model

MC requires the discretization of the processes in (48) and (49), which may
lead to the occurrence of negative variances in the sample path. Such values
are problematic since in the process of generating a stock price sample path
the square-root as well other roots of the variance need to be computed. This
issue had been studied extensively in the financial literature, mostly in the
context of the Heston model. Introducing the absolute value to the variance
process as in (49) resembles the so called reflective assumption which sets
the rule: if $v<0$, then $-v=v$. However, to implement this it is necessary to
apply this rule to all roots of $v$. Another possibility to avoid negative
variances would be to use the so called absorbing assumption which sets the
rule: if $v<0$, then $v=0$.535353See e.g. Gatheral [2006], p. 21. If a
sufficiently large number of trajectories and sample paths is chosen neither
of the two procedures is likely to have no significant influence on the
result, because the true variance is likely to be small if the discretization
error leads to a negative simulated variance.545454See Lee and Sokolinskiy
[2014], p. 8. KM do not point out which discretization scheme they use to
obtain their MC results, neither they mention the number of trajectories nor
the number of sample paths used in their simulation.
In order to discretize the processes in (48) and (49) I use the Milstein
scheme555555According to Glasserman [2003] the name Milstein scheme might be
misleading since there are several methods due to Milstein. However,
throughout this thesis I use the term to name the discretization scheme shown
below. which simply adds the next element of the Ito-Taylor expansion to the
Euler scheme565656i.e., the term
$\frac{1}{2}Diff(h)Diff(h)^{\prime}\left(Z^{2}_{h}(t_{i})-1\right)\Delta_{i}$
with $h=S~{}or~{}V$ is added. Where $Diff(h)$ denotes the diffusion term of
the respective SDE and $Diff(h)^{\prime}$ its derivative with respect to the
variable denoted by $h$. (See e.g. Glasserman [2003], p. 343)

$\displaystyle S(t_{i+1})$
$\displaystyle=S(t_{i})\left[1+r\Delta_{i}+\sqrt{|V(t_{i})|}\sqrt{\Delta_{i}}Z_{1}(t_{i})+\frac{1}{2}|V(t_{i})|\left(Z^{2}_{1}(t_{i})-1\right)\Delta_{i}\right]$
(78) $\displaystyle V(t_{i+1})$
$\displaystyle=V(t_{i})+\kappa\left(\theta-|V(t_{i})|\right)\Delta_{i}+\omega|V(t_{i})|^{\gamma}\sqrt{\Delta_{i}}Z_{2}(t_{i})$
(79)
$\displaystyle\quad\quad+\frac{1}{2}\omega^{2}\gamma|V(t_{i})|^{2\gamma-1}\left(Z^{2}_{2}(t_{i})-1\right)\Delta_{i}$
with $\displaystyle Z_{1}(t_{i})=Ran_{1},\quad Ran_{1}\sim N(0,1)$
$\displaystyle Z_{2}(t_{i})=\rho Ran_{1}+\sqrt{1-\rho^{2}}Ran_{2},\quad
Ran_{2}\sim N(0,1)$

Where $\Delta_{i}$ denotes the length of time step $i$. Note that I used the
reflective assumption to handle negative simulated variances if they occur.
Applying the Milstein discretization instead of the Euler discretization do
not increase accuracy of the MC results since both schemes are of weak order
one.575757See Glasserman [2003], p. 347. Nevertheless, using the Milstein
scheme do not significantly increase the computational burden, but reduces the
frequency of negative variances in the sample path.585858See Gatheral [2006],
p. 22. Hence, Gatheral [2006] advocates for always using the Milstein
discretization for the Heston model. Due to the similarity between the CEV and
the Heston model the positive effect of applying the Milstein scheme should
also apply in this case. This further mitigates the effect of the choice
between reflective and absorbing assumption. Additionally, I use 500 equally
sized time steps and 20,000 sample paths in all simulations below.595959The
pseudo-random numbers are obtained via Matlabs multiple recursive generator
(mrg32k3a).

###  A.5 MC simulation of the Commodity Future

It was already mentioned before that Lutz [2009] provides an analytic solution
for the model in (57) and (58). However, this solution requires the evaluation
of Kummer functions of the first and second kind.606060The type of
hypergeometric function that appears in the analytic solution depends on the
parameter values. For the special case of perfect correlation between the
underlying and the variance process the Bessel function of the first and
second kind appears instead of the Kummer function (See Lutz [2009], p. 64).
In order to implement the Kummer functions in the analytic solution Lutz
[2009] uses Matlab code provided by B. Barrows via Matlabs central file
exchange, but needs to adapt that code to handle integer arguments.616161Note
that since release R2014b of Matlab, Kummer functions of the second kind are
available as build-in functions via kummerU(a,b,z) (See Matlab online
documentation). The function is only available via the Symbolic Math Toolbox.
Alternatively the analytic solution could be implemented by using Runge-Kutta
algorithms.626262See Lutz [2009], pp. 77 - 79 for a comparison of these two
implementation approaches.
However, for simplicity I will use MC simulations instead of the analytic
solution to obtain reference values for the futures price approximated by KM’s
series expansion. Lutz [2009] also compares the analytic solution to MC
results, finding a very close match of the obtained prices. In order to
generate the MC solution I use a similar approach as in section A.4. Since the
log of the price of the underlying is simulated instead of the underlying
itself, I use a simple Euler discretization for the (57). However, by the same
arguments as in section A.4 I use a Milstein discretization for the variance
process in (58) to reduce the frequency of negative simulated variances in the
sample paths.

$\displaystyle X(t_{i+1})$
$\displaystyle=X(t_{i})+\left[\eta\left(\alpha-X(t_{i})\right)-\frac{1}{2}|V(t_{i})|\right]\Delta_{i}+\sqrt{|V(t_{i})|}\sqrt{\Delta_{i}}Z_{1}(t_{i})$
(80) $\displaystyle V(t_{i+1})$
$\displaystyle=|V(t_{i})|+\kappa\left(\theta-|V(t_{i})|\right)\Delta_{i}+\omega\sqrt{|V(t_{i})|\Delta_{i}}Z_{2}(t_{i})$
(81)
$\displaystyle\quad\quad+\frac{1}{4}\omega^{2}\left(Z^{2}_{2}(t_{i})-1\right)\Delta_{i}$
with $\displaystyle Z_{1}(t_{i})=Ran_{1},\quad Ran_{1}\sim N(0,1)$
$\displaystyle Z_{2}(t_{i})=\rho Ran_{1}+\sqrt{1-\rho^{2}}Ran_{2},\quad
Ran_{2}\sim N(0,1)$

Where $\Delta_{i}$ denotes the length of time step $i$. As the above equations
indicate I use again the reflective assumption to handle negative simulated
variances in the sample path, if they occur. For his MC simulation Lutz [2009]
uses a huge number of 2,500 time steps and 1.5 million sample paths.636363Of
which 375,000 were independent. The large number of 1.5 million paths results
from the use of antithetic sampling. Unfortunately the computational resources
that had been available to me were not sufficient to use such numbers. Hence,
I constraint myself to 1,000 time steps and 200,000 sample paths, which I
believe should be sufficient to achieve satisfying results for the reference
values. Below Figure 20(a) and 20(b) show the histogram and the density of the
simulated price of the underlying at maturity of some Future contract. The
simulation is based on the above shown discretization of the commodity model
and the following set of parameter values: $S=80$, $\bar{S}=85$ (i.e.
$\alpha=log(85)$), $\eta=1.0$, $T-t=0.5$, $\omega=0.2$, $v(t)=0.04$,
$\kappa=1.0$, $\theta=0.05$, and $\rho=-0.5$. These parameters had been
suggested by Lutz [2009], since they should yield a realistic set up.

(a) Histogram of simulated commodity prices at time $T$

(b) Density of commodity prices at time $T$

Figure 20: Simulation of commodity prices

Notes: Panel (a): 1,000 bins had been used. Panel (b): Kernel density
estimated from the data in Panel (a) by Matlab’s kdensity() function. For both
Panels: $S=80$, $\bar{S}=85$ (i.e. $\alpha=ln(85)$), $\eta=1.0$, $T-t=0.5$,
$\omega=0.2$, $v(t)=0.04$, $\kappa=1.0$, $\theta=0.05$, and $\rho=-0.5$

###  A.6 Definition of the scale measure of the CEV variance process

Generally, the scale measure $\Omega(\cdot)$ of a diffusion
$dx=\mu(x)dt+\sigma(x)dW$, with $dW$ denoting a Wiener process, is defined
as646464See Jones [2003], p. 215.

$\displaystyle\Omega(x)$ $\displaystyle=\int_{m}^{n}\Theta(x)dx$ (82)
$\displaystyle\text{with }\Theta(x)$
$\displaystyle=exp\left(-\int\frac{2\mu(x)}{\tilde{\sigma}^{2}(x)}dx\right)$
(83)

Where $\Theta(x)$ is called the scale density. Note that
$\tilde{\sigma}^{2}(x)$ denotes the variance of the process $x$ only and not
of the whole covariance matrix, if e.g. $x$ would be part of a system of SDEs.
Hence, in the case of the CEV variance process in (49) we have
$\mu(v)=\kappa\left(\theta-v\right)$ and
$\tilde{\sigma}^{2}(v)=\omega^{2}v^{2\gamma}$. Such that the scale density can
be computed as

$\displaystyle\Theta(v)$
$\displaystyle=exp\left(-\int\frac{2\mu(v)}{\tilde{\sigma}^{2}(v)}dv\right)$
$\displaystyle=exp\left(-\int\frac{2\kappa\left(\theta-v\right)}{\omega^{2}v^{2\gamma}}dv\right)$
$\displaystyle=exp\left(-\frac{2\kappa}{\omega^{2}}\left(\theta\int\frac{1}{v^{2\gamma}}dv-\int\frac{1}{v^{2\gamma-1}}dv\right)\right)$
(84)

Where the two integrals can be solved as

$\displaystyle\int\frac{1}{v^{2\gamma}}dv$
$\displaystyle=\begin{cases}-\frac{x^{1-2\gamma}}{2\gamma-1}&\quad\text{if
}\gamma\neq 1/2\\\ log(v)&\quad\text{if }\gamma=1/2\\\ \end{cases}$ (85)
$\displaystyle\int\frac{1}{v^{2\gamma-1}}dv$
$\displaystyle=\begin{cases}-\frac{x^{2-2\gamma}}{2\gamma-2}&\quad\text{if
}\gamma\neq 1\\\ log(v)&\quad\text{if }\gamma=1\\\ \end{cases}$ (86)

Plugging the first solutions from (85) and (86) into (84) and rearranging
yields the scale density of the CEV variance process as it was given in (51)

$\displaystyle\Theta(v)=exp\left(\frac{2\kappa\theta}{\omega^{2}\left(2\gamma-1\right)}\frac{1}{v^{2\gamma-1}}-\frac{\kappa}{\omega^{2}\left(\gamma-1\right)}\frac{1}{v^{2\gamma-2}}\right)$

Note that this expression is identical to the scale density for the CEV
process given in Jones [2003] with $\alpha=\kappa\theta$ and $\beta=-\kappa$.

## References

  * Ait-Sahalia and Yu [2006] Yacine Ait-Sahalia and Jialin Yu.  Saddlepoint approximations for continuous-time markov processes.  _Journal of Econometrics_ , (134):507–551, 2006\.
  * Black and Scholes [1973] Fischer Black and Myron Scholes.  The pricing of options and corporate liabilities.  _Journal of Political Economy_ , (81):637–659, 1973.
  * Bollerslev and Zhou [2002] T. Bollerslev and H. Zhou.  Estimating stochastic volatility diffusion using conditional moments of integrated volatility.  _Journal of Econometrics_ , (109):33–65, 2002\.
  * Chan et al. [1992] K.C. Chan, G.A. Karolyi, F.A. Longstaff, and A.B. Sanders.  An empirical comparison of alternative models of the short-term interest rate.  _Journal of Finance_ , (47):1209–1227, 1992.
  * Cizek et al. [2011] Pavel Cizek, Wolfgang K. Haerdle, and Rafa Weron.  _Statistical tools for finance and insurance_.  Springer, Berlin, 2011.
  * Cox et al. [1981] J.C. Cox, J.E. Ingersoll, and S.A. Ross.  The relation between forward prices and futures prices.  _Journal of Financial Economics_ , (9):321–346, 1981.
  * Crisostomo [2014] Ricardo Crisostomo.  An analysis of the heston stochastic volatility model: Implementation and calibration using matlab.  _Spanish Securities Market Commission – CNMV_ , (Working Paper No. 58), 2014.
  * Daniels [1954] H. E. Daniels.  Saddlepoint approximations in statistics.  _Annals of Mathematical Statistics_ , (25(4)):631–650, 1954.
  * Elices and Gimenez [2013] Alberto Elices and Eduard Gimenez.  Applying hedging stategies to estimate model risk and provision calculation.  _Quantitative Finance_ , (13(7)):1015–1028, 2013\.
  * Forani and Mele [2006] F. Forani and Antonnio Mele.  Approximating volatility diffusions with cev-arch models.  _Journal of Economic Dynamics and Control_ , (30):931–966, 2006.
  * Friedman [1975] Avner Friedman.  _Stochastic differential equations and applications Vol. I_.  Academic Press, New York, 1975.
  * Fries [2007] Christian Fries.  _Mathematical Finance: Theory, Modelling, Implementation_.  John Wiley & Sons, Inc., Hoboken, New Jersey, 2007.
  * Garcia [2013] Carlos A. Catalan Garcia.  Pricing equity derivatives with the kristensen-mele approach.  _Master Thesis at Universidad de Madrid, del Pais Vasco, de Valencia, and Castilla-La Mancha_ , 2013.
  * Gatheral [2006] Jim Gatheral.  _The Volatility Surface: A Practitioner’s Guide_.  John Wiley & Sons, Inc., Hoboken, New Jersey, 2006.
  * Gill et al. [1981] Philip E. Gill, Walter Murray, and Margaret H. Wright.  _Practical Optimaziation_.  Elsevir, San Diego, 1981.
  * Glasserman [2003] Paul Glasserman.  _Monte Carlo Methods in Financial Engineering_.  Springer, New York, 2003.
  * Goutis and Casella [1999] Constantino Goutis and George Casella.  Explaining the saddlepoint approximation.  _The American Statistician_ , (53(3)):216–224, 1999.
  * Hermann [2006] Martin Hermann.  _Numerische Mathematik_.  Oldenbourg Wissenschaftsverlag GmbH, München, 2nd edition, 2006.
  * Heston [1993] Steven Heston.  A closed-form solution for options with stochastic volatility with applications to bond and currency options.  _Review of Financial Studies_ , (6):327–343, 1993\.
  * Hirsa and Neftci [2014] Ali Hirsa and Salih N. Neftci.  _An Introduction to the Mathematics of Financial Derivatives_.  Elsevier Inc., London, 3rd edition, 2014.
  * Hull [2015] John C. Hull.  _Options, Futures and other Derivatives_.  Pearson, Toronto, 9th edition, 2015.
  * Hull and White [1988] John C. Hull and Alan White.  An analysis of the bias in option pricing caused by a stochastic volatility.  _Advances in futures and options research: a research annual_ , (Vol. 3):27–61, 1988.
  * Jones [2003] Christopher S. Jones.  The dynamics of stochastic volatility: evidence from underlying and options markets.  _Journal of Econometrics_ , (116):181–224, 2003\.
  * Karatzas and Shreve [1991] Ioannis Karatzas and Steven E. Shreve.  _Brownian Motion and Stochastic Calculus_.  Springer-Verlag, New York, 2nd edition, 1991.
  * Kimmel [2008] Robert L. Kimmel.  Changing times: The pricing problem in non-linear models.  _Dice Center Working Paper_ , (2008-24), 2008.
  * Kristensen and Mele [2011] Dennis Kristensen and Antonio Mele.  Adding and subtracting black-scholes: A new approach to approximate derivative prices in continuous-time models.  _Journal of Financial Economics_ , (102):390–415, 2011.
  * Lee and Sokolinskiy [2014] Cheng-Few Lee and Oleg Sokolinskiy.  R-2gam stochastic volatility model: flexibility and calibration.  _Review of Quantitative Finance and Accounting_ , (doi: 10.1007/s11156-014-0443-7), 2014.
  * Lehoczky [2012] John Lehoczky.  Simulating the greeks in finance.  Presentation slides from the conference Probability, Control, and Finance at Columbia University, NY, 2012.
  * Lewis [2000] A.L. Lewis.  _Option Valuation Under Stochastic Volatility_.  Finance Press, Newport Beach, CA, 2nd edition, 2000.
  * Lord and Kahl [2010] Rodger Lord and Christian Kahl.  Complex logarithms in heston-like models.  _Mathematical Finance_ , 20(4):671–694, 2010\.
  * Lutz [2009] Björn Lutz.  _Pricing of Derivatives on Mean-Reverting Assets_.  Doctoral Dissertation, University of Tuebingen, 2009.
  * Martins et al. [2003] Joaquim R. R. A. Martins, Peter Sturdza, and Alonso Juan J.  The complex-step derivative approximation.  _ACM Transactions on Mathematical Software_ , (29(3)):245–262, 2003.
  * Rogers and Zane [1999] L.C.G. Rogers and O. Zane.  Saddlepoint approximations to option prices.  _The Annals of Applied Probability_ , (9(2)):493–503, 1999.
  * Schöbel and Zhu [1999] Rainer Schöbel and Jianwei Zhu.  Stochastic volatility with an ornstein-uhlenbeck process: An extension.  _European Finance Review_ , (3):23–46, 1999.
  * Schwartz [1997] Eduardo S. Schwartz.  The stochastic behavior of commodity prices: Implications for valuation and hedging.  _The Journal of Finance_ , (VLII (3)):923–973, 1997.
  * Stein and Stein [1991] Elias M. Stein and Jeremy C. Stein.  Stock price distributions with stochastic volatility: An analytic approach.  _The Review of Financial Studies_ , (4(4)):727–752, 1991.
  * Yang [2006] Jin Yang.  _Stochastic Volatility Models: Option Price Approximations, Asymptotics and Maximum Likelihood Estimation_.  Doctoral Dissertation, University of Illinois at Urbana-Champaign, 2006\.
  * Younesian [2013] Houman Younesian.  Approximating option prices.  _Master Thesis at University of Rotterdam_ , 2013.
  * Zhu [2000] Jianwei Zhu.  _Modular Pricing of Options: An Application of Fourier Analysis_.  Doctoral Dissertation, University of Tuebingen, 2000.
  * Zhu [2010] Jianwei Zhu.  _Applications of Fourier Transform to Smile Modeling: Theory and Implementation_.  Springer-Verlag, Heidelberg, 2nd edition, 2010.