Closed-form approximations with respect to the mixing solution for option pricing under stochastic volatility

TL;DR

This paper develops closed-form approximations for European put options under stochastic volatility models, enabling fast calibration and accurate pricing with error bounds, by expanding around the mean of a Black-Scholes expectation.

Contribution

It introduces a second-order Taylor expansion approach for option pricing under stochastic volatility with piecewise-constant parameters, providing explicit formulas and error bounds.

Findings

Errors are within acceptable ranges for practical use

Derived bounds on Taylor expansion remainder

Fast calibration scheme demonstrated

Abstract

We consider closed-form approximations for European put option prices within the Heston and GARCH diffusion stochastic volatility models with time-dependent parameters. Our methodology involves writing the put option price as an expectation of a Black-Scholes formula and performing a second-order Taylor expansion around the mean of its argument. The difficulties then faced are simplifying a number of expectations induced by the Taylor expansion. Under the assumption of piecewise-constant parameters, we derive closed-form pricing formulas and devise a fast calibration scheme. Furthermore, we perform a numerical error and sensitivity analysis to investigate the quality of our approximation and show that the errors are well within the acceptable range for application purposes. Lastly, we derive bounds on the remainder term generated by the Taylor expansion.

Tables

Table 1. Table 7.1 . Signed error of implied volatilities in basis points, computed in the Heston model with κ 𝜅 \kappa varying from 40% to 160% of its “safe” 3-piece value.

	$\kappa$	40% 0	60% 0	80%	100%	120%	140%	160%
ATM	1M	0.77 0	1.08 0	1.23 0	1.32	1.37	1.40	1.41
	3M	-26.35	-15.43	-9.00	-4.96	-2.30	-0.45	0.87
	6M	-59.94	-31.56	-16.97	-8.82	-3.96	-0.93	1.00
	1Y	-70.41	-29.66	-12.76	-4.92	-1.05	0.95	1.98
Put 25	1M	12.81	11.01	9.65	8.61	7.79	7.13	6.59
	3M	14.05	12.89	12.08	11.43	10.89	10.41	9.98
	6M	3.83	9.08	11.35	12.14	12.17	11.82	11.29
	1Y	-7.47	6.88	11.14	11.86	11.34	10.42	9.44
Put 10	1M	34.21	27.41	22.52	18.91	16.19	14.08	12.40
	3M	86.32	62.63	47.89	38.17	31.40	26.50	22.81
	6M	115.94	78.44	57.53	44.60	35.91	29.68	25.04
	1Y	105.65	66.94	47.08	35.05	27.06	21.45	17.38

Table 2. Table 7.2 . Signed error of implied volatilities in basis points, computed in the Heston model with θ 𝜃 \theta varying from 40% to 160% of its “safe” 3-piece value.

	$\theta$	40%	60%	80%	100%	120%	140%	160%
ATM	1M	-0.12 0	0.22 0	0.32 0	0.32	0.25	0.14	0.01
	3M	-20.59	-12.98	-8.66	-6.05	-4.41	-3.35	-2.65
	6M	-32.70	-20.25	-13.51	-9.46	-6.90	-5.19	-4.02
	1Y	-21.50	-12.02	-7.09	-4.24	-2.45	-1.28	-0.47
Put 25	1M	14.04	11.50	9.69	8.35	7.31	6.47	5.77
	3M	24.69	17.26	13.10	10.49	8.68	7.33	6.26
	6M	31.34	20.93	15.91	13.07	11.24	9.94	8.93
	1Y	28.59	19.26	15.15	12.80	11.30	10.23	9.41
Put 10	1M	33.55	26.93	22.07	18.45	15.69	13.53	11.79
	3M	84.65	63.16	48.78	38.90	31.87	26.69	22.76
	6M	103.04	74.71	56.46	44.40	36.11	30.16	25.72
	1Y	81.69	58.22	43.82	34.74	28.72	24.48	21.36

Table 3. Table 7.3 . Signed error of implied volatilities in basis points, computed in the Heston model with λ 𝜆 \lambda varying from 40% to 160% of its “safe” 3-piece value.

	$\lambda$	40%	60%	80%	100%	120%	140%	160%
ATM	1M	0.04 0	0.59 0	0.78 0	0.15	-2.00	-6.32	-13.33
	3M	-0.07	0.15	-1.68	-7.00	-17.03	-32.53	-54.16
	6M	0.98	0.74	-2.26	-9.55	-22.26	-41.14	-66.66
	1Y	2.37	2.28	0.41	-4.38	-12.88	-25.70	-43.27
Put 25	1M	1.31	3.36	5.61	7.94	10.00	11.34	11.55
	3M	2.50	5.92	9.21	11.35	11.31	8.24	1.66
	6M	4.09	7.93	11.34	13.02	11.87	7.04	-2.11
	1Y	5.10	8.18	11.05	12.73	12.43	9.50	3.39
Put 10	1M	0.47	4.78	10.97	19.32	29.81	42.31	56.54
	3M	2.10	10.46	22.57	38.09	56.38	76.63	98.37
	6M	4.70	14.43	28.08	44.95	64.13	84.85	106.37
	1Y	5.13	12.44	22.53	34.81	48.54	63.19	78.41

Table 4. Table 7.4 . Signed error of implied volatilities in basis points, computed in the Heston model with ρ 𝜌 \rho varying from 160% to 40% of its “safe” 3-piece value.

	$\rho$	160%	140%	120%	100%	80%	60%	40%
ATM	1M	29.14	14.93	6.37	1.16	-2.02 0	-3.96 0	-5.10 0
	3M	27.01	10.49	-0.06	-6.97	-11.56	-14.60	-16.52
	6M	24.82	8.30	-2.34	-9.41	-14.19	-17.41	-19.49
	1Y	25.35	10.58	1.42	-4.47	-8.34	-10.91	-12.54
Put 25	1M	44.10	28.02	16.46	8.28	2.52	-1.49	-4.19
	3M	63.28	40.72	23.85	11.21	1.67	-5.51	-10.85
	6M	68.04	44.03	26.03	12.39	1.96	-6.03	-12.12
	1Y	60.63	39.51	24.27	13.08	4.73	-1.55	-6.27
Put 10	1M	28.21	28.60	24.52	18.96	13.15	7.68	2.81
	3M	64.31	58.82	49.05	38.07	27.19	17.01	7.78
	6M	76.48	68.62	56.94	44.37	32.11	20.69	10.33
	1Y	64.00	55.57	44.91	34.16	24.11	15.04	7.05

Table 5. Table 7.5 . Signed error of implied volatilities in basis points, computed in the GARCH diffusion model with κ 𝜅 \kappa varying from 40% to 160% of its “safe” 3-piece value.

	$\kappa$	40%	60%	80%	100%	120%	140%	160%
ATM	1M	-1.33 0	-1.70 0	-2.06 0	-2.35	-2.58	-2.77	-2.93
	3M	-2.69	-3.10	-3.30	-3.31	-3.22	-3.06	-2.87
	6M	-3.32	-3.16	-2.74	-2.25	-1.78	-1.36	-1.02
	1Y	-2.79	-1.80	-1.04	-0.54	-0.25	-0.08	0.00
Put 25	1M	-1.36	-1.73	-2.09	-2.38	-2.61	-2.81	-2.96
	3M	-2.71	-3.12	-3.31	-3.32	-3.23	-3.07	-2.87
	6M	-3.30	-3.13	-2.72	-2.23	-1.75	-1.34	-1.00
	1Y	-2.77	-1.77	-1.02	-0.52	-0.23	-0.07	0.02
Put 10	1M	-1.37	-1.74	-2.10	-2.39	-2.62	-2.81	-2.97
	3M	-2.75	-3.16	-3.35	-3.37	-3.27	-3.11	-2.91
	6M	-3.33	-3.17	-2.75	-2.25	-1.78	-1.36	-1.01
	1Y	-2.76	-1.76	-1.00	-0.51	-0.22	-0.06	0.02

Table 6. Table 7.6 . Signed error of implied volatilities in basis points, computed in the GARCH diffusion model with θ 𝜃 \theta varying from 40% to 160% of its “safe” 3-piece value.

	$\theta$	40%	60%	80%	100%	120%	140%	160%
ATM	1M	-1.09 0	-1.51 0	-1.95 0	-2.35	-2.71	-3.04	-3.35
	3M	-1.43	-2.14	-2.78	-3.34	-3.85	-4.31	-4.73
	6M	-0.67	-1.29	-1.82	-2.27	-2.67	-3.02	-3.35
	1Y	0.47	0.08	-0.22	-0.46	-0.66	-0.84	-1.00
Put 25	1M	-1.11	-1.54	-1.98	-2.38	-2.74	-3.07	-3.38
	3M	-1.42	-2.13	-2.77	-3.34	-3.84	-4.30	-4.73
	6M	-0.65	-1.26	-1.79	-2.24	-2.63	-2.99	-3.31
	1Y	0.43	0.04	-0.26	-0.51	-0.72	-0.90	-1.07
Put 10	1M	-1.13	-1.56	-2.00	-2.39	-2.75	-3.09	-3.40
	3M	-1.43	-2.14	-2.78	-3.35	-3.85	-4.31	-4.74
	6M	-0.69	-1.30	-1.82	-2.27	-2.66	-3.02	-3.34
	1Y	0.40	0.02	-0.29	-0.53	-0.75	-0.93	-1.10

Table 7. Table 7.7 . Signed error of implied volatilities in basis points, computed in the GARCH diffusion model with λ 𝜆 \lambda varying from 40% to 160% of its “safe” 3-piece value.

	$\lambda$	40%	60%	80%	100%	120%	140%	160%
ATM	1M	-2.51 0	-2.36 0	-2.36 0	-2.36	-2.36	-2.36	-2.35
	3M	-3.43	-3.33	-3.32	-3.32	-3.31	-3.31	-3.31
	6M	-2.33	-2.26	-2.26	-2.26	-2.26	-2.26	-2.27
	1Y	-0.55	-0.51	-0.50	-0.50	-0.50	-0.50	-0.51
Put 25	1M	-2.51	-2.36	-2.36	-2.35	-2.35	-2.35	-2.34
	3M	-3.42	-3.31	-3.31	-3.30	-3.29	-3.29	-3.28
	6M	-2.32	-2.26	-2.25	-2.25	-2.25	-2.25	-2.25
	1Y	-0.57	-0.54	-0.55	-0.56	-0.57	-0.59	-0.61
Put 10	1M	-2.51	-2.36	-2.36	-2.36	-2.36	-2.35	-2.35
	3M	-3.43	-3.33	-3.32	-3.32	-3.32	-3.31	-3.31
	6M	-2.33	-2.26	-2.26	-2.26	-2.26	-2.26	-2.25
	1Y	-0.55	-0.50	-0.50	-0.49	-0.49	-0.48	-0.48

Table 8. Table 7.8 . Runtimes for pricing a put option with maturities 6M, 1Y in the Heston model via our closed-form approximation method, where #pieces of the parameters vary. The length of each time-interval corresponding to a piece is equal. Each runtime is the average of 100 runtimes (with identical parameters).

	#Pieces		Run(ms)	Slow down($\times$)
