Closed-form approximations in derivatives pricing:
The Kristensen-Mele approach

Degree Dissertation

for the

Master Examination in Economics & Finance

at the

Faculty of Economics and Social Sciences

of the

Eberhard Karls Universität

Tübingen

Submitted by: Michael Kurz Date of submission: 14/09/2015

ABSTRACT

Kristensen and Mele (2011) developed a new approach to obtain closed-form
approximations to continuous-time derivatives pricing models. The approach
uses a power series expansion of the pricing bias between an intractable model
and some known auxiliary model. Since the resulting approximation formula has
closed-form it is straightforward to obtain approximations of greeks. In this
thesis I will introduce Kristensen and Mele’s methods and apply it to a
variety of stochastic volatility models of European style options as well as a
model for commodity futures. The focus of this thesis is the effect of
different model choices and different model parameter values on the numerical
stability of Kristensen and Mele’s approximation.

Keywords: Closed-form approximations, Option pricing theory, Stochastic
volatility, Continuous-time models, Commodity futures, Greeks

’As far as the laws of mathematics refer to reality, they are not certain; and
as far as they are certain, they do not refer to reality.’

\- Albert Einstein

###### Contents

  1. 1 Introduction
  2. 2 Closed-form approximations to asset prices and greeks
    1. 2.1 Preliminaries: Infinitesimal generators
    2. 2.2 Asset price approximations
    3. 2.3 Approximating greeks
    4. 2.4 Connection to other approaches
      1. 2.4.1 Risk-neutral densities and saddlepoint approximations
      2. 2.4.2 A similar expansion: Yang’s approach
  3. 3 Approximating the Heston model
    1. 3.1 KM expansion for the Heston model
      1. 3.1.1 Digression on optimal nuisance parameters
    2. 3.2 Accuracy of the approximation
      1. 3.2.1 Approximating Put Prices
      2. 3.2.2 Approximating Greeks in the Heston model
  4. 4 Approximating the CEV model
    1. 4.1 KM’s expansion for the CEV model
    2. 4.2 Numerical accuracy
  5. 5 Stochastic volatility with an OU - process
    1. 5.1 Approximation of the SZ model
    2. 5.2 Accuracy of the approximation
      1. 5.2.1 Approximating Greeks in the SZ model
  6. 6 Commodity futures with stochastic volatility
    1. 6.1 A baseline model for commodity futures
    2. 6.2 KM expansion
    3. 6.3 Numerical accuracy
  7. 7 Conclusion
  8. A Reference solutions
    1. A.1 Fourier transform solution of the Heston model
      1. A.1.1 Analytic Greeks in the Heston model
    2. A.2 Fourier transform solution of the SZ model
    3. A.3 Digression on complex logarithms
    4. A.4 MC simulation for the CEV model
    5. A.5 MC simulation of the Commodity Future
    6. A.6 Definition of the scale measure of the CEV variance process

###### List of Figures

  1. 1 Simulated vs. B-S implied volatility
  2. 2 KM approximation of the Heston [1993] model
  3. (a) Using $\sigma_{0}=\sqrt{v(t)}$
  4. (b) Using $\sigma_{0}=\sqrt{\theta}$
  5. 3 KM approximation of the Heston [1993] model
  6. (a) Time to maturity 1.5 years
  7. (b) Time to maturity 2.0 years
  8. (c) Time to maturity 4.0 years
  9. 4 Effect of maturity and moneyness on KM’s expansion
  10. 5 Influence of correlation on KM’s approximation
  11. (a) Time to maturity one year
  12. (b) Time to maturity 0.5 year
  13. 6 Effect of maturity and correlation on KM’s expansion
  14. 7 KM approximation of European Put option in the Heston [1993] model
  15. (a) Time to maturity one year
  16. (b) Time to maturity two year
  17. (c) Time to maturity 0.25 year
  18. 8 Different CEV parameter values
  19. (a) CEV with $\gamma=0.6$
  20. (b) CEV with $\gamma=1.33$
  21. 9 CEV: Correlation and Moneyness for $\gamma=1.33$
  22. (a) $T-t=0.5$
  23. (b) $T-t=1.0$
  24. 10 CEV: Correlation and Moneyness for $\gamma=0.6$
  25. (a) $T-t=0.5$
  26. (b) $T-t=1.0$
  27. 11 Level of elasticity of variance and moneyness
  28. 12 Approximation of the Stein and Stein [1991] model.
  29. (a) Time to maturity $T-t=0.25$
  30. (b) Time to maturity $T-t=0.50$
  31. (c) Time to maturity $T-t=1.00$
  32. 13 Approximation of the Schöbel and Zhu [1999] model.
  33. (a) Time to maturity $T-t=0.25$
  34. (b) Time to maturity $T-t=0.50$
  35. (c) Time to maturity $T-t=1.00$
  36. 14 Sensitivity towards increasing time to maturity
  37. 15 SZ-model with six corrective terms
  38. (a) Time to maturity equals one month
  39. (b) Time to maturity equals three months
  40. 16 MC result vs. KM approximation - maturity 0.5 years
  41. 17 MC result vs. KM approximation - maturity 0.25 years
  42. 18 MC result vs. KM approximation - maturity 1.0 years
  43. 19 $z_{0}$ path around the origin
  44. 20 Simulation of commodity prices
  45. (a) Histogram of simulated commodity prices at time $T$
  46. (b) Density of commodity prices at time $T$

###### List of Tables

  1. 1 Iterations of the pricing error for the Heston [1993] model
  2. 2 Approximating the Heston [1993] model
  3. 3 Approximating Heston’s $\Delta$
  4. 4 Approximating Heston’s $\Gamma$
  5. 5 Approximating Heston’s $\mathscr{V}$
  6. 6 Iterations of the pricing error for the CEV model
  7. 7 Approximating the CEV model - $\gamma=0.6$
  8. 8 Approximating the CEV model - $\gamma=1.33$
  9. 9 Iterations of the pricing error for the SZ/S&S-model
  10. 10 Approximating $\Delta$ for the SZ-model
  11. 11 Approximating $\Gamma$ for the SZ-model
  12. 12 Approximating SZ’s $\mathscr{V}$
  13. 13 Iterations of the pricing error for the Commodity Futures model

##  1 Introduction

In the field of financial engineering, the continuous-time framework had
become the golden standard in the modeling of derivative securities.
Continuous-time models for the prices of these securities are build on the
notion of stochastic differential equations (SDEs) to describe the dynamics of
the market of the assets or the risks on which the derivative is written on.
The price of a specific derivative security then usually arises from the
solution of a partial differential equation (PDE) which can be constructed
from the markets SDEs. One of the most celebrated of such solutions is the
Black and Scholes [1973] model, which resulted from the pioneering works of
Robert C. Merton, Fischer Black and Myron Scholes, and yields the price of an
European Call option written on a financial asset. While the solution to the
Back-Scholes PDE can be derived in closed-form, this is not possible for the
vast majority of the more realistic and thus more complex models nowadays used
by practitioners and in academia. The continuous-time framework, in fact, owes
its success at least partially to this very issue, as it opens convenient ways
to cope with security prices that are not available in closed-form.
One widely used approach are Fourier transform methods, which allow for so
called semi closed-form solutions. Where the word semi is attached since the
solution is only given in closed-form for the models characteristic functions,
which are complex valued functions of the underlying and time. The price of
the derivative security is then computed by numerically integrating over the
characteristic functions, what might be a mathematically non-trivial task and
may introduce sources of numerical instability which are not easy to
spot.111See the example on complex valued logarithms in Appendix A.3 of this
thesis. Two classical approaches to obtain approximations in case
(semi-)closed-form solutions are unavailable, are Finite Difference schemes
and Monte Carlo simulations. While both approaches are widely applicable they
might be computationally burdensome, as Finite Differences requires the
repeated solution of linear systems of equations, and Monte Carlo requires the
simulation of huge numbers of sample paths of the underlying based on
discretized SDEs over a fine grid of time increments.
Even though these approaches can be described as standards to cope with PDEs
for which closed-form solutions are not available, there appeared also one
branch of academic literature attempting to develop approaches to derive
closed-form approximations to PDEs. These closed-form approximations are
usually built on Taylor series expansions of conditional expectations. The
approach suggested by Kristensen and Mele [2011] (henceforth KM) fits into
this branch of literature. However, KM do not focus on approximating the asset
price directly, but target on developing an approximation of the expected bias
which arises if an oversimplifying baseline model (or "auxiliary" in KM’s
terms) would be used to value some derivative, e.g. when the Black-Scholes
model would be used while the true market dynamics are assumed to obey the
Heston model. The derivatives price is then computed as the sum of the price
resulting from the baseline model and the expected bias. Where the bias is
approximated by a power series expansion around the baseline model, what is
among the main features distinguishing KM’s approach from other closed-form
approximations available in the literature. In a nutshell, KM’s approximation
consists of two steps: (i) Firstly the Feynman-Kac representation of the
solution to a PDE which the bias obeys is derived. This Feynman-Kac
representation consists of two conditional expectations, one describing the
difference in the payoffs and one describing the difference in the market
dynamics between the true and the baseline model. (ii) The second step is to
approximate these expectations through a series expansion around the baseline
model. One important restriction regarding the choice of a suitable model
arises with the differentiability of the payoff function. KM’s approach can
only accommodate non-differentiable payoffs if one chooses a baseline model
that has a payoff identical to the one of the model of interest.
The general idea of approximating the pricing bias instead of the asset price
directly is not entirely new to the financial engineering literature. Hull and
White [1988] already suggested an approach to approximate the price of a
derivative in the Heston model through a series expansion of the bias between
the Heston and the Black-Scholes model.222Note that in 1988 there were no
solutions available for such models for non-zero correlation cases. (See Hull
and White [1988], pp. 30 - 31.) The model received the name Heston model after
Steve Heston published his influential work in 1991, showing how such models
can be solved in (semi-)closed-form through Fourier transform methods.
However, while Hull/White and KM are using similar PDEs to determine the
pricing bias their series expansions are very different. Specifically,
Hull/White use a series of the bias in the model’s volatility of variance
parameter, whereby the elements of which the series consists are functions of
the state variables and time which are defined through a system of
differential equations. Hull/White provide closed-form solutions for up to
three elements in their series.333See Hull and White [1988], equations (10)
and (16) to (18). Note that Hull/White’s and KM’s PDEs for the bias are only
consistent if one assumes that the Black-Scholes constant volatility parameter
is identical to Heston’s spot volatility parameter.
Despite its merit in approximating derivative prices, the approach of focusing
on computing a pricing bias is also interesting in its own right. Using the
same illustrative example as KM, assume that the true market dynamics would be
given by the Heston model but a trader uses the Black-Scholes delta to
determine the amount of stocks he needs to buy to hedge his short position in
an option. KM claim that the pricing bias obtained through their approach can
be interpreted as the expected total hedging costs from using a wrong model,
such as in the described situation. Such an interpretation of the bias is
confirmed by Elices and Gimenez [2013] who show that the expectation of the
costs of a hedging strategy is equal to the difference between the price
obtained through the true model and some simpler models used for valuation and
portfolio management.
Younesian [2013] tests the empirical option valuation performance of the CEV
and the Heston model both approximated by KM’s approach, whereby he uses data
on prices of S&P 500 call index options for the period of years from 2002 to
2011. Since the approximation has closed-form, Younesian [2013] is able to
estimate the model’s parameters directly from cross-sectional option data.
Nevertheless, Younesian [2013] reports that the obtained parameter estimates
are quite imprecise. Especially the correlation parameter $\rho$ in the CEV
model seems to be difficult to estimate if the model is approximated via KM’s
approach. As he notes this may also seriously flaw the estimates of the other
model parameters.
Garcia [2013] investigates how KM’s approximation performs when applied to
jump processes. Specifically, he compares two alternative approaches to
compute KM’s pricing bias for the CEV-Merton model using the Black-Scholes
model as baseline. Garcia [2013] reports to be able to only compute the first
two elements in the series expansion such that the obtain approximations are
comparatively imprecise and it is difficult to infer whether the series
expansion converges.
The implementation of KM’s series expansion requires the computation of
derivatives of the baseline model up to higher orders. If e.g. one wishes to
implement KM’s approximation for the Heston model, while using the Black-
Scholes model as baseline and using five terms in the series expansion, this
requires taking derivatives of the Black-Scholes model up to an order of 10 in
stock price direction as well as a number cross-derivatives in stock price and
time direction. While KM themselves do report how they implemented their
approximation, Garcia [2013] reports to have all baseline model derivatives
implemented through finite difference approximations. Whereas this appears to
be a convenient and flexible solution, it substantially increases the
computation times and also creates an approximation that does not have a fully
closed-form. Hence, I believe that this approach undermines the main
advantages offered by KM’s method. Younesian [2013] uses Maple to implement
the whole approximations, but reports difficulties in saving the obtained
expressions in Matlab files if expansions using more than two terms are
computed. Implementing KM’s approach by using Maple calls inside Matlab avoids
the need of saving the expressions to some file but slows down computation
time substantially.
In this thesis I will compute most parts of the series expansion manually. I
only use Maple to compute time-stock price cross-derivatives of the Black-
Scholes model, which I could successfully save to Matlab files. By using this
approach all of my KM approximations are given in closed-form completely.
Additionally, I am able to compute expansion using up to five or six terms
without any loss in computation speed.
In terms of models I will focus on stochastic volatility models, which provide
a good chance to test the performance and convergence behavior of KM’s
approximation when varying different model parameters. In particular, I will
apply the method to approximate option prices in the Heston and CEV model as
well as the Schöbel/Zhu model and to approximate Commodity Future prices in a
stochastic volatility model which is due to Lutz [2009]. Technically it is
straightforward to define a suitable baseline model for each stochastic
volatility model by simply taking the same price process but replacing
stochastic by constant volatility. In case of the option pricing examples this
leads to the Black-Scholes model as baseline, while I use model 1 from
Schwartz [1997] as baseline for the commodity future example.
The remainder of this thesis develops as follows: In section 2 KM’s asset
price approximation as well as KM’s approximation for hedge ratios are
introduced and the links to similar approximation approaches available in the
literature are shown. Section 3 shows the application to KM’s approximation
for option prices as well as different hedge ratios in the Heston model. This
is followed by a discussion on the performance in approximating a slight
generalization the Heston model - Jones [2003]’s CEV model - in section 4.
Section 5 continues the discussion by showing an application of KM’s approach
to option prices and hedge ratios under the non-affine Schöbel and Zhu [1999]
model. In section 6 the focus is shifted away from option prices towards
prices of commodity futures, where the price process itself is assumed to be
mean-reverting. Section 7 summarizes and concludes the discussion.

##  2 Closed-form approximations to asset prices and greeks

In the next two subsections I will, closely following Kristensen and Mele
[2011, sect. 3.1 & 3.2], state the KM asset price representation, its
approximation formula, as well as an approximation to partial derivatives of
the asset price.

###  2.1 Preliminaries: Infinitesimal generators

The approximation approach of KM uses heavily the notion of infinitesimal
operators. Hence, I will start the description of KM’s approximation by
briefly introducing these operator.
Let $z(t)$ be a $(d\times 1)$ vector of state variables and assume that under
the risk-neutral probability measure $z(t)$ satisfies the following vector
SDE,

$\displaystyle\begin{bmatrix}dz(t)_{1}\\\ \vdots\\\ dz(t)_{d}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\mu_{1}\\\ \vdots\\\
\mu_{d}\end{bmatrix}dt+\begin{bmatrix}\sigma_{11}&\cdots&\sigma_{1d}\\\
\vdots&\ddots&\vdots\\\
\sigma_{d1}&\cdots&\sigma_{dd}\end{bmatrix}\begin{bmatrix}dW(t)_{1}\\\
\vdots\\\ dW(t)_{d}\end{bmatrix}$ $\displaystyle dz(t)$
$\displaystyle=\mu\left(z(t)\right)dt+\sigma\left(z(t)\right)dW(t),~{}~{}t\in\left[0,\infty\right)$
(1)

where $dW(t)_{i},~{}i\in\left[1,d\right]$ are the increments of arbitrarily
correlated standard Brownian motion processes and $\mu\left(z(t)\right)$ and
$\sigma\left(z(t)\right)$ define some time-homogenous drift and diffusion
terms. Further, let $f(z(t))$ be a continuously twice differentiable function
of $z(t)$. Then e.g. Hirsa and Neftci [2014] define the expected rate of
change in $f(z(t))$ through the operator $\mathscr{A}$, such that

$\displaystyle\mathscr{A}f(z(t))=\lim_{\Delta\rightarrow
0}\frac{\mathbb{E}_{z,t}\left[f(z(t+\Delta))\right]-f(z(t))}{\Delta}$ (2)

