1 Introduction
We consider a discrete-time approximation for the price of European call option in the model of financial market with stochastic volatility driven by the Ornstein–Uhlenbeck process. An analytic expression for the price of the option is derived in [9]; however, the resulting formula is complicated and difficult to apply in most of available software. The discrete-time approximation is ready to be modeled even in the nonspecific software.
The problem of construction of discrete-time analogues for stochastic volatility models of financial markets is studied in a series of works including [5, 7, 2, 16, 1, 6, 18]. Various techniques are implemented, for example, multilevel Monte Carlo [5], conditional Monte Carlo [2, 18], exact simulation [2, 16], and Itô–Taylor approximations [7].
In most of the works, authors construct discrete-time approximations both for processes that describe the evolution of the price of asset and for processes driving the volatility of asset price. The model considered in this paper allows us to apply another approach: we only discretize the volatility process. The resulting discrete-time volatility process is then averaged in a special way and substituted into the option pricing formula. The option price is determined conditionally on the path of volatility process, and thus the conditional Monte Carlo approach is used. The rate of convergence of the option price calculated using the discrete-time volatility to the true option price for a given trajectory of volatility process is estimated.
Discretization of the model is naturally connected with the problem of discrete-time approximations to the solutions of stochastic differential equations. These matters are widely investigated and systematized in [8, 14, 17]. The simplest discrete-time approximation is the stochastic generalization of Euler approximation for deterministic differential equations proposed in [11], which is also referred to as the Euler–Maruyama scheme. Another suitable for implementation and effective method is the Milstein scheme [12]. Since the model under consideration is a diffusion with additive noise, both schemes coincide which is referred to below. It is worth noticing that Euler and Milstein schemes both belong to the class of Itô–Taylor approximations and have orders of convergence 0.5 and 1, respectively. For some diffusions, the approximation schemes can be enhanced to provide higher-order convergence, but this usually results in great increase in computation time.
Although exact simulation provides more precision compared to the Euler approximation, in this paper, we use the latter. This is motivated by the fact that the Euler approximation is cheaper in terms of computation time and by our desire to assess the rate of convergence of conditional option prices when the volatility is discretized using the Euler scheme.
This paper is structured as follows. We begin with the definition of the model under consideration and the discretization scheme used. In Section 3, the prices of the European call option are compared for discrete-time and continuous volatility processes to derive the estimate of strong convergence order. Section 4 provides numeric results of the simulation. In Section 5, we demonstrate the precision of discrete-time approximation for the case of deterministic volatility. Appendix A contains definitions and auxiliary results on discretization schemes and orders of their convergence mostly coming from [8].
2 The model and discrete approximation of volatility process
Let $\{\varOmega ,\mathcal{F},\mathbf{F}=\{{\mathcal{F}_{t}^{(B,Z)}},t\ge 0\},\mathbb{P}\}$ be a complete probability space with filtration generated by Wiener processes $\{B_{t},\hspace{0.1667em}Z_{t},\hspace{0.1667em}0\le t\le T\}$. We consider the model of the market where one risky asset is traded, its price evolves according to the geometric Brownian motion $\{S_{t},\hspace{0.2778em}0\le t\le T\},$ and its volatility is driven by a stochastic process. More precisely, the market is described by the pair of stochastic differential equations
We denote by $S_{0}=S$ and $Y_{0}=Y$ the deterministic initial values of the processes specified by Eqs. (1)–(2).
We further impose the following assumptions:
For example, the conditions mentioned in assumption (C2) are satisfied for the measurable function $\sigma (x)$ such that $c\le {\sigma }^{2}(x)\le C$ for $0<x<T$ and some constants $0<c<C$. Moreover, given the square integrability of $\sigma (Y_{s})$, the solution of differential equation (1) is given by
which yields that $S_{t}$ is continuous. Hence, the product $\sigma (Y_{s})S_{t}$ is square integrable: ${\int _{0}^{T}}{\sigma }^{2}(Y_{t}){S_{t}^{2}}dt<\infty $ a.s.
(3)
\[S_{t}=S_{0}\exp \Bigg(\mu t-\frac{1}{2}{\int _{0}^{t}}{\sigma }^{2}(Y_{s})ds+{\int _{0}^{t}}\sigma (Y_{s})dB_{s}\Bigg),\]The unique solution of the Langevin equation (2) $Y_{t}$ is the so called Ornstein–Uhlenbeck (OU) process. Its properties make it a suitable tool for modeling volatility in financial markets. One of the most important of the features is the mean-reversion property. The OU process is Gaussian with the following characteristics:
Moreover, the OU process is Markov and admits the explicit representation
Following [9], we proceed to the risk-neutral setting characterized by the minimal martingale measure $\mathbb{Q}$. With r being the interest rate, Eqs. (1)–(2) are now in the following form (see Section 5 in [9]):
where
(5)
\[\begin{array}{r@{\hskip0pt}l}\displaystyle dS_{t}& \displaystyle =rS_{t}dt+\sigma (Y_{t})S_{t}d{B_{t}^{\mathbb{Q}}},\\{} \displaystyle dY_{t}& \displaystyle =-\alpha Y_{t}dt+kd{Z_{t}^{\mathbb{Q}}},\end{array}\]
\[\begin{array}{r@{\hskip0pt}l}\displaystyle {B_{t}^{\mathbb{Q}}}& \displaystyle =B_{t}+{\int _{0}^{t}}\frac{\mu -r}{\sigma (Y_{s})}ds\hspace{1em}\text{and}\hspace{1em}{Z_{t}^{\mathbb{Q}}}=Z_{t}\end{array}\]
are independent Wiener processes w.r.t. $\mathbb{Q}$.This continuous-time model admits a variety of discrete-time approximations. In this paper, we apply the familiar Euler–Maruyama scheme, also referred to as the Euler scheme. The Euler–Maruyama approximation to the true solution of the Langevin equation (2) is the Markov chain ${Y}^{(m)}$ defined as follows:
The continuous-time process ${Y_{t}^{(m)}}$ is a step-type process defined by
where $[x]$ denotes an integer part of x.
-
• the partition of the interval $[0,T]$ into m equal subintervals of width $\Delta t=T/m$ is considered;
-
• the initial value of the scheme is set: ${Y_{0}^{(m)}}=Y_{0}$;
-
• ${Y_{l+1}^{(m)}}$, which we will use as a shorthand for ${Y_{(l+1)T/m}^{(m)}}$, $0\le l\le m-1$, is recursively defined by where $\Delta {Z_{l}^{\mathbb{Q}}}={Z_{(l+1)T/m}^{\mathbb{Q}}}-{Z_{lT/m}^{\mathbb{Q}}}$.
3 The price of European call option
The price of European call option V in the initial time moment of in model (5) is provided by
The inner expectation is conditional on the path of $Y_{s}$, $0\le s\le T$, and therefore, it actually is the Black–Scholes price for a model with deterministic time-dependent volatility. According to Lemma 2.1 in [13], the inner expectation in (7) has the following representation:
where $\bar{\sigma }:=\sqrt{\frac{1}{T}{\int _{0}^{T}}{\sigma }^{2}(Y_{s})ds}\ge 0$, Φ is the standard normal distribution function. The function $\bar{\sigma }$ may be viewed as the volatility averaged from the initial moment of time to maturity. The arguments of Φ in (8) are denoted as $d_{1}$ and $d_{2}$.
(8)
\[\begin{array}{r@{\hskip0pt}l}\displaystyle P& \displaystyle :={\mathbb{E}}^{\mathbb{Q}}\big\{{\big({S_{T}^{\mathbb{Q}}}-K\big)}^{+}\hspace{0.1667em}\big|\hspace{0.1667em}Y_{s},0\le s\le T\big\}={\operatorname{e}}^{\ln S+rT}\varPhi (d_{1})-K\varPhi (d_{2})\\{} & \displaystyle :={\operatorname{e}}^{\ln S+rT}\varPhi \bigg(\frac{\ln S+(r+\frac{1}{2}{\bar{\sigma }}^{2})T-\ln K}{\bar{\sigma }\sqrt{T}}\bigg)\\{} & \displaystyle \hspace{1em}-K\varPhi \bigg(\frac{\ln S+(r-\frac{1}{2}{\bar{\sigma }}^{2})T-\ln K}{\bar{\sigma }\sqrt{T}}\bigg),\end{array}\]Our aim is to estimate the error arising as a result of approximation of the exact formula (7) by application of the Euler approximation to the process that drives volatility. Thus, we need to assess the expectation of R given by
where m is the number of discretization points dividing the time interval $[0,T]$ into equal intervals, $\hat{P}_{m}$ denotes the price of the option in discrete setting calculated using the formula similar to (8):
where
with
where ${Y_{l}^{(m)}}$ is defined in (6).
(10)
\[\hat{P}_{m}={\operatorname{e}}^{\ln S+rT}\varPhi \big({d_{1}^{(m)}}\big)-K\varPhi \big({d_{2}^{(m)}}\big),\](13)
\[\bar{\sigma }_{m}=\sqrt{\frac{1}{T}\sum \limits_{l=1}^{m}{\sigma }^{2}\big({Y_{l}^{(m)}}\big)\frac{T}{m}}=\sqrt{\frac{1}{m}\sum \limits_{l=1}^{m}{\sigma }^{2}\big({Y_{l}^{(m)}}\big)},\]It is unlikely that we are able to find an exact or even approximate value for R. However, what really makes interest for investigation of the above bundle of models is the rate of convergence of the discrete setting to the continuous one. In order to assess the rate of convergence, the expression for an upper bound of R in terms of m needs to be derived.
Comparing (8) and (10), we can see that the approximation error arises solely due to the difference between $\bar{\sigma }$ and $\bar{\sigma }_{m}$. So, the first step would assessing the upper bound of expectation of absolute value of this difference w.r.t. m. After that, R might be expressed in terms of $R_{\sigma }:=\mathbb{E}|\bar{\sigma }-\bar{\sigma }_{m}|$.
Proof.
Since $\bar{\sigma }_{m}$ and $\bar{\sigma }$ are both square root functions, it is be more convenient to work with ${\bar{\sigma }_{m}^{2}}$ and ${\bar{\sigma }}^{2}$. To this end, we will use Hölder’s inequality:
\[\begin{array}{r@{\hskip0pt}l}\displaystyle \mathbb{E}|\bar{\sigma }_{m}-\bar{\sigma }|& \displaystyle =\mathbb{E}\Bigg|\sqrt{\frac{1}{T}{\int _{0}^{T}}{\sigma }^{2}(Y_{s})ds}-\sqrt{\frac{1}{m}\sum \limits_{i=1}^{m}{\sigma }^{2}\big({Y_{i}^{(m)}}\big)}\Bigg|\\{} & \displaystyle \le \mathbb{E}\Bigg|\sqrt{\Bigg|\frac{1}{T}{\int _{0}^{T}}{\sigma }^{2}(Y_{s})ds-\frac{1}{m}\sum \limits_{i=1}^{m}{\sigma }^{2}\big({Y_{i}^{(m)}}\big)\Bigg|}\Bigg|\\{} & \displaystyle \le {\Bigg(\mathbb{E}\Bigg|\frac{1}{T}{\int _{0}^{T}}{\sigma }^{2}(Y_{s})ds-\frac{1}{m}\sum \limits_{i=1}^{m}{\sigma }^{2}\big({Y_{i}^{(m)}}\big)\Bigg|\Bigg)}^{1/2}.\end{array}\]
Now we represent the integral as a sum of integrals over shorter intervals. Since the second summand does not depend on s, we may move it inside the integral sign, multiplying it by the inverse to the interval length:
\[\begin{array}{r@{\hskip0pt}l}\displaystyle \mathbb{E}|\bar{\sigma }_{m}-\bar{\sigma }|& \displaystyle \hspace{0.1667em}\le \hspace{0.1667em}{\Bigg(\mathbb{E}\Bigg|\sum \limits_{i=0}^{m-1}\Bigg(\frac{1}{T}{\int _{iT/m}^{(i+1)T/m}}{\sigma }^{2}(Y_{s})ds-\frac{1}{m}{\sigma }^{2}\big({Y_{i+1}^{(m)}}\big)\Bigg)\Bigg|\Bigg)}^{1/2}\\{} & \displaystyle \hspace{0.1667em}=\hspace{0.1667em}{\Bigg(\mathbb{E}\Bigg|\sum \limits_{i=0}^{m-1}\Bigg(\frac{1}{T}{\int _{iT/m}^{(i+1)T/m}}\hspace{-0.1667em}{\sigma }^{2}(Y_{s})ds\hspace{0.1667em}-\hspace{0.1667em}\frac{1}{m}\frac{m}{T}{\int _{iT/m}^{(i+1)T/m}}\hspace{-0.1667em}{\sigma }^{2}\big({Y_{i+1}^{(m)}}\big)ds\Bigg)\Bigg|\Bigg)}^{1/2}\\{} & \displaystyle \hspace{0.1667em}=\hspace{0.1667em}{\Bigg(\mathbb{E}\Bigg|\frac{1}{T}\sum \limits_{i=0}^{m-1}{\int _{iT/m}^{(i+1)T/m}}\big({\sigma }^{2}(Y_{s})-{\sigma }^{2}\big({Y_{i+1}^{(m)}}\big)\big)ds\Bigg|\Bigg)}^{1/2}.\end{array}\]
We apply the Hölder property of ${\sigma }^{2}(x)$:
\[\begin{array}{r@{\hskip0pt}l}& \displaystyle {\Bigg(\mathbb{E}\Bigg|\frac{1}{T}\sum \limits_{i=0}^{m-1}{\int _{iT/m}^{(i+1)T/m}}\big({\sigma }^{2}(Y_{s})-{\sigma }^{2}\big({Y_{i+1}^{(m)}}\big)\big)ds\Bigg|\Bigg)}^{1/2}\\{} & \displaystyle \hspace{1em}\le {\Bigg(\frac{L}{T}\mathbb{E}\Bigg(\sum \limits_{i=0}^{m-1}{\int _{iT/m}^{(i+1)T/m}}{\big|Y_{s}-{Y_{i+1}^{(m)}}\big|}^{\gamma }ds\Bigg)\Bigg)}^{1/2}\\{} & \displaystyle \hspace{1em}={\Bigg(\frac{L}{T}\sum \limits_{i=0}^{m-1}{\int _{iT/m}^{(i+1)T/m}}\mathbb{E}{\big|Y_{s}-{Y_{i+1}^{(m)}}\big|}^{\gamma }ds\Bigg)}^{1/2}.\end{array}\]
Recall that ${Y_{i}^{(m)}}$ is a shorthand for ${Y_{iT/m}^{(m)}}={Y_{s}^{(m)}}$, $s\in [iT/m,(i+1)T/m)$, and Proposition from the Appendix A yields that $\mathbb{E}|Y_{s}-{Y_{i+1}^{(m)}}|\le C_{1}{m}^{-1}$, where $C_{1}$ is some positive constant. We use Hölder’s inequality to derive that $\mathbb{E}|Y_{s}-{Y_{i+1}^{(m)}}{|}^{\gamma }\le {C_{1}^{\gamma }}{m}^{-\gamma }$ and arrive at
\[\mathbb{E}|\bar{\sigma }_{m}-\bar{\sigma }|\le {\bigg(\frac{L}{T}m\frac{T}{m}{C_{1}^{\gamma }}{m}^{-\gamma }\bigg)}^{1/2}=C{m}^{-\gamma /2}\]
for $C:=\sqrt{L{C_{1}^{\gamma }}}$, which proves the lemma. □The above lemma enables us to prove the main result of this section.
Theorem 3.1.
Let ${\sigma }^{2}(x)$ satisfy Hölder condition (14). Then $\mathbb{E}R\le D{m}^{-\gamma /2}$, where D is some positive constant.
Proof.
The function $\varPhi (x)$ has a continuous bounded derivative on $\mathbb{R}$; hence, we can use its Lipschitz property:
where c is a positive constant, and the last inequality is due to the assumption that $\sigma (x)$ is bounded away from zero for any $x\in \mathbb{R}$ (see assumption (C2)). Hence, using Lemma 3.1, we get
\[\begin{array}{r@{\hskip0pt}l}\displaystyle \mathbb{E}R& \displaystyle =\mathbb{E}|P-\hat{P}_{m}|\\{} & \displaystyle \le \mathbb{E}\big(S{\operatorname{e}}^{rT}\big|\varPhi (d_{1})-\varPhi \big({d_{1}^{(m)}}\big)\big|+K\big|\varPhi (d_{2})-\varPhi \big({d_{2}^{(m)}}\big)\big|\big)\\{} & \displaystyle \le \mathbb{E}\Big(\underset{x}{\sup }\big|f(x)\big|\big(S{\operatorname{e}}^{rT}\big|d_{1}-{d_{1}^{(m)}}\big|+K\big|d_{2}-{d_{2}^{(m)}}\big|\big)\Big),\end{array}\]
where $f(x)=\frac{1}{\sqrt{2\pi }}{\operatorname{e}}^{-\frac{{x}^{2}}{2}}$ is the density of the standard normal distribution. In the above representation, (15)
\[\begin{array}{r@{\hskip0pt}l}\displaystyle \big|d_{1}-{d_{1}^{(m)}}\big|& \displaystyle =\bigg|\bigg(\frac{1}{\bar{\sigma }}-\frac{1}{\bar{\sigma }_{m}}\bigg)\frac{\ln S-\ln K+rT}{\sqrt{T}}+\frac{1}{2}\sqrt{T}(\bar{\sigma }-\bar{\sigma }_{m})\bigg|\end{array}\]
\[\mathbb{E}\big|d_{1}-{d_{1}^{(m)}}\big|\le C_{1}\mathbb{E}|\bar{\sigma }-\bar{\sigma }_{m}|\le D_{1}{m}^{-\gamma /2},\]
where $C_{1}:=|\frac{\ln (S/K)+rT}{{c}^{2}\sqrt{T}}+\frac{\sqrt{T}}{2}|$ and $D_{1}$ are positive constants.Similarly, $\mathbb{E}|d_{2}-{d_{2}^{(m)}}|\le D_{2}{m}^{-\gamma /2}$, $D_{2}=const>0$, and we arrive at
for a positive constant D.
The theorem is proved. □
4 Numeric examples
Theorem 6.1 in [9] provides an analytic representation for the price of European call option for the stochastic volatility model under consideration. However, using it to calculate the price of an option is rather difficult and time-consuming. We further present the results of calculation of the price of European call option using simulation techniques.
The calculation process is performed in Matlab 7.9.0 and is structured as follows:
-
1. The choice of discrete ranges of values of input parameters;
-
2. The choice of the function $\sigma (Y_{s})$;
-
3. For each combination of input parameters we generate 1000 trajectories of an Ornstein–Uhlenbeck process by splitting the time interval into subintervals of length $\Delta t=0.001$ and modeling values of the OU process at these points (that is, generating normally distributed variables with known mean and standard deviation using relationship (6)). For each trajectory, (10) is applied to calculate ${\bar{\sigma }_{m}^{2}}$ and the price of an option. The results for all trajectories are then averaged and discounted to provide the sample average of the price denoted by $\hat{\mathbb{E}}\hat{P}_{m}$. The average volatility over all trajectories and time interval is denoted by $\hat{\mathbb{E}}{\bar{\sigma }_{m}^{2}}$.
Table 1.
${\sigma }^{2}(Y_{s})=a|Y_{s}|+b$
| T | k | r | K | a | b | $\hat{\mathbb{E}}{\bar{\sigma }_{m}^{2}}$ | $\hat{\mathbb{E}}\hat{P}_{m}$ | $\hat{\mathbb{E}}{\bar{\sigma }_{m}^{2}}$ | $\hat{\mathbb{E}}\hat{P}_{m}$ |
| $\alpha =1$ | $\alpha =100$ | ||||||||
| 0.25 | 0.1 | 0 | 0.8 | 1 | 0 | 0.088 | 0.204 | 0.009 | 0.200 |
| 0.5 | 0.1 | 0 | 0.8 | 1 | 0 | 0.082 | 0.213 | 0.007 | 0.200 |
| 1 | 0.1 | 0 | 0.8 | 1 | 0 | 0.073 | 0.227 | 0.007 | 0.200 |
| 0.25 | 0.5 | 0 | 0.8 | 1 | 0 | 0.147 | 0.211 | 0.031 | 0.200 |
| 0.5 | 0.5 | 0 | 0.8 | 1 | 0 | 0.185 | 0.235 | 0.030 | 0.201 |
| 1 | 0.5 | 0 | 0.8 | 1 | 0 | 0.216 | 0.280 | 0.029 | 0.207 |
| 0.25 | 1 | 0 | 0.8 | 1 | 0 | 0.264 | 0.224 | 0.059 | 0.201 |
| 0.5 | 1 | 0 | 0.8 | 1 | 0 | 0.338 | 0.264 | 0.058 | 0.207 |
| 1 | 1 | 0 | 0.8 | 1 | 0 | 0.412 | 0.334 | 0.058 | 0.221 |
| 0.25 | 0.1 | 0.01 | 1 | 1 | 0.2 | 0.289 | 0.108 | 0.209 | 0.092 |
| 0.5 | 0.1 | 0.01 | 1 | 1 | 0.2 | 0.281 | 0.151 | 0.207 | 0.130 |
| 1 | 0.1 | 0.01 | 1 | 1 | 0.2 | 0.273 | 0.210 | 0.207 | 0.184 |
| 0.25 | 0.5 | 0.01 | 1 | 1 | 0.2 | 0.346 | 0.117 | 0.231 | 0.097 |
| 0.5 | 0.5 | 0.01 | 1 | 1 | 0.2 | 0.375 | 0.172 | 0.230 | 0.137 |
| 1 | 0.5 | 0.01 | 1 | 1 | 0.2 | 0.414 | 0.254 | 0.229 | 0.193 |
| 0.25 | 1 | 0.01 | 1 | 1 | 0.2 | 0.459 | 0.134 | 0.259 | 0.102 |
| 0.5 | 1 | 0.01 | 1 | 1 | 0.2 | 0.532 | 0.203 | 0.258 | 0.145 |
| 1 | 1 | 0.01 | 1 | 1 | 0.2 | 0.617 | 0.305 | 0.258 | 0.204 |
| 0.25 | 0.1 | 0.02 | 1.2 | 1 | 1 | 1.089 | 0.141 | 1.009 | 0.134 |
| 0.5 | 0.1 | 0.02 | 1.2 | 1 | 1 | 1.079 | 0.228 | 1.007 | 0.218 |
| 1 | 0.1 | 0.02 | 1.2 | 1 | 1 | 1.073 | 0.347 | 1.007 | 0.335 |
| 0.25 | 0.5 | 0.02 | 1.2 | 1 | 1 | 1.148 | 0.147 | 1.031 | 0.136 |
| 0.5 | 0.5 | 0.02 | 1.2 | 1 | 1 | 1.178 | 0.240 | 1.030 | 0.221 |
| 1 | 0.5 | 0.02 | 1.2 | 1 | 1 | 1.216 | 0.371 | 1.029 | 0.339 |
| 0.25 | 1 | 0.02 | 1.2 | 1 | 1 | 1.262 | 0.157 | 1.059 | 0.138 |
| 0.5 | 1 | 0.02 | 1.2 | 1 | 1 | 1.341 | 0.260 | 1.058 | 0.225 |
| 1 | 1 | 0.02 | 1.2 | 1 | 1 | 1.414 | 0.402 | 1.058 | 0.344 |
To begin with, let us recall the notation of input parameters along with ranges of values assigned to them in the process of simulation:
T – time to maturity, $T=0.25;\hspace{2.5pt}0.5;\hspace{2.5pt}1$;
k – volatility of OU process, $k=0.1;\hspace{2.5pt}0.5;\hspace{2.5pt}1$;
α – mean-reversion rate, $\alpha =1;\hspace{2.5pt}100$;
r – interest rate, $r=0;\hspace{2.5pt}0.01;\hspace{2.5pt}0.02$;
K – strike price, $K=0.8;\hspace{2.5pt}1;\hspace{2.5pt}1.2$;
$S_{0}$ – initial price of stock, $S_{0}=1$;
$Y_{0}$ – initial value of OU process, $Y_{0}=0.1$.
In order to produce numerical results, we choose the following options for the function $\sigma (Y_{s})$:
The results of simulations are split into groups by the mean-reversion rate α and function $\sigma (Y_{s})$. Meaningless and uninteresting results provided by some distinct combinations of inputs are ignored.
Mean-reversion of 1 corresponds to slow reverting models, and fast mean-reverting models are characterized by $\alpha =100$. Matters of speed of mean-reversion are addressed, for example, in [4].
Table 2.
${\sigma }^{2}(Y_{s})={\operatorname{e}}^{Y_{s}}+c$
| T | k | r | K | $\hat{\mathbb{E}}{\bar{\sigma }_{m}^{2}}$ | $\hat{\mathbb{E}}\hat{P}_{m}$ | $\hat{\mathbb{E}}{\bar{\sigma }_{m}^{2}}$ | $\hat{\mathbb{E}}\hat{P}_{m}$ |
| $\alpha =1$ | $\alpha =100$ | ||||||
| 0.25 | 0.1 | 0 | 0.8 | 1.113 | 0.303 | 1.024 | 0.297 |
| 0.5 | 0.1 | 0 | 0.8 | 1.103 | 0.372 | 1.022 | 0.363 |
| 1 | 0.1 | 0 | 0.8 | 1.088 | 0.465 | 1.021 | 0.456 |
| 0.25 | 0.5 | 0 | 0.8 | 1.135 | 0.305 | 1.025 | 0.297 |
| 0.5 | 0.5 | 0 | 0.8 | 1.131 | 0.374 | 1.023 | 0.363 |
| 1 | 0.5 | 0 | 0.8 | 1.119 | 0.468 | 1.022 | 0.456 |
| 0.25 | 1 | 0 | 0.8 | 1.184 | 0.307 | 1.027 | 0.297 |
| 0.5 | 1 | 0 | 0.8 | 1.212 | 0.380 | 1.025 | 0.363 |
| 1 | 1 | 0 | 0.8 | 1.238 | 0.478 | 1.024 | 0.456 |
| 0.25 | 0.1 | 0.01 | 1 | 1.112 | 0.209 | 1.024 | 0.201 |
| 0.5 | 0.1 | 0.01 | 1 | 1.103 | 0.291 | 1.022 | 0.281 |
| 1 | 0.1 | 0.01 | 1 | 1.086 | 0.401 | 1.021 | 0.390 |
| 0.25 | 0.5 | 0.01 | 1 | 1.121 | 0.209 | 1.025 | 0.201 |
| 0.5 | 0.5 | 0.01 | 1 | 1.129 | 0.294 | 1.023 | 0.281 |
| 1 | 0.5 | 0.01 | 1 | 1.128 | 0.405 | 1.022 | 0.390 |
| 0.25 | 1 | 0.01 | 1 | 1.178 | 0.213 | 1.026 | 0.201 |
| 0.5 | 1 | 0.01 | 1 | 1.206 | 0.299 | 1.025 | 0.281 |
| 1 | 1 | 0.01 | 1 | 1.216 | 0.412 | 1.023 | 0.390 |
| 0.25 | 0.1 | 0.02 | 1.2 | 1.110 | 0.143 | 1.024 | 0.135 |
| 0.5 | 0.1 | 0.02 | 1.2 | 1.103 | 0.231 | 1.022 | 0.220 |
| 1 | 0.1 | 0.02 | 1.2 | 1.087 | 0.349 | 1.021 | 0.338 |
| 0.25 | 0.5 | 0.02 | 1.2 | 1.133 | 0.145 | 1.025 | 0.135 |
| 0.5 | 0.5 | 0.02 | 1.2 | 1.128 | 0.233 | 1.023 | 0.220 |
| 1 | 0.5 | 0.02 | 1.2 | 1.115 | 0.352 | 1.021 | 0.338 |
| 0.25 | 1 | 0.02 | 1.2 | 1.162 | 0.147 | 1.027 | 0.135 |
| 0.5 | 1 | 0.02 | 1.2 | 1.201 | 0.239 | 1.025 | 0.220 |
| 1 | 1 | 0.02 | 1.2 | 1.255 | 0.367 | 1.023 | 0.338 |
We may observe that, under faster mean-reversion, the average volatility $\hat{\mathbb{E}}{\bar{\sigma }_{m}^{2}}$ and, consequently, the price of the option are lower, which is exactly what is expected from the model.
Tables 3 and 4 illustrate how the price of the option changes with the decrease of time step in discrete model.
Table 3.
${\sigma }^{2}(Y_{s})=|Y_{s}|+0.2$, $K=1$, $r=0.02$, $k=0.1$, $T=1$. Convergence
| $\Delta t$ | α | $\hat{\mathbb{E}}{\bar{\sigma }_{m}^{2}}$ | ${\bar{d}_{1}^{(m)}}$ | ${\bar{d}_{2}^{(m)}}$ | $\hat{\mathbb{E}}\hat{P}_{m}$ |
| ${10}^{-2}$ | 1 | 0.272367 | 0.299043 | −0.222056 | 0.213552 |
| ${10}^{-3}$ | 1 | 0.271043 | 0.298506 | −0.221338 | 0.213073 |
| ${10}^{-4}$ | 1 | 0.272534 | 0.299123 | −0.222179 | 0.213631 |
| ${10}^{-5}$ | 1 | 0.271837 | 0.298822 | −0.221753 | 0.213351 |
| ${10}^{-6}$ | 1 | 0.271421 | 0.298667 | −0.221560 | 0.213220 |
| ${10}^{-2}$ | 100 | 0.208910 | 0.272291 | −0.184776 | 0.189047 |
| ${10}^{-3}$ | 100 | 0.206599 | 0.271267 | −0.183264 | 0.188073 |
| ${10}^{-4}$ | 100 | 0.206439 | 0.271196 | −0.183159 | 0.188005 |
| ${10}^{-5}$ | 100 | 0.206413 | 0.271184 | −0.183142 | 0.187994 |
| ${10}^{-6}$ | 100 | 0.206443 | 0.271198 | −0.183162 | 0.188007 |
Table 4.
${\sigma }^{2}(Y_{s})={\operatorname{e}}^{Y_{s}}+0.2$, $K=1$, $r=0.02$, $k=0.1$, $T=1$. Convergence
| $\Delta t$ | α | $\hat{\mathbb{E}}{\bar{\sigma }_{m}^{2}}$ | ${\bar{d}_{1}^{(m)}}$ | ${\bar{d}_{2}^{(m)}}$ | $\hat{\mathbb{E}}\hat{P}_{m}$ |
| ${10}^{-2}$ | 1 | 1.265414 | 0.556279 | −0.519082 | 0.431865 |
| ${10}^{-3}$ | 1 | 1.269504 | 0.579243 | −0.543620 | 0.432472 |
| ${10}^{-4}$ | 1 | 1.266274 | 0.584925 | −0.549670 | 0.431990 |
| ${10}^{-5}$ | 1 | 1.266030 | 0.566934 | −0.530485 | 0.431948 |
| ${10}^{-6}$ | 1 | 1.265635 | 0.576169 | −0.540343 | 0.431892 |
| ${10}^{-2}$ | 100 | 1.201083 | 0.566092 | −0.529585 | 0.422128 |
| ${10}^{-3}$ | 100 | 1.201047 | 0.566500 | −0.530021 | 0.422123 |
| ${10}^{-4}$ | 100 | 1.201026 | 0.566203 | −0.529703 | 0.422120 |
| ${10}^{-5}$ | 100 | 1.201036 | 0.566693 | −0.530228 | 0.422121 |
| ${10}^{-6}$ | 100 | 1.201023 | 0.566052 | −0.529542 | 0.422119 |
In view of Section 3, it is also of certain interest to compare calculations obtained over one trajectory but under different discretization steps. We constructed 2000 trajectories with time-step size of ${10}^{-6}$: 1000 for the case $\alpha =1$ and 1000 for the case $\alpha =100$. These trajectories are considered to be “true” continuous-time trajectories of the Ornstein–Uhlenbeck process $Y_{t}$. The corresponding values of ${\bar{\sigma }_{m}^{2}}$ are considered to be “true” continuous-time values of ${\bar{\sigma }}^{2}$. The calculations were then performed for wider discretization intervals using the points of constructed trajectories. Thus, the samples of discretization errors for ${\bar{\sigma }_{m}^{2}}$ were derived. Probably, the estimate of ${\bar{\sigma }_{m}^{2}}$ is more valuable in such context since one would not usually calculate the price of an option over one trajectory. However, the estimate of volatility is usually derived from past data, which is in essence one distinct realization of the space of all possible scenarios.
Tables 5 and 6 provide characteristics of the samples of discretization errors. Errors are measured as a percentage of the “true” value.
Table 5.
${\sigma }^{2}(Y_{s})=|Y_{s}|+0.2$, $K=1$, $r=0.02$, $k=0.1$, $T=1$. Characteristics of sample of errors
| ${10}^{-2}$ | ${10}^{-3}$ | ${10}^{-4}$ | ${10}^{-5}$ | |
| $\alpha =1$ | ||||
| Average | 0.08710% | 0.00834% | 0.00081% | 0.00008% |
| St. error | 0.0000427 | 0.0000042 | 0.0000004 | 0 |
| Median | 0.0009575 | 0.0000834 | 0.000008 | 0.0000007 |
| St. deviation | 0.0013517 | 0.0001334 | 0.0000137 | 0.0000013 |
| Excess | −0.217306 | −0.191189 | −0.143295 | −0.021156 |
| Skewness | 0.0492335 | −0.002248 | 0.023124 | 0.0577173 |
| Min | −0.29706% | −0.03669% | −0.00303% | −0.00036% |
| Max | 0.52352% | 0.04766% | 0.00502% | 0.00044% |
| Count | 1000 | 1000 | 1000 | 1000 |
| $\alpha =100$ | ||||
| Average | 0.07790% | 0.00742% | 0.00083% | 0.00007% |
| St. error | 0.000043 | 0.0000044 | 0.0000004 | 0 |
| Median | 0.0008379 | 0.0000728 | 0.0000083 | 0.0000007 |
| St. deviation | 0.0013602 | 0.0001379 | 0.0000136 | 0.0000014 |
| Excess | −0.234452 | −0.302723 | −0.352995 | −0.054568 |
| Skewness | −0.024765 | 0.0922374 | 0.0055451 | 0.0229423 |
| Min | −0.30504% | −0.03231% | −0.00323% | −0.00037% |
| Max | 0.46265% | 0.04974% | 0.00454% | 0.00050% |
| Count | 1000 | 1000 | 1000 | 1000 |
It can be seen from the tables that approximation results do not differ significantly for various time-steps. Even the widest investigated discretization interval provides acceptable precision for most applications.
Table 6.
${\sigma }^{2}(Y_{s})={\operatorname{e}}^{Y_{s}}+0.2$, $K=1$, $r=0.02$, $k=0.1$, $T=1$. Characteristics of sample of errors
| ${10}^{-2}$ | ${10}^{-3}$ | ${10}^{-4}$ | ${10}^{-5}$ | |
| $\alpha =1$ | ||||
| Average | 0.02496% | 0.00268% | 0.00026% | 0.00002% |
| St. error | 0.0000113 | 0.0000011 | 0.0000001 | 0.00000001 |
| Median | 0.0002559 | 0.0000266 | 0.0000027 | 0.0000002 |
| St. deviation | 0.0003584 | 0.0000354 | 0.0000035 | 0.0000003 |
| Excess | 0.1947561 | 0.1687356 | −0.0576859 | 0.0700827 |
| Skewness | −0.1691937 | −0.0097185 | −0.1643507 | −0.0522007 |
| Min | −0.09961% | −0.00861% | −0.00088% | −0.00011% |
| Max | 0.12871% | 0.01464% | 0.00126% | 0.00013% |
| Count | 1000 | 1000 | 1000 | 1000 |
| $\alpha =100$ | ||||
| Average | 0.02692% | 0.00268% | 0.00025% | 0.00002% |
| St. error | 0.0000118 | 0.0000012 | 0.0000001 | 0 |
| Median | 0.0002712 | 0.0000265 | 0.0000027 | 0.0000002 |
| St. deviation | 0.0003735 | 0.0000377 | 0.0000036 | 0.0000003 |
| Excess | 0.17242 | 0.070383 | 0.3383414 | 0.0853763 |
| Skewness | −0.0299531 | −0.0195205 | −0.1914745 | −0.0371876 |
| Min | −0.09174% | −0.01068% | −0.00112% | −0.00011% |
| Max | 0.16291% | 0.01411% | 0.00139% | 0.00014% |
| Count | 1000 | 1000 | 1000 | 1000 |
5 Checking approximation precision in the case of deterministic volatility
In this section, we compare the option prices obtained for the Euler scheme (6) with the true prices of European call option for different sets of parameters for the case of deterministic time-dependent volatility.
The models with deterministic time-dependent volatility are the natural extension of the Black–Scholes model. The expression for the price of the option is the same as in the classical model except for the fact that, instead of constant volatility, it operates with average (or root mean square) volatility over the time interval to maturity (see, e.g., [10, 19]). Thus, the formula remains similar to (8) and (10).
It has been shown that deterministic volatility does not reflect the real-world stochastic dynamics correctly [3, 15], and such models have begun falling out of favor in the mid-1980s. The shift to stochastic volatility models was boosted by rapid development of computational tools.
Nevertheless, deterministic volatility is suitable for the purpose of our investigation since we can calculate the exact price of the option for the continuous time model.
In order to analyze the deterministic time-dependent volatility case, it looks natural to let the Brownian noise term in the definition of $Y_{t}$ vanish. Thus, we get
which is a familiar linear differential equation solved by
For the same transformation functions σ and sets of parameters as in the previous section, we calculate the prices of European call option in the continuous case using (8) and compare it with the prices of the same option calculated using (10)–(13) with
We use the time step of 0.01 and only 10 simulations per combination of inputs. As before, all calculations are performed in Matlab 7.9.0.
Table 7 presents the results of calculations. Comparison of two approaches reveals that the Euler–Maruyama scheme provides a good approximation for the exact option price. In the case of fast mean-reversion, the results coincide when rounded to sixth digit.
Table 7.
Approximate option prices versus true option prices for deterministic volatility
| T | α | r | K | a | b | $\hat{\mathbb{E}}\hat{P}_{m}$ | $\mathbb{E}V$ | $\hat{\mathbb{E}}\hat{P}_{m}$ | $\mathbb{E}V$ |
| ${\sigma }^{2}(Y_{s})=a\| Y_{s}\| +b$ | ${\sigma }^{2}(Y_{s})={\operatorname{e}}^{Y_{s}}+0.2$ | ||||||||
| 0.25 | 1 | 0 | 0.8 | 1 | 0 | 0.203891 | 0.203888 | 0.316223 | 0.316220 |
| 0.5 | 1 | 0 | 0.8 | 1 | 0 | 0.211556 | 0.211549 | 0.390150 | 0.390147 |
| 1 | 1 | 0 | 0.8 | 1 | 0 | 0.223003 | 0.222994 | 0.490305 | 0.490302 |
| 0.25 | 1 | 0.01 | 1 | 1 | 0.2 | 0.107942 | 0.107935 | 0.224736 | 0.224733 |
| 0.5 | 1 | 0.01 | 1 | 1 | 0.2 | 0.150207 | 0.150199 | 0.312794 | 0.312791 |
| 1 | 1 | 0.01 | 1 | 1 | 0.2 | 0.206464 | 0.206457 | 0.429067 | 0.429064 |
| 0.25 | 1 | 0.02 | 1.2 | 1 | 1 | 0.141317 | 0.141313 | 0.159958 | 0.159954 |
| 0.5 | 1 | 0.02 | 1.2 | 1 | 1 | 0.227633 | 0.227629 | 0.253710 | 0.253706 |
| 1 | 1 | 0.02 | 1.2 | 1 | 1 | 0.345261 | 0.345257 | 0.379955 | 0.379952 |
| 0.25 | 100 | 0 | 0.8 | 1 | 0 | 0.200000 | 0.200000 | 0.309950 | 0.309950 |
| 0.5 | 100 | 0 | 0.8 | 1 | 0 | 0.200000 | 0.200000 | 0.382107 | 0.382106 |
| 1 | 100 | 0 | 0.8 | 1 | 0 | 0.200000 | 0.200000 | 0.481610 | 0.481610 |
| 0.25 | 100 | 0.01 | 1 | 1 | 0.2 | 0.091044 | 0.091044 | 0.217149 | 0.217149 |
| 0.5 | 100 | 0.01 | 1 | 1 | 0.2 | 0.128449 | 0.128449 | 0.303457 | 0.303457 |
| 1 | 100 | 0.01 | 1 | 1 | 0.2 | 0.181507 | 0.181507 | 0.419198 | 0.419198 |
| 0.25 | 100 | 0.02 | 1.2 | 1 | 1 | 0.133108 | 0.133108 | 0.152065 | 0.152065 |
| 0.5 | 100 | 0.02 | 1.2 | 1 | 1 | 0.217100 | 0.217100 | 0.243748 | 0.243748 |
| 1 | 100 | 0.02 | 1.2 | 1 | 1 | 0.333759 | 0.333759 | 0.369312 | 0.369312 |
Remark 5.1.
In this paper, we consider the price of the option at the initial time moment. However, all the above considerations are applicable for any valuation date t between the initial time moment and maturity. Some obvious changes need to be made, for example, the function $\bar{\sigma }_{t}:=\sqrt{\frac{1}{T-t}{\int _{t}^{T}}{\sigma }^{2}(Y_{s})ds}\ge 0$ needs to be introduced instead of $\bar{\sigma }$, and T needs to be substituted by $T-t$ in (8)–(13).