6M	1		0.77	1.00
	2		1.43	1.87
	3		2.48	3.24
	4		4.70	6.13
	5		7.87	10.27
	6		11.05	14.43
	7		16.98	22.18
	8		22.21	29.00
	9		31.53	41.17
	10		40.73	53.19
1Y	1		0.80	1.00
	2		1.71	2.13
	3		2.65	3.30
	4		4.60	5.74
	5		6.59	8.21
	6		10.94	13.64
	7		15.80	19.69
	8		22.17	27.65
	9		29.75	37.09
	10		40.08	49.97

Table 9. Table 7.9 . Runtimes for pricing a put option with maturities 6M, 1Y in the GARCH diffusion model via our closed-form approximation method, where the #pieces of parameters vary. The length of each time-interval corresponding to a piece is equal. Each runtime is the average of 100 runtimes (with identical parameters).

	#Pieces		Run(ms)	Slow down($\times$)
6M	1		0.69	1.00
	2		1.37	1.98
	3		2.81	4.05
	4		4.56	6.58
	5		8.30	11.96
	6		13.07	18.84
	7		20.28	29.22
	8		30.56	44.03
	9		44.38	63.94
	10		62.86	90.57
1Y	1		0.67	1.00
	2		1.45	2.17
	3		2.63	3.93
	4		5.08	7.59
	5		8.75	13.07
	6		12.85	19.19
	7		20.79	31.04
	8		29.50	44.05
	9		44.24	66.07
	10		65.43	97.71

Equations

$$\mathrm{d}{S}_t$$

$$\mathrm{d}{V}_t$$

$$\mathrm{d}⟨W,B{⟩}_t$$

$$\mathrm{d}{S}_t$$

$$\mathrm{d}{V}_t$$

$$\mathrm{d}⟨W,B{⟩}_t$$

$$\mathrm{d}{S}_t$$

$$\mathrm{d}{\sigma }_t$$

$$\mathrm{d}⟨W,B{⟩}_t$$

$$\mathrm{d}{Y}_t={\mu }_t{Y}_t\mathrm{d}t+{\nu }_t{Y}_t\mathrm{d}{\tilde{B}}_t,Y_0=y_0,$$

$$S_T$$

$${\xi }_t$$

$$\mathrm{d}{\xi }_t$$

$${\xi }_t^n=\exp \left({\int }_0^t{\rho }_u\sqrt{{\sigma }_u}{\text{1}}_{\{u\le {\tau }_n\}}\mathrm{d}{B}_u-\frac{1}{2}{\int }_0^t{\rho }_u^2{\sigma }_u{\text{1}}_{\{u\le {\tau }_n\}}\mathrm{d}u\right).$$

$$\underset{n}{\sup }\mathbb{E}({\xi }_T^n\log ({\xi }_T^n))<\infty$$

$${\xi }_T^n\log ({\xi }_T^n)={\xi }_T^n\left({\int }_0^T{\rho }_t\sqrt{{\sigma }_t}{\text{1}}_{\{t\le {\tau }_n\}}\mathrm{d}{B}_t-\frac{1}{2}{\int }_0^T{\rho }_t^2{\sigma }_t{\text{1}}_{\{t\le {\tau }_n\}}\mathrm{d}t\right).$$

$$\mathbb{E}({\xi }_T^n\log ({\xi }_T^n))$$

$$=\hat{\mathbb{E}}\left({\int }_0^T{\rho }_t\sqrt{{\sigma }_t}{\text{1}}_{\{t\le {\tau }_n\}}\mathrm{d}{\hat{B}}_t+\frac{1}{2}{\int }_0^T{\rho }_t^2{\sigma }_t{\text{1}}_{\{t\le {\tau }_n\}}\mathrm{d}t\right)$$

$$=\frac{1}{2}{\int }_0^T\hat{\mathbb{E}}({\rho }_t^2{\sigma }_t{\text{1}}_{\{t\le {\tau }_n\}})\mathrm{d}t$$

$$\le \frac{1}{2}{\int }_0^T\hat{\mathbb{E}}({\sigma }_t)\mathrm{d}t.$$

$$\displaystyle\mathrm{d}\sigma_{t}=\Big{[}\rho_{t}\beta(t,\sigma_{t})\sqrt{\sigma_{t}}\textbf{1}_{\{t\leq\tau_{n}\}}+\alpha(t,\sigma_{t})\Big{]}\mathrm{d}t+\beta(t,\sigma_{t})\mathrm{d}\hat{B}_{t}.$$

$${\rho }_t\beta (t,x)\sqrt{x}+\alpha (t,x)\le M(1+x),$$

$$\mathrm{d}{\sigma }_t\le M(1+{\sigma }_t)\mathrm{d}t+\beta (t,{\sigma }_t)\mathrm{d}{\hat{B}}_t.$$

$$\hat{\mathbb{E}}({\sigma }_t)$$

$$={\sigma }_0+M{\int }_0^t(1+\hat{\mathbb{E}}({\sigma }_u))\mathrm{d}u.$$

$$m_t\le c_0(1+t)+c{\int }_0^t{m}_u\mathrm{d}u.$$

$$m_t\le c_0(1+t)e^{ct}.$$

$$\text{Put}:=e^{-{\int }_0^T{r}_t^d\mathrm{d}t}\mathbb{E}(K-S_T{)}_{+}$$

$$\displaystyle\begin{split}\text{Put}&=\mathbb{E}\left(e^{-\int_{0}^{T}r_{t}^{d}\mathrm{d}t}\mathbb{E}\big{[}(K-S_{T})_{+}|\mathcal{F}_{T}^{B}\big{]}\right)\\ &=\mathbb{E}\left(\text{Put}_{\text{BS}}\left(S_{0}\xi_{T},\int_{0}^{T}\sigma_{t}(1-\rho_{t}^{2})\mathrm{d}t\right)\right),\end{split}$$

$$\begin{array}{rl}{\text{Put}}_{\text{BS}}(x,y)& :=K{e}^{-{\int }_0^T{r}_t^d\mathrm{d}t}\mathcal{N}(-d_{-})-x{e}^{-{\int }_0^T{r}_t^f\mathrm{d}t}\mathcal{N}(-d_{+}),\\ d_{\pm }(x,y):=d_{\pm }& :=\frac{\ln (x/K)+{\int }_0^T(r_t^d-r_t^f)\mathrm{d}t}{\sqrt{y}}\pm \frac{1}{2}\sqrt{y}.\end{array}$$

$$\begin{array}{rl}{\text{Put}}^{(2)}& ={\text{Put}}_{\text{BS}}(\hat{x},\hat{y})\\ & +\frac{1}{2}{\partial }_{xx}{\text{Put}}_{\text{BS}}(\hat{x},\hat{y})S_0^2\mathbb{E}({\xi }_T-1{)}^2+\frac{1}{2}{\partial }_{yy}{\text{Put}}_{\text{BS}}(\hat{x},\hat{y})\mathbb{E}{\left({\int }_0^T(1-{\rho }_t^2)({\sigma }_t-\mathbb{E}({\sigma }_t))\mathrm{d}t\right)}^2\\ & +{\partial }_{xy}{\text{Put}}_{\text{BS}}(\hat{x},\hat{y})S_0\mathbb{E}\left\{({\xi }_T-1)\left({\int }_0^T(1-{\rho }_t^2)({\sigma }_t-\mathbb{E}({\sigma }_t))\mathrm{d}t\right)\right\},\end{array}$$

$${\text{Put}}_{\text{BS}}\left(S_0{\xi }_T,{\int }_0^T{\sigma }_t(1-{\rho }_t^2)\mathrm{d}t\right)\approx {\text{Put}}_{\text{BS}}(\hat{x},\hat{y})$$

$$+{\partial }_x{\text{Put}}_{\text{BS}}(\hat{x},\hat{y})S_0({\xi }_T-1)+{\partial }_y{\text{Put}}_{\text{BS}}(\hat{x},\hat{y})\left({\int }_0^T(1-{\rho }_t^2)({\sigma }_t-\mathbb{E}({\sigma }_t))\mathrm{d}t\right)$$

$$+\frac{1}{2}{\partial }_{xx}{\text{Put}}_{\text{BS}}(\hat{x},\hat{y})S_0^2({\xi }_T-1{)}^2+\frac{1}{2}{\partial }_{yy}{\text{Put}}_{\text{BS}}(\hat{x},\hat{y}){\left({\int }_0^T(1-{\rho }_t^2)({\sigma }_t-\mathbb{E}({\sigma }_t))\mathrm{d}t\right)}^2$$

$$+{\partial }_{xy}{\text{Put}}_{\text{BS}}(\hat{x},\hat{y})S_0({\xi }_T-1)\left({\int }_0^T(1-{\rho }_t^2)({\sigma }_t-\mathbb{E}({\sigma }_t))\mathrm{d}t\right).$$

$$\mathbb{E}({\xi }_T-1{)}^2,$$

$$\mathbb{E}{\left({\int }_0^T(1-{\rho }_t^2)({\sigma }_t-\mathbb{E}({\sigma }_t))\mathrm{d}t\right)}^2,$$

Full text

Keywords: Stochastic volatility; Closed-form expansion; Closed-form approximation; Heston; GARCH

††Research supported by an Australian Government Research Training Program (RTP) Scholarship.

1. Introduction

We consider the European put option pricing problem within the Heston and GARCH111Generalised AutoRegressive Conditional Heteroskedasticity. diffusion stochastic volatility models, see [18] and [30, 36] respectively. Our goal is to study how European put option prices can be approximated in each of these frameworks via expansion of the so-called mixing solution, which will be detailed later in this article. We find that through a second-order Taylor expansion of the mixing solution, we can derive accurate approximations to European put option prices in the aforementioned models. In addition, our methodology works naturally with time-dependent parameters. This is seen as a major advantage compared to other methodologies, for example transform methods, which cannot handle time-dependent parameters well. Furthermore, the approximation formulas are written solely in terms of iterated integrals, which we show obey convenient recursive properties when parameters are assumed to be piecewise-constant. This results in the approximation being essentially closed-form, and moreover leads to a fast calibration scheme. Our method builds upon [14], in which the Heston model is considered with constant parameters. A sensitivity analysis is performed in order to assess our approximation numerically. Furthermore, we present mathematical bounds on the error term in terms of higher order moments of the underlying variance process.

It has been well established that volatility is highly dependent on the strike and maturity of European option contracts. This phenomenon is called the volatility smile or skew, an attribute the well known Black-Scholes model [6] fails to take into account, see for example [16] for more on this. In response, there have been a number of frameworks proposed to explain the volatility smiles and skews observed in the market, for example, local volatility models, stochastic volatility models, and local-stochastic volatility models. In particular, stochastic volatility models have been introduced, where the volatility itself is a stochastic process possibly correlated with the spot. However with this added complexity, often option prices cannot be computed in a closed-form fashion. This is detrimental, as closed-form solutions lead to rapid option pricing, an attribute necessary for fast calibration of financial models. Without a closed-form solution, option pricing must be done numerically via Monte-Carlo or PDE methods, both methods being computationally costly.