Hirsa and Neftci [2014] state that since the expected change in $f(z(t))$ is a
smoother function than $f(z(t))$ itself, the limit in (2) can be defined even
as a Brownian motion process itself is non-differentiable. This derivative-
like notion of $\mathscr{A}$ relates the rate change directly to Itô’s Lemma.
Applying Itô’s Lemma to $f(z(t))$ and replacing the $dW(t)$ terms in resulting
expression by it’s drift, which is zero for $dW(t)$, and finally dividing the
remainder by $dt$ yields the same result as taking the limit in (2).444See
e.g. Hirsa and Neftci [2014], p. 356. Hence, the expected rate of change can
be expressed as

$\displaystyle\mathscr{A}f(z(t))=\sum^{d}_{i=1}\mu_{i}(z)\dfrac{\partial
f(z(t))}{\partial
z_{i}}+\dfrac{1}{2}\sum^{d}_{i=1}\sum^{d}_{j=1}\sigma^{2}_{ij}(z)\dfrac{\partial^{2}f(z(t))}{\partial
z_{i}\partial z_{j}}$ (3)

where
$\sigma^{2}(z(t))=\sigma(z(t))\sigma(z(t))^{\prime}\in\mathbb{R}^{d\times d}$
is the covariance matrix of (1).
The representation in (3) also shows clearly the connection to PDEs. Using the
definition of $\mathscr{A}$ the PDE for the derivative $f(z(t))$ can be
constructed as555Note that this is not a derivation of the PDE, but simply
shows the technical connection between the concepts of an infinitesimal
generator and the PDE. While $\mathscr{A}f(z(t))$ can be interpreted, as
discussed before, as an expected rate of change, the PDE determines the
intertemporal change in the derivatives value (See e.g. Fries [2007], pp. 46 -
47). It is easy to see that both concepts are mathematically related. However,
generally PDEs are derived by arbitrage arguments, such that their economic
meaning becomes clear. Throughout this thesis I will not derive the PDEs of
the considered models by arbitrage arguments in order to keep the focus of
this work on the method of developing closed-form approximations.

$\displaystyle\underbrace{\dfrac{\partial f(z,t)}{\partial
t}+\mathscr{A}f(z,t)}_{\equiv~{}\mathscr{L}f(z,t)}+c(z,t)=R(z,t)f(z,t)$ (4)

Where $R(z,t)$ is the instantaneous short-term rate and $c(z,t)$ denotes the
instantaneous coupon rate. Following KM it will be assumed that $c(z,t)$ is
identical zero and hence will be dropped in the following. Based on this
notion of the infinitesimal generator, KM define a new operator $\mathscr{L}$
as shown in the above expression.
In the subsequent section the operator $\mathscr{L}$ will be repeatedly
applied to different functions. Hence, I found it useful to also define the
operator independent from any specific function as

$\displaystyle\mathscr{L}$ $\displaystyle=\dfrac{\partial}{\partial
t}+\mathscr{A}$ $\displaystyle=\dfrac{\partial}{\partial
t}+\sum^{d}_{i=1}\mu_{i}(z)\dfrac{\partial}{\partial
z_{i}}+\dfrac{1}{2}\sum^{d}_{i=1}\sum^{d}_{j=1}\sigma^{2}_{ij}(z)\dfrac{\partial^{2}}{\partial
z_{i}\partial z_{j}}$ (5)

The notation in (5) emphasizes that applying the operator $\mathscr{L}$ to any
function practically means constructing a sum of the derivatives and cross-
derivatives of that function in direction of all of its state variables as
well as in time direction and weigh these derivatives with the elements of the
drift vector and the covariance matrix of the underlying system of SDEs.
Note that KM name the operator $\mathscr{L}$ an infinitesimal generator.
However, in the literature an infinitesimal generator often is defined as
$\mathscr{A}$ only.

###  2.2 Asset price approximations

KM’s approximation approach essentially builds on a power series expansion
around a ’baseline’ model for which an exact closed-form solution is
available. Hence, the approximation requires the formulation of two models.
First, assume a d-dimensional vector SDE as in (1) and assume that $f(z,t)$
denotes the price of a derivative written on $z(t)$.666For notational
convenience I will write $z(t)$ simply as $z$ and $f(z(t))$ as $f(z,t)$ to
emphasis that $f$ is also a function of time. Using the definition of the
operator $\mathscr{L}$ in (5), $f(z,t)$ satisfies the partial differential
equation

$\displaystyle\mathscr{L}f(z,t)-R(z,t)f(z,t)$ $\displaystyle=0$ (6)
$\displaystyle\text{s.t.}~{}~{}f(z,T)$ $\displaystyle=b(z)~{}$

Where $b(z)$ defines the payoff function at some maturity date $T>t$ and
$R(z,t)$ again denotes the instantaneous short-term interest rate. Further,
assume that there is no closed-form solution to (6) available.
As baseline model consider the price of another derivative $f_{0}(y,t)$, which
is written on a $(m\times 1)$ vector of state variables $y(t)$, where $m\leq
d$ and $y(t)$ follows the (vector) SDE,

$\displaystyle
dy(t)=\mu_{y}\left(y(t)\right)dt+\sigma_{y}\left(y(t)\right)dW_{y}(t),~{}~{}t\in\left[0,\infty\right)$
(7)

Where $\mu_{y}$, $\sigma_{y}$ are time-homogenous drift and diffusion terms
and $dW_{y}(t)$ denotes the increments of a m-dimensional standard Brownian
motion. If the baseline model is of lower dimension as the model we desire to
approximate, it is necessary to perform a correction in order to fit the
dimensions of the two models.
KM suggest to perform this correction in the following way: Let $z_{0}(t)$ be
a $(d\times 1)$ vector process, such that $z_{0}(t)=y(t)$ and
$dz_{0}(t)=dy(t)$ if $m=d$ and
$z_{0}(t)=\left[y(t)^{\prime}~{}\aleph\right]^{\prime}$ if $m<d$. Where
$\aleph$ is a $(d-m\times 1)$ vector of zeros. Further, if $m<d$ define the
SDE for $z_{0}(t)$ such that,

$\displaystyle dz_{0}(t)$
$\displaystyle=\mu_{0}\left(z_{0}(t)\right)dt+\sigma_{0}\left(z_{0}(t)\right)dW_{y}(t),~{}~{}t\in\left[0,\infty\right)$
(8) $\displaystyle\mu_{0,i}\left(z_{0}(t)\right)$
$\displaystyle=\begin{cases}\mu_{y,i}\left(y(t)\right)&\quad\text{if }1\leq
i\leq m\\\ 0&\quad\text{otherwise}\\\ \end{cases}$ (9)
$\displaystyle\sigma_{0,ij}\left(z_{0}(t)\right)$
$\displaystyle=\begin{cases}\sigma_{y,ij}\left(y(t)\right)&\quad\text{if
}1\leq i,j\leq m\\\ 0&\quad\text{otherwise}\\\ \end{cases}$ (10)

Hence, the drift and diffusion terms of the SDEs in (1) and (8) will be of
same dimension, what later on, will be necessary to obtain an expression for
the pricing difference between the two models. Note that $f_{0}(\cdot)$
satisfies the partial differential equation

$\displaystyle\mathscr{L}_{0}f_{0}(z_{0},t)$
$\displaystyle=R(z_{0},t)f(z_{0},t)$ (11)
$\displaystyle\text{s.t.}~{}~{}f(z_{0},T)$ $\displaystyle=b_{0}(z_{0})$

Where $b_{0}(z_{0})$ denotes the payoff of the derivative $f_{0}(\cdot)$ at
maturity. $\mathscr{L}_{0}$ is essentially the same operator as in (5), but
defined for the m-dimensional system (8).
If a closed-form solution to (11) exists, then $f_{0}(\cdot)$ can be used as a
baseline model to derive an approximate closed-form expression for $f(z,t)$.
Such a closed-form approximation can be derived in two steps.
First, denote the difference in the price between the unknown model $f(z,t)$
and the auxiliary model $f_{0}(\cdot)$ as $\Delta f(z,t)=f(z,t)-f_{0}(z,t)$.
If $f(z,t)$ and $f_{0}(z,t)$ are well-defined solutions to their respective
PDEs given in (6) and (11), then the pricing difference $\Delta f(z,t)$ will
satisfy the PDE777See Kristensen and Mele [2011], p. 394.

$\displaystyle\mathscr{L}\Delta f(z,t)+\delta(z,t)$
$\displaystyle=R(z,t)\Delta f(z,t)$ (12) $\displaystyle\text{s.t.}~{}~{}\Delta
f(z,T)$ $\displaystyle=d(z)$

Where $d(z)=b(z)-b_{0}(z)$ adjusts the boundary condition for possible
mismatches in the payoff functions of the two models. The second adjustment
term $\delta(z,t)$ corrects for the differences in the underlying driving
forces of the two models, i.e. for the differences between the SDEs given in
(1) and (8) which KM define as

$\displaystyle\delta(z,t)=(\mathscr{L}-\mathscr{L}_{0})f_{0}(z,t)$ (13)

Recalling the definition of the operator $\mathscr{L}$ (and $\mathscr{L}_{0}$)
it is straightforward to derive an expression for $\delta(z,t)$.

$\displaystyle\left(\mathscr{L}-\mathscr{L}_{0}\right)f_{0}(z,t)=$
$\displaystyle~{}\dfrac{\partial f_{0}}{\partial
t}+\mathscr{A}f_{0}-\dfrac{\partial f_{0}}{\partial t}-\mathscr{A}_{0}f_{0}$
$\displaystyle=$ $\displaystyle~{}\sum^{d}_{i=1}\mu_{i}\dfrac{\partial
f_{0}}{\partial
z_{i}}+\dfrac{1}{2}\sum^{d}_{i=1}\sum^{d}_{j=1}\sigma^{2}_{ij}(z)\dfrac{\partial^{2}f_{0}}{\partial
z_{i}\partial z_{j}}-\sum^{d}_{i=1}\mu_{0,i}\dfrac{\partial f_{0}}{\partial
z_{i}}+\dfrac{1}{2}\sum^{d}_{i=1}\sum^{d}_{j=1}\sigma^{2}_{0,ij}(z)\dfrac{\partial^{2}f_{0}}{\partial
z_{i}\partial z_{j}}$ $\displaystyle=$
$\displaystyle~{}\sum^{d}_{i=1}\left(\mu_{i}(z)-\mu_{0,i}(z)\right)\dfrac{\partial
f_{0}}{\partial
z_{i}}+\dfrac{1}{2}\sum^{d}_{i=1}\sum^{d}_{j=1}\left(\sigma^{2}_{ij}(z)-\sigma^{2}_{0,ij}(z)\right)\dfrac{\partial^{2}f_{0}}{\partial
z_{i}\partial z_{j}}$ $\displaystyle=$
$\displaystyle~{}\sum^{d}_{i=1}\Delta\mu_{i}\dfrac{\partial
f_{0}(z_{t},t)}{\partial
z_{i}}+\dfrac{1}{2}\sum^{d}_{i=1}\sum^{d}_{j=1}\Delta\sigma^{2}_{ij}(z)\dfrac{\partial^{2}f(z,t)}{\partial
z_{i}\partial z_{j}}$ (14)

Where $\Delta\mu_{i}(z)=\mu_{i}(z)-\mu_{0,i}(z)$ and
$\Delta\sigma^{2}_{ij}(z)=\sigma^{2}_{ij}(z)-\sigma^{2}_{0,ij}(z)$. Recalling
(2) and (3) the first line in (14) suggests that KM’s $\delta(z,t)$ represents
the difference in expected rate of change of $f_{0}$ under the true and the
baseline dynamics.
Since $f_{0}(\cdot)$ and $d(\cdot)$ are known, both adjustment terms are
straightforward to obtain and a solution to the PDE in (12) can be derived as
a function of the two adjustments $d(z)$ and $\delta(z,t)$. KM do this by
applying Theorem 7.6 in Karatzas and Shreve [1991], which states the Feynman-
Kac representation of the solution to (12)

$\displaystyle\Delta
f(z,t)=\mathbb{E}_{z,t}\left[e^{\left(-\int_{t}^{T}R\left(z(s),s\right)ds\right)}d(z(T))+\int_{t}^{T}e^{\left(-\int_{t}^{T}R\left(z(u),u\right)du\right)}\delta(z(s),s)ds\right]$
(15)

Recalling that $\Delta f(z,t)=f(z_{t},t)-f_{0}(z_{t},t)$ and rearranging
yields KM’s Theorem 1:

$\displaystyle f(z,t)=f_{0}(z,t)$
$\displaystyle+\mathbb{E}_{z,t}\left[e^{\left(-\int_{t}^{T}R\left(z(s),s\right)ds\right)}d(z(T))\right]$
$\displaystyle+\int_{t}^{T}\mathbb{E}_{z,t}\left[e^{\left(-\int_{t}^{T}R\left(z(u),u\right)du\right)}\delta(z(s),s)\right]ds$
(16)

The expression in (15) yields an exact formulation of the pricing bias between
the baseline and the true model. Hence, if $f(z,t)$ represents e.g. the Heston
[1993] model and this is assumed to reflect the true market dynamics while
$f_{0}(z,t)$ represents e.g. the Black-Scholes model and this is used to
manage a portfolio of derivative contracts, then (15) yields the bias from
stochastic volatility. KM note that this bias can be interpreted as the
expected hedging costs arising from the use of the wrong model. By using a
different approach Elices and Gimenez [2013] show that this interpretation of
the bias or pricing difference $\Delta f(\cdot)$ indeed applies. Note that the
right hand side in (15) depends only on the solution known solution to the
baseline.
Since the aim of KM’s method is to derive asset price approximations, it is
more useful to consider the formulation in (16). However, the expression in
(16) relies on the evaluation of two conditional moments, which have no
trivial solution. Generally there would be several possibilities to evaluate
(16). As KM note, it would be possible to compute the two conditional moments
by Monte Carlo integration. In such a case the chosen baseline model would
have a similar role as a control variate. Hence, KM conclude that there might
be no obvious advantage to directly estimating $f(z,t)$ using Monte Carlo with
$f_{0}(z,t)$ as control variate.
This directly leads to the second step of KM’s approach, which is to
approximate the conditional moments in (16) by a series expansions. The
expansion suggested for this purpose is given in KM’s Definition 1:

$\displaystyle
f_{N}(z,t)=f_{0}(z,t)+\sum_{n=0}^{N}\dfrac{1}{n!}d_{n}(z,t)\left(T-t\right)^{n}+\sum_{n=0}^{N}\dfrac{1}{(n+1)!}\delta_{n}(z,t)\left(T-t\right)^{n+1}$
(17)

Where $\delta_{0}(z,s)$ is defined as in (14) and $d_{0}(z,t)=b(z)-b_{0}(z)$.
While for $n>0$

$\displaystyle\delta_{n}(z,s)=$
$\displaystyle\mathscr{L}\delta_{n-1}(z,t)-R(z,t)\delta_{n-1}(z,t)$ (18)
$\displaystyle d_{n}(z,s)=$
$\displaystyle\mathscr{L}d_{n-1}(z,t)-R(z,t)d_{n-1}(z,t)$ (19)

The number $N$ denotes the order of the approximation, such that $f_{N}(z,t)$
denotes the $N^{th}$ order approximation to $f(z,t)$. Under some regulatory
conditions, which I will discuss in a moment, $f_{N}(z,t)\rightarrow f(z,t)$
as $N\rightarrow\infty$. KM note that, in order to compute the iterations on
the adjustment terms, it is must be ensured that both, $\delta(z,t)$ and
$d(z,t)$, are $2N$ times differentiable with respect to $z$ and additionally
that $\delta(z,t)$ is $N$ times differentiable with respect to $t$. Hence, if
the payoff function $b(z)$ is non-differentiable, as it is the case e.g. for a
plain vanilla options, the choice of a baseline model is restricted to models
with $b_{0}(z)=b(z)$, such that $d_{n}(z,t)=0~{}\forall~{}n$.
The formula in (17) provides a method to obtain a closed-form approximation of
the asset price or the expected hedging costs even for cases in which the true
model is not known in closed-form. Hence, computation will be fast. This
feature provides the greatest advantage of KM’s method over approaches like
Elices and Gimenez [2013], which is computationally costly as it is based on
simulations of hedging strategies.888See Elices and Gimenez [2013], p. 1016.
Both steps of KM’s approximation approach rely on a number of assumption.
These are given in KM’s Proposition A.1 & A.2, which I summarize as follows:

  1. 1.

