We consider a discrete-time approximation of paths of an Ornstein–Uhlenbeck process as a mean for estimation of a price of European call option in the model of financial market with stochastic volatility. The Euler–Maruyama approximation scheme is implemented. We determine the estimates for the option price for predetermined sets of parameters. The rate of convergence of the price and an average volatility when discretization intervals tighten are determined. Discretization precision is analyzed for the case where the exact value of the price can be derived.
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].
The model and discrete approximation of volatility process
Let {Ω,F,F={Ft(B,Z),t≥0},P} be a complete probability space with filtration generated by Wiener processes {Bt,Zt,0≤t≤T}. We consider the model of the market where one risky asset is traded, its price evolves according to the geometric Brownian motion {St,0≤t≤T}, and its volatility is driven by a stochastic process. More precisely, the market is described by the pair of stochastic differential equations
dSt=μStdt+σ(Yt)StdBt,dYt=−αYtdt+kdZt.
We denote by S0=S and Y0=Y the deterministic initial values of the processes specified by Eqs. (1)–(2).
We further impose the following assumptions:
The Wiener processes B and Z are uncorrelated;
the volatility function σ:R→R+ is measurable, bounded away from zero by a constant c:
σ(x)≥c>0,x∈R,0,\hspace{1em}x\in \mathbb{R},\]]]>
and satisfies the condition ∫0Tσ2(Yt)dt<∞ a.s.;
the coefficients α, μ, and k are positive.
For example, the conditions mentioned in assumption (C2) are satisfied for the measurable function σ(x) such that c≤σ2(x)≤C for 0<x<T and some constants 0<c<C. Moreover, given the square integrability of σ(Ys), the solution of differential equation (1) is given by
St=S0exp(μt−12∫0tσ2(Ys)ds+∫0tσ(Ys)dBs),
which yields that St is continuous. Hence, the product σ(Ys)St is square integrable: ∫0Tσ2(Yt)St2dt<∞ a.s.
The unique solution of the Langevin equation (2) Yt 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:
E[Yt]=Y0e−αt,Var[Yt]=k22α(1−e−2αt).
Moreover, the OU process is Markov and admits the explicit representation
Yt=Y0e−αt+k∫0te−α(t−s)dZs.
Following [9], we proceed to the risk-neutral setting characterized by the minimal martingale measure Q. With r being the interest rate, Eqs. (1)–(2) are now in the following form (see Section 5 in [9]):
dSt=rStdt+σ(Yt)StdBtQ,dYt=−αYtdt+kdZtQ,
where
BtQ=Bt+∫0tμ−rσ(Ys)dsandZtQ=Zt
are independent Wiener processes w.r.t. 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 partition of the interval [0,T] into m equal subintervals of width Δt=T/m is considered;
the initial value of the scheme is set: Y0(m)=Y0;
Yl+1(m), which we will use as a shorthand for Y(l+1)T/m(m), 0≤l≤m−1, is recursively defined by
Yl+1(m)=(1−αΔt)Yl(m)+kΔZlQ,
where ΔZlQ=Z(l+1)T/mQ−ZlT/mQ.
The continuous-time process Yt(m) is a step-type process defined by
Yt(m)=Y[tm/T]T/m(m),t∈[0,T],
where [x] denotes an integer part of x.
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
V=e−rTEQ{EQ{(STQ−K)+|Ys,0≤s≤T}}.
The inner expectation is conditional on the path of Ys, 0≤s≤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:
P:=EQ{(STQ−K)+|Ys,0≤s≤T}=elnS+rTΦ(d1)−KΦ(d2):=elnS+rTΦ(lnS+(r+12σ¯2)T−lnKσ¯T)−KΦ(lnS+(r−12σ¯2)T−lnKσ¯T),
where σ¯:=1T∫0Tσ2(Ys)ds≥0, Φ is the standard normal distribution function. The function σ¯ may be viewed as the volatility averaged from the initial moment of time to maturity. The arguments of Φ in (8) are denoted as d1 and d2.
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
R:=|P−Pˆm|,
where m is the number of discretization points dividing the time interval [0,T] into equal intervals, Pˆm denotes the price of the option in discrete setting calculated using the formula similar to (8):
Pˆm=elnS+rTΦ(d1(m))−KΦ(d2(m)),
where d1(m)=lnS+(r+12σ¯m2)T−lnKσ¯mT,d2(m)=lnS+(r−12σ¯m2)T−lnKσ¯mT, with
σ¯m=1T∑l=1mσ2(Yl(m))Tm=1m∑l=1mσ2(Yl(m)),
where Yl(m) is defined in (6).
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 σ¯ and σ¯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σ:=E|σ¯−σ¯m|.
Letσ2(x)satisfy the Hölder condition|σ2(x)−σ2(y)|≤L|x−y|γ,where0<γ≤1, and L is some positive constant. ThenERσ≤Cm−0.5γ, where C is some positive constant.
Since σ¯m and σ¯ are both square root functions, it is be more convenient to work with σ¯m2 and σ¯2. To this end, we will use Hölder’s inequality:
E|σ¯m−σ¯|=E|1T∫0Tσ2(Ys)ds−1m∑i=1mσ2(Yi(m))|≤E||1T∫0Tσ2(Ys)ds−1m∑i=1mσ2(Yi(m))||≤(E|1T∫0Tσ2(Ys)ds−1m∑i=1mσ2(Yi(m))|)1/2.
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:
E|σ¯m−σ¯|≤(E|∑i=0m−1(1T∫iT/m(i+1)T/mσ2(Ys)ds−1mσ2(Yi+1(m)))|)1/2=(E|∑i=0m−1(1T∫iT/m(i+1)T/mσ2(Ys)ds−1mmT∫iT/m(i+1)T/mσ2(Yi+1(m))ds)|)1/2=(E|1T∑i=0m−1∫iT/m(i+1)T/m(σ2(Ys)−σ2(Yi+1(m)))ds|)1/2.
We apply the Hölder property of σ2(x):
(E|1T∑i=0m−1∫iT/m(i+1)T/m(σ2(Ys)−σ2(Yi+1(m)))ds|)1/2≤(LTE(∑i=0m−1∫iT/m(i+1)T/m|Ys−Yi+1(m)|γds))1/2=(LT∑i=0m−1∫iT/m(i+1)T/mE|Ys−Yi+1(m)|γds)1/2.
Recall that Yi(m) is a shorthand for YiT/m(m)=Ys(m), s∈[iT/m,(i+1)T/m), and Proposition from the Appendix A yields that E|Ys−Yi+1(m)|≤C1m−1, where C1 is some positive constant. We use Hölder’s inequality to derive that E|Ys−Yi+1(m)|γ≤C1γm−γ and arrive at
E|σ¯m−σ¯|≤(LTmTmC1γm−γ)1/2=Cm−γ/2
for C:=LC1γ, which proves the lemma. □
The above lemma enables us to prove the main result of this section.
Letσ2(x)satisfy Hölder condition (14). ThenER≤Dm−γ/2, where D is some positive constant.
The function Φ(x) has a continuous bounded derivative on R; hence, we can use its Lipschitz property:
ER=E|P−Pˆm|≤E(SerT|Φ(d1)−Φ(d1(m))|+K|Φ(d2)−Φ(d2(m))|)≤E(supx|f(x)|(SerT|d1−d1(m)|+K|d2−d2(m)|)),
where f(x)=12πe−x22 is the density of the standard normal distribution. In the above representation, |d1−d1(m)|=|(1σ¯−1σ¯m)lnS−lnK+rTT+12T(σ¯−σ¯m)|≤|σ¯−σ¯m||1σ¯σ¯mln(S/K)+rTT+T2|≤|σ¯−σ¯m||ln(S/K)+rTc2T+T2|, where c is a positive constant, and the last inequality is due to the assumption that σ(x) is bounded away from zero for any x∈R (see assumption (C2)). Hence, using Lemma 3.1, we get
E|d1−d1(m)|≤C1E|σ¯−σ¯m|≤D1m−γ/2,
where C1:=|ln(S/K)+rTc2T+T2| and D1 are positive constants.
Similarly, E|d2−d2(m)|≤D2m−γ/2, D2=const>00$]]>, and we arrive at
ER=E|P−Pˆm|≤12π(D1SerTm−γ/2+D2Km−γ/2)=Dm−γ/2
for a positive constant D.
The theorem is proved. □
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:
The choice of discrete ranges of values of input parameters;
The choice of the function σ(Ys);
For each combination of input parameters we generate 1000 trajectories of an Ornstein–Uhlenbeck process by splitting the time interval into subintervals of length Δ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 σ¯m2 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 EˆPˆm. The average volatility over all trajectories and time interval is denoted by Eˆσ¯m2.
σ2(Ys)=a|Ys|+b
T
k
r
K
a
b
Eˆσ¯m2
EˆPˆm
Eˆσ¯m2
EˆPˆm
α=1
α=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;0.5;1;
k – volatility of OU process, k=0.1;0.5;1;
α – mean-reversion rate, α=1;100;
r – interest rate, r=0;0.01;0.02;
K – strike price, K=0.8;1;1.2;
S0 – initial price of stock, S0=1;
Y0 – initial value of OU process, Y0=0.1.
In order to produce numerical results, we choose the following options for the function σ(Ys):
σ2(Ys)=a|Ys|+b, where a={0,1}, b={0,0.2,1} (Table 1);
σ2(Ys)=eYs+c, c=0.02 (Table 2).
The results of simulations are split into groups by the mean-reversion rate α and function σ(Ys). 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 α=100. Matters of speed of mean-reversion are addressed, for example, in [4].
σ2(Ys)=eYs+c
T
k
r
K
Eˆσ¯m2
EˆPˆm
Eˆσ¯m2
EˆPˆm
α=1
α=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 Eˆσ¯m2 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.
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 α=1 and 1000 for the case α=100. These trajectories are considered to be “true” continuous-time trajectories of the Ornstein–Uhlenbeck process Yt. The corresponding values of σ¯m2 are considered to be “true” continuous-time values of σ¯2. The calculations were then performed for wider discretization intervals using the points of constructed trajectories. Thus, the samples of discretization errors for σ¯m2 were derived. Probably, the estimate of σ¯m2 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.
σ2(Ys)=|Ys|+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
α=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
α=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.
σ2(Ys)=eYs+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
α=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
α=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
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 Yt vanish. Thus, we get
dYt=−αYtdt,
which is a familiar linear differential equation solved by
Yt=Y0e−αt.
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
Yl+1(m)=(1−αΔt)Yl(m).
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.
Approximate option prices versus true option prices for deterministic volatility
T
α
r
K
a
b
EˆPˆm
EV
EˆPˆm
EV
σ2(Ys)=a‖Ys‖+b
σ2(Ys)=eYs+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
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 σ¯t:=1T−t∫tTσ2(Ys)ds≥0 needs to be introduced instead of σ¯, and T needs to be substituted by T−t in (8)–(13).
The Euler scheme: definitions and auxiliary results
The reader is advised to refer to [8], which provides in-depth study of numerical approximations of stochastic differential equations.
Consider the stochastic differential equation
dXt=a(t,Xt)dt+b(t,Xt)dWt,t∈[t0,T],
and assume that there is a unique strong solution X(t) with X(t0)=X0. In order for this to be the case, certain assumptions need to be made about the functions a and b. Namely, refer to the following assumptions (assumptions (A1)–(A4) in [8], pp. 128–129):
a=a(t,x) and b=b(t,x) are jointly L2-measurable in (t,x)∈[t0,T]×R;
the functions a and b satisfy the Lipschitz condition w.r.t. x, that is, there exists a constant K>00$]]> such that
|a(t,x)−a(t,y)|≤K|x−y|
and
|b(t,x)−b(t,y)|≤K|x−y|
for all t∈[t0,T] and x,y∈R;
there exists a constant K>00$]]> such that
|a(t,x)|2≤K|1+|x|2|
and
|b(t,x)|2≤K|1+|x|2|
for all t∈[t0,T] and x,y∈R;
Xt0 is Ft0-measurable with E|Xt0|2<∞.
Let Xt(m) be a discretization scheme of the process Xt.
(See [8].) We shall say that an approximating process Xt(m) converges in the strong sense with order γ∈(0,∞] to the true process Xt if there exists a finite constant K such that
E(|Xt−Xt(m)|)≤Km−γ.
The same terminology will be applied to the functions of approximating processes.
(See [8].) We shall say that a discrete time approximation scheme Xt(m) is strongly consistent if there exists a nonnegative function c=c(m) with
limm→∞c(m)=0
such that
E(|E(Xi+1(m)−Xi(m)T/m|FiT/m)−a(iTm,Xi(m))|2)≤c(m)
and
E(mT|Xi+1(m)−Xi(m)−E(Xi+1(m)−Xi(m)|FiT/m)−b(iTm,Xi(m))ΔWi|2)≤c(m)
for all fixed values Xi(m)=y and i=0,1,…,m.
(See [8], 9.6.2, p. 324.) Let assumptions (A1)–(A4) hold for (22). Then a strongly consistent equidistant-time discrete approximationX(m)of the process X on[t0,T], withXt0(m)=Xt0, converges strongly to X.
Evidently, the Euler scheme Y(m) introduced to approximate Y in Section 2 satisfies all the above requirements and hence converges strongly. Moreover, it is a well-known fact that, in general, the convergence of the Euler approximation is of order 0.5. One can check these propositions using the estimates of the rate of convergence provided in [8] by the proof of Theorem 9.6.2 and Exercise 9.6.3.
However, our case is more specific since Y(m) approximates the diffusion process with additive noise, that is, b(t,x)=k is constant.Hence, the following proposition holds.
Y(m)is the Milstein scheme and thus converges strongly with order 1.
Really, the only difference in representation of Y(m) as the Milstein scheme compared to the Euler one is in the additional summand of the form
12bb′((ΔZQ)2−T/m),
which is identically zero for the constant function b. The Milstein scheme is known to converge with order 1 (see, e.g., [8], Theorem 10.6.3, p. 361).
ReferencesAndersen, L.: Simple and efficient simulation of the Heston stochastic volatility model. J. Comput. Finance11(3), 1–42 (2008) Broadie, M., Kaya, O.: Exact simulation of stochastic volatility and other affine jump diffusion processes. Oper. Res.54(2), 217–231 (2006). MR2222897. doi:10.1287/opre.1050.0247Brockman, P., Chowdhury, M.: Deterministic versus stochastic volatility: Implications for option pricing models. Appl. Financ. Econ.7, 499–505 (1997) Fouque, J.-P., Papanicolaou, G., Sircar, K.R.: Derivatives in Financial Markets with Stochastic Volatility. Cambridge University Press, USA (2000). MR1768877Giles, M.B.: Multilevel Monte Carlo path simulation. Oper. Res.56(3), 607–617 (2008). MR2436856. doi:10.1287/opre.1070.0496Hull, J., White, A.: The pricing of options on assets with stochastic volatilities. J. Finance42(2), 281–300 (1987) Jourdain, B., Sbai, M.: High order discretization schemes for stochastic volatility models. arXiv:0908.1926 (2009) Kloeden, P., Platen, E.: Numerical Solution of Stochastic Differential Equations. Springer, New York (1992). MR1214374. doi:10.1007/978-3-662-12616-5Kuchuk-Iatsenko, S., Mishura, Y.: Pricing the European call option in the model with stochastic volatility driven by Ornstein–Uhlenbeck process. Exact formulas. Modern Stoch., Theory Appl.2(3), 233–249 (2015) Lewis, A.: Option Valuation Under Stochastic Volatility: With Mathematica Code. Finance Press, Newport Beach, CA (2000). MR1742310Maruyama, G.: Continuous Markov processes and stochastic equations. Rend. Circ. Mat. Palermo4, 48–90 (1955). MR0071666Milstein, G.: Approximate integration of stochastic differential equations. Theory Probab. Appl.19(3), 557–562 (1975). MR0356225Mishura, Y., Rizhniak, G., Zubchenko, V.: European call option issued on a bond governed by a geometric or a fractional geometric Ornstein–Uhlenbeck process. Modern Stoch., Theory Appl.1, 95–108 (2014). MR3314796. doi:10.15559/vmsta-2014.1.1.2Platen, E., N., B.-L.: Numerical Solution of Stochastic Differential Equations with Jumps in Finance. Springer, New York (2010). MR2723480. doi:10.1007/978-3-642-13694-8Sheikh, A.M.: Transaction data tests of S&P 100 call option pricing. J. Financ. Quant. Anal.26, 459–475 (1991) van Haastrecht, A., Pelsser, A.: Efficient, almost exact simulation of the Heston stochastic volatility model. Int. J. Theor. Appl. Finance31(1), 1–42 (2010). MR2646972. doi:10.1142/S0219024910005668Wang, H.: Monte Carlo Simulation with Applications to Finance. Chapman and Hall/CRC, Boca Raton, FL (2012). MR2962599Willard, G.: Calculating prices and sensitivities for path-independent derivative securities in multifactor models. J. Deriv.5(1), 45–61 (1997) Wilmott, P., Howinson, S., Dewynne, I.: The Mathematics of Financial Derivatives. Cambridge University Press (1995). MR1357666. doi:10.1017/CBO9780511812545