If we assume that the characteristic function of the log-spot is known explicitly, then the option price can be computed quasi-explicitly (meaning in terms of at most one-dimensional complex integrals, where the integrands are explicit function), albeit under the restrictive assumption of constant or piecewise-constant parameters [18, 8, 29]. One class of models where this occurs is the class of affine models such as the Heston and Schöbel-Zhu model. However, for non-affine models, the characteristic function of the log-spot is rarely known explicitly, and such a procedure will not be effective. Non-affine models, although usually intractable compared to their affine counterparts, are often far more realistic. This has been shown in a number of studies, see for example [15, 10, 21]. For these reasons, numerical procedures such as PDE and Monte-Carlo methods have been substantially developed in the literature [34, 2, 22].

Closed-form approximations are an alternative methodology for option pricing, where the option price is approximated by a closed-form expression. The main advantages are that the option price can be computed rapidly and since transform methods are not used, time-dependent parameters can usually be handled well. One motivation for quick option pricing formulas is calibration, where the option price must be computed several times within an optimisation procedure.

There have been a plethora of results on closed-form expansions in the literature. For example, [28] derive a general closed-form expression for the price of an option via a PDE approach, as well as its corresponding implied volatility. [17] use singular perturbation techniques to obtain an explicit approximation for the option price and implied volatility in their SABR model. [1] show that from the mixing solution, one can approximate the put option price by decomposing it into a sum of two terms, one being completely correlation independent and the other dependent on correlation. However, neither terms are explicit. Furthermore, similar to our work, [4, 3] show that under the assumption of small correlation, an expansion can be performed with respect to the mixing solution, where the resulting expectations can be computed using Malliavin calculus techniques. Similarly, in the case of the time-dependent Heston model, [5] consider the mixing solution and expand around vol-of-vol, performing a combination of Taylor expansions and computing the resulting terms via Malliavin calculus techniques.

Stochastic volatility models usually either model the volatility directly, or indirectly via the variance process. A critical assumption is that volatility or variance has some sort of mean reversion behaviour, and this is supported by empirical evidence, see for example [16]. Specifically, for modelling the variance, a large class of stochastic volatility models is given by

[TABLE]

whereas for modelling the volatility, this class is of the form

[TABLE]

for some $\tilde{\mu},\hat{\mu}$ and $\mu\in\mathbb{R}$.222There exist other classes of stochastic volatility models, for example the exponential Ornstein-Uhlenbeck model [35] is not a part of either of these classes. Our model formulation here is for FX market purposes, but can be easily adjusted for equity and fixed income markets. In this article, we will focus on this class of stochastic volatility models. Some popular models in the literature include:

Model Variance/Volatility Dynamics of $V$ $\hat{\mu}$ $\tilde{\mu}$ $\mu$

Heston [18] Variance $\mathrm{d}V_{t}=\kappa_{t}(\theta_{t}-V_{t})\mathrm{d}t+\lambda_{t}\sqrt{V_{t}}\mathrm{d}B_{t}$ 0 1 1/2

Schöbel/Zhu [33] Volatility $\mathrm{d}V_{t}=\kappa_{t}(\theta_{t}-V_{t})\mathrm{d}t+\lambda_{t}\mathrm{d}B_{t}$ 0 1 0

GARCH [30, 36] Variance $\mathrm{d}V_{t}=\kappa_{t}(\theta_{t}-V_{t})\mathrm{d}t+\lambda_{t}V_{t}\mathrm{d}B_{t}$ 0 1 1

Inverse-Gamma [24] Volatility $\mathrm{d}V_{t}=\kappa_{t}(\theta_{t}-V_{t})\mathrm{d}t+\lambda_{t}V_{t}\mathrm{d}B_{t}$ 0 1 1

3/2 Model [25] Variance $\mathrm{d}V_{t}=\kappa_{t}(\theta_{t}V_{t}-V_{t}^{2})\mathrm{d}t+\lambda_{t}V_{t}^{3/2}\mathrm{d}B_{t}$ 1 2 3/2

Verhulst [26, 9] Volatility $\mathrm{d}V_{t}=\kappa_{t}(\theta_{t}V_{t}-V_{t}^{2})\mathrm{d}t+\lambda_{t}V_{t}\mathrm{d}B_{t}$ 1 2 1

Note that the Verhulst model also goes by the names Logistic model and XGBM model.

This article is dedicated to detailing how a second-order expansion of the so-called mixing solution for stochastic volatility models with time-dependent parameters can result in an approximation for the price of a European put option. This approximation is always in terms of some specific iterated integrals. We then show that these iterated integrals obey a convenient recursive scheme. Not only does this lead to the approximation being essentially closed-form, but this also results in a fast calibration scheme. The tractability of our methodology relies largely on the dynamics of the underlying variance process. Our methodology extends that of [14], in which the Heston model is considered with constant parameters. Our contribution is the following. We apply the expansion methodology in the Heston and GARCH diffusion stochastic volatility models with time-dependent parameters. In particular, the GARCH diffusion model is an example of a non-affine model, which practitioners prefer, see for example [15, 10, 21]. We include a robust error analysis, design a fast calibration scheme and give extensive numerical results. The sections are organised as follows:

• Section 2 details preliminary calculations, where we express the put option price as the mixing solution. Once done, a second-order Taylor expansion is performed, expressing the approximation formula in terms of a number of expectations of functionals of the underlying variance process.
• Section 3 details how to derive more convenient expressions for these expectations obtained in Section 2 through change of measure techniques. Specifically, we rewrite the spot $S_{T}$ as a convenient expression so as to construct a term which is a Doléans-Dade exponential, thereby defining a Radon-Nikodym derivative. This term allows us to change measure, allowing for more convenient calculations.
• Section 4 introduces specific models. As precise dynamics are assumed, the objective is to derive expressions in terms of specific iterated integrals for the expectations from Section 3. In particular, we consider the Heston and GARCH diffusion models.
• In Section 5 we perform an error analysis, bounding the error in the expansion in terms of higher order moments of the underlying variance process.
• In Section 6, under the assumption of piecewise-constant parameters, we obtain closed-form expressions for our pricing formulas. Moreover, we design a fast calibration scheme. Specifically, we show that the iterated integrals in which the pricing formulas in Section 4 are in terms of satisfy some convenient recursive properties.
• Section 7 is dedicated to a numerical error and sensitivity analysis for the Heston and GARCH diffusion models.

2. Preliminaries

Suppose the spot $S$ with variance $\sigma$ follows the dynamics

[TABLE]

where $W$ and $B$ are Brownian motions with deterministic, time-dependent instantaneous correlation $(\rho_{t})_{0\leq t\leq T}$, defined on the filtered probability space $(\Omega,\mathcal{F},(\mathcal{F}_{t})_{0\leq t\leq T},\mathbb{Q})$. Here $T$ is a finite time horizon, where $(r_{t}^{d})_{0\leq t\leq T}$ and $(r_{t}^{f})_{0\leq t\leq T}$ are the deterministic, time-dependent domestic and foreign interest rates respectively. In addition, $\rho_{t}\in[-1,1]$ for any $t\in[0,T]$. Furthermore, $(\mathcal{F}_{t})_{0\leq t\leq T}$ is the filtration generated by $(W,B)$ which satisfies the usual assumptions. In the following, $\mathbb{E}(\cdot)$ denotes the expectation under $\mathbb{Q}$, where $\mathbb{Q}$ is a domestic risk-neutral measure which we assume to be chosen.

Assumption 2.1**.**

$\sigma$ possesses a non-negative pathwise unique strong solution which does not blow up in finite time.

We will enforce 2.1 for the rest of this article. We briefly remark that sufficient conditions which ensure pathwise uniqueness of the solution are the usual Itô conditions (that is, Lipschitz continuity on the drift and diffusion coefficients) or the Yamada-Watanabe condition (Theorem 1 in [37]). Moreover, in order to ensure that the solution does not blow up in finite time, linear growth conditions on the drift and diffusion coefficients are sufficient.

Remark 2.1** (Geometric Brownian motion process).**

Let $\tilde{B}$ be an arbitrary Brownian motion process. We call $Y$ a GBM$(y_{0};\mu_{t},\nu_{t})$ if it solves the SDE

[TABLE]

where $\mu$ and $\nu$ are adapted to the Brownian filtration and satisfy some regularity conditions (for example, $\mu$ and $\nu$ bounded on $[0,T]$ is sufficient).

We decompose the Brownian motion $W$ as $W_{t}=\int_{0}^{t}\rho_{u}\mathrm{d}B_{u}+\int_{0}^{t}\sqrt{1-\rho_{u}^{2}}\mathrm{d}Z_{u}$, where $Z$ is a Brownian motion under $\mathbb{Q}$ such that $B$ and $Z$ are independent. Then, noticing $S$ is a GBM$(S_{0};r_{t}^{d}-r_{t}^{f},\sqrt{\sigma_{t}})$, we obtain the expression

[TABLE]

where

[TABLE]

Assumption 2.2**.**

$\sup_{t\in[0,T]}\limsup_{x\to\infty}\frac{\rho_{t}\beta(t,x)\sqrt{x}+\alpha(t,x)}{x}<\infty.$

We will enforce 2.2 for the rest of this article.

Lemma 2.1**.**

$\xi$ is a martingale.

Proof.

The proof is inspired by arguments made in Theorem 2.4 (i) of [27]. First notice that

[TABLE]

Let $\tau_{n}:=\inf\{t\geq 0:\sigma_{t}>n\}$ be the first time $\sigma$ crosses the level $n$. Then $\tau_{n}\uparrow\infty$ $\mathbb{Q}$ a.s. Define $\xi_{t}^{n}:=\xi_{t\wedge\tau_{n}}$. Then by solving eq. 2.1 up to $\tau_{n}$, this yields the pathwise unique strong solution

[TABLE]

Since the integrand of the above Itô integral is bounded, $(\xi_{t}^{n})_{t}$ is a martingale for each $n\in\mathbb{N}$. Furthermore, $\xi$ is a non-negative local-martingale. Hence by Fatou’s lemma, we have that $\xi$ is a non-negative supermartingale. Our goal now is to show that

[TABLE]

since then, by the Vallée Poussin theorem, $(\xi_{T}^{n})_{n}$ is uniformly integrable and thus ${\mathbb{E}(\xi_{T})=1}$. Combined with the fact that $\xi$ is a non-negative supermartingale, this will then imply that $\xi$ is a martingale. Now,

[TABLE]

Clearly $\xi_{T}^{n}$ is a Radon-Nikodym derivative which defines a measure $\hat{\mathbb{Q}}$. Hence by Girsanov’s theorem, $\hat{B}_{t}=B_{t}-\int_{0}^{t}\rho_{u}\sqrt{\sigma_{u}}\textbf{1}_{\{u\leq\tau_{n}\}}\mathrm{d}u$ is a Brownian motion under $\hat{\mathbb{Q}}$. Denote the expectation under $\hat{\mathbb{Q}}$ by $\hat{\mathbb{E}}$. Then

[TABLE]

So it suffices to determine a bound on $\hat{\mathbb{E}}(\sigma_{t})$ for $t\leq T$. Under $\hat{\mathbb{Q}}$, we have

[TABLE]

Now as a consequence of 2.2, this implies that there exists some constant $M>0$ such that