Since (16) is basically just an application of the Feynman-Kac formula, the
expression holds under the same regulatory conditions under which this formula
holds. These condition are given in Karatzas and Shreve [1991, Theorem 5.7.6]
or in Kristensen and Mele [2011, condition (A.3)-(A.5), p. 412]. Among these
conditions only the linear growth condition which is imposed on the drift and
diffusion terms of the SDEs may not hold in some models considered in finance.
However, this condition only serves to ensure the existence of the solution to
the model’s SDEs. Hence, if existence can be ensured for both the baseline and
the true model by other means then the linear growth condition can be
ignored.999For the last statement see the footnote in KM, p. 394.

  2. 2.

In order to ensure that the expansion in (17) converges assumptions on the
operator $\mathscr{L}$ need to be imposed. Specifically, $\mathscr{L}$ must
have a transition density $p_{t}(y|x)$ with respect to Lebesgue measure.
Hence, the considered diffusion model must have a transition density, what is
generally true for most diffusion models. Further, there must be a measure
$\pi$, such that $\pi(x)p_{t}(y|x)=\pi(y)p_{t}(x|y)$. This can be understood
as a generalization of time-reversibility, which e.g. is always satisfied for
stationary, univariate processes.

  3. 3.

The last set of assumption concerns the corrective terms $\delta(z,t)$ and
$d(z,t)$. Assume that for a function space $H$, which contains some function
norm $||\cdot||_{H}$, there exists some $\hat{\tau}>0$ and some functions
$\psi_{\delta},\psi_{d}\in H$ such that
$\delta:\mathbb{R}^{n}\times\mathbb{R}_{+}\rightarrowtail\mathbb{R}$ and
$d:\mathbb{R}^{n}\times\mathbb{R}_{+}\rightarrowtail\mathbb{R}$ satisfy
$\mathbb{E}\left[\psi_{\delta}\left(x(\hat{\tau})\right)\mid
x(0)=x\right]=\delta(x,\hat{\tau})$ and
$\mathbb{E}\left[\psi_{d}\left(x(\hat{\tau})\right)\mid
x(0)=x\right]=d(x,\hat{\tau})$, where additionally
$t\rightarrowtail\delta(x,t)$ is analytic uniformly in $||\cdot||_{H}$. In
words this means that it must be possible to match the two adjustment terms
$\delta$ and $d$ by conditional moments. Further, it is required that the
instantaneous short-term rate is time homogeneous and $sup_{x}\mid
R(x,t)=R(x)\mid<\infty$.

These assumptions are difficult to verify for a specific model.101010See
Kristensen and Mele [2011], p. 411. In fact, KM do not verify them in their
experiments but rather show empirically that their expansion converges for the
models they consider.

###  2.3 Approximating greeks

So far the KM approximation was introduced as a method to obtain prices of
financial derivatives or hedging costs. Besides this the computation of hedge
ratios for a derivative contract, the so called greeks, is of equal importance
in many applications, such as risk management. The greeks are defined as
partial derivatives of the respective pricing formula with respect to relevant
parameters and hence can be interpreted as the sensitivity of the derivatives
price to changes in these parameters. I will consider the following greeks in
the numerical examples below:

  1. 1.

$\Delta~{}=~{}\dfrac{\partial C}{\partial S}$ : Sensitivity of the derivatives
price to changes in the underlying

  2. 2.

$\Gamma~{}=~{}\dfrac{\partial^{2}C}{\partial S^{2}}$ : Sensitivity of the
derivatives price to changes in $\Delta$

  3. 3.

$\mathscr{V}~{}=~{}\dfrac{\partial C}{\partial v}$ : Sensitivity of the
derivatives price to changes in spot variance

From the definition of the greeks it is obvious that their computation is more
or less trivial as long as a closed-form solution for the derivatives price is
available. For cases in which a closed-form solution is not available the
literature developed a number of numerical approaches which are often based on
Monte Carlo simulation. However, these approach often increase estimation risk
and, if the payoff function is non-differentiable, straightforward procedures
such as finite differences perform poorly.111111See Lehoczky [2012], p. 4 - 6,
for an example on the estimation of the $\Delta$ of a Knock-in option via
finite differences. Lehoczky [2012] provides a good overview of publications
on the computation of greeks.
KM argue that their approximation as given in (17) can be used to also
estimate the greeks of the derivative by differentiating the approximation
formula with respect to the parameter of interest. Hence, KM state the
following expression to compute the above mentioned greeks

$\displaystyle\dfrac{\partial^{k}f_{N}(z,t)}{\partial
z^{k}_{i}}=\dfrac{\partial^{k}f_{0}(z,t)}{\partial
z^{k}_{i}}+\sum_{n=0}^{N}\dfrac{\left(T-t\right)^{n}}{n!}\dfrac{\partial^{k}d_{n}(z,t)}{\partial
z^{k}_{i}}+\sum_{n=0}^{N}\dfrac{\left(T-t\right)^{n+1}}{(n+1)!}\dfrac{\partial^{k}\delta_{n}(z,t)}{\partial
z^{k}_{i}}$ (20)

where the first and second partial derivatives of the adjustment terms are
given by $k=1$ and $k=2$. KM provide the results for the analytic derivatives
of first and second order for these adjustment terms.
Alternatively, KM suggest that these derivatives could be evaluated
numerically, e.g. by finite-differences. Note, that the before mentioned poor
performance of finite difference methods in the case of non-differentiable
payoffs does not apply here. If the payoff function is non-differentiable,
then a baseline model must be chosen such that $d_{n}(z,t)=0~{}\forall~{}n$,
whereby the function $\delta_{n}(z,t)$ necessarily needs to be differentiable
in any region relevant for the computation of greeks. Since non-
differentiability of the payoff function is given for many derivatives I will
restrict the attention to that case.
Since the approximation expressions obtained from (17) becomes quite lengthy
even for small numbers of corrective terms analytic differentiation might be
cumbersome. Hence, numerical differentiation seems to be a more efficient
choice.
In order to obtain numerical first and second derivatives of in (20) I apply
central differences to each of the corrective terms, while using the analytic
derivative of the baseline model. Hence, the derivatives of the corrective
terms in (20) are estimated as

$\displaystyle\dfrac{\partial\delta_{n}(z,t)}{\partial z_{i}}~{}\approx$
$\displaystyle~{}\frac{\delta_{n}(z_{i}+h,t)-\delta_{n}(z_{i}-h,t)}{2h}$ (21)
$\displaystyle\dfrac{\partial^{2}\delta_{n}(z,t)}{\partial
z_{i}^{2}}~{}\approx$
$\displaystyle~{}\frac{\delta_{n}(z_{i}+h,t)+2\delta_{n}(z_{i},t)-\delta_{n}(z_{i}-h,t)}{h^{2}}$
(22) $\displaystyle\text{for }n=$ $\displaystyle 1,\ldots,N$

The numerical derivatives of $d_{n}(z,t)$ could be computed analogously. Using
central differences yields an error of order $\mathscr{O}\left(h^{2}\right)$
for the approximation of the derivatives. This error comes additionally to the
overall approximation error of KM’s method. However, my numerical results
shown in the later sections suggest that this approach nevertheless achieves
sufficient accuracy. Following the suggestion in Hermann [2006, p. 257 & Ch.
8.1.3] the step size $h$ is chosen to be a function of the differentiation
variable and the machine precision of the computer used to perform the
computations, such that

$\displaystyle h=\left(z_{i}+1\right)10^{\frac{log_{10}(eps)}{3}-1}$ (23)

Where $eps$ denotes the machine precession. On the system I used to run all
Matlab scripts related to this thesis, $eps$ had a value of $2.2204\cdot
10^{-16}$. Note that, while in general the step size should be chosen as small
as possible to minimize truncation errors, choosing too small steps would lead
to round-off errors.121212See e.g. Gill et al. [1981], p. 11 - 12. Hence, (23)
attempts to balance these two effects. Note that the specific functional form
in (23) does not necessarily lead to an optimal step size. Better solutions
for $h$ might be available in the literature. Also, one might try the complex-
step approach suggested by Martins et al. [2003]. Therein the real function of
which we want to estimate the derivatives is transformed into the complex
plane. Since this yields a situation in which no difference operation is
required to estimate the derivative the previously mentioned trade-off between
truncation and round-off error is circumvented.131313See Martins et al. [2003,
Eq. 6] However, as I will show in section 3.2.2, the specification in (23)
turned out to yield precisely the same results as were reported in KM’s
original publication, so I decided to keep this definition. Hermann [2006,
Programm 5.2, p. 257] provides Matlab code to estimate the Jacobi-Matrix of a
function based on forward differences and the step size in (23). I adapted
this code to estimate the four greeks previously defined.

###  2.4 Connection to other approaches

In the finance literature there are two closed-form approximation approaches
bearing an especially close relationship to KM’s approach, that have not been
mentioned in the introduction. In the following I will provide a short
overview of these approaches by summarizing KM’s treatment of the two methods.

####  2.4.1 Risk-neutral densities and saddlepoint approximations

The closed-form approximation approach from KM closely relates to saddlepoint
approximations of risk-neutral transition densities and cumulative
distribution functions, which also can be used to approximate prices of
derivatives and, as claimed by KM, arise as a special case of their expansion.
The notion of saddlepoint approximations to probability distribution dates
back to work of Daniels [1954] and had been applied in the context of option
pricing e.g. by Rogers and Zane [1999] and Ait-Sahalia and Yu [2006]. Whereas
Goutis and Casella [1999] provide a good introduction on how this method can
be applied. As noted by Ait-Sahalia and Yu [2006], saddlepoint approximations
can be applied to models for which a characteristic function of the transition
density can found in closed-form, but the transition density itself is
unknown. Hence, the method can be seen as an alternative to Fourier inversion,
but comes with the advantage of avoiding numerical integration.141414See
Rogers and Zane [1999], p. 494. To see the connection to KM’s approach recall
from the definition of the Feynmann-Kac theorem that the price of a derivative
equals the discounted exception of its payoff, i.e.
$f(z,t)=\mathbb{E}\left[e^{\left(-\int_{t}^{T}R(z(s),s)\right)ds}b(z(T))\right]$
. Where $b(z(T))$, as before, denotes the assets payoff. Following KM, the
instantaneous short-term rate is assumed to be identical zero within this
section, such the exponential term in the expectation drops out. Again
considering the model of interest $f(z,t)$ and a suitable baseline model
$f_{0}(z,t)$, both asset prices can be expressed through the Feynmann-Kac
formula as151515See Kristensen and Mele [2011], p. 396.

$\displaystyle f(z,t)=$
$\displaystyle\int_{\mathbb{R}^{d}}b(z(T))p\left(z(T),T\mid z,t\right)dz$ (24)
$\displaystyle f_{0}(z,t)=$
$\displaystyle\int_{\mathbb{R}^{d}}b_{0}(z(T))p_{0}\left(z(T),T\mid
z,t\right)dz$ (25)

Where the expectation is expressed as an Riemann-integral and $p$ and $p_{0}$
denote the risk-neutral transition densities of the two models. If the
baseline model is chosen to match the payoff of the actual model at maturity,
i.e. $b(z(T))=b_{0}(z(T))$, the pricing bias between the models is

$\displaystyle f(z,t)-f_{0}(z,t)=\int_{\mathbb{R}^{d}}b(z(T))\Delta
p\left(z(T),T\mid z,t\right)dz$ (26)

Where $\Delta p=p-p_{0}$ denotes the transition discrepancy, the difference
between the two transition densities. Comparing (26) with (15) and recalling
that the short-term rate is assumed to be zero suggests that161616See also
Kristensen and Mele [2011], p. 397.

$\displaystyle\int_{\mathbb{R}^{d}}b(z(T))\Delta p\left(z(T),T\mid
z,t\right)dz=\int_{t}^{T}\mathbb{E}\left[\delta(z(s),s)\mid z,t\right]ds$ (27)

Such that applying KM’s approximation would be equivalent to computing a
derivatives price through an estimate of its transition discrepancy. Indeed,
Kristensen and Mele [2011, Appendix B] provide a proof that the above equation
holds. Nevertheless, the mere equality in (27) is not of practical use until
it is specified how the integrals on the left and the right hand side should
be evaluated. For the left-hand side, saddlepoint approximations can be used
to first estimate the transition discrepancy and then the pricing error can be
computed through the evaluation of the Riemann-integral. KM’s approach, as
described at the beginning, directly approximates the pricing error by means
of a series expansion around the baseline model. However, KM note that the
equality in (27) hinges on the assumption that the instantaneous short-term
rate is identical zero. In section 6, I will show that KM’s approximation
applied to a future indeed yields a representation of the pricing bias as on
right hand side of (27). Such that the use of saddlepoint approximations for
transition discrepancies and KM’s series expansion of the pricing bias indeed
are equivalent in this case. However, this is not true for e.g. option pricing
models.

####  2.4.2 A similar expansion: Yang’s approach

Yang [2006] suggests a closed-form approximation for option prices in the
context of stochastic volatility models, which appears to be very similar to
the approach suggested by KM. They provide a good comparison between their
approximation and the one in Yang [2006], which I attempt to summarize in this
section.171717The approach is a part of Yang’s doctoral dissertation, which
unfortunately was not accessible. Hence, while I still include the original
reference to Yang’s dissertation it should be noted that this subsection
solely summarizes the treatment of the method as described in KM.
Yang’s approximation is also built on the notion of a pricing bias between
some baseline model $f^{0}(z,t)$ known in closed-form and some unknown model
$f(z,t)$. Following the treatment of the method in KM the instantaneous short-
term rate $R(z,t)$ is assumed to be identical zero, and the true model and the
baseline model are assumed to have the same payoff at maturity. The difference
between the two models, denoted by $\Delta f(z,t)=f(z,t)-f^{(0)}(z,t)$,
satisfies the PDE

$\displaystyle\mathscr{L}_{0}\Delta f(z,t)+\tilde{\delta}(z,t)=0$ (28)
$\displaystyle\text{with
}\tilde{\delta}(z,t)=\left(\mathscr{L}-\mathscr{L}_{0}\right)f(z,t)$ (29)

The two operators $\mathscr{L}$ and $\mathscr{L}_{0}$ are defined the same as
for KM’s approximation. However, KM note two important differences between
their and Yang’s approximation. First, in (28) the operator associated with
the baseline model, instead the one of the true model, is applied to the
difference in the two asset prices. Second, the operator
$\left(\mathscr{L}-\mathscr{L}_{0}\right)$, even though it is identical to the
one used in KM’s approximation, is applied to the true model instead of the
baseline model. Hence, $\tilde{\delta}$ can not be expressed in closed-form,
as directly as KM’s $\delta$. By applying the Feynmann-Kac formula to (28) a
solution similar to the one in (16) is obtained

$\displaystyle
f(z,t)=f^{(0)}(z,t)+\int_{t}^{T}\mathbb{E}^{0}\left[\tilde{\delta}(z,t)\mid
z,t\right]$ (30)

Where $\mathbb{E}^{0}$ indicates that the expectation is taken under the
probability measure of the baseline model. Note that in the case of KM’s
approach this expectation is taken under the probability of the true model.
Again similar to KM, Yang suggests a series expansion to approximate the
integral expression

$\displaystyle f(z,t)=f^{(0)}(z,t)+\sum_{m=1}^{M}f^{(m)}(z,t)$ (31)

Where each $f^{(m)}(z,t)$ is a solution to the PDE
$\mathscr{L}_{0}f^{(m)}+\left(\mathscr{L}-\mathscr{L}_{0}\right)f^{(m-1)}=0$
with boundary condition $f^{(m)}(z,T)=0$. Such that the $f^{(m)}(z,t)$ can be
computed recursively starting from the known $f^{(0)}(z,t)$ until
$f^{(M)}(z,t)$. Yang suggests to do so by applying the Feynmann-Kac formula to
each of the PDEs

$\displaystyle
f^{(m)}=\int_{t}^{T}\mathbb{E}^{0}\left[\left(\mathscr{L}-\mathscr{L}_{0}\right)f^{(m-1)}\mid
z,t\right]$ (32)

Note that for $m=1$ the expression is, besides the expectation operator, the
same as in (16) of KM’s approximation and thus the expression inside the
expectation can be expressed in closed-form. The whole recursion can be solved
e.g. using standard symbolic software. KM compare the performance of their
method to the performance of Yang’s expansion for the Heston [1993] as well as
the more general CEV model. They find that, while both models perform
reasonably well, Yang’s expansion is slightly more accurate in approximating a
non-affine formulation of the CEV model, but is slightly less accurate for the
Heston model.

##  3 Approximating the Heston model

The abstract description of KM’s expansion in section 2.2 already showed that
the implementation of an approximation can be cumbersome if a large number of
corrective terms would be required to achieve precise approximations. Hence,
for practical applications, the approximation should converge fast to the true
price of the derivative contract of interest, such that only a small number of
corrective terms need to be computed. In order to assess the convergence
behavior of the approximation I follow KM and apply their expansion to the
well-known Heston [1993] model. Since for this model (semi-)closed form
solutions are available it is straightforward to analyze the accuracy of the
approximation. Hence, in this section I first attempt to replicate the results
already shown in KM, but also to go further and analyze in greater depth the
effect of different parameter values on the convergence behavior of the
approximation.
The Heston model is defined by the following system of SDEs

$\displaystyle dS(t)~{}=$ $\displaystyle~{}rSdt+\sqrt{v(t)}SdW_{1}(t)$ (33)
$\displaystyle dv(t)~{}=$
$\displaystyle~{}\kappa\left(\theta-v(t)\right)dt+\omega\sqrt{v(t)}dW_{2}(t)$
(34) $\displaystyle dW_{1}(t)dW_{2}(t)~{}=$ $\displaystyle~{}\rho dt$

Where $\kappa$ denotes the speed of mean-reversion parameter, $\theta$ and
$\omega$ the long-run variance and the instantaneous volatility of variance
respectively. The instantaneous short-term rate is constant over time and all
states, such that $R(S,t)=r~{}\forall~{}S,t$ in the notation of section 2.2.
This implies, i.e. for equation (16), that
$\int_{t}^{T}R(S,u)du=r\left(T-t\right)$.

###  3.1 KM expansion for the Heston model

In this section I will use the principles outlined in section 2.2 to derive
the KM approximation of the Heston model. As a first step one needs to derive
the vector $\mu(S(t),v(t))$ and the matrix $\sigma(S(t),v(t))$ as described in
(1).
Trivially,
$\mu(S(t),v(t))=\left[rS~{}~{}\kappa\left(\theta-v(t)\right)\right]^{\prime}$.
$\sigma(S(t),v(t))$ is derived by decomposing $dW_{1}(t),~{}dW_{2}(t)$ into
independent Wiener processes denoted by $\hat{dW_{1}(t)},~{}\hat{dW_{2}(t)}$,
such that

$\displaystyle\begin{bmatrix}dW(t)_{1}\\\
dW(t)_{2}\end{bmatrix}=\begin{bmatrix}1&0\\\
\rho&\sqrt{1-\rho^{2}}\end{bmatrix}\times\begin{bmatrix}\hat{dW(t)_{1}}\\\
\hat{dW(t)_{2}}\end{bmatrix}$ (35)

Where $\hat{dW_{1}(t)}\hat{dW_{2}(t)}=0$. From this $\sigma(S(t),v(t))$ as in
(1) is given by

$\displaystyle\begin{bmatrix}\sqrt{v(t)}S&0\\\
0&\omega\sqrt{v(t)}\end{bmatrix}\times\begin{bmatrix}1&0\\\
\rho&\sqrt{1-\rho^{2}}\end{bmatrix}=\begin{bmatrix}\sqrt{v(t)}S&0\\\
\rho\omega\sqrt{v(t)}&\sqrt{1-\rho^{2}}\omega\sqrt{v(t)}\end{bmatrix}$ (36)

Hence, the covariance matrix
$\sigma^{2}(\cdot)=\sigma(\cdot)\sigma(\cdot)^{\prime}$ of the Heston [1993]
model is

$\displaystyle\begin{bmatrix}\sqrt{v(t)}S(t)&0\\\
\rho\omega\sqrt{v(t)}&\sqrt{1-\rho^{2}}\omega\sqrt{v(t)}\end{bmatrix}$
$\displaystyle\times\begin{bmatrix}\sqrt{v(t)}S(t)&\rho\omega\sqrt{v(t)}\\\
0&\sqrt{1-\rho^{2}}\omega\sqrt{v(t)}\end{bmatrix}=$
$\displaystyle=\begin{bmatrix}v(t)S(t)^{2}&\rho\omega v(t)S(t)\\\ \rho\omega
v(t)S(t)&\omega^{2}v(t)\end{bmatrix}$ (37)

Denoting the price of a European style call option written on $S(t)$ by
$C\left(S(t),t\right)$ and applying the operator $\mathscr{L}$ to this price
one obtains,

$\displaystyle\mathscr{L}C\left(S(t),t\right)=\dfrac{\partial C}{\partial
t}+rS\dfrac{\partial C}{\partial S}+$
$\displaystyle+\kappa\left(\theta-v\right)\dfrac{\partial C}{\partial v}+$
$\displaystyle\dfrac{1}{2}\left(vS^{2}\dfrac{\partial^{2}C}{\partial
S^{2}}+\omega^{2}v\dfrac{\partial^{2}C}{\partial v^{2}}\right)+\rho\omega
vS\dfrac{\partial^{2}C}{\partial S\partial v}$ (38)

Such that the PDE of a European style call in the Heston model can be
constructed from this expression as

$\displaystyle rC\left(S,t\right)=\dfrac{\partial C}{\partial
t}+rS\dfrac{\partial C}{\partial S}+$
$\displaystyle+\kappa\left(\theta-v\right)\dfrac{\partial C}{\partial v}+$
$\displaystyle\dfrac{1}{2}\left(vS^{2}\dfrac{\partial^{2}C}{\partial
S^{2}}+\omega^{2}v\dfrac{\partial^{2}C}{\partial v^{2}}\right)+\rho\omega
vS\dfrac{\partial^{2}C}{\partial S\partial v}$ (39) $\displaystyle\text{s.t.
}b(S,T)~{}=$ $\displaystyle~{}max\left[S-K,0\right]$

In a next step one needs to define a suitable baseline model at which KM’s
series expansion can be developed. As mentioned earlier the Black-Scholes
model would be a convenient choice in the set up of stochastic volatility
models. Assume that the underlying follows the same process as in (33), but
replacing $\sqrt{v(t)}$ with a constant $\eta_{0}$, i.e.

$\displaystyle dS(t)~{}=$ $\displaystyle~{}rSdt+\eta_{0}SdW_{1}(t)$ (40)

By applying the principles for defining a baseline model described in section
2.2 the drift vector and the covariance matrix of the Black-Scholes model are
given by

$\displaystyle\mu_{0}=\left[rS~{}~{}0\right]^{\prime}\quad\text{,}\quad\sigma^{2}_{0}(\cdot)=\begin{bmatrix}\eta_{0}^{2}S^{2}&0\\\
0&0\end{bmatrix}$ (41)

Denoting by $C^{BS}(S,t)$ the price of a European call option under the Black-
Scholes model the models PDE can be again written by using the definition of
the operator $\mathscr{L}$

$\displaystyle rC^{BS}\left(S,t\right)$ $\displaystyle=\dfrac{\partial
C^{BS}}{\partial t}+rS\dfrac{\partial C^{BS}}{\partial
S}+\dfrac{1}{2}\eta_{0}^{2}S^{2}\dfrac{\partial^{2}C^{BS}}{\partial S^{2}}$
(42) $\displaystyle\text{s.t. }b(S,T)~{}=$
$\displaystyle~{}max\left[S-K,0\right]$

The solution to this PDE is the well-known Black-Scholes formula

$\displaystyle C^{BS}(S,t)$
$\displaystyle=\mathbb{N}(d_{1})S-\mathbb{N}(d_{2})Ke^{-r(T-t)}$ (43)
$\displaystyle d_{1}$
$\displaystyle=\frac{log(S/K)+\left(r+\sigma^{2}_{0}/2\right)\left(T-t\right)}{\eta_{0}\sqrt{T-t}}$
$\displaystyle d_{2}$ $\displaystyle=d_{1}-\sigma_{0}\sqrt{T-t}$

Where $\mathbb{N}(\cdot)$ denotes a standard normal distribution.
Using the PDEs in (39) and (42) as well as the definition in (12) the PDE of
the pricing difference $\Delta
C(S,t)=C\left(S,t\right)-C^{BS}\left(S,t\right)$ is defined as

$\displaystyle\mathscr{L}\Delta C(S,t)+\delta(S,t)=r\Delta C(S,t)$ (44)

As both model specifications, Heston and Black-Scholes, share identical
boundary conditions $d_{n}(S)=0~{}\forall~{}n$ in (17). As previously
mentioned it is necessary to match the boundary condition if the approximated
model does have a non-differentiable boundary.
Applying (14) yields the initial pricing error

$\displaystyle\delta_{0}(S,t)=$ $\displaystyle\left(rS-
rS\right)\dfrac{\partial C^{BS}}{\partial
S}+\left(\kappa\left(\theta-v(t)\right)-0\right)\dfrac{\partial
C^{BS}}{\partial v}+$
$\displaystyle\dfrac{1}{2}\left(\left(v(t)S^{2}-\eta_{0}^{2}S^{2}\right)\dfrac{\partial^{2}C^{BS}}{\partial
S^{2}}+\left(\omega^{2}v(t)-0\right)\dfrac{\partial^{2}C^{BS}}{\partial
v^{2}}\right)+\left(\rho\omega
Sv(t)^{-}0\right)\dfrac{\partial^{2}C^{BS}}{\partial S\partial v}$
$\displaystyle=$
$\displaystyle\dfrac{1}{2}\left(v(t)-\eta_{0}\right)S^{2}\dfrac{\partial^{2}C^{BS}}{\partial
S^{2}}$ (45)

Where the last equality in (45) results because $\partial C^{BS}/\partial
v=\partial^{2}C^{BS}/\partial v^{2}=\partial^{2}C^{BS}/\partial S\partial
v=0$. Below Table 1 shows how the series of corrective terms develops with the
order of the approximation, starting from (45) up to any positive integer $N$.
Each of the rows in Table 1 is developed according to the rule in (18). KM
note that, through the initial pricing error only a convexity adjustment is
added to Black-Scholes price, while for $N=1$ and $N=2$ (see also Table 1)
more information about the variance process is included.
Using the equations from Table 1 in (17) yields the $N^{th}$ order KM
approximation for a Heston call option

$\displaystyle C_{N}(S,t;\eta_{0})=C^{BS}(S,t,\eta_{0})+$
$\displaystyle(T-t)\delta_{0}+\frac{(T-t)^{2}}{2}\delta_{1}+$
$\displaystyle+\frac{(T-t)^{3}}{6}\delta_{2}+\ldots+\frac{(T-t)^{N+1}}{(N+1)!}\delta_{N}$
(46)

Since $C_{N}(S,t;\sigma_{0})\rightarrow C(S,t)$ as $N\rightarrow\infty$ each
additional corrective term $\delta_{n}$ should increase the accuracy of the
approximation. However, such convergence should occur already after only a
view corrective terms since otherwise the approach would loose it tangibility.
A feature of the KM approximation is that the approximated price $C_{N}$
depends on the constant volatility parameter of the Black-Scholes model while
the Heston model does not. Such a dependence on one or more nuisance
parameters originating from the chosen baseline model will usually arise in
the application of the KM approximation. Hence, it is necessary to make a
choice regarding the Black-Scholes volatility.

Table 1: Iterations of the pricing error for the Heston [1993] model

n |  Pricing error $\delta_{n}(S,t)$
---|---
0 |  $\dfrac{1}{2}\left(v-\eta_{0}\right)S^{2}\dfrac{\partial^{2}C^{BS}}{\partial S^{2}}$
1 |  $\dfrac{\partial\delta_{0}}{\partial t}+rS\dfrac{\partial\delta_{0}}{\partial S}+\kappa\left(\theta-v\right)\dfrac{\partial\delta_{0}}{\partial v}+\frac{1}{2}vS^{2}\dfrac{\partial^{2}\delta_{0}}{\partial S^{2}}+\rho\omega Sv\dfrac{\partial^{2}\delta_{0}}{\partial S\partial v}-r\delta_{0}$ (*)
2 |  $\dfrac{\partial\delta_{1}}{\partial t}+rS\dfrac{\partial\delta_{1}}{\partial S}+\kappa\left(\theta-v\right)\dfrac{\partial\delta_{1}}{\partial v}+\frac{1}{2}\left(vS^{2}\dfrac{\partial^{2}\delta_{1}}{\partial S^{2}}+\omega^{2}v\dfrac{\partial^{2}\delta_{1}}{\partial v^{2}}\right)+\rho\omega Sv\dfrac{\partial^{2}\delta_{1}}{\partial S\partial v}-r\delta_{1}$
$\vdots$ |  $\vdots$
N |  $\dfrac{\partial\delta_{N-1}}{\partial t}+rS\dfrac{\partial\delta_{N-1}}{\partial S}+\kappa\left(\theta-v\right)\dfrac{\partial\delta_{N-1}}{\partial v}+\frac{1}{2}\left(vS^{2}\dfrac{\partial^{2}\delta_{N-1}}{\partial S^{2}}+\omega^{2}v\dfrac{\partial^{2}\delta_{N-1}}{\partial v^{2}}\right)+\rho\omega Sv\dfrac{\partial^{2}\delta_{N-1}}{\partial S\partial v}-r\delta_{N-1}$

(*) Note that $\partial^{2}\delta_{0}/\partial v^{2}=0$.

####  3.1.1 Digression on optimal nuisance parameters

For the present case of the Heston model, KM suggest to choose either
$\eta_{0}=\sqrt{v(t)}$ or $\eta_{0}=\sqrt{\theta}$. As I will show in the next
sections, the specific choice of $\eta_{0}$ has no big influence on the result
of the approximation. However, setting $\eta_{0}=\sqrt{v(t)}$ appears to be
superior regarding accuracy compared to other choices. KM themselves also use
this specification in all their experiments based on the Heston model181818And
its generalization the CEV model., but do not further explain the reason for
this particular choice. However, since $v(t)$ denotes the current spot
volatility in the market, this choice also makes intuitive sense. Recall that
(15) showed clearly that the KM approximation mainly provides a correction
term if one e.g. would compute an option price by the Black-Scholes model when
the true market dynamics are determined by the Heston model. However, both
models require the definition of a "spot volatility" parameter. The only
difference is that in the Black-Scholes case this value is assumed to be a
constant over time, while in the Heston case this value is interpreted as the
initial value of the stochastic variance process. Since $v(t)$ denotes a
variance in (34), the choice $\eta_{0}=\sqrt{v(t)}$ implies that one would use
the same spot volatility estimate to compute Black-Scholes and Heston option
prices. Since such spot volatilities are usually estimated from observed
prices of options and their respective underlying this seems to be reasonable.
KM suggest that one might also use an optimal nuisance parameter, defined via
the optimization problem:

$\displaystyle\eta^{Opt}_{N}(S,v,t)=\text{arg}~{}\underset{\eta_{0}}{\text{min}}\left(C_{N}(S,t;\eta_{0})-C^{BS}(S,t,\eta_{0})\right)^{2}$
(47)

Note that this is equivalent to minimizing the square of the pricing error. As
KM point out, $\eta^{Opt}_{N}(S,v,t)$ convergences to the implied volatility
of the Black-Scholes model as $N\rightarrow\infty$. Consider the experiment
shown in Figure 1, there the blue line shows one possible sample path of the
variance process as defined in (34). The red line in the same figure shows
Black-Scholes implied volatilities, where the option prices had been computed
under the Heston model using the respective value $v(t)$ from the blue line in
each time step. Hence, the red line shows the implied volatilities one would
estimate at each point in time if working under the assumption of the Black-
Scholes model, while the true market dynamics are governed by the Heston
model. Figure 1 shows clearly that the spot volatility of the Heston model and
the implied volatility of the Black-Scholes model are highly correlated, where
the implied Black-Scholes volatility always slightly underestimates the true
volatility. However, this admittedly quite heuristic example shows that the
choice $\eta_{0}=\sqrt{v(t)}$ is almost equivalent to using the optimization
criterion in (47), explaining why this specific choice is superior to the
suggested alternative regarding accuracy of the approximation if a finite and
low number of corrective terms is used.

Figure 1: Simulated vs. B-S implied volatility

Notes: The variances series had been a single sample path generated by an
Euler discretazation of (48) & (49) with 500 time increments. Option prices
had been computed by Fourier inversion using numerical integration, implied
Black-Scholes volatilities by Matlabs bsimpl() function. Parameters:
$\kappa=2.00$, $\theta=0.04$, $\omega=0.10$, $\rho=-0.5$ and $r=0.10$.

###  3.2 Accuracy of the approximation

In order to investigate the accuracy of the approximations I will follow KM
and express approximation errors relative to the analytic solution as

$\displaystyle\%Diff=\frac{C_{Approximation}-C_{Analytic}}{C_{Analytic}}\cdot
100$