[TABLE]

for all $x\geq 0$, uniformly in $t\leq T$. Utilising eq. 2.2 yields

[TABLE]

Hence

[TABLE]

Define $m_{t}:=\hat{\mathbb{E}}(\sigma_{t})$. After redefining constants, we obtain the following integral inequality,

[TABLE]

Then by Gronwall’s inequality, we obtain

[TABLE]

∎

Let

[TABLE]

be the price of a put option on $S$. Denote by $(\mathcal{F}_{t}^{B})_{0\leq t\leq T}$ the filtration generated by the Brownian motion $B$, as well as $\mathcal{N}(\cdot)$ and $\phi(\cdot)$ the standard Normal distribution and density functions respectively.

Proposition 2.1**.**

Put can be expressed as

[TABLE]

where

[TABLE]

Proof.

This is a consequence of the mixing solution methodology, which is detailed in Appendix A. ∎

Proposition 2.2** (Second-order put option price approximation).**

The second-order put option price approximation, denoted by $\text{Put}^{(2)}$, is given by

[TABLE]

where $(\hat{x},\hat{y}):=(S_{0},\int_{0}^{T}(1-\rho_{t}^{2})\mathbb{E}(\sigma_{t})\mathrm{d}t)$.

Proof.

Notice that $\text{Put}_{\text{BS}}$ is a composition of smooth functions. Hence, $\text{Put}_{\text{BS}}$ is smooth on $(\mathbb{R}^{2}_{+};\mathbb{R}_{+})$. We Taylor expand $\text{Put}_{\text{BS}}$ around the mean of $\left(S_{0}\xi_{T},\int_{0}^{T}\sigma_{t}(1-\rho_{t}^{2})\mathrm{d}t\right)$ up to second-order. By Lemma 2.1, the expansion point is $(\hat{x},\hat{y})=(S_{0},\int_{0}^{T}(1-\rho_{t}^{2})\mathbb{E}(\sigma_{t})\mathrm{d}t)$. Thus

[TABLE]

Taking expectation yields a second-order approximation to the put option price, that is, $\text{Put}^{(2)}$. Notice that $\text{Put}_{\text{BS}}(\hat{x},\hat{y})$ is a deterministic quantity, thus the first-order terms will vanish. ∎

The partial derivatives of $\text{Put}_{\text{BS}}$ are analogous to the Black-Scholes Greeks, which are explicit and provided in Appendix B. In order to simplify $\text{Put}^{(2)}$, what remains to be done is the calculation of each of the expectations from Proposition 2.2, which are

[TABLE]

Remark 2.2** (Call options).**

When considering the pricing of call options, the approximation methodology is essentially the same as that of put options. It is clear that the mixing solution methodology from Appendix A can be easily adapted to the case of call options, in which case we obtain

[TABLE]

where

[TABLE]

Moreover, put-call parity states that

[TABLE]

Thus the second-order Call partial derivatives will be the same as the second-order Put partial derivatives. Comparing the mixing solution expressions for Call and Put, one can see that the only difference between the call and put option second-order approximation is the zero-th order term. For this reason, in this article, we need only consider the pricing of put options.

Remark 2.3** (Greeks).**

One can obtain second-order approximations of put option Greeks by simply differentiating the expression for $\text{Put}^{(2)}$. These expressions are provided in Appendix C.

3. Calculation of expectations

Lemma 3.1**.**

Let $(\mathbb{Q}_{n})_{n\geq 0}$ be a sequence of probability measures equivalent to $\mathbb{Q}$, defined by the Radon-Nikodym derivatives

[TABLE]

where $\mathbb{Q}_{0}:=\mathbb{Q}$ and $B^{0}:=B$. Under $\mathbb{Q}_{n}$, $B_{t}^{n}:=B_{t}^{n-1}-\int_{0}^{t}\rho_{u}\sqrt{\sigma_{u}}\mathrm{d}u$ is a Brownian motion.

Proof.

This is a clear consequence of Girsanov’s theorem. ∎

Remark 3.1**.**

Through Lemma 3.1, we have the following relationships:

[TABLE]

and

[TABLE]

These relationships allow for alternative and often more convenient calculations of expectations under $\mathbb{Q}$.

Using Remark 3.1, we can now give alternative expressions for the expectations seen in eqs. 2.6 to 2.8.

3.1. $\mathbb{E}(\xi_{T}-1)^{2}$

First, expanding eq. 2.6 gives

[TABLE]

This second moment can be dealt with a number of changes of measures.

[TABLE]

Under the assumption of constant parameters333That is, $\alpha(t,\sigma_{t})=\alpha(\sigma_{t})$, $\beta(t,\sigma_{t})=\beta(\sigma_{t})$ and $\rho_{t}=\rho$. we may calculate eq. 3.3 explicitly via the Laplace transform for certain processes $\sigma$. However to our knowledge, there exists no explicit solution when parameters are time-dependent, see [20]. Instead, we approximate eq. 3.3 by expanding the exponential around the mean of $\int_{0}^{T}\rho_{t}^{2}\sigma_{t}\mathrm{d}t$ to second-order.

[TABLE]

where we have used the fact that $\left(\int_{0}^{T}f(t)\mathrm{d}t\right)^{2}=2\int_{0}^{T}f(t)\left(\int_{0}^{t}f(s)\mathrm{d}s\right)\mathrm{d}t$.

3.2. $\mathbb{E}\left(\int_{0}^{T}(1-\rho_{t}^{2})(\sigma_{t}-\mathbb{E}(\sigma_{t}))\mathrm{d}t\right)^{2}$

To calculate eq. 2.7, we use the same approach from Section 3.1.

[TABLE]

3.3. $\mathbb{E}\left\{(\xi_{T}-1)\left(\int_{0}^{T}(1-\rho_{t}^{2})(\sigma_{t}-\mathbb{E}(\sigma_{t}))\mathrm{d}t\right)\right\}$

Calculation of the mixed expectation eq. 2.8 gives

[TABLE]

4. Pricing for specific models

We now introduce specific dynamics for both the spot and its underlying variance process. From Section 3, it is apparent that simplifying the expression for $\text{Put}^{(2)}$ will largely depend on the tractability of the variance process $\sigma$ under the original measure $\mathbb{Q}$, as well as the artificial measures $\mathbb{Q}_{1}$ and $\mathbb{Q}_{2}$.

Remark 4.1** (Notation for key processes).**

Let $(\tilde{\kappa}_{t})_{0\leq t\leq T},(\tilde{\theta}_{t})_{0\leq t\leq T}$ and $(\tilde{\lambda}_{t})_{0\leq t\leq T}$ be time-dependent, deterministic, strictly positive and bounded on $[0,T]$. Let $\tilde{B}$ be an arbitrary Brownian motion. We will utilise the following terminology:

1. (1)

If $\tilde{V}$ solves

[TABLE]

then we call $\tilde{V}$ a CIR$(\tilde{v}_{0};\tilde{\kappa}_{t},\tilde{\theta}_{t},\tilde{\lambda}_{t})$. 2. (2)

If $\tilde{V}$ solves

[TABLE]

then we call $\tilde{V}$ an IGa$(\tilde{v}_{0};\tilde{\kappa}_{t},\tilde{\theta}_{t},\tilde{\lambda}_{t})$.

We point the reader towards Appendix D for further information on these processes.

4.1. Heston model

Suppose the spot $S$ with variance $V$ follows the Heston dynamics

[TABLE]

where $(\kappa_{t})_{0\leq t\leq T},(\theta_{t})_{0\leq t\leq T}$ and $(\lambda_{t})_{0\leq t\leq T}$ are time-dependent, deterministic, strictly positive and bounded on $[0,T]$.

It is clear that the variance process $V$ in eq. 4.1 is a CIR$(v_{0};\kappa_{t},\theta_{t},\lambda_{t})$.

Remark 4.2**.**

Note that the Heston model satisfies 2.2, since

[TABLE]

Hence by Lemma 2.1, $\xi$ is a martingale.

Lemma 4.1**.**

Let $(\mathbb{Q}_{n})_{n\geq 0}$ be a sequence of probability measures equivalent to $\mathbb{Q}$, defined by the Radon-Nikodym derivatives

[TABLE]

where $\mathbb{Q}_{0}:=\mathbb{Q}$ and $B^{0}:=B$. Under $\mathbb{Q}_{n}$, $B_{t}^{n}:=B_{t}^{n-1}-\int_{0}^{t}\rho_{u}\sqrt{V_{u}}\mathrm{d}u$ is a Brownian motion. For $n\geq 0$, the dynamics of $V$ under the measure $\mathbb{Q}_{n}$ are

[TABLE]

which is a CIR$(v_{0};\kappa_{t}-n\lambda_{t}\rho_{t},\frac{\theta_{t}\kappa_{t}}{\kappa_{t}-n\lambda_{t}\rho_{t}},\lambda_{t})$.

Proof.

This lemma is simply obtained through Lemma 3.1, then expressing $V$ under the new measures $\mathbb{Q}_{n}$. ∎

Thus the variance process $V$ is a CIR process under all measures considered.

4.1.1. Pricing under the Heston framework

The second-order approximation of the put option price in the Heston framework is given by the following theorem.

Theorem 4.1** (Second-order Heston put option price).**

Let $\hat{x}=S_{0}$ and

[TABLE]

The second-order approximation to the put option price in the Heston model, denoted by $\text{Put}_{\text{H}}^{(2)}$, is

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

Proof.

Use Proposition 2.2 and adapt Section 3 to the Heston framework. The moments of $V$ correspond to the moments of a CIR process, which are obtained from Section D.1. ∎

4.2. GARCH diffusion model

Suppose the spot $S$ with variance $V$ follows the GARCH diffusion dynamics

[TABLE]

where $(\kappa_{t})_{0\leq t\leq T},(\theta_{t})_{0\leq t\leq T}$ and $(\lambda_{t})_{0\leq t\leq T}$ are time-dependent, deterministic, strictly positive and bounded on $[0,T]$. It is evident that the variance process $V$ in eq. 4.3 is an IGa$(v_{0};\kappa_{t},\theta_{t},\lambda_{t})$.

Remark 4.3**.**

Unlike the Heston model, the GARCH diffusion model does not satisfy 2.2, since

[TABLE]

Thus, there is no guarantee that $\xi$ is indeed a martingale. Instead, we will now assume that $\xi$ is a martingale and show that, even if it is indeed a martingale, it would provide no benefit in our expansion procedure.

Lemma 4.2**.**

Let $(\mathbb{Q}_{n})_{n\geq 0}$ be a sequence of probability measures equivalent to $\mathbb{Q}$, defined by the Radon-Nikodym derivatives

[TABLE]

where $\mathbb{Q}_{0}:=\mathbb{Q}$ and $B^{0}:=B$. Under $\mathbb{Q}_{n}$, $B_{t}^{n}:=B_{t}^{n-1}-\int_{0}^{t}\rho_{u}\sqrt{V_{u}}\mathrm{d}u$ is a Brownian motion. For $n\geq 0$, the dynamics of $V$ under the measure $\mathbb{Q}_{n}$ are

[TABLE]

Proof.

This lemma is simply obtained through Lemma 3.1, then expressing $V$ under the new measures $\mathbb{Q}_{n}$. ∎

Remark 4.4**.**

Let $\mathbb{Q}_{n}$ be defined as in Lemma 4.2. Under the measures $\mathbb{Q}_{n}$, $n\geq 1$, $V$ has no known expression for its solution, nor known expression for its moments.

Validity of Remark 4.4.

The SDE in Lemma 4.2 is a linear diffusion type SDE. From Appendix E, it is known that if an explicit solution exists, it is given by

[TABLE]

where $F$ is a GBM$(1;\lambda_{t}^{2},-\lambda_{t})$ and $Y$ is the solution to the integral equation (written in differential form)

[TABLE]

Define $A_{t}:=\kappa_{t}\theta_{t}F_{t}$ and $C_{t}:=\frac{n\lambda_{t}\rho_{t}}{\kappa_{t}}F_{t}^{-1/2}$. Note that $A_{t}$ and $C_{t}$ are both non-differentiable in $t$. Thus

[TABLE]

As far as we know, there is no explicit solution to these types of integral equations in the literature, even when $A$ and $C$ are constants. As for explicit moments, it is unclear how to approach this problem. There seems to be no approach to this problem in the literature, especially in the case of time-dependent parameters, see for example [23] Chapter 4.4 for a comprehensive list of explicitly solvable SDEs. Furthermore, as an explicit solution does not exist, we cannot use the method of approximating moments via the SDE’s solution. ∎

4.2.1. Pricing under the GARCH diffusion framework: $\rho=0$

In the GARCH diffusion model, 2.2 is violated, and thus there is no guarantee $\xi$ is a martingale. Additionally, assuming $\xi$ is a martingale is not helpful, as the change of measure technique gives an intractable dynamic for $V$; we cannot appeal to it for calculating expectations. However, in the case of $\rho=0$ a.e., this implies $\xi_{T}=1$ $\mathbb{Q}$ a.s., and one will notice that the terms in the expansion requiring a change of measure will disappear. Of course, the cost is the additional restrictive assumption that spot and volatility movements are uncorrelated. We hope to mitigate this issue in future work by combining this approach with small correlation expansion methods, see [4, 3].

Theorem 4.2** (Second-order GARCH put option price).**

Assume $\rho=0$ a.e. Let $\hat{x}=S_{0}$ and $\hat{y}=\int_{0}^{T}\mathbb{E}(V_{t})\mathrm{d}t$. Then the second-order put option price in the GARCH diffusion model, denoted by $\text{Put}^{(2)}_{\text{GARCH}}$, is

[TABLE]

Both $\mathrm{Cov}(V_{s},V_{t})$ and $\mathbb{E}(V_{t})$ are given in Section D.2.

Proof.

Use Proposition 2.2 under the assumption of $\rho=0$ a.e. Then notice that

[TABLE]

∎

5. Error analysis

We present an explicit bound on the error term in our expansion in terms of higher order moments of the corresponding variance process. Specifically, this means bounding the remainder term in the second-order expansion of the function Put${}_{\text{BS}}$, and for the case when $\rho\neq 0$, the error term associated with the expansion of $e^{\int_{0}^{T}\rho_{u}^{2}\sigma_{u}\mathrm{d}u}$.

We will need explicit expressions for the error terms. These are given by Taylor’s theorem, which we will present here to fix notation. We only consider the results up to second-order.

Theorem 5.1** (Taylor’s theorem for $f:\mathbb{R}\to\mathbb{R}$).**

Let $A\subseteq\mathbb{R},B\subseteq\mathbb{R}$ and $f:A\to B$ be a $C^{3}$ function in a closed interval about the point $a\in A$. Then

[TABLE]

where the remainder $R$ is given by

[TABLE]

We will prefer the integration bounds to be from [math] to $1$ rather than $a$ to $x$, as $a$ and $x$ will correspond to random variables.

Theorem 5.2** (Taylor’s theorem for $g:\mathbb{R}^{2}\to\mathbb{R}$).**

Let $A\subseteq\mathbb{R}^{2},B\subseteq\mathbb{R}$ and $g:A\to B$ be a $C^{3}$ function in a closed ball about the point $(a,b)\in A$. Then

[TABLE]

where the remainder $R$ is given by

[TABLE]

with $\alpha:=(\alpha_{1},\alpha_{2})$ and $|\alpha|:=\alpha_{1}+\alpha_{2}$.

5.1. Expression for error term

The representation for the total error generated by the expansion can be summarised by the following theorem.

Theorem 5.3** (Total expansion error).**

As a functional of the underlying variance process $\sigma$, let $\mathcal{E}_{\text{BS}}(\sigma)$ and $\mathcal{\tilde{E}}(\sigma)$ correspond to the error induced by the expansion of $\text{Put}_{\text{BS}}$ and $e^{\int_{0}^{T}\rho_{u}^{2}\sigma_{u}\mathrm{d}u}$ respectively. The error generated by Taylor expansions for a general variance process $\sigma$ is given by

[TABLE]

where

[TABLE]

and

[TABLE]

Proof.

First, we deal with the error term associated with the function Put${}_{\text{BS}}$, that is, $\mathcal{E}_{\text{BS}}(\sigma)$. Recall the expansion of Put${}_{\text{BS}}$ around the point $(\hat{x},\hat{y}):=(S_{0},\int_{0}^{T}\mathbb{E}(\sigma_{t})(1-\rho_{t}^{2})\mathrm{d}t)$ evaluated at $(S_{0}\xi_{T},\int_{0}^{T}\sigma_{t}(1-\rho_{t}^{2})\mathrm{d}t)$ for a general variance process $\sigma$:

[TABLE]

Using Theorem 5.2 for the function Put${}_{\text{BS}}$, this gives the error term as

[TABLE]

We now investigate the error term associated with the calculation of $\mathbb{E}\xi_{T}^{2}$, that is, $\mathcal{\tilde{E}}(\sigma)$. Let us look at this term without the expectation.

[TABLE]

We expand $e^{\int_{0}^{T}\rho_{u}^{2}\sigma_{u}\mathrm{d}u}$ around the expectation of the exponential’s argument under $\mathbb{Q}_{2}$. Note that $\xi_{T}^{2}e^{-\int_{0}^{T}\rho_{u}^{2}\sigma_{u}\mathrm{d}u}$ is the Radon-Nikodym derivative which changes measure from $\mathbb{Q}$ to $\mathbb{Q}_{2}$. Expanding to second-order gives

[TABLE]

Finally, the coefficient in front of $\xi_{T}^{2}$ in the pricing formula is $\frac{1}{2}\partial_{xx}\text{Put}_{\text{BS}}(\hat{x},\hat{y})S_{0}^{2}$. Thus, the error term $\mathcal{\tilde{E}}(\sigma)$ can be written as

[TABLE]

∎

Corollary 5.1** (Total expansion error: $\rho=0$).**

The error due to Taylor expansions for a general variance process $\sigma$ with $\rho=0$ a.e., denoted by $\mathcal{E}_{0}(\sigma)$, is given by

[TABLE]

which is just $\mathcal{E}_{\text{BS}}(\sigma)$ when $\rho=0$ a.e.

5.2. Bounding error term

The hope now is to be able to bound $\mathbb{E}(\mathcal{E}(\sigma))$ in terms of the higher order moments of the variance process $\sigma$. To do this, we will show that the partial derivatives $\frac{\partial^{3}\text{Put}_{\text{BS}}}{\partial x^{\alpha_{1}}\partial y^{\alpha_{2}}}(F(u),G(u))$ for $u\in(0,1)$ where $\alpha_{1}+\alpha_{2}=3$, are functions of $T$ and $K$ which are bounded on $\mathbb{R}^{2}_{+}$.

Lemma 5.1**.**

Consider the third-order partial derivatives of $\text{Put}_{\text{BS}}$, $\frac{\partial^{3}\text{Put}_{\text{BS}}}{\partial x^{\alpha_{1}}\partial y^{\alpha_{2}}}$, where $\alpha_{1}+\alpha_{2}=3$. Let $\hat{f}(x):=\ln(x/K)+\int_{0}^{T}\left(r_{t}^{d}-r_{t}^{f}\right)\mathrm{d}t$. Then

[TABLE]

Furthermore, this is the only case where the partial derivatives explode.

Proof.

In the following, we will repeatedly let $F$ be an arbitrary polynomial of some degree, as well as $A$ to be an arbitrary constant. That is, they may be different on each use. From Appendix B, it can be seen that as a function of $x$ and $y$, the third-order partial derivatives are of the form

[TABLE]

where we recall

[TABLE]

Written in this form eq. 5.1, it is evident that the partial derivatives could only blow up if either $x$ or $y$ tend to [math] or infinity. We need only look at these limits independently of the other variable. In the following, we shall say $f=o(g)$ if and only if $\lim_{x\to\infty}\frac{f(x)}{g(x)}=0$ and $f=o_{0}(g)$ if and only if $\lim_{x\downarrow 0}\frac{f(x)}{g(x)}=0$.

1. (1)

For fixed $x$: From eq. 5.1 the partial derivatives are of the form $A\frac{\phi(d_{+})F(d_{+},d_{-},\sqrt{y})}{y^{m/2}}$. It can be shown that

[TABLE]

where $D_{2}=\frac{1}{2}(\hat{f}(x))^{2}$ and $D_{1}=1/8$. Hence both $D_{2}$ and $D_{1}$ are non-negative. However, there will be two cases to consider, when $D_{2}>0$ or $D_{2}=0$.

1. (a)

Suppose $D_{2}>0$, then $\hat{f}(x)\neq 0$. As $F$ is a polynomial in $d_{+},d_{-}$, and $\sqrt{y}$, we can say that $F(d_{+},d_{-},\sqrt{y})=o_{0}(1/y^{M_{0}/2})$ and $F(d_{+},d_{-},\sqrt{y})=o(y^{M/2})$ for some $M,M_{0}\in\mathbb{N}$. Thus

[TABLE]

and also

[TABLE]

Then as $y\downarrow 0$ or $y\to\infty$, the partial derivatives tend to [math]. 2. (b)

Suppose $D_{2}=0$, then $\hat{f}(x)=0$. Thus $d_{+}=\frac{1}{2}\sqrt{y}$ and $\phi(d_{+})=Ae^{-D_{1}y}$. Evidently, $F(d_{+},d_{-},\sqrt{y})=\sum_{i=0}^{N}C_{i}y^{i/2}$ for some $N\in\mathbb{N}$ and constants $C_{0},\dots,C_{N}$. Thus

[TABLE]

This quantity tends to 0 as $y\to\infty$, as the exponential decay makes the polynomial growth/decay irrelevant. However, when $y\downarrow 0$, then this limit depends on the polynomial $F$. If $N>m$ and if one of the $C_{0},C_{1},\dots,C_{m}$ are non-zero then this quantity tends to $\infty$ as $y\downarrow 0$. If $N<m$, then this quantity tends to $\infty$ as $y\downarrow 0$. For each of the partial derivatives, it can be shown that either of these cases are satisfied. Thus the partial derivatives tend to $\infty$ when $y\downarrow 0$.