I will use this definition of the approximation error for all models
considered in this thesis and not only for the Heston model. If an analytic
solution is not available I use prices obtained via Monte Carlo Simulation
like analytic prices in the above expression.
Firstly I use the same parameters as KM did in their initial example, which
are also roughly the same as in Heston [1993]. Such that $\kappa=2.00$,
$\theta=0.04$, $\omega=0.10$, $\rho=-0.5$ and $r=0.10$. The spot volatility is
$v(t)=0.04$, the strike price is set to $100$ and time to maturity is one
year. In Figure 2(a) the Black-Scholes volatility is set equal to
$\eta_{0}=\sqrt{v(t)}$, such that the resulting figure replicates fig. 3 in
KM. For Figure 2(b) the nuisance parameter is set equal to
$\eta_{0}=\sqrt{\theta}$. In both cases the approximation is convergent,
achieving sufficient precession after only five corrective terms. As already
indicated before choosing the nuisance parameter as in Figure 2(a) appears to
be superior to the alternative. Note that setting $\eta_{0}=\sqrt{v(t)}$
implies for the initial corrective term $\delta_{0}=0$. Hence, in Figure 2(a)
the $N=0$ line effectively shows just the percentage difference between the
Black-Scholes and the Heston price. This is not the case in Figure 2(b), where
$\eta_{0}=\sqrt{\theta}\neq\sqrt{v(t)}$. Since in this case the approximation
looses accuracy, this loss seems to be related to the occurrence of the
$\delta_{0}$ term in the series expansion.

(a) Using $\sigma_{0}=\sqrt{v(t)}$

(b) Using $\sigma_{0}=\sqrt{\theta}$

Figure 2: KM approximation of the Heston [1993] model

Notes: Analytic prices had been computed via Fourier transforms. $N$ = numer
of corrective terms. Strike price = 100, time to maturity = 1 year,
$\kappa=2.00$, $\theta=0.04$, $\omega=0.10$, $\rho=-0.5$, $r=0.10$,
$v(t)=0.04$. All B-S parameters that coincide with the Heston [1993]
parameters are set equal.

Note that in both cases there is a range of stock prices for which the
approximation involving four corrective terms performs slightly better than
the one using five. However, the maximum absolute error is lowest for $N=4$.
Hence, I will use $\eta_{0}=\sqrt{v(t)}$ and $N=4$ throughout this section.
Figures 3(a) to 3(c) show the errors of the same approximation for different
maturities. These examples show an interesting pattern. If time to maturity is
increased the accuracy of all approximations except the $N=0$ approximation
reduces. For a very long maturity of four years the convergence even seems to
turn into divergence, since the $N=0$ approximation yields the most precise
result. Recall that the nuisance parameter of the approximation was set such
that the $N=0$ approximation consists only of the Black-Scholes price. When
analyzing volatility smiles in options data it is a well documented fact that
the smile is far more pronounced for short maturities than for long
maturities, i.e. the volatility smile tends to flatten with increasing
volatility. Since the Black-Scholes model predicts a completely flat "smile"
the true option price will be closer to the Black-Scholes price for long
maturities compared to short maturities. This suggests that KM’s approach,
which simply adds non-zero corrective terms to the Black-Scholes price in
Figure 3(c), gets increasingly imprecise with increasing time to maturities.
Hence, convergence of KM’s series expansion might not be uniform over time.
This issue had been also detected in KM for the approximation of prices of
long-run bonds. However, in case of equity or FX options such long maturities
are rare.

(a) Time to maturity 1.5 years

(b) Time to maturity 2.0 years

(c) Time to maturity 4.0 years

Figure 3: KM approximation of the Heston [1993] model

Analytic prices had been computed via Fourier transforms. $N$ = numer of
corrective terms. Strike price = 100, $\kappa=2.00$, $\theta=0.04$,
$\omega=0.10$, $\rho=-0.5$, $r=0.10$, $v(t)=0.04$. All B-S parameters that
coincide with the Heston [1993] parameters are set equal.

Figure 4 below shows how the accuarcy of the approximation changes when time
to maturity and moneyness are varied, whereas now also maturities below one
year are considered. As expected the performance of the approximation
increases for lower maturities. The pattern that the approximation of ITM (in-
the-money, moneyness > 1) is better than the approximation of OTM (out-of-the-
money, moneyness < 1) calls remains over all maturities. Interestingly for
very long maturities the approach overestimates the call price while for the
shortest maturity the price is underestimated. It should be noted that, while
for maturities below one year the approximation error of far-ITM calls
virtually becomes zero, it remains slightly above zero for all longer
maturities. Overall Figure 4 suggests a high accuracy of the approximation for
most maturities.

Figure 4: Effect of maturity and moneyness on KM’s expansion

Notes: Analytic prices had been computed via Fourier transforms. KM’s
approximation uses five corrective terms. $\kappa=2.00$, $\theta=0.04$,
$\omega=0.10$, $\rho=-0.5$, $r=0.10$, $v(t)=0.04$. All B-S parameters that
coincide with the Heston [1993] parameters are set equal.

Next I analyze the effect of correlation between the stock price and variance
process. In all previous figures the coefficient of correlation $\rho$ was set
equal to $-0.5$. Figure 5 shows KM’s approximation of the Heston model when
varying $\rho$ as well as the stock price simultaneously, whereas the number
of corrective terms is fixed to five (i.e. $N=4$). Panel (a) shows the result
for a maturity of one year, Panel (b) for a maturity of 0.5 year. Both plots
show that the approximation loses precision fast for OTM calls as correlation
increases. With increasing correlation accuracy of the approximation of OTM
calls in Figure 5(b) is significantly reduced. However, the precision of the
approximation of ITM calls appears to be almost unaffected by correlation.
Figure 5(a) clearly shows a lower overall accuracy of the approximation for a
longer time to maturity. There correlation always affects accuracy over the
range of different moneyness. Nevertheless, a similar pattern as in Figure
5(b) applies. For low levels of correlation the difference in accuracy between
ITM and OTM calls is far less pronounced than for high correlations. Overall,
Figures 4 and 5 show that KM’s approximation is strongly affected by
correlation, moneyness and time to maturity. Since Figure 5 also indicates a
connection between time to maturity and correlation, Figure 6 shows explicitly
the effect on the approximation error of KM’s expansion for an ATM call (at-
the-money, moneyness = 1), if these two quantities are varied. For a time to
maturity blow 6 months the level of correlation has no significant effect on
the accuracy. However, for longer times to maturity correlations between $0.0$
and $-0.5$ lead to an overestimation of the call price, while correlations
between $-0.5$ and $-1.0$ lead to an underestimation. The highest
approximation errors in either case are reached at the extreme ends of the
considered spectrum of correlation and the longest considered time to maturity
of two years.

(a) Time to maturity one year

(b) Time to maturity 0.5 year

Figure 5: Influence of correlation on KM’s approximation

Notes: Analytic prices had been computed via Fourier transforms. KM’s
approximation uses five corrective terms. Time to maturity = 1 year (Panel A)
or = 0.5 year (Panel B), $\kappa=2.00$, $\theta=0.04$, $\omega=0.10$,
$r=0.10$, $v(t)=0.04$. All B-S parameters that coincide with the Heston [1993]
parameters are set equal.

Figure 6: Effect of maturity and correlation on KM’s expansion

Notes: Analytic prices were obtained via Fourier inversion. KM’s approximation
uses five corrective terms. $\kappa=2.00$, $\theta=0.04$, $\omega=0.10$,
$r=0.10$, $v(t)=0.04$. All Black-Scholes parameters that coincide with the
Heston [1993] parameters are set equal.

All figures that appeared in this section so far had in common that the
accuracy of KM’s approximation had it’s maximum absolute error in the range of
OTM options. The errors increase fast the deeper the option is in the OTM
region and decreased as the option price moves in the ITM region. Hence, it
seems like the moneyness of the option has a significant influence on the
accuracy of the approximation. For short maturities and/or low levels of
correlation the approximation appears to be quite accurate over a wide range
of moneyness. However, if maturity and/or correlation increases precision
deteriorates strongest for OTM options.
KM also analyze their approximation of the Heston model for another set of
parameters. Specifically KM used the following values: Strike $=1000$,
$T-t=1/12$, $\kappa=0.1465$, $\theta=0.5172$, $\omega=0.5786$, $\rho=-0.0243$,
$r=0.00$, $v(t)=\theta$ and $\gamma=0.5$. These parameters had been estimated
by Bollerslev and Zhou [2002] in the context of FX options.191919As also
observed by Younesian [2013] this parameter choice violates the Feller
condition, since $4\kappa\theta/v=0.5860<2$. However, Cizek et al. [2011] note
that this commonly happens when calibrating the Heston model to market data in
the FX context. It merely implies that the variance process hits the zero
boundary recurrently but leaves it immediately, i.e. the time spent at the
boundary is zero. (See Cizek et al. [2011], p. 144). Hence, the violation of
the Feller condition should not affect the results of KM’s approximation. The
results of KM’s approximation can be seen in Table 2, which replicates the
results already shown in KM. Panel A shows how spot moneyness influences the
accuracy of the approximation, whereas Panel B shows the influence of
different values for the spot variance $v(t)$ for an ATM call. Again accuracy
of the approximation is very high in both Panels. The extreme high accuracy in
this case might also be related to the very low level of correlation and the
short time to maturity. The previous examples indicated that the both is
likely to improve the approximation.

Table 2: Approximating the Heston [1993] model

Panel A:
---
Stock price | Fourier Transform | KM approx. | %Diff
950 | 57.8425 | 57.8449 | 0.00418
960 | 62.3711 | 62.3738 | 0.0042574
970 | 67.1005 | 67.1033 | 0.0042447
980 | 72.0291 | 72.0321 | 0.0041553
990 | 77.1553 | 77.1584 | 0.0040021
1000 | 82.4766 | 82.4797 | 0.003797
1010 | 87.9903 | 87.9934 | 0.0035513
1020 | 93.6933 | 93.6964 | 0.003275
1030 | 99.5822 | 99.5852 | 0.0029773
1040 | 105.6532 | 105.656 | 0.0026663
1050 | 111.9021 | 111.9048 | 0.0023492

Panel B:
---
$v(t=0)$ | Fourier Transform | KM approx. | %Diff
0.1 | 36.4488 | 36.4854 | 0.10045
0.2 | 51.4125 | 51.4255 | 0.025319
0.3 | 62.8997 | 62.9068 | 0.011276
0.4 | 72.5792 | 72.5838 | 0.0063472
0.5 | 81.1007 | 81.104 | 0.0040628
0.6 | 88.7981 | 88.8006 | 0.002821
0.7 | 95.8702 | 95.8721 | 0.002072
0.8 | 102.4465 | 102.4481 | 0.0015857
0.9 | 108.6171 | 108.6184 | 0.0012524
1 | 114.4477 | 114.4488 | 0.001014
1.1 | 119.9878 | 119.9888 | 0.00083766

Comparsion of European Call option prices under the dynamics of the Heston
[1993] model. The Fourier transform prices had been obtained via numerical
integration using Matlab’s integral() function. KM’s approximation uses the
Black-Scholes model as baseline model. Panel A: Accuracy for differing spot
moneyness of the option. Strike price = 1,000. Panel B: At-the-money options
(S(t) = strike = 1,000) for differing spot variance. In both Panles: Time to
maturity $=1/12$, $\kappa=0.1465$, $\theta=0.5172$, $\omega=0.5786$,
$\rho=-0.0243$, $r=0.00$, $v(t)=\theta$ and $\gamma=0.5$. The volatility for
the Black-Scholes model is $\sigma_{0}=\sqrt{v(t)}$. All B-S parameters that
coincide with the Heston [1993] parameters are set equal.

####  3.2.1 Approximating Put Prices

Independent from the chosen parameter values the approximation always
performed comparatively poor for OTM calls and comparatively good for ITM
calls. In the following I will analyze whether this effect also appears in the
case of European Put options. KM do not provide an example for the
approximation of a put option. However, the series expansion of a put is
straightforward to derive. Recall the definitions of the corrective terms
$d(\cdot)$ and $\delta_{n}(\cdot)$ in the series expansion in (46). $d(\cdot)$
was defined as the difference between the final payoffs of the true model and
the baseline model. If one uses a Black-Scholes put option as baseline model
to approximate a European put option in the Heston framework, the term
$d(\cdot)$ remains equal to zero. Recall that it is necessary to choose a
baseline model with an identical final payoff, since the payoff function of a
plain vanilla option is non-differentiable. The other corrective term
$\delta_{n}$ accounts for the difference in the driving forces of the true and
the baseline market. Technically this correction is built on the difference of
the coefficients in the two model’s systems of SDEs and the derivatives of the
baseline model (see Table 1).
The initial pricing error and all of the higher order corrective terms given
in Table 1 remain unchanged if one switches from call to put options. In order
to see this recall the following: The coefficients are the same in either case
and the derivatives of the baseline model appearing in the corrective terms
are also the same. In the Black-Scholes model only the first derivative in
stock price direction differs between call and put, while the second
derivative is identical in both cases. In the initial pricing error
$\delta_{0}$ only the second derivative of the Black-Scholes formula appears
and all other derivatives of the baseline model that are used in the
approximation are taken from this second derivative. Hence, the whole series
of corrective terms $\delta_{n}$ is identical for a call and put option.
Hence, one can switch from call to put prices in the KM approximation simply
by using the Black-Scholes price of a Put option as the first element in the
series expansion in (46) and leaving the rest unchanged. This result can be
checked by applying the put-call parity relationship202020See e.g. Hull
[2015], p. 242. $C(S,t)+Ke^{-r(T-t)}=Put(S,t)+S_{0}$ to the call prices
approximated via KM’s series expansion, i.e. Put prices obtained from Call
prices through put-call parity should be identical to the Put prices obtained
through the series expansion using the Black-Scholes Put prices as its first
element. Applying both approaches I could confirm that this is indeed true.
Figure 7(a) shows the convergence of the KM approximation for a Put option.
Regarding the parameter values I again use the values from the first part of
the last section, i.e. $\kappa=2.00$, $\theta=0.04$, $\omega=0.10$,
$\rho=-0.5$, $r=0.10$, $v(t)=0.04$, with a strike of 100 and a time to
maturity of one year. Overall, each of the KM approximations show lower
absolute maximum errors than the approximation of the call price with the same
number of corrective terms. Each additional corrective term appears to
increase the accuracy of the approximation. However, as before there is a
range of option prices for which the expansion including four terms performs
slightly better than the one including five terms. The difference in accuracy
in the approximation of the OTM and ITM options is less clear in the case of a
put option. While the effect, that the precision of the approximation of OTM
options is higher than for ITM options, is present in all the approximations
except for the $N=2$ expansion, it is far less pronounced than for the Call
option. Figures 7(b) and 7(c) repeat the analysis for a time to maturity of
two and $0.25$ years respectively. As expected the accuracy increases
significantly for the short maturity and decreases for the long maturity. Also
the effect of better approximations of ITM options becomes more pronounced. As
before, in the long maturity case the convergence pattern starts to dissolve.

(a) Time to maturity one year

(b) Time to maturity two year

(c) Time to maturity 0.25 year

Figure 7: KM approximation of European Put option in the Heston [1993] model

Notes: Convergance behavior of the KM approximation for a European Put. True
model is Heston [1993], baseline model is the Black-Scholes model. Percentage
errors are in pct. of the analytic Heston price. Analytic prices had been
computed via Fourier transforms and numerical integration using Matlab’s
integral() function. $N$ = numer of corrective terms. Strike price = 100, time
to maturity = 1 year in Panel (a) and 0.25 year in Panel (b), $\kappa=2.00$,
$\theta=0.04$, $\omega=0.10$, $\rho=-0.5$, $r=0.10$, $v(t)=0.04$ and
$\gamma=0.5$. The nuisance parameter is set to $\sigma_{0}=\sqrt{v(t)}$.

####  3.2.2 Approximating Greeks in the Heston model

So far I analyzed the performance of KM’s series expansion in approximating
the price of European call and put options in the framework of the Heston
model. Next, I will assess the accuracy of the KM approximation for the
computation of hedge ratios in the Heston framework.
In order to compute the KM approximation of the hedge ratios I use the series
expansion in (20) and the finite difference approximation to compute the
derivatives therein. I calculate theses greeks using the same set of
parameters as for the computation of Table 2. This example is identical to the
one conducted in KM. However, KM do not explain whether they used numerical or
analytic approximations to obtain their series expansion values. By using
finite differences I am able to obtain the same values for all Greeks as
reported in KM for the case of varying monyness. For the case of varying spot
volatility, both my approximations as well as my analytic hedge ratios deviate
slightly from the values reported in KM. The results can been seen in Table 3
and 4. As I represent the options $\Delta$ in $\Gamma$ in percent the
approximation errors are reported in percentage points.
Table 5 shows the result when approximating $\mathscr{V}$ by KM’s series
expansion. The approximation errors are given as absolute deviations here. For
$\mathscr{V}$ I obtain the most significant deviation of my results from KM’s,
when the spot volatility is varied (see Panel B). This is somewhat surprising
since I obtained exactly KM’s results for varying moneyness, using the same
method of computation in both cases.
Note that in the case of $\mathscr{V}$ there is no corresponding Greek of the
baseline model. Even though the value of the Black-Scholes volatility
parameter is set equal to spot volatility parameter in the Heston model, the
variables are not identical. This was also emphasized when deriving the
initial pricing error in (45), where the corresponding derivatives of the
Black-Scholes model were considered to be equal to zero. Hence, when
considering the series expansion in (20) the first element representing the
derivative of the baseline model equals zero in the approximation of
$\mathscr{V}$.

Table 3: Approximating Heston’s $\Delta$

Panel A:
---
Stock price | Fourier Transform | KM approx. | Diff., pp.
950 | 44.2794 | 44.2819 | -0.0025864
960 | 46.2918 | 46.294 | -0.0021588
970 | 48.2928 | 48.2945 | -0.0016931
980 | 50.2776 | 50.2788 | -0.001201
990 | 52.2414 | 52.2421 | -0.00069418
1000 | 54.18 | 54.1801 | -0.00018398
1010 | 56.0893 | 56.089 | 0.00031874
1020 | 57.9657 | 57.9649 | 0.00080396
1030 | 59.8058 | 59.8046 | 0.0012626
1040 | 61.6066 | 61.6049 | 0.0016869
1050 | 63.3654 | 63.3633 | 0.00207

Panel B:
---
$v(t=0)$ | Fourier Transform | KM approx. | Diff., pp.
0.1 | 51.9512 | 51.9516 | -0.00044321
0.2 | 52.6614 | 52.6509 | 0.010469
0.3 | 53.2189 | 53.2149 | 0.0040575
0.4 | 53.6929 | 53.6917 | 0.0012105
0.5 | 54.1121 | 54.1121 | 0.0000348
0.6 | 54.492 | 54.4932 | -0.0012002
0.7 | 54.8416 | 54.8419 | -0.00021102
0.8 | 55.1673 | 55.1588 | 0.0085098
0.9 | 55.4732 | 55.4419 | 0.031343
1 | 55.7625 | 55.69 | 0.072555
1.1 | 56.0376 | 55.9066 | 0.13101

Notes: Comparsion of $\Delta_{S}$ for European Call option under the dynamics
of the Heston [1993] model. All Greek are in pct, differences in percentage
points. The Fourier transform prices had been obtained via numerical
integration using Matlab’s integral() function. Greeks for the KM approx. are
obtained via finite differences. Panel A: Differing spot moneyness of the
option. Strike price = 1,000. Panel B: At-the-money options (S(t) = strike =
1,000) for differing spot variance. In both Panles: Time to maturity $=1/12$,
$\kappa=0.1465$, $\theta=0.5172$, $\omega=0.5786$, $\rho=-0.0243$, $r=0.00$,
$v(t)=\theta$. The volatility for the Black-Scholes model is
$\sigma_{0}=\sqrt{v(t)}$. All B-S parameters that coincide with the Heston
[1993] parameters are set equal.

Table 4: Approximating Heston’s $\Gamma$

Panel A:
---
Stock price | Fourier Transform | KM approx. | Diff., pp.
950 | 0.20165 | 0.20161 | 4.0526e-05
960 | 0.20076 | 0.20071 | 4.4782e-05
970 | 0.19937 | 0.19932 | 4.8141e-05
980 | 0.1975 | 0.19745 | 5.0149e-05
990 | 0.19519 | 0.19514 | 5.0908e-05
1000 | 0.19246 | 0.19241 | 5.0701e-05
1010 | 0.18935 | 0.1893 | 4.9582e-05
1020 | 0.18588 | 0.18583 | 4.7187e-05
1030 | 0.18209 | 0.18205 | 4.4327e-05
1040 | 0.17802 | 0.17798 | 4.0407e-05
1050 | 0.1737 | 0.17366 | 3.6158e-05

Panel B:
---
$v(t=0)$ | Fourier Transform | KM approx. | Diff., pp.
0.1 | 0.44642 | 0.38533 | 0.061096
0.2 | 0.31234 | 0.30392 | 0.0084156
0.3 | 0.25395 | 0.25273 | 0.0012173
0.4 | 0.21935 | 0.21918 | 0.0001669
0.5 | 0.1958 | 0.19575 | 5.1236e-05
0.6 | 0.17844 | 0.17843 | 1.3833e-05
0.7 | 0.16496 | 0.1652 | -0.00023702
0.8 | 0.15409 | 0.15448 | -0.00038083
0.9 | 0.14509 | 0.1436 | 0.0014905
1 | 0.13748 | 0.1273 | 0.010179
1.1 | 0.13092 | 0.096136 | 0.034784

Notes: Comparsion of $\Gamma_{S}$ for European Call option under the dynamics
of the Heston [1993] model. All Greek are in pct, differences in percentage
points. The Fourier transform prices had been obtained via numerical
integration using Matlab’s integral() function. Greeks for the KM approx. are
obtained via finite differences. Panel A: Differing spot moneyness of the
option. Strike price = 1,000. Panel B: At-the-money options (S(t) = strike =
1,000) for differing spot variance. In both Panles: Time to maturity $=1/12$,
$\kappa=0.1465$, $\theta=0.5172$, $\omega=0.5786$, $\rho=-0.0243$, $r=0.00$,
$v(t)=\theta$. The volatility for the Black-Scholes model is
$\sigma_{0}=\sqrt{v(t)}$. All B-S parameters that coincide with the Heston
[1993] parameters are set equal.

Table 5: Approximating Heston’s $\mathscr{V}$

Panel A:
---
Stock price | Analytic $\mathscr{V}$ | KM approx. | Diff, abs.
950 | 74.9687 | 74.9679 | 0.00076581
960 | 76.221 | 76.2212 | -0.00023367
970 | 77.2834 | 77.2847 | -0.0013377
980 | 78.1538 | 78.1563 | -0.0025286
990 | 78.8316 | 78.8354 | -0.0037851
1000 | 79.3178 | 79.3229 | -0.005083
1010 | 79.6148 | 79.6212 | -0.0063961
1020 | 79.7259 | 79.7336 | -0.0076966
1030 | 79.6561 | 79.6651 | -0.0089566
1040 | 79.4111 | 79.4213 | -0.010148
1050 | 78.9977 | 79.0090 | -0.0112

Panel B:
---
$v(t=0)$ | Fourier Transform | KM approx. | Diff., abs.
0.1 | 180.4329 | 151.6085 | 28.8244
0.2 | 127.7884 | 122.9584 | 4.8299
0.3 | 104.3134 | 103.4803 | 0.83308
0.4 | 90.279 | 90.1519 | 0.12705
0.5 | 80.6826 | 80.6861 | -0.0035699
0.6 | 73.5874 | 73.5304 | 0.056988
0.7 | 68.0654 | 67.8674 | 0.19792
0.8 | 63.6084 | 63.6146 | -0.0061472
0.9 | 59.9122 | 61.4241 | -1.512
1 | 56.7816 | 62.6834 | -5.9018
1.1 | 54.0853 | 69.5145 | -15.4293

Notes: Comparsion of $\mathscr{V}$ for European Call option under the dynamics
of the Heston [1993] model. The Fourier transform prices had been obtained via
numerical integration using Matlab’s integral() function. Greeks for the KM
approx. are obtained via finite differences. Panel A: Differing spot moneyness
of the option. Strike price = 1,000. Panel B: At-the-money options (S(t) =
strike = 1,000) for differing spot variance. In both Panles: Time to maturity
$=1/12$, $\kappa=0.1465$, $\theta=0.5172$, $\omega=0.5786$, $\rho=-0.0243$,
$r=0.00$, $v(t)=\theta$. The volatility for the Black-Scholes model is
$\sigma_{0}=\sqrt{v(t)}$. All B-S parameters that coincide with the Heston
[1993] parameters are set equal.

##  4 Approximating the CEV model

In this section I will continue to follow the path of KM and apply their
series expansion method to another model of stochastic volatility. This also
had been done by KM themselves, such that I also in this section intend to
start my numerical analysis by replicating KM’s results. However, as in the
previous case of the Heston model, I will also consider some own examples.
The constant elasticity of variance (CEV) model212121I use the same
nomenclature as in KM. However, in the literature the term CEV is often used
to describe models in which variance is given by a deterministic function e.g.
of the stock price. was first introduced by Chan et al. [1992] to model
stochastic interest rates and had been further analyzed in the context of
stochastic volatility, among others, by Jones [2003] or Lewis [2000] and is
defined by the following system of SDEs

$\displaystyle dS(t)~{}=$ $\displaystyle~{}rSdt+\sqrt{v(t)}SdW_{1}(t)$ (48)
$\displaystyle dv(t)~{}=$
$\displaystyle~{}\kappa\left(\theta-v(t)\right)~{}dt+\omega|v(t)|^{\gamma}~{}dW_{2}(t)$
(49) $\displaystyle dW_{1}(t)dW_{2}(t)~{}=$ $\displaystyle~{}\rho dt$

Where $\kappa$ denotes the speed of mean-reversion parameter, $\theta$ and
$\omega$ the long-run variance and the diffusion coefficient respectively and
$r$ the constant instantaneous short-term rate. $\gamma$ denotes the CEV
parameter. The representation in (48) and (49) is the same as in KM. However,
usually the variance process includes $v(t)^{\gamma}$ instead of
$|v(t)|^{\gamma}$. KM do not point out their reasons to choose this version of
the CEV model, but it does not influence the series expansion in any way. The
CEV model appears to be a generalization of the Heston model. By setting
$\gamma=0.5$ one obtains the square root process in (34).
Despite the popularity of the Heston model, cases of the CEV model with
$\gamma\neq\frac{1}{2}$ are of empirical relevance as shown in Jones [2003].
While in the Heston model the volatility of instantaneous volatility is not
level dependent, the volatility of instantaneous of variance appears to be
level dependent.222222See Jones [2003] Jones [2003] concludes that, level
dependence appears to be a feature of the data and hence the CEV model might
perform better in replicating the observed volatility smiles. Hence, the
process in (49) with choices for $\gamma\neq 1/2$ might be a sensible
alternative to Heston’s specification. Indeed, by analyzing large samples of
daily return data from the S&P 500 index with 3537 observations, Jones [2003]
estimated values of $\gamma$ of $1.33$. As reported therein such values,
indicating a high elasticity of variance, are consistent with other findings
in empirical finance.232323See the discussion in Jones [2003], p. 196. Using a
more current data set on S&P 500 call index options Younesian [2013] estimates
a CEV parameter of approximately $0.6$. Note that whereas Jones [2003] relies
on Monte Carlo simulations and a Baysian estimation approach to obtain his
estimates, Younesian [2013] estimates the parameters directly from cross-
sections of option prices using non-linear least square estimation together
with KM’s closed-form approximation formula. However, Younesian [2013] notes
that the estimate of the CEV parameter might seriously flawed through
difficulties in the correct estimation of the correlation coefficient which he
linked to KM’s series approximation.
Since there are no analytic solutions available for any case of the CEV model
other than $\gamma=\frac{1}{2}$, I will rely on Monte Carlo simulations to
obtain reference values for the option prices estimated by the KM approach.
Hence, all approximation errors reported in this section can be understood as
percentage or absolute deviation from the price obtained via Monte Carlo (MC).

###  4.1 KM’s expansion for the CEV model

Since the CEV model is a simple generalization of the Heston model, it is
straightforward to derive KM’s series expansion for this model. One can simply
repeat the steps taken in the case of the Heston model, just replacing the
$\sqrt{v(t)}$ term by $v(t)^{\gamma}$ in all equations in section 3.1.
Following the approach of KM I again use the Black-Scholes model as baseline
model. As I again intend to approximate European call options under the CEV
model and also use calls in the Black-Scholes model, the corrective terms
adjusting the boundary are equal to zero.
Applying (14) yields an initial pricing error that is identical to the one
obtained for the Heston model

$\displaystyle\delta_{0}(S,t)=\dfrac{1}{2}\left(v(t)-\eta_{0}\right)S^{2}\dfrac{\partial^{2}C^{BS}}{\partial
S^{2}}$ (50)

Where (50) again results because $\partial C^{BS}/\partial
v=\partial^{2}C^{BS}/\partial v^{2}=\partial^{2}C^{BS}/\partial S\partial
v=0$. Below Table 6 shows how the series of pricing errors develop with the
order of the approximation, starting from (50) up to any positive integer $N$.
Comparing Table 6 with Table 1 shows that the only difference in the
corrective terms lies in the variance term $v(t)$, which is taken to the power
of $\gamma$ instead of $1/2$. The resulting series expansion is the same as in
(46), just using the corrective terms form Table 6 instead.
With respect to the assumptions underlying KM’s approximation it might be
interesting to investigate one property of the CEV process more carefully.
Specifically, if $\gamma>1$ then $v\rightarrow\infty$ and thus the growth
condition fails to hold.242424See Jones [2003], p. 187. Recalling the brief
summary of the assumptions on which KM’s approach is built on, the growth
condition was mentioned in the context of the Feynman-Kac representation of
the solution to the PDE which determines the difference between the true and
the baseline model (see section 2.2). However, this condition was only
necessary to ensure existence of a solution to the SDEs of the true model.
Hence, it would be sufficient to prove existence differently if the growth
condition is violated. In order to show the existence of a solution to the CEV
SDEs I follow closely the prove in Jones [2003, Appendix C].252525Jones [2003]
shows a proof of the existence of a solution but uses a different formulation
of the CEV process. Hence, I am basically redoing the steps in Jones [2003]
using the formulation of the CEV model as it is used in this thesis.
Note that the price process of the CEV model essentially is a geometric
Brownian motion. For the case of constant variance the existence of a solution
to this kind of SDE had been proven elsewhere in the literature. For the
present case of stochastic volatility Jones [2003] shows that a solution
exists if there exists a solution to the SDE determining the stochastic
behavior of the variance. Hence, in order to show that a solution to the CEV
model exists and thus the Feynman-Kac representation of the solution to the
PDE determining the pricing bias in KM’s approximation is well defined, it is
sufficient to show the existence of a solution to (49). To do so it is
sufficient to show that the scale measure of (49) is unattainable in 0 and
$+\infty$.262626See Jones [2003], p. 217 or the sources cited therein. Whereby
the scale measure of (49) is

$\displaystyle\Omega(v)$ $\displaystyle=\int_{m}^{n}\Theta(v)dv$ (51)
$\displaystyle\text{with }\Theta(v)$
$\displaystyle=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)$

Appendix A.6 includes a brief derivation of the above scale measure.
Jones [2003] states that the upper bound $+\infty$ is unattainable if
$\Omega(v)=\int_{m}^{+\infty}\Theta(v)dv=+\infty$. For $\gamma>1$ this is
true, since then $\frac{1}{v^{2\gamma-1}}\rightarrow 0$ and
$\frac{1}{v^{2\gamma-2}}\rightarrow 0$ such that $\Theta(v)\rightarrow 1$ as
$v\rightarrow\infty$. Hence, the integral fails to converge what causes
$+\infty$ to be unattainable.
Accordingly Jones [2003] denotes the lower bound 0 as unattainable if
$\Omega(v)=\int_{0}^{n}\Theta(v)dv=+\infty$. Note that for $\gamma>1$,
$\Theta(v)\approx
exp\left(\frac{2\kappa\theta}{\omega^{2}\left(2\gamma-1\right)}\frac{1}{v^{2\gamma-1}}\right)$
as $v\rightarrow 0$.272727See Jones [2003], p. 216. If additionally to
$\gamma>1$ also $\kappa>0$ and naturally $\theta>0$, then
$\Xi\equiv\frac{2\kappa\theta}{\omega^{2}\left(2\gamma-1\right)}>0$ holds
true. For small values of $v$ this implies the following relation282828See
Jones [2003], p. 216.

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

Since $1/v\rightarrow+\infty$ as $v\rightarrow 0$ the above inequality
indicates that also $\Theta(v)\rightarrow+\infty$ as $v\rightarrow 0$. Hence,
for the scale measure of the CEV variance process both 0 and $+\infty$ are not
attainable and thus a solution to (49) exists, such that the existence of a
solution to the SDEs determining the CEV model can be ensured even though the
growth condition is violated for $\gamma>1$. Hence, at least the first step of
KM’s approximation, i.e. the Feynman-Kac representation of the solution of the
PDE which the pricing error obeys, should be well defined.