To conclude, for fixed $x$, if $\hat{f}(x)=0$, then the partial derivatives tend to [math] if $y\to\infty$, and to $\infty$ if $y\downarrow 0$. When $\hat{f}(x)\neq 0$, the partial derivatives tend to [math] if $y\downarrow 0$ or $y\to\infty$. 2. (2)

For fixed $y$: From eq. 5.1, the partial derivatives are of the form $\frac{\phi(d_{+})F(d_{+},d_{-})}{x^{n}}$. It can be shown that

[TABLE]

where $E_{2}>0$ and $E_{1}\in\mathbb{R}$. Furthermore, as $F$ is a polynomial in $d_{+}$ and $d_{-}$, then

[TABLE]

for some $N\in\mathbb{N}$ and constants $C_{0},\dots C_{N}$. It is clear $F=o(x)$ and $F=o_{0}(\ln^{N+1}(x))$. Thus

[TABLE]

Then as $x\to\infty$, this quantity tends to [math]. In addition

[TABLE]

Then as $x\downarrow 0$ this quantity tends to 0. So for fixed $y$, the partial derivatives tend to 0 as $x\downarrow 0$ or $x\to\infty$.

∎

We will however, be concerned with the behaviour of $u\mapsto\frac{\partial^{3}}{\partial x^{\alpha_{1}}\partial y^{\alpha_{2}}}\text{Put}_{\text{BS}}(F(u),G(u))$, meaning we will have to consider both arguments simultaneously, as they are both linear functions of $u$.

Lemma 5.2**.**

Consider the third-order partial derivatives of $\text{Put}_{\text{BS}}$, $\frac{\partial^{3}\text{Put}_{\text{BS}}}{\partial x^{\alpha_{1}}\partial y^{\alpha_{2}}}$, where $\alpha_{1}+\alpha_{2}=3$ as well as the linear functions $h_{1},h_{2}:[0,1]\to\mathbb{R}_{+}$ such that $h_{1}(u)=u(d_{1}-c_{1})+c_{1}$ and $h_{2}(u)=u(d_{2}-c_{2})+c_{2}$. Assume there exists no point $a\in(0,1)$ such that

[TABLE]

Then there exists functions $M_{\alpha}$ bounded on $\mathbb{R}_{+}^{2}$ such that

[TABLE]

Furthermore, the behaviour of $M_{\alpha}$ for fixed $K$ and $T$ is characterised by the functions $\zeta$ and $\eta$ respectively, where

[TABLE]

with $\tilde{r}(T):=\int_{0}^{T}(r_{t}^{d}-r_{t}^{f})\mathrm{d}t$ and $E_{2}>0$, $E_{1}\in\mathbb{R},\hat{A}\in\mathbb{R}$, $n\in\mathbb{N}$ and $c_{0},\dots,c_{n}$ are constants, and

[TABLE]

with $D_{2}>0,D_{1}\in\mathbb{R},\tilde{A}\in\mathbb{R},N\in\mathbb{N}$ and $C_{0},\dots,C_{N}$ are constants.

Proof.

Under the assumptions of $h_{1}$ and $h_{2}$, this supremum will be bounded by a direct application of Lemma 5.1. Next, we need to show that $M_{\alpha}$ are bounded on $\mathbb{R}_{+}^{2}$ and behave as $\zeta$ and $\eta$ for fixed $K$ and $T$ respectively.

In the following, $A$ denotes an arbitrary constant and $F$ an arbitrary polynomial of some degree. They may be different on each use.

1. (1)

Behaviour in $T$: Fix all variables as constant except $T$. Then we can write the partial derivatives as

[TABLE]

Expanding and collecting terms with $T$, we can write the partial derivatives in the form

[TABLE]

where $E_{2}>0$, $E_{1}\in\mathbb{R}$, $n\in\mathbb{N}$ and $c_{0},\dots,c_{n}$ are constants. As $\zeta$ is a composition of polynomials and exponentials of $\tilde{r}(T)$, then it is bounded for any closed interval not containing [math]. Now since $\sup_{t\in[0,T]}(|r_{t}^{d}-r_{t}^{f}|)=:R<1$ then $|\tilde{r}(T)|<RT$. Thus $\tilde{r}^{i}(T)=o(T^{i})$ and $\tilde{r}^{i}(T)=o_{0}(T^{i})$. Hence $\zeta$ tends to [math] as $T\downarrow 0$ or $T\to\infty$. Thus $\zeta$ is bounded on $\mathbb{R}_{+}$. 2. (2)

Behaviour in $K$: Now fix all variables as constant except $K$. Then the partial derivatives can be written as

[TABLE]

Expanding and collecting terms with $K$, the partial derivatives can be written in the form

[TABLE]

where $D_{2}>0$ and $D_{1}\in\mathbb{R}$. Then writing out the polynomial explicitly

[TABLE]

where $N\in\mathbb{N}$ and $C_{0},\dots,C_{N}$ are constants. Thus

[TABLE]

$\eta$ is bounded for any closed interval not containing [math] since it is a composition of exponentials and logarithms. Then as $\ln^{i}(K)=o(K)$ and $\ln^{i}(K)=o_{0}(\ln^{N+1}(K))$, $\eta$ tends to [math] as $K\downarrow 0$ or $K\to\infty$. Thus $\eta$ is bounded on $\mathbb{R}_{+}$.

∎

Proposition 5.1**.**

There exists functions $M_{\alpha}$ as in Lemma 5.2 such that

[TABLE]

Proof.

Since $F$ and $G$ are linear functions, then from Lemma 5.2, this claim is immediately true if we can show that $G$ is strictly greater than [math] for any $u\in(0,1)$ $\mathbb{Q}$ a.s. Recall

[TABLE]

$G$ corresponds to the linear interpolation of $\int_{0}^{T}(1-\rho_{t}^{2})\mathbb{E}(\sigma_{t})\mathrm{d}t$ and $\int_{0}^{T}(1-\rho_{t}^{2})\sigma_{t}\mathrm{d}t$. It is clear $\sup_{t\in[0,T]}(1-\rho_{t}^{2})>0$. As $\sigma$ corresponds to the variance process, by 2.1 it is a non-negative process. Hence, the set $\{t\in[0,T]:\sigma_{t}>0\}$ has non-zero Lebesgue measure. Thus, these integrals are strictly positive and hence $G$ is strictly larger than [math] for any $u\in(0,1)$ $\mathbb{Q}$ a.s. ∎

We obtain the following corollary.

Corollary 5.2**.**

There exists a function $M$ as in Lemma 5.2 such that

[TABLE]

Proof.

Recall

[TABLE]

Then by the same argument in the proof of Proposition 5.1, $\tilde{G}$ is strictly greater than [math] $\mathbb{Q}$ a.s. Hence by Lemma 5.2, the claim is true. ∎

Theorem 5.4** (Error bounds for general $\sigma$).**

The error term in the pricing formula is bounded as

[TABLE]

where $\tilde{M}(T,K)=\partial_{xx}\text{Put}_{\text{BS}}(\hat{x},\hat{y})$ is bounded on $\mathbb{R}_{+}^{2}$ and $C=1/12,C_{\alpha}=\frac{1}{3}\frac{|\alpha|}{\alpha_{1}!\alpha_{2}!}$ are constants, the latter depending on $\alpha$.

Proof.

First, by Proposition 5.1, we have that $E_{\alpha}(S_{0}\xi_{T},\int_{0}^{T}\sigma_{t}(1-\rho_{t}^{2}))\leq\frac{1}{3}M_{\alpha}(T,K)$. By Theorem 5.3, the error is decomposed as $\mathcal{E}(\sigma)=\mathcal{E}_{\text{BS}}(\sigma)+\mathcal{\tilde{E}}(\sigma)$. We will make use of the well known integral inequality

[TABLE]

For the term $\mathcal{E}_{\text{BS}}(\sigma)$, we have

[TABLE]

where we have used the integral inequality eq. 5.2. Applying the Cauchy-Schwarz inequality, we obtain

[TABLE]

where we have used the integral inequality eq. 5.2 for the second inequality.

For the term $\mathcal{\tilde{E}}(\sigma)$, notice that

[TABLE]

Now for $u\in(0,1)$, $e^{u\int_{0}^{T}\rho_{m}^{2}(\sigma_{m}-\mathbb{E}_{\mathbb{Q}_{2}}(\sigma_{m}))\mathrm{d}m}\leq e^{u\int_{0}^{T}\rho_{m}^{2}|\sigma_{m}-\mathbb{E}_{\mathbb{Q}_{2}}(\sigma_{m})|\mathrm{d}m}$. Thus

[TABLE]

Hence

[TABLE]

Thus

[TABLE]

Finally, using the Cauchy-Schwarz inequality and the integral inequality eq. 5.2, we obtain

[TABLE]

Furthermore, notice that $\partial_{xx}\text{Put}_{\text{BS}}(\hat{x},\hat{y})=\tilde{M}(T,K)$, where $\tilde{M}(T,K)$ is a function which behaves like $M(T,K)$. ∎

Corollary 5.3** (Error bounds for general $\sigma$: $\rho=0$).**

For $\rho=0$ a.e., the error term in the pricing formula is bounded as

[TABLE]

where $C=1/6$ is a constant.

Proof.

From Corollary 5.1, we have

[TABLE]

Notice that $(\int_{0}^{T}|f(t)|\mathrm{d}t)^{3}\leq T^{2}\int_{0}^{T}|f(t)|^{3}\mathrm{d}t$ and use Corollary 5.2. Then the result is immediate. ∎

Remark 5.1**.**

Theorem 5.4 and Corollary 5.3 present bounds on the error generated by the expansion procedure in terms of higher-order moments of the variance process. Note that there are two components in the expansion procedure that mainly contribute to the error, the third-order partial derivatives of $\text{Put}_{\text{BS}}$ and the higher-order moments of the variance process. This is a clear consequence of how the remainder term behaves in Taylor’s theorem; namely, the remainder is governed by the chosen function and the difference in the expansion and evaluation point. Proposition 5.1 and Corollary 5.2 demonstrate that $\frac{\partial^{3}\text{Put}_{\text{BS}}}{\partial x^{\alpha_{1}}\partial y^{\alpha_{2}}}(F(u),G(u))$ for $u\in(0,1)$ are dominated by functions such as $M(T,K)$, which itself are bounded on $\mathbb{R}^{2}_{+}$. The behaviour of $M(T,K)$ for a fixed $K$ and fixed $T$ is governed by the functions $\zeta(T)$ and $\eta(K)$ respectively. Interestingly, the error bound as a function of strike $K$ is purely governed by the function $\eta(K)$, and thus tends to zero for very small or very large strikes. However, the total behaviour in $T$ of the error bound depends not just on $\zeta(T)$, but on other terms as well, most notably the moments of the underlying variance process.

Indeed, in order to go a step further, one can consider a specific example of a variance process, and then attempt to bound the moments that appear in Theorem 5.4 and Corollary 5.3. However, for the purposes of application, it is perhaps more beneficial to instead consider a numerical sensitivity analysis in order to investigate the behaviour of the error with respect to parameter adjustments. Indeed, we have performed this in Section 7 for the Heston and GARCH diffusion models.

6. Closed-form formulas and fast calibration

In this section, under the assumption of piecewise-constant parameters, we provide closed-form formulas for the price of a put option as well as a fast calibration scheme. To do this, we appeal to the results from Section 4 on the approximation of prices of put options in the various models. We recognise that these prices are in terms of specific iterated integrals. Our goal is to show that these iterated integrals obey a convenient recursive property when parameters are piecewise-constant. To this end, define the integral operator

[TABLE]

In addition, we define the $n$-fold integral operator using the following recurrence:

[TABLE]

The rest of the section is dedicated to making these integral operators closed-form when parameters are piecewise-constant.

Let $\mathcal{T}=\{0=T_{0},T_{1},\dots,T_{N-1},T_{N}=T\}$, where $T_{i}<T_{i+1}$ be a collection of maturity dates on $[0,T]$, with $\Delta T_{i}:=T_{i+1}-T_{i}$ and $\Delta T_{0}\equiv 1$. When the dummy functions are piecewise-constant, that is, $l^{(n)}_{t}=l^{(n)}_{i}$ on $t\in[T_{i},T_{i+1})$ and similarly for $k^{(n)}$, then we can recursively calculate the integral operators eq. 6.1 and eq. 6.2. Define

[TABLE]

where $\gamma_{i}(u):=(u-T_{i})/\Delta T_{i}$ and $p\in\mathbb{N}\cup\{0\}$. In addition, define the $n$-fold extension of $\varphi_{T_{i},\cdot}^{(\cdot,\cdot)}$ recursively by

[TABLE]

where $p_{n}\in\mathbb{N}\cup\{0\}$. With the assumption that the dummy functions are piecewise-constant, we can obtain the integral operator at time $T_{i+1}$ expressed by terms at time $T_{i}$.

[TABLE]

The only terms here that are not explicit are the functions $e_{\cdot}^{(\cdot,\dots,\cdot)}$ and $\varphi_{T_{i},\cdot}^{(\cdot,\cdot),\dots,(\cdot,\cdot)}$. For $t\in(T_{i},T_{i+1}]$, we can derive the following:

[TABLE]

where $e_{0}^{(k^{(n)},\dots,k^{(1)})}=1$. Using integration by parts and basic integration properties, we obtain

[TABLE]

In addition, for $n\geq 2$,

[TABLE]

6.1. Closed-form prices

We now express the second-order put option prices under the Heston and GARCH diffusion models in terms of the integral operators eq. 6.1 and eq. 6.2. These are given by the following propositions. By doing so, we essentially obtain a closed-form expression for the put option prices when parameters are piecewise-constant. This is a consequence of the relationships derived in this section. Since we are essentially expressing the results from Section 4 with new notation, we omit the proofs for these propositions.

Proposition 6.1** (Second-order Heston explicit put option price).**

The second-order put option price in the Heston model can be written as

[TABLE]

where

[TABLE]

[TABLE]

Furthermore $\hat{x}=S_{0}$ and $\hat{y}=v_{0}\omega_{T}^{(-\kappa,1-\rho^{2})}+\omega_{T}^{(-\kappa,1-\rho^{2}),(\kappa,\kappa\theta)}.$

Proposition 6.2** (Second-order GARCH explicit put option price).**

The second-order put option price in the GARCH diffusion model can be written as

[TABLE]

where

[TABLE]

Furthermore $\hat{x}=S_{0}$ and

[TABLE]

Remark 6.1** (Fast calibration).**

Since the integral operators obey a recursion, we can exploit dynamic programming to greatly speed up the calibration procedure. To implement our fast calibration scheme, one executes the following algorithm. Let $\mu_{t}\equiv\mu=(\mu^{(1)},\mu^{(2)},\dots,\mu^{(n)})$ be an arbitrary set of parameters and denote by $\omega_{t}$ an arbitrary integral operator.

• Calibrate $\mu$ over $[0,T_{1})$ to obtain $\mu_{0}$. This involves computing $\omega_{T_{1}}$.
• Calibrate $\mu$ over $[T_{1},T_{2})$ to obtain $\mu_{1}$. This involves computing $\omega_{T_{2}}$ which is in terms of $\omega_{T_{1}}$, the latter already being computed in the previous step.
• Repeat until time $T_{N}$.

7. Numerical tests and sensitivity analysis

In this section, we will numerically investigate the accuracy of our closed-form approximation formula in the Heston and GARCH diffusion models (Proposition 6.1 and Proposition 6.2 respectively) by considering the sensitivity of our approximation formula whilst varying one parameter at a time, with the others all being fixed. Namely, for an arbitrary set of piecewise-constant parameters $(\kappa_{t},\theta_{t},\lambda_{t},\rho_{t})\equiv(\kappa,\theta,\lambda,\rho)=:\mu$, we vary only one of $\kappa$, $\theta$, $\lambda$, $\rho$ at a time and keep the rest fixed. Then, we will compute the difference (signed error) in implied volatilities via our approximation formula as well as the Monte-Carlo for maturity times $T\in\{1/12,3/12,6/12,1\}\equiv\{1\text{M},3\text{M},6\text{M},1\text{Y}\}$ and strikes corresponding to Put 10, 25, and ATM deltas. Thus, the signed error of the implied volatility for a given parameter set $\mu$, maturity $T$, and strike $K$ is

[TABLE]

The Monte-Carlo simulation is implemented by taking advantage of a variation of the mixing solution methodology from Appendix A. Namely, consider the log-spot $X_{t}=\ln S_{t}$ and log-strike $k=\ln K$. Let $\sigma$ be the variance process. Then

[TABLE]

where

[TABLE]

We will not provide a proof of this result, as it can be derived by adapting the proof of Theorem A.1 to the case of modelling the log-spot in a straightforward manner. Using this relationship for Monte-Carlo simulation, there is no need to simulate $S$, one need only simulate $\sigma$. This reduces the runtime as well as standard error of the procedure.

We will employ the usual Euler-Maruyama method to simulate the variance process. For all our simulations, we use 2,000,000 Monte Carlo paths, and 24 time steps per day, where a year is comprised of 252 trading days. This is to reduce the Monte-Carlo and discretisation errors sufficiently well.

Remark 7.1**.**

The code utilised to obtain the numerical results in this section is available at GitHub [13]. In particular, what is provided is:

• A routine which computes our closed-form approximation of put option prices for both the Heston model and GARCH diffusion model (the GARCH diffusion model with $\rho=0$ a.e.) with piecewise-constant parameter inputs.
• A routine which implements the Monte-Carlo simulation via the mixing solution methodology for the pricing of put option prices in the Heston model and GARCH diffusion model with piecewise-constant parameter inputs.
• A routine which compares the accuracy and runtimes of the aforementioned methods.

7.1. Heston sensitivity analysis

Remark 7.2**.**

A piecewise-constant parameter which is piecewise-constant over $n$ intervals will be called an $n$-piece parameter or a parameter with $n$ pieces.

We start from a “safe” parameter set, which are constant parameters calibrated by Bloomberg USD/JPY FX option price data [7]. The safe parameter set is

[TABLE]

with

[TABLE]

However, we would like to consider piecewise-constant parameters, as our closed-form approximation method is designed in order to take advantage of piecewise-constant parameter inputs. To do so, we will make our parameters piecewise-constant over three time intervals, the length of the first, second and third time interval having proportions $1/4$, $1/4$, $1/2$ of the maturity time $T$ respectively. For example, if $T=1/12$ then a 3-piece parameter is piecewise-constant on the time intervals $[0,1/48),[1/48,2/48),[2/48,4/48)$. To choose the “safe” values of the 3-piece parameters, we will simply (albeit somewhat artificially) perturb the value of each safe parameter over each time interval, yielding the following table of “safe” 3-piece parameter sets:

[TABLE]

In our numerical analysis, we vary one of the 3-piece parameters $(\kappa,\theta,\lambda,\rho)$ with the rest fixed, and then compute implied volatilities via both the closed-form approximation formula as well as the Monte-Carlo method as described above. Specifically, we select a 3-piece parameter from $(\kappa,\theta,\lambda,\rho)$, and start at 40% of its safe 3-piece parameter value, then increase the value of each piece in increments of 20%, all the way up to 160%, whilst keeping the other three 3-piece parameters fixed at their safe value. We then repeat this process with each of the other three 3-piece parameters. The relevant tables are Table 7.1, Table 7.2, Table 7.3, and Table 7.4, for the analysis of $\kappa$, $\theta$, $\lambda$, and $\rho$ respectively. Note that the $\rho$ table values are reversed, as $\rho$ is negative, and thus 160% of the safe 3-piece $\rho$ is approximately $(-0.594,-0.658,-0.626)$. This ensures that the parameter values are increasing for all tables.

Remark 7.3**.**

The Feller condition is

[TABLE]

Tables 7.1 to 7.4 use red text to indicate when the Feller condition is not satisfied. Note that in application, this condition is almost always violated. That is, parameters calibrated from market data almost always violate the Feller condition, see for example [11, 12, 31]. In fact, this is a reason why practitioners now favour non-affine models.

0

0

0

The sensitivity analysis is consistent with what we expect. For example, for large maturity $T$, it can be easily shown that the component-wise variance of the difference in the expansion and evaluation point increases. Moreover, for large vol-of-vol $\lambda$ or large correlation $|\rho|$, one would expect the variance of the difference in the expansion and evaluation point to increase. Thus, when these parameters are large, we expect our closed-form approximation formula to be less accurate. Indeed, the numerical results confirm this behaviour. However, there seems to be a peculiarity when comparing the numerical results pertaining to the maturities $T=1/2$ and $T=1$. Specifically, it seems that the approximation formula for $T=1$ is often more accurate than for $T=1/2$. By studying the error bound from Theorem 5.4, one explanation for this behaviour is that the terms that are exponentially decaying in $T$ from the function $\zeta(T)$ (recall this is the function $M(T,K)$ for fixed $K$, see Lemma 5.2) may start to contribute more around this region of $T$. Of course, this is just conjecture, as such a claim is dependent on all the other terms in the error bound Theorem 5.4 (most notably, the behaviour of the moments of the underlying variance process). However, it seems like a reasonable claim, as the other terms in the error bound are all increasing functions of $T$. A precise quantification of this claim would require further investigation of the moments of the underlying variance process. In terms of application, for realistic parameter values (that is, values around the “safe” 3-piece parameter sets) we see that the magnitude of error is around 10-50bps, which would be deemed reasonable for use in industry.

7.2. GARCH diffusion model sensitivity analysis

We start from the same “safe” 3-piece parameter set from Remark 7.2, albeit with $\rho=0$ on each piece. In our analysis, we vary one of the 3-piece parameters $(\kappa,\theta,\lambda)$ at a time and keep the rest fixed. We employ the same strategy as in the Heston framework, we start at 40% of the safe 3-piece parameter value, and increase it in increments of 20%, all the way up to 160%, whilst the rest are fixed at their safe value. The relevant tables are Table 7.5, Table 7.6, and Table 7.7 for the analysis of $\kappa$, $\theta$, and $\rho$ respectively.

0

0