Table 6: Iterations of the pricing error for the CEV model

n |  Pricing error $\delta_{n}(S,t)$
---|---
0 |  $\dfrac{1}{2}\left(v-\sigma_{0}\right)S^{2}\dfrac{\partial^{2}C^{BS}}{\partial S^{2}}$
1 |  $\dfrac{\partial\delta_{0}}{\partial t}+rS\dfrac{\partial\delta_{0}}{\partial S}+\kappa\left(\theta-v\right)\dfrac{\partial\delta_{0}}{\partial v}+\frac{1}{2}vS^{2}\dfrac{\partial^{2}\delta_{0}}{\partial S^{2}}+\rho\omega Sv(t)^{\gamma+\frac{1}{2}}\dfrac{\partial^{2}\delta_{0}}{\partial S\partial v}-r\delta_{0}$ (*)
2 |  $\dfrac{\partial\delta_{1}}{\partial t}+rS\dfrac{\partial\delta_{1}}{\partial S}+\kappa\left(\theta-v\right)\dfrac{\partial\delta_{1}}{\partial v}+\frac{1}{2}\left(vS^{2}\dfrac{\partial^{2}\delta_{1}}{\partial S^{2}}+\omega^{2}v^{2\gamma}\dfrac{\partial^{2}\delta_{1}}{\partial v^{2}}\right)+\rho\omega Sv^{\gamma+\frac{1}{2}}\dfrac{\partial^{2}\delta_{1}}{\partial S\partial v}-r\delta_{1}$
$\vdots$ |  $\vdots$
N |  $\dfrac{\partial\delta_{N-1}}{\partial t}+rS\dfrac{\partial\delta_{N-1}}{\partial S}+\kappa\left(\theta-v\right)\dfrac{\partial\delta_{N-1}}{\partial v}+\frac{1}{2}\left(vS^{2}\dfrac{\partial^{2}\delta_{N-1}}{\partial S^{2}}+\omega^{2}v^{2\gamma}\dfrac{\partial^{2}\delta_{N-1}}{\partial v^{2}}\right)+\rho\omega Sv^{\gamma+\frac{1}{2}}\dfrac{\partial^{2}\delta_{N-1}}{\partial S\partial v}-r\delta_{N-1}$

(*) Note that $\partial^{2}\delta_{0}/\partial v^{2}=0$.

###  4.2 Numerical accuracy

Figure 8 shows KM’s approximation of a European call in the CEV model for the
two CEV parameter choices $\gamma=0.6$ (Panel a) and $\gamma=1.33$ (Panel b).
Time to maturity is one year and the strike price is set to 100. The other
parameters are: $\theta=0.04$, $\omega=0.10$, $\rho=-0.5$, $r=0.1$ and the
spot variance $v=0.05$, time to maturity is one year and the strike price is
100. The Black-Scholes volatility is $\sqrt{v}$. As in the previous section,
$N$ denotes the order of the approximation. Since there is no analytic
solution to the CEV model for these two choices of the CEV parameter, the
reference values are obtained through Monte Carlo simulation. Both panels in
Figure 8 clearly show that KM’s series expansion converges to the Monte Carlo
solution. The approximation for $\gamma=1.33$ appears to be slightly more
precise than for $\gamma=0.6$. Overall the convergence pattern pretty much
resembles the pattern already observed for the Heston model. For both cases of
the CEV parameter the approximation improves with each new corrective term,
such that for both values of $\gamma$ clearly approximations of order $N=4$
should be used. The results indicate a slightly higher accuracy for
$\gamma=1.33$. However, the error in both panels seems to not completely die
out in the ITM region as this was the case for the Heston model.

(a) CEV with $\gamma=0.6$

(b) CEV with $\gamma=1.33$

Figure 8: Different CEV parameter values

Notes: Approximation of the CEV model. Panel (a): $\gamma=0.6$, Panel (b):
$\gamma=1.33$. In both panels: $\theta=0.04$, $\omega=0.10$, $\rho=-0.5$,
$r=0.1$ and the spot variance $v=0.05$, time to maturity is one year and the
strike price is 100. Black-Scholes volatility is $\sqrt{v}$. Reference values
are obtained via Monte Carlo simulation using the specification in Appendix
A.4.

KM analyze the accuracy of their approximation for the CEV model for
$\gamma=0.6$ and use the set of parameters estimated by Bollerslev and Zhou
[2002] for FX options. Table 7 below attempts to replicate parts of Kristensen
and Mele [2011, Table 2]. Time to maturity is one month, $\kappa=0.1465$,
$\theta=0.5172$, $\omega=0.5786$, $\rho=-0.0243$, $r=0.00$, and $v(t)=\theta$.
The strike of the Call option is set to 1,000. I also include 95% confidence
intervals in the table which are not reported in KM. Both prices, the MC as
well as the KM prices, differ from the ones reported therein. KM mention that
they design their MC simulation in a way such that the MC price obtained with
$\gamma=0.6$ can not deviate more than 0.5% from the analytic Heston price,
i.e. from the price obtained via Fourier inversion with the same set of
parameters but $\gamma=0.5$.292929See Kristensen and Mele [2011], p. 401.
Nevertheless, when comparing the option prices reported in Kristensen and Mele
[2011, Table 1] and Kristensen and Mele [2011, Table 2], it turns out that the
MC prices in their Table 2 deviate from the Fourier prices in their Table 1 by
more than 0.5%. Additionally, the procedure of fixing the MC prices to the
Fourier results might be questionable. Even though the parameter values are
close, there are structural consequence on the resulting model. As mentioned
already before, while the Heston model (i.e. $\gamma=0.5$) is affine, the
general CEV case with $\gamma=0.6$ yields a non-affine model. Further, the
barrier for the deviation of 0.5% seems arbitrary. Hence, I will not follow
KM’s suggestion to restrict the MC prices. Consequently, the prices based on
the MC simulation may deviate from the Heston prices in another way than
reported in KM. However, KM’s MC prices lie inside the confidence intervals of
my own MC results.
The results obtained via KM’s series expansion deviate from the values
reported in Kristensen and Mele [2011, Table 2]. The program which I employed
to obtain the approximation in Table 7 is the same as for Table 2. Since there
the results matched the values reported in KM exactly, the deviation is
surprising. While KM’s original results overestimate their MC prices
constantly by around 1%, my own results constantly underestimate my MC prices.
However, my MC and KM approximation prices differ by only about -0.55%.
Nevertheless, the conclusion would remain unchanged if I would use KM’s
original results. For the non-affine CEV case KM’s series expansion seems to
yield less precise approximations than it did for the affine Heston model.
Table 8 uses the same parameter values as Table 7, but now setting $\gamma$ to
$1.33$ in order to analyze the case of high elasticity of variance. The
difference to the MC price remains below 1%, but now the series expansion
constantly yields prices that are higher than the MC prices. The approximation
errors also show now an unexpected pattern. They increase when the option
moves from OTM to ITM. This seems to be related to specific parameter values
and the high elasticity of variance.

Table 7: Approximating the CEV model - $\gamma=0.6$

Panel A:
---
Stock price | Monte Carlo | Conf. Lower | Conf. Upper | KM approx. | % Diff
950 | 58.178 | 56.5888 | 59.7672 | 57.8674 | -0.53393
960 | 62.7283 | 61.0752 | 64.3814 | 62.3967 | -0.52869
970 | 67.4865 | 65.7693 | 69.2036 | 67.1266 | -0.5333
980 | 72.4486 | 70.6673 | 74.2298 | 72.0555 | -0.54255
990 | 77.6149 | 75.7696 | 79.4602 | 77.1817 | -0.55814
1000 | 82.9715 | 81.0622 | 84.8809 | 82.5029 | -0.56483
1010 | 88.5223 | 86.5491 | 90.4955 | 88.0163 | -0.57164
1020 | 94.2615 | 92.2246 | 96.2983 | 93.7188 | -0.57575
1030 | 100.1746 | 98.0743 | 102.2748 | 99.6069 | -0.56664
1040 | 106.2693 | 104.106 | 108.4326 | 105.677 | -0.55739
1050 | 112.5434 | 110.3175 | 114.7693 | 111.9249 | -0.54959

Panel B:
---
$v(t=0)$ | Monte Carlo | Conf. Lower | Conf. Upper | KM approx. | % Diff
0.1 | 36.7795 | 35.9868 | 37.5721 | 36.6167 | -0.44249
0.2 | 51.7656 | 50.6285 | 52.9026 | 51.5021 | -0.50903
0.3 | 63.2958 | 61.8824 | 64.7093 | 62.9573 | -0.53486
0.4 | 73.0207 | 71.3661 | 74.6754 | 72.6188 | -0.55043
0.5 | 81.5878 | 79.7144 | 83.4613 | 81.1286 | -0.56283
0.6 | 89.3302 | 87.2538 | 91.4066 | 88.8177 | -0.57367
0.7 | 96.4461 | 94.1786 | 98.7137 | 95.8836 | -0.58329
0.8 | 103.0649 | 100.6155 | 105.5143 | 102.455 | -0.59176
0.9 | 109.2765 | 106.6527 | 111.9003 | 108.6217 | -0.59925
1 | 115.1459 | 112.3539 | 117.9379 | 114.449 | -0.60525
1.1 | 120.7237 | 117.7688 | 123.6786 | 119.9864 | -0.61074

Notes: Comparsion of European Call option prices under the dynamics of the CEV
model with $\gamma=0.6$. MC prices are obtained using the Milstein scheme, 500
time steps and 20,000 sample paths. Confidence intervals are computed on the
95% level. KM’s approximation uses the Black-Scholes model as baseline model.
Panel A: Accuracy for differing spot moneyness of the option. Strike price =
1,000. Panel B: At-the-money options (S(t) = strike = 1,000) for differing
spot variance. In both Panles: Time to maturity $=1/12$, $\kappa=0.1465$,
$\theta=0.5172$, $\omega=0.5786$, $\rho=-0.0243$, $r=0.00$, and $v(t)=\theta$.
The volatility for the Black-Scholes model is $\sigma_{0}=\sqrt{v(t)}$.

Table 8: Approximating the CEV model - $\gamma=1.33$

Panel A:
---
Stock price | Monte Carlo | Conf. Lower | Conf. Upper | KM approx. | % Diff
950 | 57.7911 | 56.2022 | 59.3801 | 57.9685 | 0.30688
960 | 62.2887 | 60.6359 | 63.9416 | 62.4995 | 0.33837
970 | 66.974 | 65.257 | 68.6911 | 67.2303 | 0.38266
980 | 71.8629 | 70.0816 | 73.6442 | 72.1595 | 0.4128
990 | 76.9577 | 75.1121 | 78.8032 | 77.2853 | 0.42574
1000 | 82.2442 | 80.3345 | 84.1539 | 82.6053 | 0.43906
1010 | 87.7206 | 85.7469 | 89.6943 | 88.1168 | 0.45168
1020 | 93.3737 | 91.3361 | 95.4113 | 93.8168 | 0.47455
1030 | 99.2136 | 97.1124 | 101.3148 | 99.7018 | 0.49215
1040 | 105.2331 | 103.0687 | 107.3975 | 105.7682 | 0.50848
1050 | 111.4358 | 109.2087 | 113.6629 | 112.0119 | 0.51693

Panel B:
---
$v(t=0)$ | Monte Carlo | Conf. Lower | Conf. Upper | KM approx. | % Diff
0.1 | 36.6629 | 35.8701 | 37.4558 | 36.8541 | 0.52153
0.2 | 51.4396 | 50.3018 | 52.5775 | 51.6922 | 0.49096
0.3 | 62.818 | 61.4036 | 64.2324 | 63.1147 | 0.47232
0.4 | 72.4188 | 70.7634 | 74.0743 | 72.7493 | 0.45633
0.5 | 80.8778 | 79.0039 | 82.7517 | 81.235 | 0.44158
0.6 | 88.5227 | 86.4464 | 90.5991 | 88.9015 | 0.42796
0.7 | 95.5489 | 93.2821 | 97.8157 | 95.9457 | 0.41536
0.8 | 102.0848 | 99.6369 | 104.5326 | 102.4961 | 0.40293
0.9 | 108.2186 | 105.5973 | 110.8399 | 108.642 | 0.39128
1 | 114.0145 | 111.2261 | 116.8029 | 114.4488 | 0.38094
1.1 | 119.5232 | 116.573 | 122.4734 | 119.9658 | 0.37029

Notes: Comparsion of European Call option prices under the dynamics of the CEV
model with $\gamma=1.33$. MC prices are obtained using the Milstein scheme,
500 time steps and 20,000 sample paths. Confidence intervals are computed on
the 95% level. KM’s approximation uses the Black-Scholes model as baseline
model. Panel A: Accuracy for differing spot moneyness of the option. Strike
price = 1,000. Panel B: At-the-money options (S(t) = strike = 1,000) for
differing spot variance. In both Panles: Time to maturity $=1/12$,
$\kappa=0.1465$, $\theta=0.5172$, $\omega=0.5786$, $\rho=-0.0243$, $r=0.00$,
and $v(t)=\theta$. The volatility for the Black-Scholes model is
$\sigma_{0}=\sqrt{v(t)}$.

Focusing on the $\gamma=1.33$ case, Figure 9 shows for different maturities
that correlation seems to have no significant influence on the accuracy of the
approximation. All parameters are the same as for Figure 8(b), despite the
level correlation and the time to maturity. For the shorter time to maturity
of six months in Figure 9(a) KM’s approximation achieves exactly the same
accuracy for each level of correlation. For the longer time to maturity of one
year in Figure 9(b) the lines indicating the accuracy of the approximations
move slightly apart with an increasing level of correlation, however the
differences are extremely small. Hence, it seems that the effect of
correlation on the accuracy of KM’s approximation for the CEV model with
$\gamma=1.33$ is not significant. Since the level of correlation previously
had an effect on the accuracy when approximating the Heston model, this effect
now seems to be offset by the high elasticity of variance. Figure 10 performs
the same analyzes for a CEV parameter of $\gamma=0.6$, leaving all other
parameters equal. The result supports the view that effect of correlation on
the approximation is offset by a high elasticity of variance. For $\gamma=0.6$
correlation does have a visible effect, however the effect vanishes as the
option moves deeper in the ITM region. As expected the influence of
correlation is also increasing with increasing time to maturity.

(a) $T-t=0.5$

(b) $T-t=1.0$

Figure 9: CEV: Correlation and Moneyness for $\gamma=1.33$

Notes: Effect mooneyness and correlation on the KM approximation of the CEV
model, using five corrective terms. Panel (a): Time to maturity is three
months. Panel (b): Time to maturity is six months. $\gamma=1.33$,
$\theta=0.04$, $\omega=0.10$, $\rho=-0.5$, $r=0.1$ and the spot variance
$v=0.05$. Black-Scholes volatility is set to $\sqrt{v}$. Reference values are
obtained via Monte Carlo simulation using the specification in Appendix A.4.

(a) $T-t=0.5$

(b) $T-t=1.0$

Figure 10: CEV: Correlation and Moneyness for $\gamma=0.6$

Notes: Effect mooneyness and correlation on the KM approximation of the CEV
model, using five corrective terms. Panel (a): Time to maturity is three
months. Panel (b): Time to maturity is six months. $\gamma=0.6$,
$\theta=0.04$, $\omega=0.10$, $\rho=-0.5$, $r=0.1$ and the spot variance
$v=0.05$. Black-Scholes volatility is set to $\sqrt{v}$. Reference values are
obtained via Monte Carlo simulation using the specification in Appendix A.4.

Finally, it might be interesting to analyze the effect of $\gamma$ itself in a
more focused way. Figure 11 attempts to do this. The parameters are the same
as for Figure 8, but now letting $\gamma$ vary from $0.6$ to $1.6$. As in the
previous graphs I use an order of $N=4$ for the KM approximation. The figure
suggests that there is no significant influence of the elasticity parameter on
the accuracy of the approximation for $\gamma>0.6$. For all values the series
expansion shows the expected pattern of being less precise for OTM options.
Nevertheless, the absolute maximum difference to the MC price never exceeds
1.5%.

Figure 11: Level of elasticity of variance and moneyness

Notes: Effect of different elasticity of variance parameters and moneyness on
the KM approximation of the CEV model, using five corrective terms. Time to
maturity is one year, $\theta=0.04$, $\omega=0.10$, $\rho=-0.5$, $r=0.1$ and
the spot variance $v=0.05$. Black-Scholes volatility is set to $\sqrt{v}$.
Reference values are obtained via Monte Carlo simulation using the
specification in Appendix A.4.

##  5 Stochastic volatility with an OU - process

In this section I will change the focus to another class of stochastic
volatility models for option prices, in which stochastic volatility is modeled
directly and follows an Ornstein - Uhlenbeck process. Consider the following
system of SDEs

$\displaystyle dS(t)~{}=$ $\displaystyle~{}r~{}dt+\sigma(t)~{}dW_{1}(t)$ (53)
$\displaystyle d\sigma(t)~{}=$
$\displaystyle~{}\kappa\left(\theta-\sigma(t)\right)~{}dt+\omega~{}dW_{2}(t)$
(54) $\displaystyle dW_{1}(t)dW_{2}(t)~{}=$ $\displaystyle~{}\rho dt$

Where $\kappa$ denotes the speed of mean-reversion parameter, $\theta$ the
long-run mean and $\omega$ the instantaneous volatility of volatility. This
model had been originally suggested by Stein and Stein [1991] (henceforth S&S)
for the case of zero correlation between the stock price and the volatility
process, such that in the original formulation $\rho\equiv 0$. Since
correlation between the volatility and the stock price process plays an
important role in replicating volatility smiles observed in the data, zero
correlation seems to be rather restrictive. However, the model had been
extended by Schöbel and Zhu [1999](henceforth SZ) to also cover arbitrary
correlation between the two processes. I consider this model since it falls
outside of the class of affine models and hence yields a more common non-
affine model than the previous non-affine CEV example. Schöbel and Zhu [1999]
provide a (semi-)closed-form solution via Fourier inversion techniques for
this model, such that the accuracy and convergence behavior of KM’s
approximation can be checked against an analytic solution.
Note that the SZ model is not a special case of the Heston model or vice
versa. For the parameter restriction
$\kappa=\kappa_{Heston}/2,~{}\omega=\omega_{Heston}/2,~{}\theta=0$ and
$\theta_{Heston}=\omega^{2}/\kappa_{Heston}$ the Schöbel and Zhu [1999] model
is equivalent to the Heston model. However, these restrictions over-determine
the Heston model such that the two models are not consistent for a wide range
of parameter values.303030See Schöbel and Zhu [1999], p. 30. This provides
another advantage in assessing the robustness of KM’s approximation. The non-
affine and affine specifications of the CEV model shared the same general
structure of the variance process, with the one specification arising as a
special case of the other. While KM do not consider any non-affine option
pricing model other than the CEV case, they do so in the context of interest
rate models. Specifically, KM apply their approximation to the non-affine two-
factor model suggested by Forani and Mele [2006].313131See Kristensen and Mele
[2011], p. 407. In this context KM found that the approximation errors only
decreased with the number of corrective terms for maturities less than three
years and conclude that for longer maturities probably many corrective terms
are needed to achieve accurate approximations. Since the accuracy of KM’s
method also deteriorated with increasing time to maturity in the Heston, it
might be worthwhile to assess the accuracy of the approximation for different
times to maturity in the context of the Stein and Stein [1991] and Schöbel and
Zhu [1999] models.

###  5.1 Approximation of the SZ model

The steps to derive KM’s approximation for the SZ/S&S - model are essentially
the same as for the Heston model. Again the Black-Scholes formula serves as
the baseline model. The first step is to derive the covariance matrix of the
SZ - model.323232Since the SZ- and S&S- model are equivalent, except for the
choice of $\rho$ I will base the whole expansion on the SZ specification. An
approximation of the S&S-model is simply obtained by setting $\rho$ to zero.
By applying the same principle as in (36) in section 3.1 the covariance matrix
is given by

$\displaystyle\begin{bmatrix}\sigma(t)S(t)&0\\\
\rho\omega&\sqrt{1-\rho^{2}}\omega\end{bmatrix}\times\begin{bmatrix}\sigma(t)S(t)&\rho\omega\\\
0&\sqrt{1-\rho^{2}}\omega\end{bmatrix}=\begin{bmatrix}\sigma(t)^{2}S(t)^{2}&\rho\omega\sigma(t)S(t)\\\
\rho\omega\sigma(t)S(t)&\omega^{2}\end{bmatrix}$ (55)

Applying (14) with the right-hand side in the above matrix equation yields the
initial pricing error

$\displaystyle\delta_{0}=\frac{1}{2}\left(\sigma(t)^{2}-\eta_{0}^{2}\right)S(t)^{2}\dfrac{\partial^{2}C^{BS}}{\partial
S(t)^{2}}$ (56)

Note that $\delta_{0}$ is the same as in (50) by just replacing $v(t)$ with
$\sigma(t)^{2}$. Based on this initial pricing error all further corrective
terms can be computed iteratively according to the rule in (18). The iteration
of the corrective terms is shown in Table 9 below. Note that for the case of
the S&S-model, i.e. with $\rho=0$, all terms including the cross-derivatives
$\partial^{2}\delta_{n}/\partial S\partial\sigma(t)$ drop out. Using (46)
together with the corrective terms in Table 9 the price of a plain vanilla
call under the dynamics of the SZ-model can be expressed. For this model I
developed the expansion for the SZ-model until $N=5$, since this was the
optimal order of the series expansion used by KM for their approximation of
the Forani and Mele [2006] model.

Table 9: Iterations of the pricing error for the SZ/S&S-model n |  Pricing error $\delta_{n}(S,t)$
---|---
1 |  $\dfrac{\partial\delta_{0}}{\partial t}+rS\dfrac{\partial\delta_{0}}{\partial S}+\kappa\left(\theta-\sigma(t)\right)\dfrac{\partial\delta_{0}}{\partial\sigma(t)}+\frac{1}{2}\left(\sigma(t)^{2}S^{2}\dfrac{\partial^{2}\delta_{0}}{\partial S^{2}}+\omega^{2}\dfrac{\partial^{2}\delta_{0}}{\partial\sigma(t)^{2}}\right)+\rho\omega S\sigma(t)\dfrac{\partial^{2}\delta_{0}}{\partial S\partial\sigma(t)}-r\delta_{0}$
2 |  $\dfrac{\partial\delta_{1}}{\partial t}+rS\dfrac{\partial\delta_{1}}{\partial S}+\kappa\left(\theta-\sigma(t)\right)\dfrac{\partial\delta_{1}}{\partial\sigma(t)}+\frac{1}{2}\left(\sigma(t)^{2}S^{2}\dfrac{\partial^{2}\delta_{1}}{\partial S^{2}}+\omega^{2}\dfrac{\partial^{2}\delta_{1}}{\partial\sigma(t)^{2}}\right)+\rho\omega S\sigma(t)\dfrac{\partial^{2}\delta_{1}}{\partial S\partial\sigma(t)}-r\delta_{1}$
$\vdots$ |  $\vdots$
N |  $\dfrac{\partial\delta_{N-1}}{\partial t}+rS\dfrac{\partial\delta_{N-1}}{\partial S}+\kappa\left(\theta-\sigma(t)\right)\dfrac{\partial\delta_{N-1}}{\partial\sigma(t)}+\frac{1}{2}\left(\sigma(t)^{2}S^{2}\dfrac{\partial^{2}\delta_{N-1}}{\partial S^{2}}+\omega^{2}\dfrac{\partial^{2}\delta_{N-1}}{\partial\sigma(t)^{2}}\right)+\rho\omega S\sigma(t)\dfrac{\partial^{2}\delta_{N-1}}{\partial S\partial\sigma(t)}-r\delta_{N-1}$

###  5.2 Accuracy of the approximation

Figure 12 follows the same idea as the initial figures in the previous
sections, showing the behavior of the approximation errors for different stock
prices, when approximating the S&S-model, i.e. the zero correlation case of
the model. The parameter values are taken from Schöbel and Zhu [1999]. Figure
12(a), with a maturity of $0.25$ years, shows the expected convergence
pattern. For $N=4$ or $N=5$ corrective terms KM’s approximation yields quite
precise results with percentage errors not exceeding $0.6458\%$ for $N=4$ and
$1.1282\%$ for $N=5$. That the maximum absolute error for $N=4$ is slightly
less then for $N=5$ indicates that a higher number of corrective terms do not
necessarily for all models leads to greater precession. A feature that makes
it difficult to determine the optimal order of the approximation if an
analytic solution is indeed unknown.
With increasing time to maturity accuracy and also the convergence behavior
deteriorates fast. Figures 12(b) and 12(c) show the same approximation for
times to maturity of $0.5$ years and one year. While in Figure 12(b) no
obvious convergence pattern seems to be present, the last two Panels even seem
to hint towards divergence instead of convergences, at least for the
considered number of corrective terms. This indicates that the approximation
becomes instable much faster with increasing time to maturity when using an
Ornstein-Uhlenbeck instead of a square-root like process as in the previous
sections. There instability occurred for maturities well above one year only.
Figure 13 shows the result of approximating the more general SZ-model by KM’s
approximation. I assume a correlation coefficient of $\rho=-0.5$, while all
other parameters are the same as before. The convergence pattern appears to be
the same as for the S&S model. The approximation converges for a short time to
maturity of $0.25$ years as can be seen in Panel 13(a). However, the Panels
13(b) and 13(c) show that for times to maturity of half a year and one year
the convergence seems to turn into divergence. Note that even for the short
maturity case the approximation error for $N=5$ remains substantial for OTM
calls. The approximation error dies out fast and more uniformly with an
increasing stock price compared to the zero correlation case. Additionally, in
this case the approximation indeed unambiguously gets more precise with each
new corrective term, while in the zero correlation case for the same maturity
the sixth corrective term slightly reduced the accuracy. When time to maturity
increases, the approximation of the SZ models shows the same behavior as the
S&S model. For a maturity of half a year the convergence pattern starts to
dissolve and turns into divergence for a maturity of one year. Overall the
compatibility of KM’s approximation approach seems to be lower for this class
of models compared to the CEV or the Heston model. For the Heston model a
pattern of convergence remained for maturities of up to two years. Kimmel
[2008] finds that such divergence behavior for some maturities occurs if the
series expansions used for pricing has singularities at some complex
maturities. In such a case the series expansion will start to diverge for
maturities of the same absolute magnitude as the complex maturity at which the
singularity occurred.333333See Kimmel [2008], p. 1. More specifically Kimmel
[2008] argues that there exists a circle in the complex plane, which radius is
defined by distance from zero to singularity closest to zero, such the series
expansion converges if the maturity lies within the circle and diverges if it
lies outside. This suggests that for the SZ model maturities over 0.5 years
already lie outside this circle.
Note that both, Figure 12 and 13, show that the most precise KM expansion is
the one only including the first corrective terms, i.e. $N=0$, when time to
maturity is $0.5$ years or longer. For long maturities of one year these
expansions seem to be strikingly precise. However, this can simply be
explained by the specific choice of the nuisance parameter and the general
behavior of the volatility smile. By choosing $\eta_{0}=\sigma(t)$ one
effectively sets the first corrective term equal to zero, such that the $N=0$
expansion just uses the Black-Scholes price as an approximation. As I argued
already in the discussion of the approximation of the Heston model, for long
maturities the volatility smile tends to flatten and thus the Black-Scholes
price will be closer to the true price.

(a) Time to maturity $T-t=0.25$

(b) Time to maturity $T-t=0.50$

(c) Time to maturity $T-t=1.00$

Figure 12: Approximation of the Stein and Stein [1991] model.

Notes: For all Panels: Strike Price = 100, $\theta=0.2$, $\kappa=4.0$,
$\omega=0.1$, $r=0.0953$, spot volatility $\sigma=0.2$. Black-Scholes
volatility $\sigma_{0}=\sigma=0.2$. Analytic prices via Fourier inversion
solution of the SZ-model with $\rho=0$.

(a) Time to maturity $T-t=0.25$

(b) Time to maturity $T-t=0.50$

(c) Time to maturity $T-t=1.00$

Figure 13: Approximation of the Schöbel and Zhu [1999] model.

Notes: For all Panels: Strike Price = 100, $\theta=0.2$, $\kappa=4.0$,
$\omega=0.1$, $r=0.0953$, spot volatility $\sigma=0.2$, $\rho=-0.5$. Black-
Scholes volatility $\sigma_{0}=\sigma=0.2$. Analytic prices via Fourier
inversion.

It appears like the approximation only yields results with sufficient accuracy
if time to maturity is below six month. This is summarized by Figure 14, which
clearly shows that the accuracy of the approximation deteriorates fast with
increasing time to maturity. Hence, I will first focus on short maturities of
three month and one month to investigate the effect of different levels of
correlation on the performance of the approximation. Figures 12(a) and 13(a)
indicate significantly lower accuracy for OTM calls than for ITM calls when
correlation is at a level of $-0.5$.

Figure 14: Sensitivity towards increasing time to maturity

Notes: Instability of KM’s approximation of the SZ model for increasing time
to maturity and six corrective terms. Strike Price = 100, $\theta=0.2$,
$\kappa=4.0$, $\omega=0.1$, $r=0.0953$, spot volatility $\sigma=0.2$,
$\rho=-0.5$. Black-Scholes volatility $\sigma_{0}=\sigma=0.2$. Analytic prices
via Fourier inversion.

Figure 15 shows the behavior of the approximation for maturities of one month
and three months when using six corrective terms (i.e. $N=5$) for different
levels of correlation and moneyness. The figure indicates a behavior similar
as in the case of the Heston model (see Figure 5 for a comparison). For the
case of a time to maturity of one month in Figure 15(a) only levels
correlation of $-0.5$ to $-1$ have an influence on the accuracy of the
approximation, whereby this influence is also restricted to OTM options. ITM
options are approximated with extremely high accuracy for all levels of
correlation. If a longer time to maturity of three month, as in Figure 15(b),
is considered the influence of correlation increases significantly. For the
zero correlation case the approximation remains completely stable over the
whole moneyness range. An increasing level of correlation reduces the accuracy
of the approximation visibly, whereby the effect is now not only restricted to
OTM options but also appears for ITM options. Even though it is less
pronounced for the latter. Longer maturities are not considered for the SZ
model due to the instability of the approximation in these cases.

(a) Time to maturity equals one month

(b) Time to maturity equals three months

Figure 15: SZ-model with six corrective terms

Note: Analytic prices via Fourier transform. Strike $=~{}100$, $\kappa=4.0$,
$\theta=0.2$, $\omega=0.1$, $r=0.0953$, and spot volatility
$\sigma(t)=\theta$. Black-Scholes volatility equal spot volatility. All
approximations are of order $N=5$.

####  5.2.1 Approximating Greeks in the SZ model

In order to obtain KM’s approximation of hedge ratios for a Call option under
the dynamics of the SZ-model I again use equation (20). As in the case of the
Heston model I will also again use central finite differences to obtain
derivatives of the corrective terms. Analytic reference values can be obtained
via Fourier transforms following the same principle as for the Heston model,
which I outlined in appendix A.1.1.343434For a detailed description of the
Greeks in the SZ-model see Schöbel and Zhu [1999], p. 29.
Tables 10 and 11 show KM’s approximation of $\Delta$ and $\Gamma$ of a call
option under the dynamics of the SZ-model. As before the greeks are reported
in percent and thus the approximation error is given in percentage points. The
chosen parameters are the same as before and time to maturity is $0.25$. The
accuracy of the approximation is, as expected, lower as it was for the case of
Heston model. Table 12 shows the approximation of the options $\mathscr{V}$.
Approximation errors are reported as absolute deviations. Note that also here
the Black-Scholes $\mathscr{V}$ is equal to zero. Panel A in Table 12 shows
the approximation for varying moneyness. The approximation errors remain below
one percent and show magnitudes similar to the other approximations. However,
note that the spot volatility is set to $0.2$. Panel B of the same table
indicates that the approximation is highly sensitive to the choice of the spot
volatility parameter. The approximation yields accurate results for spot
volatilities of $0.1$ and $0.2$. With increasing spot volatility the
approximation of $\mathscr{V}$ becomes instable, yielding extremely high
negative values for the hedge ratio. Despite the time to maturity this reveals
another source of instability of the approximation procedure in case of the SZ
model.

Table 10: Approximating $\Delta$ for the SZ-model

Panel A:
---
Stock price | Fourier Transform, % | KM approx., % | Diff., pp.
80 | 1.5781 | 1.9609 | -0.38278
85 | 7.6057 | 8.6632 | -1.0575
90 | 22.6168 | 22.4057 | 0.21117
95 | 44.2531 | 42.1098 | 2.1433
100 | 64.6821 | 62.9238 | 1.7583
105 | 79.5248 | 79.2665 | 0.2583
110 | 88.791 | 89.2885 | -0.4975
115 | 94.0791 | 94.6592 | -0.58008
120 | 96.9379 | 97.473 | -0.53503

Panel B:
---
$\sigma(t=0)$ | Fourier Transform, % | KM approx., % | Diff., pp.