The error in implied volatility for the closed-form approximation formula in the GARCH diffusion model behaves well, with most being less than 3bp in magnitude. In contrast to the Heston model analysis, this is most likely due to the fact that the correlation $\rho$ is assumed to be 0 always, and hence one would expect the variance of the difference in the expansion and evaluation point to be smaller. Otherwise, the approximation behaves as we expect, with errors being larger for large maturity $T$ and vol-of-vol $\lambda$, as the variance of the difference in the expansion and evaluation point will grow with these parameters. Again, there is an anomaly when comparing the numerical results pertaining to the maturities $T=1/2$ and $T=1$, as the latter seems to be more accurate. For similar reasons to the Heston model case, we attribute this behaviour to the function $\zeta(T)$ (see Lemma 5.2 and Corollary 5.3).

7.3. Computational runtime analysis

The use of our closed-form approximation formula yields approximations of put option prices essentially instantaneously. However, it is worth investigating the speed up/slow down when considering $n$-piece parameters for different $n$. We will consider pricing a put option with maturities 6M and 1Y, whilst varying the number of pieces the parameters possess in both the Heston and GARCH diffusion models va our closed-form approximation formulas. The values of the parameters do not affect the runtime, and thus we have not provided them in our numerical tests. Each runtime we display will be the average of 100 runtimes with identical parameters.

The Heston model runtime analysis is given in Table 7.8, whilst the GARCH diffusion model runtime analysis is given in Table 7.9. Unsurprisingly, the maturity value does not really affect the runtime, as varying the maturity does not contribute to any computational burden in the formula. For 1-piece parameters (that is, constant parameters), the runtimes of the approximation formulas for both the Heston and GARCH diffusion models are almost the same (just less than a millisecond). However, for parameters with a large number of pieces, the runtime of the GARCH diffusion model formula is longer than for the Heston model. For example, for 10-piece parameters, running the Heston model closed-form approximation formula takes around 40 milliseconds, whereas the GARCH diffusion model one takes around 65 milliseconds. This may seem surprising, as one may expect that the GARCH diffusion model closed-form approximation formula would be less costly to run, since here we are in the simpler setting where $\rho=0$ a.e.. However, the reason why the GARCH diffusion model’s closed-form approximation formula is computationally more costly is because it possesses a 5-fold integral operator term, whereas the Heston model closed-form approximation formula has only 4-fold and below integral operators.

8. Conclusion

We have derived closed-form approximation formulas for European put option prices in the context of stochastic volatility models with time-dependent parameters. Our method involves a second-order Taylor expansion of the mixing solution, followed by simplification of a number of expectations via the use of change of measure techniques. Such a method has been considered by Drimus with respect to the Heston model with constant parameters [14]. We extend this method as we consider time-dependent parameters, which requires an additional approximation of one of the expectations. Furthermore, we obtain the expression for the error term induced by the expansion. We give general bounds on the induced error term in terms of higher order moments of the underlying variance process. Under the Heston framework, we find that all terms induced by the expansion are able to be expressed as iterated integrals. Furthermore, we attempt to generalise this method to the GARCH diffusion model. We show that obtaining an expression for the solution or moments of the resulting variance process after the change of measure is non-trivial to achieve. By assuming $\rho=0$ a.e. in the GARCH diffusion model, we are able to address this problem, albeit with the added assumption of uncorrelated spot and volatility movements. Additionally, we show that the iterated integrals obey a convenient recursive property when parameters are piecewise-constant. By doing so, the approximation formulas become closed-form. In addition, we devise a fast calibration scheme which exploits the recursive property of our integral operators. Lastly, we perform a numerical error and sensitivity analysis to investigate the quality of our approximation in the Heston and GARCH diffusion models. We find that the error is well within the range appropriate for application purposes and behaves as we expect for certain parameter values, such as long maturity, large vol-of-vol and large correlation.

The purely probabilistic mixing solution approach, which is the backbone of our expansion method, is very appealing due to its generality and ability to handle time-dependent parameters. Further research would be needed to combine it, with no correlation restriction, with the type of non-affine stochastic volatility models favoured by practitioners.

Appendix A Mixing solution

In this appendix, we give a derivation of the result referred to as the mixing solution by [19]. This result is crucial for the expansion methodology implemented in Section 2. Hull and White first established the expression for the case of independent Brownian motions $W$ and $B$. Later on, this was extended for the correlated Brownian motions case, see [36, 32].

Under a domestic risk-neutral measure $\mathbb{Q}$, suppose that the spot $S$ with variance $\sigma$ follows the dynamics

[TABLE]

Theorem A.1** (Mixing solution).**

Let $\text{Put}=e^{-\int_{0}^{T}r_{t}^{d}\mathrm{d}t}\mathbb{E}(K-S_{T})_{+}$. Then

[TABLE]

where $\text{Put}_{\text{BS}}$ is given in eq. 2.4.

Proof.

By writing the driving Brownian motion of the spot as $W_{t}=\int_{0}^{t}\rho_{u}\mathrm{d}B_{u}+\int_{0}^{t}\sqrt{1-\rho_{u}^{2}}\mathrm{d}Z_{u}$, where $Z$ is a Brownian motion under $\mathbb{Q}$ which is independent of $B$, this gives the explicit pathwise unique strong solution of $S$ as

[TABLE]

First, notice that both $\sigma$ and $\xi$ are adapted to the filtration $(\mathcal{F}_{t}^{B})_{0\leq t\leq T}$. Thus, it is evident that $S_{T}|\mathcal{F}_{T}^{B}$ will have a log-normal distribution, namely

[TABLE]

Hence, the calculation of e^{-\int_{0}^{T}r_{t}^{d}\mathrm{d}t}\mathbb{E}\big{[}(K-S_{T})_{+}|\mathcal{F}_{T}^{B}\big{]} will result in a Black-Scholes like formula.

[TABLE]

It is immediate that e^{-\int_{0}^{T}r_{t}^{d}\mathrm{d}t}\mathbb{E}\big{[}(K-S_{T})_{+}|\mathcal{F}_{T}^{B}\big{]}=\text{Put}_{\text{BS}}\left(S_{0}\xi_{T},\tilde{\sigma}^{2}(T)\right). ∎

Appendix B $\text{Put}_{\text{BS}}$ partial derivatives

This appendix contains some partial derivatives for the Black-Scholes put option formula $\text{Put}_{\text{BS}}$, given in eq. 2.4. One can think of these partial derivatives as being analogous to Black-Scholes Greeks. However, these are slightly different as our Black-Scholes formulas are parameterised with respect to integrated variance rather than volatility.

B.1. First-order $\textnormal{Put}_{\textnormal{BS}}$

[TABLE]

B.2. Second-order $\textnormal{Put}_{\textnormal{BS}}$

[TABLE]

B.3. Third-order $\textnormal{Put}_{\textnormal{BS}}$

[TABLE]

B.4. Fourth-order $\textnormal{Put}_{\textnormal{BS}}$

[TABLE]

Appendix C Greeks

Expressions for second-order approximations of put option Greeks can be obtained via simple partial differentiation of $\text{Put}^{(2)}$ (eq. 2.5).

The Put Delta approximation is obtained via partial differentiation of $\text{Put}^{(2)}$ with respect to the underlying $S_{0}$.

[TABLE]

The Put Gamma approximation is obtained via partial differentiation of the Put Delta approximation with respect to the underlying $S_{0}$.

[TABLE]

Notice that the above expectations are the same as those from eqs. 2.6 to 2.8. These expectations are made explicit in Section 3.

Appendix D Calculation of moments

In this appendix, we derive expressions for some of the moments, mixed moments and covariances of the CIR and IGa processes utilised in this article. The results for the CIR process are well known, however, we have not found a source for the IGa moments with time-dependent parameters in the literature.

D.1. Moments of the CIR process

Let $V$ be a CIR$(v_{0};\kappa_{t},\theta_{t},\lambda_{t})$. It satisfies the SDE

[TABLE]

where we assume $(\kappa_{t})_{0\leq t\leq T},(\theta_{t})_{0\leq t\leq T}$ and $(\lambda_{t})_{0\leq t\leq T}$ are time-dependent, deterministic, strictly positive and bounded on $[0,T]$. For $s<t$, it can be integrated to obtain

[TABLE]

In particular, for $s=0$,

[TABLE]

Proposition D.1**.**

$V$ has the following moments:

[TABLE]

all for $m,n\geq 1$ and $s<t$.

Proof.

We give an outline for obtaining $\mathrm{Var}(V_{t})$ and $\mathrm{Cov}(V_{s},V_{t})$. The other terms follow a similar methodology. Notice that $\mathrm{Var}(V_{t})=\mathbb{E}(V_{t}-\mathbb{E}(V_{t}))^{2}$. Then using eq. D.2 and $\mathbb{E}(V_{t})$,

[TABLE]

Assume $s<t$. Using the representation of $V_{t}$ in terms of $V_{s}$ eq. D.1, we have

[TABLE]

where we have used that $V_{s}$ is independent of the Itô integral $\int_{s}^{t}e^{-\int_{u}^{t}\kappa_{z}\mathrm{d}z}\lambda_{u}\sqrt{V_{u}}\mathrm{d}B_{u}$. ∎

D.2. Moments of the IGa process

Let $V$ be an IGa$(v_{0};\kappa_{t},\theta_{t},\lambda_{t})$. It satisfies the SDE

[TABLE]

where we assume $(\kappa_{t})_{0\leq t\leq T},(\theta_{t})_{0\leq t\leq T}$ and $(\lambda_{t})_{0\leq t\leq T}$ are time-dependent, deterministic, strictly positive and bounded on $[0,T]$. Then for $s<t$, $V$ has the explicit strong solution

[TABLE]

In particular, for $s=0$,

[TABLE]

Proposition D.2**.**

$V$ has the following moments:

[TABLE]

all for $m,n\geq 1$ and $s<t$.

Proof.

We show how to obtain $\mathbb{E}(V_{s}^{n}V_{t}^{m})$. The other terms follow a similar methodology. We consider the differential of $V^{n}$.

[TABLE]

Multiplying both sides by $V_{s}^{m}$ and taking expectation yields

[TABLE]

Differentiating both sides in $t$ and letting $M^{m,n}_{s}(t):=\mathbb{E}(V_{s}^{m}V_{t}^{n})$, then

[TABLE]

This is a first order ODE, which can be solved with the integrating factor method by integrating from $s$ to $t$. ∎

Appendix E Solutions to SDEs with linear diffusion

Suppose the diffusion $U$ solves the SDE

[TABLE]

where $(\nu_{t})_{0\leq t\leq T}$ is adapted to the Brownian filtration and $f$ and $\nu$ satisfy some regularity conditions so that a pathwise unique strong solution for $U$ exists. For example, $f$ Lipschitz in $x$ uniformly in $t$, and $\nu$ bounded on $[0,T]$ is sufficient.

Proposition E.1**.**

The solution to eq. E.1 can be expressed as

[TABLE]

where $F$ is a GBM$(1;\nu_{t}^{2},-\nu_{t})$, That is,

[TABLE]

and $Y$ solves the integral equation (written in differential form)

[TABLE]

Proof.

We essentially verify that this form of $U$ satisfies the SDE eq. E.1.

[TABLE]

∎